1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
|
C SUBROUTINE INTERP.F
C
C SOURCE
C BEVINGTON, PAGES 266-267.
C
C PURPOSE
C INTERPOLATE BETWEEN DATA POINTS TO EVALUATE A FUNCTION
C
C USAGE
C CALL INTERP (X, Y, NPTS, NTERMS, XIN, YOUT)
C
C DESCRIPTION OF PARAMETERS
C X - ARRAY OF DATA POINTS FOR INDEPENDENT VARIABLE
C Y - ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE
C NPTS - NUMBER OF PAIRS OF DATA POINTS
C NTERMS - NUMBER OF TERMS IN FITTING POLYNOMIAL
C XIN - INPUT VALUE OF X
C YOUT - INTERPOLATED VALUE OF Y
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C NONE
C
C COMMENTS
C DIMENSION STATEMENT VALID FOR NTERMS UP TO 10
C VALUE OF NTERMS MAY BE MODIFIED BY THE PROGRAM
C
SUBROUTINE INTERP (X,Y,NPTS,NTERMS,XIN,YOUT)
DOUBLE PRECISION DELTAX,DELTA,A,PROD,SUM
DIMENSION X(1),Y(1)
DIMENSION DELTA(10),A(10)
C
C SEARCH FOR APPROPRIATE VALUE OF X(1)
C
11 DO 19 I=1,NPTS
IF (XIN-X(I)) 13,17,19
13 I1=I-NTERMS/2
IF (I1) 15,15,21
15 I1=1
GOTO 21
17 YOUT=Y(I)
18 GOTO 61
19 CONTINUE
I1=NPTS-NTERMS+1
21 I2=I1+NTERMS-1
IF (NPTS-I2) 23,31,31
23 I2=NPTS
I1=I2-NTERMS+1
25 IF (I1) 26,26,31
26 I1=1
27 NTERMS=I2-I1+1
C
C EVALUATE DEVIATIONS DELTA
C
31 DENOM=X(I1+1)-X(I1)
DELTAX=(XIN-X(I1))/DENOM
DO 35 I=1,NTERMS
IX=I1+I-1
35 DELTA(I)=(X(IX)-X(I1))/DENOM
C
C ACCUMULATE COEFFICIENTS A
C
40 A(1)=Y(I1)
41 DO 50 K=2,NTERMS
PROD=1.
SUM=0.
IMAX=K-1
IXMAX=I1+IMAX
DO 49 I=1,IMAX
J=K-I
PROD=PROD*(DELTA(K)-DELTA(J))
49 SUM=SUM-A(J)/PROD
50 A(K)=SUM+Y(IXMAX)/PROD
C
C ACCUMULATE SUM OF EXPANSION
C
51 SUM=A(1)
DO 57 J=2,NTERMS
PROD=1.
IMAX=J-1
DO 56 I=1,IMAX
56 PROD=PROD*(DELTAX-DELTA(I))
57 SUM=SUM+A(J)*PROD
60 YOUT=SUM
61 RETURN
END
|
