ddatrp.f

Go to the documentation of this file.

      SUBROUTINE ddatrp (X, XOUT, YOUT, YPOUT, NEQ, KOLD, PHI, PSI)
C***BEGIN PROLOGUE  DDATRP
C***SUBSIDIARY
C***PURPOSE  Interpolation routine for DDASSL.
C***LIBRARY   SLATEC (DASSL)
C***TYPE      DOUBLE PRECISION (SDATRP-S, DDATRP-D)
C***AUTHOR  PETZOLD, LINDA R., (LLNL)
C***DESCRIPTION
C-----------------------------------------------------------------------
C     THE METHODS IN SUBROUTINE DDASTP USE POLYNOMIALS
C     TO APPROXIMATE THE SOLUTION. DDATRP APPROXIMATES THE
C     SOLUTION AND ITS DERIVATIVE AT TIME XOUT BY EVALUATING
C     ONE OF THESE POLYNOMIALS,AND ITS DERIVATIVE,THERE.
C     INFORMATION DEFINING THIS POLYNOMIAL IS PASSED FROM
C     DDASTP, SO DDATRP CANNOT BE USED ALONE.
C
C     THE PARAMETERS ARE:
C     X     THE CURRENT TIME IN THE INTEGRATION.
C     XOUT  THE TIME AT WHICH THE SOLUTION IS DESIRED
C     YOUT  THE INTERPOLATED APPROXIMATION TO Y AT XOUT
C           (THIS IS OUTPUT)
C     YPOUT THE INTERPOLATED APPROXIMATION TO YPRIME AT XOUT
C           (THIS IS OUTPUT)
C     NEQ   NUMBER OF EQUATIONS
C     KOLD  ORDER USED ON LAST SUCCESSFUL STEP
C     PHI   ARRAY OF SCALED DIVIDED DIFFERENCES OF Y
C     PSI   ARRAY OF PAST STEPSIZE HISTORY
C-----------------------------------------------------------------------
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   830315  DATE WRITTEN
C   901009  Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
C   901019  Merged changes made by C. Ulrich with SLATEC 4.0 format.
C   901026  Added explicit declarations for all variables and minor
C           cosmetic changes to prologue.  (FNF)
C***END PROLOGUE  DDATRP
C
      INTEGER  neq, kold
      DOUBLE PRECISION  x, xout, yout(*), ypout(*), phi(neq,*), psi(*)
C
      INTEGER  i, j, koldp1
      DOUBLE PRECISION  c, d, gamma, temp1
C
C***FIRST EXECUTABLE STATEMENT  DDATRP
      koldp1=kold+1
      temp1=xout-x
      DO 10 i=1,neq
         yout(i)=phi(i,1)
10       ypout(i)=0.0d0
      c=1.0d0
      d=0.0d0
      gamma=temp1/psi(1)
      DO 30 j=2,koldp1
         d=d*gamma+c/psi(j-1)
         c=c*gamma
         gamma=(temp1+psi(j-1))/psi(j)
         DO 20 i=1,neq
            yout(i)=yout(i)+c*phi(i,j)
20          ypout(i)=ypout(i)+d*phi(i,j)
30       CONTINUE
      RETURN
C
C------END OF SUBROUTINE DDATRP------
      END