polint.F

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file polint.F

C

      SUBROUTINE polint (XA,YA,N,X,Y,DY)
 
      IMPLICIT DOUBLE PRECISION (a-h, o-z)
C                                        Adapted from "Numerical Recipes" 
      parameter(nmax=10)
      dimension xa(n),ya(n),c(nmax),d(nmax)
      ns=1
      dif=abs(x-xa(1))
      DO 11 i=1,n
        dift=abs(x-xa(i))
        IF (dift.LT.dif) THEN
          ns=i
          dif=dift
        ENDIF
        c(i)=ya(i)
        d(i)=ya(i)
11    CONTINUE
      y=ya(ns)
      ns=ns-1
      DO 13 m=1,n-1
        DO 12 i=1,n-m
          ho=xa(i)-x
          hp=xa(i+m)-x
          w=c(i+1)-d(i)
          den=ho-hp
          IF(den.EQ.0.)pause
          den=w/den
          d(i)=hp*den
          c(i)=ho*den
12      CONTINUE
        IF (2*ns.LT.n-m)THEN
          dy=c(ns+1)
        ELSE
          dy=d(ns)
          ns=ns-1
        ENDIF
        y=y+dy
13    CONTINUE
      RETURN
      END