← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to scipy/integrate/odepack/intdy.f

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
scipy/integrate/odepack/intdy.f
 
1
      subroutine intdy (t, k, yh, nyh, dky, iflag)
 
2
clll. optimize
 
3
      integer k, nyh, iflag
 
4
      integer iownd, iowns,
 
5
     1   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
 
6
     2   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
 
7
      integer i, ic, j, jb, jb2, jj, jj1, jp1
 
8
      double precision t, yh, dky
 
9
      double precision rowns,
 
10
     1   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
 
11
      double precision c, r, s, tp
 
12
      dimension yh(nyh,1), dky(1)
 
13
      common /ls0001/ rowns(209),
 
14
     2   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
 
15
     3   iownd(14), iowns(6),
 
16
     4   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
 
17
     5   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
 
18
c-----------------------------------------------------------------------
 
19
c intdy computes interpolated values of the k-th derivative of the
 
20
c dependent variable vector y, and stores it in dky.  this routine
 
21
c is called within the package with k = 0 and t = tout, but may
 
22
c also be called by the user for any k up to the current order.
 
23
c (see detailed instructions in the usage documentation.)
 
24
c-----------------------------------------------------------------------
 
25
c the computed values in dky are gotten by interpolation using the
 
26
c nordsieck history array yh.  this array corresponds uniquely to a
 
27
c vector-valued polynomial of degree nqcur or less, and dky is set
 
28
c to the k-th derivative of this polynomial at t.
 
29
c the formula for dky is..
 
30
c              q
 
31
c  dky(i)  =  sum  c(j,k) * (t - tn)**(j-k) * h**(-j) * yh(i,j+1)
 
32
c             j=k
 
33
c where  c(j,k) = j*(j-1)*...*(j-k+1), q = nqcur, tn = tcur, h = hcur.
 
34
c the quantities  nq = nqcur, l = nq+1, n = neq, tn, and h are
 
35
c communicated by common.  the above sum is done in reverse order.
 
36
c iflag is returned negative if either k or t is out of bounds.
 
37
c-----------------------------------------------------------------------
 
38
      iflag = 0
 
39
      if (k .lt. 0 .or. k .gt. nq) go to 80
 
40
      tp = tn - hu -  100.0d0*uround*(tn + hu)
 
41
      if ((t-tp)*(t-tn) .gt. 0.0d0) go to 90
 
42
c
 
43
      s = (t - tn)/h
 
44
      ic = 1
 
45
      if (k .eq. 0) go to 15
 
46
      jj1 = l - k
 
47
      do 10 jj = jj1,nq
 
48
 10     ic = ic*jj
 
49
 15   c = dfloat(ic)
 
50
      do 20 i = 1,n
 
51
 20     dky(i) = c*yh(i,l)
 
52
      if (k .eq. nq) go to 55
 
53
      jb2 = nq - k
 
54
      do 50 jb = 1,jb2
 
55
        j = nq - jb
 
56
        jp1 = j + 1
 
57
        ic = 1
 
58
        if (k .eq. 0) go to 35
 
59
        jj1 = jp1 - k
 
60
        do 30 jj = jj1,j
 
61
 30       ic = ic*jj
 
62
 35     c = dfloat(ic)
 
63
        do 40 i = 1,n
 
64
 40       dky(i) = c*yh(i,jp1) + s*dky(i)
 
65
 50     continue
 
66
      if (k .eq. 0) return
 
67
 55   r = h**(-k)
 
68
      do 60 i = 1,n
 
69
 60     dky(i) = r*dky(i)
 
70
      return
 
71
c
 
72
 80   call xerrwv('intdy--  k (=i1) illegal      ',
 
73
     1   30, 51, 0, 1, k, 0, 0, 0.0d0, 0.0d0)
 
74
      iflag = -1
 
75
      return
 
76
 90   call xerrwv('intdy--  t (=r1) illegal      ',
 
77
     1   30, 52, 0, 0, 0, 0, 1, t, 0.0d0)
 
78
      call xerrwv(
 
79
     1  '     t not in interval tcur - hu (= r1) to tcur (=r2)       ',
 
80
     1   60, 52, 0, 0, 0, 0, 2, tp, tn)
 
81
      iflag = -2
 
82
      return
 
83
c----------------------- end of subroutine intdy -----------------------
 
84
      end