]>
Commit | Line | Data |
---|---|---|
21886bb6 | 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 |