home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / odepack / intdy.f < prev    next >
Text File  |  1996-09-28  |  3KB  |  85 lines
  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 = DBLE(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 = DBLE(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(30HINTDY--  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(30HINTDY--  T (=R1) ILLEGAL      ,
  77.      1   30, 52, 0, 0, 0, 0, 1, T, 0.0D0)
  78.       CALL XERRWV(
  79.      1  60H      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 
  85.