-
Notifications
You must be signed in to change notification settings - Fork 0
/
interpolation.f
executable file
·153 lines (142 loc) · 3.86 KB
/
interpolation.f
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
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
module interpolation
logical :: ffc_control
private
public fival,cfival,FFC,FFR4,ffc_control
contains
!************************************************************************
!* REAL 4-point lagrange interpolation routine.
!* interpolates thr FUNCTION value fival at point r from an
!* array of points stored in fdis(ndm). this array is assumed
!* to be defined such that the first element fdis(1) CONTAINS
!* the FUNCTION value at r=xv(1) and xv(2 .. ndm) are monotonically
!* increasing.
!************************************************************************
FUNCTION cfival(r,xv,fdis,ndm,alpha)
IMPLICIT REAL*8(A-H,O-Z)
COMPLEX*16 cfival,fdis(ndm),y1,y2,y3,y4
DIMENSION xv(ndm)
IF(r.GT.xv(ndm)) go to 9
DO 5 k=1,ndm-2
5 IF(r.LT.xv(k)) go to 6
k=ndm-2
6 nst=MAX(k-1,1)
x1=xv(nst)
x2=xv(nst+1)
x3=xv(nst+2)
x4=xv(nst+3)
y1=fdis(nst+0)
y2=fdis(nst+1)
y3=fdis(nst+2)
y4=fdis(nst+3)
pii1=(x1-x2)*(x1-x3)*(x1-x4)
pii2=(x2-x1)*(x2-x3)*(x2-x4)
pii3=(x3-x1)*(x3-x2)*(x3-x4)
pii4=(x4-x1)*(x4-x2)*(x4-x3)
xd1=r-x1
xd2=r-x2
xd3=r-x3
xd4=r-x4
pi1=xd2*xd3*xd4
pi2=xd1*xd3*xd4
pi3=xd1*xd2*xd4
pi4=xd1*xd2*xd3
cfival=y1*pi1/pii1+y2*pi2/pii2+y3*pi3/pii3+y4*pi4/pii4
RETURN
9 cfival=fdis(ndm) * EXP(alpha*(xv(ndm)-r))
RETURN
END
c-----------------------------------------------------------------------
************************************************************************
* REAL 4-point lagrange interpolation routine.
* interpolates thr FUNCTION value fival at point r from an
* array of points stored in fdis(ndm). this array is assumed
* to be defined such that the first element fdis(1) CONTAINS
* the FUNCTION value at r=xv(1) and xv(2 .. ndm) are monotonically
* increasing.
************************************************************************
FUNCTION fival(r,xv,fdis,ndm,alpha)
IMPLICIT REAL*8(A-H,O-Z)
REAL*8 fdis(ndm),y1,y2,y3,y4
DIMENSION xv(ndm)
IF(r.GT.xv(ndm)) go to 9
DO 5 k=1,ndm-2
5 IF(r.LT.xv(k)) go to 6
k=ndm-2
6 nst=MAX(k-1,1)
x1=xv(nst)
x2=xv(nst+1)
x3=xv(nst+2)
x4=xv(nst+3)
y1=fdis(nst+0)
y2=fdis(nst+1)
y3=fdis(nst+2)
y4=fdis(nst+3)
pii1=(x1-x2)*(x1-x3)*(x1-x4)
pii2=(x2-x1)*(x2-x3)*(x2-x4)
pii3=(x3-x1)*(x3-x2)*(x3-x4)
pii4=(x4-x1)*(x4-x2)*(x4-x3)
xd1=r-x1
xd2=r-x2
xd3=r-x3
xd4=r-x4
pi1=xd2*xd3*xd4
pi2=xd1*xd3*xd4
pi3=xd1*xd2*xd4
pi4=xd1*xd2*xd3
fival=y1*pi1/pii1+y2*pi2/pii2+y3*pi3/pii3+y4*pi4/pii4
RETURN
9 fival=fdis(ndm) * EXP(alpha*(xv(ndm)-r))
RETURN
END
c-----------------------------------------------------------------------
FUNCTION FFC(PP,F,N)
COMPLEX*16 FFC,F(N)
REAL*8 PP
PARAMETER(X=.16666666666667)
I=PP
IF(I.LE.0) GO TO 2
IF(I.GE.N-2) GO TO 4
1 P=PP-I
P1=P-1.
P2=P-2.
Q=P+1.
FFC=(-P2*F(I)+Q*F(I+3))*(P*P1*X)+(P1*F(I+1)-P*F(I+2))*(Q*P2*.5)
RETURN
2 IF(I.LT.0) GO TO 3
I=1
GO TO 1
3 FFC=F(1)
RETURN
4 IF(I.GT.N-2) GO TO 5
I=N-3
GO TO 1
5 FFC=F(N)
RETURN
END
FUNCTION FFR4(Y,F,N)
IMPLICIT REAL*8(A-H,O-Z)
REAL*8 F(N),P,P1,P2,Q,X,FFR4
REAL*8 Y
PARAMETER(X=.16666666666667)
P=Y
I=P
IF(I.LE.0) GO TO 2
IF(I.GE.N-2) GO TO 4
1 P=P-I
P1=P-1.
P2=P-2.
Q=P+1.
FFR4=(-P2*F(I)+Q*F(I+3))*P*P1*X+(P1*F(I+1)-P*F(I+2))*Q*P2*.5
RETURN
2 IF(I.LT.0) GO TO 3
I=1
GO TO 1
3 FFR4=F(1)
RETURN
4 IF(I.GT.N-2) GO TO 5
I=N-3
GO TO 1
5 FFR4=F(N)
RETURN
END
end module