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
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
|
subroutine intrp (itab, xtab, ytab, ntab, x, y, ierr)
c
c Interpolator using CODIM1 algorithm which is admittedly
c obscure but works well.
c
c itab - a label between 1 and 20 to identify the table and its
c most recent search index
c xtab - array of length ntab containing the x-values
c ytab - y-values
c ntab - number of x,y pairs in the table
c x - independent for which a y-value is desired
c y - returned interpolated (or extrapolated) value
c ierr - =0 for ok, -1 for extrapolation
c
double precision xtab(ntab), ytab(ntab), x, y
integer itab, ierr, index
double precision t(4), u(4)
c integer savind
c data savind/-1/
c
c----- Only 1 pt in table
if (ntab .eq. 1) then
y = ytab(1)
ierr = 0
return
endif
c
c-----
c Locate search index
call srch (itab, x, xtab, ntab, index, ierr)
c if (index .eq. savind) go to 2000
c savind = index
c
c-----
c Set interpolator index flags
i1 = 2
i2 = 3
iload = max0 (index-2, 1)
c
if (ntab .gt. 2) then
if (index.eq. 2) i2 = 4
c
if (index.eq.ntab) i1 = 1
endif
c
if (index.gt.2 .and. index.lt.ntab) then
i1 = 1
i2 = 4
endif
c-----
c Load interpolation arrays
do 1000 i = i1, i2
j = iload + (i-i1)
t(i) = xtab (j)
u(i) = ytab (j)
1000 continue
c
c-----
c Get interpolated value
2000 call codim1 (x, t, u, i1, i2, y)
return
end
c
c--------------------------------------------------------------
c
subroutine srch (itab, x, xtab, ntab, index, ierr)
c
c Search table of x-values to bracket the desired interpolant, x
c
c The returned search index will be:
c 2 - if extrapolation below the table is required
c ntab - above
c index - points to value just above x in the table if bounded.
c
c The index is saved as a starting point for subsequent entries
c in an array indexed through 'itab' which serves to label the
c set of saved search indices. Itab may be between 1 and 20.
c
c itab - The table identifier (1-20)
c x - The value for which an index is desired
c xtab - The table containing the x-values (array of length ntab)
c ntab - number of elements in the table
c index - returned index into the table (points just above x)
c ierr - 0 for ok, -1 for extrapolation
c
c Modified to remove entry points. 3/20/86 Valdes
c
integer ntab, index, init
double precision xtab(ntab), x
c
c common for subroutines intrp0 and intrpi
c
common /insvcm/ insave(20)
c
c initialize
data init/0/
c
c-----
c Initialize
if (init.eq.0) then
do 1110 i = 1, 20
1110 insave(i) = 0
init = 1
endif
c
c Determine direction of table, ascending or descending
idir = sign (1.0d0, xtab(ntab) - xtab(1))
c
c-----
c Reset error flag
ierr = 0
c
c-----
c Check for previous insaved index
last = insave(itab)
if (last .eq. 0 .or. last .gt. ntab) then
c
c-----
c no previous entry
isrch = 1
c check for extrapolation
if ((x-xtab( 1)) * idir .lt. 0.0d0) go to 2000
if ((x-xtab(ntab)) * idir .gt. 0.0d0) go to 2100
else
c
c-----
c previous entry left a valid index
isrch = last
c
c check for still wihin bounds - difference from above should be opposite
c sign of difference from below
c
if ((xtab(last)-x) * (xtab(last-1)-x) .lt. 0.0d0) then
index = last
return
endif
endif
c
c -----
c Begin searching - first determine direction
c
c This change made because x = xtab(1) was considered extrapolation.
c if ((x - xtab(isrch)) * idir .gt. 0.0d0) then
if ((x - xtab(isrch)) * idir .ge. 0.0d0) then
c forward
do 1100 i = isrch+1, ntab
if ((x-xtab(i)) * idir .gt. 0.0d0) go to 1100
go to 1500
1100 continue
c fall thru implies extrapolation required at high end
go to 2100
else
c
c-----
c negative direction search
do 1200 i = isrch-1,1,-1
if ((x-xtab(i)) * idir .lt. 0.0d0) go to 1200
go to 1400
1200 continue
c fall through implies extrapolation at low end
go to 2000
endif
c
c-----
c point has been bounded
1400 index = i + 1
go to 3000
1500 index = i
go to 3000
c
c-----
c extrapolations
2000 index = 2
ierr = -1
go to 3000
2100 index = ntab
ierr = -1
go to 3000
c
c-----
c insave index
3000 insave(itab) = index
end
c
c------
c Subroutine to reset saved index
subroutine intrp0 (itab)
integer itab
common /insvcm/ insave(20)
c
insave(itab) = 0
end
c
c-----
c Subroutine to return current index
subroutine intrpi (itab, ind)
integer itab, ind
common /insvcm/ insave(20)
c
ind = insave(itab)
end
c
c-------------------------------------------------------------------
c
subroutine codim1 (x, t, u, i1, i2, y)
c
c this subroutine performs an interposlation in a fashion
c not really understandable, but it works well.
c
c x - input independent variable
c t - array of 4 table independents surrounding x if possible
c u - array of 4 table dependents corresponding to the t array
c
c i1, i2 - indicators as follows:
c
c i1 = 1, i2 = 4 : 4 pts available in t and u arrays
c i1 = 1, i2 = 3 : 3 pts available (x near right edge of table)
c i1 = 2, i2 = 4 : (x near left edge of table)
c i1 = 2, i2 = 3 : 2 pts available
c i1 = 3, i3 = 3 : 1 pt available
c
c y - output interpolated (or extrapolated) dependent value
c
double precision t(4), u(4), x, y
integer i1, i2
double precision s, v, z, a1, a2, a3, c1, c2, c3, a4, c4, c5, c6
double precision e1, e2, p1, p2, slope1, slope2, al, bt, xe
c
c variable xk affects the extrapolation procedure. a value of -1.0
c appears to be a reliable value.
c
data xk/-1.0d0/
c
v = x
c the following code is extracted from an original source
c
a2=v-t(2)
al=a2/(t(3)-t(2))
s=al*u(3)+(1.-al)*u(2)
if(i1.gt.1.and.i2.lt.4)goto1530
a3=v-t(3)
if(i1.gt.1)goto1185
1180 a1=v-t(1)
c1=a2/(t(1)-t(2))*a3/(t(1)-t(3))
c2=a1/(t(2)-t(1))*a3/(t(2)-t(3))
c3=a1/(t(3)-t(1))*a2/(t(3)-t(2))
p1=c1*u(1)+c2*u(2)+c3*u(3)
if(i2.lt.4)goto1400
1185 a4=v-t(4)
c4=a3/(t(2)-t(3))*a4/(t(2)-t(4))
c5=a2/(t(3)-t(2))*a4/(t(3)-t(4))
c6=a2/(t(4)-t(2))*a3/(t(4)-t(3))
p2=c4*u(2)+c5*u(3)+c6*u(4)
if(i1.eq.1)goto1500
1200 if(xk.lt.0.)goto1230
xe=xk
goto1260
1230 slope1=abs((u(4)-u(3))/(t(4)-t(3)))
slope2=abs((u(3)-u(2))/(t(3)-t(2)))
xe=1.0d0
if(slope1+slope2.ne.0.)xe=1.-abs(slope1-slope2)/(slope1+slope2)
1260 p1=s+xe*(p2-s)
goto1500
1400 if(xk.lt.0.)goto1430
xe=xk
goto1460
1430 slope1=abs((u(2)-u(1))/(t(2)-t(1)))
slope2=abs((u(3)-u(2))/(t(3)-t(2)))
xe=1.0d0
if(slope1+slope2.ne.0.)xe=1.-abs(slope1-slope2)/(slope1+slope2)
1460 p2=s+xe*(p1-s)
1500 e1=abs(p1-s)
e2=abs(p2-s)
if(e1+e2.gt.0.)goto1560
1530 z=s
goto1700
1560 bt=(e1*al)/(e1*al+(1.-al)*e2)
z=bt*p2+(1.-bt)*p1
c
1700 y = z
return
end
c
c----------------------------------------------------------------------
c
subroutine lintrp (itab, xtab, ytab, ntab, x, y, ierr)
c
c Linear interpolator with last index save
c
c Arguments are identical to INTRP, and uses the same index search
c scheme so that values for ITAB should not clash with calls
c to INTRP and LINTRP.
c
double precision xtab(ntab), ytab(ntab), x , y
integer itab, ierr
c
c----- Only 1 pt in table
if (ntab .eq. 1) then
y = ytab (1)
ierr = 0
return
endif
c
c-----locate search index
call srch (itab, x, xtab, ntab, index, ierr)
c
c----- index points just above x
y = ytab(index-1) + (x - xtab(index-1)) *
1 (ytab(index) - ytab(index-1)) / (xtab(index) - xtab(index-1))
c
return
end
|
