PDF/spdf/dpolin.F

   1 #include "pdf/pilot.h"
   2       SUBROUTINE DPOLIN(XA,YA,N,X,Y,DY)
   3 #include "pdf/impdp.inc"
   4       PARAMETER (NMAX=10)
   5       DIMENSION XA(NMAX),YA(NMAX),C(NMAX),D(NMAX)
   6       Y=0.
   7       IF(N.GT.NMAX) RETURN
   8       NS=1
   9       DIF=ABS(X-XA(1))
  10       DO 11 I=1,N
  11         DIFT=ABS(X-XA(I))
  12         IF (DIFT.LT.DIF) THEN
  13           NS=I
  14           DIF=DIFT
  15         ENDIF
  16         C(I)=YA(I)
  17         D(I)=YA(I)
  18 11    CONTINUE
  19       Y=YA(NS)
  20       NS=NS-1
  21       DO 13 M=1,N-1
  22         DO 12 I=1,N-M
  23           HO=XA(I)-X
  24           HP=XA(I+M)-X
  25           W=C(I+1)-D(I)
  26           DEN=HO-HP
  27 C         IF(DEN.EQ.0.)PAUSE
  28           DEN=W/DEN
  29           D(I)=HP*DEN
  30           C(I)=HO*DEN
  31 12      CONTINUE
  32         IF (2*NS.LT.N-M)THEN
  33           DY=C(NS+1)
  34         ELSE
  35           DY=D(NS)
  36           NS=NS-1
  37         ENDIF
  38         Y=Y+DY
  39 13    CONTINUE
  40       RETURN
  41       END