1 SUBROUTINE SINTDY (T, K, YH, NYH, DKY, IFLAG) 2C***BEGIN PROLOGUE SINTDY 3C***SUBSIDIARY 4C***PURPOSE Interpolate solution derivatives. 5C***TYPE SINGLE PRECISION (SINTDY-S, DINTDY-D) 6C***AUTHOR Hindmarsh, Alan C., (LLNL) 7C***DESCRIPTION 8C 9C SINTDY computes interpolated values of the K-th derivative of the 10C dependent variable vector y, and stores it in DKY. This routine 11C is called within the package with K = 0 and T = TOUT, but may 12C also be called by the user for any K up to the current order. 13C (See detailed instructions in the usage documentation.) 14C 15C The computed values in DKY are gotten by interpolation using the 16C Nordsieck history array YH. This array corresponds uniquely to a 17C vector-valued polynomial of degree NQCUR or less, and DKY is set 18C to the K-th derivative of this polynomial at T. 19C The formula for DKY is: 20C q 21C DKY(i) = sum c(j,K) * (T - tn)**(j-K) * h**(-j) * YH(i,j+1) 22C j=K 23C where c(j,K) = j*(j-1)*...*(j-K+1), q = NQCUR, tn = TCUR, h = HCUR. 24C The quantities nq = NQCUR, l = nq+1, N = NEQ, tn, and h are 25C communicated by COMMON. The above sum is done in reverse order. 26C IFLAG is returned negative if either K or T is out of bounds. 27C 28C***SEE ALSO SLSODE 29C***ROUTINES CALLED XERRWV 30C***COMMON BLOCKS SLS001 31C***REVISION HISTORY (YYMMDD) 32C 791129 DATE WRITTEN 33C 890501 Modified prologue to SLATEC/LDOC format. (FNF) 34C 890503 Minor cosmetic changes. (FNF) 35C 930809 Renamed to allow single/double precision versions. (ACH) 36C 010412 Reduced size of Common block /SLS001/. (ACH) 37C 031105 Restored 'own' variables to Common block /SLS001/, to 38C enable interrupt/restart feature. (ACH) 39C 050427 Corrected roundoff decrement in TP. (ACH) 40C***END PROLOGUE SINTDY 41C**End 42 INTEGER K, NYH, IFLAG 43 REAL T, YH, DKY 44 DIMENSION YH(NYH,*), DKY(*) 45 INTEGER INIT, MXSTEP, MXHNIL, NHNIL, NSLAST, CNYH, 46 1 IALTH, IPUP, LMAX, MEO, NQNYH, NSLP, 47 1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 48 2 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 49 3 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU 50 REAL CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO, 51 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND 52 COMMON /SLS001/ CONIT, CRATE, EL(13), ELCO(13,12), 53 1 HOLD, RMAX, TESCO(3,12), 54 1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, 55 2 INIT, MXSTEP, MXHNIL, NHNIL, NSLAST, CNYH, 56 3 IALTH, IPUP, LMAX, MEO, NQNYH, NSLP, 57 3 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, 58 4 LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, 59 5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU 60 INTEGER I, IC, J, JB, JB2, JJ, JJ1, JP1 61 REAL C, R, S, TP 62 CHARACTER*80 MSG 63C 64C***FIRST EXECUTABLE STATEMENT SINTDY 65 IFLAG = 0 66 IF (K .LT. 0 .OR. K .GT. NQ) GO TO 80 67 TP = TN - HU - 100.0E0*UROUND*SIGN(ABS(TN) + ABS(HU), HU) 68 IF ((T-TP)*(T-TN) .GT. 0.0E0) GO TO 90 69C 70 S = (T - TN)/H 71 IC = 1 72 IF (K .EQ. 0) GO TO 15 73 JJ1 = L - K 74 DO 10 JJ = JJ1,NQ 75 10 IC = IC*JJ 76 15 C = IC 77 DO 20 I = 1,N 78 20 DKY(I) = C*YH(I,L) 79 IF (K .EQ. NQ) GO TO 55 80 JB2 = NQ - K 81 DO 50 JB = 1,JB2 82 J = NQ - JB 83 JP1 = J + 1 84 IC = 1 85 IF (K .EQ. 0) GO TO 35 86 JJ1 = JP1 - K 87 DO 30 JJ = JJ1,J 88 30 IC = IC*JJ 89 35 C = IC 90 DO 40 I = 1,N 91 40 DKY(I) = C*YH(I,JP1) + S*DKY(I) 92 50 CONTINUE 93 IF (K .EQ. 0) RETURN 94 55 R = H**(-K) 95 DO 60 I = 1,N 96 60 DKY(I) = R*DKY(I) 97 RETURN 98C 99 80 CALL XERRWD('SINTDY- K (=I1) illegal ', 100 1 30, 51, 0, 1, K, 0, 0, 0.0D0, 0.0D0) 101 IFLAG = -1 102 RETURN 103 90 CALL XERRWD('SINTDY- T (=R1) illegal ', 104 1 30, 52, 0, 0, 0, 0, 1, DBLE (T), 0.0D0) 105 CALL XERRWD( 106 1 ' T not in interval TCUR - HU (= R1) to TCUR (=R2) ', 107 1 60, 52, 0, 0, 0, 0, 2, DBLE (TP), DBLE (TN)) 108 IFLAG = -2 109 RETURN 110C----------------------- END OF SUBROUTINE SINTDY ---------------------- 111 END 112