intdy.f

Go to the documentation of this file.

      SUBROUTINE intdy (T, K, YH, NYH, DKY, IFLAG)
CLLL. OPTIMIZE
      INTEGER k, nyh, iflag
      INTEGER iownd, iowns,
     1   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
     2   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
      INTEGER i, ic, j, jb, jb2, jj, jj1, jp1
      DOUBLE PRECISION t, yh, dky
      DOUBLE PRECISION rowns, 
     1   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
      DOUBLE PRECISION c, r, s, tp
      dimension yh(nyh,*), dky(*)
      COMMON /ls0001/ rowns(209),
     2   ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
     3   iownd(14), iowns(6), 
     4   icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
     5   maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
C-----------------------------------------------------------------------
C INTDY COMPUTES INTERPOLATED VALUES OF THE K-TH DERIVATIVE OF THE
C DEPENDENT VARIABLE VECTOR Y, AND STORES IT IN DKY.  THIS ROUTINE
C IS CALLED WITHIN THE PACKAGE WITH K = 0 AND T = TOUT, BUT MAY
C ALSO BE CALLED BY THE USER FOR ANY K UP TO THE CURRENT ORDER.
C (SEE DETAILED INSTRUCTIONS IN THE USAGE DOCUMENTATION.)
C-----------------------------------------------------------------------
C THE COMPUTED VALUES IN DKY ARE GOTTEN BY INTERPOLATION USING THE
C NORDSIECK HISTORY ARRAY YH.  THIS ARRAY CORRESPONDS UNIQUELY TO A
C VECTOR-VALUED POLYNOMIAL OF DEGREE NQCUR OR LESS, AND DKY IS SET
C TO THE K-TH DERIVATIVE OF THIS POLYNOMIAL AT T. 
C THE FORMULA FOR DKY IS..
C              Q
C  DKY(I)  =  SUM  C(J,K) * (T - TN)**(J-K) * H**(-J) * YH(I,J+1)
C             J=K
C WHERE  C(J,K) = J*(J-1)*...*(J-K+1), Q = NQCUR, TN = TCUR, H = HCUR.
C THE QUANTITIES  NQ = NQCUR, L = NQ+1, N = NEQ, TN, AND H ARE
C COMMUNICATED BY COMMON.  THE ABOVE SUM IS DONE IN REVERSE ORDER.
C IFLAG IS RETURNED NEGATIVE IF EITHER K OR T IS OUT OF BOUNDS.
C-----------------------------------------------------------------------
      iflag = 0
      IF (k .LT. 0 .OR. k .GT. nq) go to 80
      tp = tn - hu -  100.0d0*uround*(tn + hu)
      IF ((t-tp)*(t-tn) .GT. 0.0d0) go to 90
C
      s = (t - tn)/h
      ic = 1
      IF (k .EQ. 0) go to 15
      jj1 = l - k
      DO 10 jj = jj1,nq
 10     ic = ic*jj
 15   c = dble(ic)
      DO 20 i = 1,n 
 20     dky(i) = c*yh(i,l)
      IF (k .EQ. nq) go to 55 
      jb2 = nq - k
      DO 50 jb = 1,jb2
        j = nq - jb 
        jp1 = j + 1 
        ic = 1
        IF (k .EQ. 0) go to 35
        jj1 = jp1 - k
        DO 30 jj = jj1,j
 30       ic = ic*jj
 35     c = dble(ic)
        DO 40 i = 1,n
 40       dky(i) = c*yh(i,jp1) + s*dky(i)
 50     CONTINUE
      IF (k .EQ. 0) RETURN
 55   r = h**(-k)
      DO 60 i = 1,n 
 60     dky(i) = r*dky(i)
      RETURN
C
 80   CALL xerrwd('INTDY--  K (=I1) ILLEGAL      ',
     1   30, 51, 0, 1, k, 0, 0, 0.0d0, 0.0d0)
      iflag = -1
      RETURN
 90   CALL xerrwd('INTDY--  T (=R1) ILLEGAL      ',
     1   30, 52, 0, 0, 0, 0, 1, t, 0.0d0)
      CALL xerrwd(
     1  '      T NOT IN INTERVAL TCUR - HU (= R1) TO TCUR (=R2)      ',
     1   60, 52, 0, 0, 0, 0, 2, tp, tn) 
      iflag = -2
      RETURN
C----------------------- END OF SUBROUTINE INTDY -----------------------
      END