Line data Source code
1 0 : subroutine do_RGBI3P_db(MD,NXD,NYD,XD,YD,ZD,NIP,XI,YI, ZI,IER, WK)
2 : !
3 : ! Rectangular-grid bivariate interpolation
4 : ! (a master subroutine of the RGBI3P/RGSF3P_db subroutine package)
5 : !
6 : ! Hiroshi Akima
7 : ! U.S. Department of Commerce, NTIA/ITS
8 : ! Version of 1995/08
9 : !
10 : ! This subroutine performs interpolation of a bivariate function,
11 : ! z(x,y), on a rectangular grid in the x-y plane. It is based on
12 : ! the revised Akima method.
13 : !
14 : ! In this subroutine, the interpolating function is a piecewise
15 : ! function composed of a set of bicubic (bivariate third-degree)
16 : ! polynomials, each applicable to a rectangle of the input grid
17 : ! in the x-y plane. Each polynomial is determined locally.
18 : !
19 : ! This subroutine has the accuracy of a bicubic polynomial, i.e.,
20 : ! it interpolates accurately when all data points lie on a
21 : ! surface of a bicubic polynomial.
22 : !
23 : ! The grid lines can be unevenly spaced.
24 : !
25 : ! The input arguments are
26 : ! MD = mode of computation
27 : ! = 1 for new XD, YD, or ZD data (default)
28 : ! = 2 for old XD, YD, and ZD data,
29 : ! NXD = number of the input-grid data points in the x
30 : ! coordinate (must be 2 or greater),
31 : ! NYD = number of the input-grid data points in the y
32 : ! coordinate (must be 2 or greater),
33 : ! XD = array of dimension NXD containing the x coordinates
34 : ! of the input-grid data points (must be in a
35 : ! monotonic increasing order),
36 : ! YD = array of dimension NYD containing the y coordinates
37 : ! of the input-grid data points (must be in a
38 : ! monotonic increasing order),
39 : ! ZD = two-dimensional array of dimension NXD*NYD
40 : ! containing the z(x,y) values at the input-grid data
41 : ! points,
42 : ! NIP = number of the output points at which interpolation
43 : ! of the z value is desired (must be 1 or greater),
44 : ! XI = array of dimension NIP containing the x coordinates
45 : ! of the output points,
46 : ! YI = array of dimension NIP containing the y coordinates
47 : ! of the output points.
48 : !
49 : ! The output arguments are
50 : ! ZI = array of dimension NIP where the interpolated z
51 : ! values at the output points are to be stored,
52 : ! IER = error flag
53 : ! = 0 for no errors
54 : ! = 1 for NXD = 1 or less
55 : ! = 2 for NYD = 1 or less
56 : ! = 3 for identical XD values or
57 : ! XD values out of sequence
58 : ! = 4 for identical YD values or
59 : ! YD values out of sequence
60 : ! = 5 for NIP = 0 or less.
61 : !
62 : ! The other argument is
63 : ! WK = three dimensional array of dimension 3*NXD*NYD used
64 : ! internally as a work area.
65 : !
66 : ! The very first call to this subroutine and the call with a new
67 : ! XD, YD, and ZD array must be made with MD=1. The call with MD=2
68 : ! must be preceded by another call with the same XD, YD, and ZD
69 : ! arrays. Between the call with MD=2 and its preceding call, the
70 : ! WK array must not be disturbed.
71 : !
72 : ! The constant in the parameter statement below is
73 : ! NIPIMX = maximum number of output points to be processed
74 : ! at a time.
75 : ! The constant value has been selected empirically.
76 : !
77 : ! This subroutine calls the RGPD3P_db, RGLCTN_db, and RGPLNL_db subroutines.
78 : !
79 :
80 : ! Specification statements
81 : ! .. Parameters ..
82 : integer :: NIPIMX
83 : parameter (NIPIMX=51)
84 :
85 : ! .. Scalar Arguments ..
86 : integer :: IER,MD,NIP,NXD,NYD
87 :
88 : ! .. Array Arguments ..
89 : double precision :: WK(3,NXD,NYD),XD(NXD),XI(NIP),YD(NYD),YI(NIP),ZD(NXD,NYD),ZI(NIP)
90 :
91 : ! .. Local Scalars ..
92 : integer :: IIP,IX,IY,NIPI
93 :
94 : ! .. Local Arrays ..
95 : integer :: INXI(NIPIMX),INYI(NIPIMX)
96 :
97 : ! .. External Subroutines ..
98 : EXTERNAL RGLCTN_db,RGPD3P_db,RGPLNL_db
99 :
100 : ! .. Intrinsic Functions ..
101 : INTRINSIC MIN
102 :
103 : ! Preliminary processing
104 : ! Error check
105 0 : if (NXD <= 1) GOTO 40
106 0 : if (NYD <= 1) GOTO 50
107 0 : do IX = 2,NXD
108 0 : if (XD(IX) <= XD(IX-1)) GOTO 60
109 : end do
110 0 : do IY = 2,NYD
111 0 : if (YD(IY) <= YD(IY-1)) GOTO 70
112 : end do
113 0 : if (NIP <= 0) GOTO 80
114 0 : IER = 0
115 : ! Calculation
116 : ! Estimates partial derivatives at all input-grid data points
117 : ! (for MD=1).
118 0 : if (MD /= 2) then
119 0 : call RGPD3P_db(NXD,NYD,XD,YD,ZD, WK)
120 : ! call RGPD3P_db(NXD,NYD,XD,YD,ZD, PDD)
121 : end if
122 : ! DO-loop with respect to the output point
123 : ! Processes NIPIMX output points, at most, at a time.
124 0 : do IIP = 1,NIP,NIPIMX
125 0 : NIPI = MIN(NIP- (IIP-1),NIPIMX)
126 : ! Locates the output points.
127 0 : call RGLCTN_db(NXD,NYD,XD,YD,NIPI,XI(IIP),YI(IIP), INXI,INYI)
128 : ! call RGLCTN_db(NXD,NYD,XD,YD,NIP,XI,YI, INXI,INYI)
129 : ! Calculates the z values at the output points.
130 0 : call RGPLNL_db(NXD,NYD,XD,YD,ZD,WK,NIPI,XI(IIP),YI(IIP),INXI,INYI, ZI(IIP))
131 : ! call RGPLNL_db(NXD,NYD,XD,YD,ZD,PDD,NIP,XI,YI,INXI,INYI, ZI)
132 : end do
133 0 : return
134 : ! Error exit
135 0 : 40 WRITE (*,FMT=9000)
136 0 : IER = 1
137 0 : GOTO 90
138 0 : 50 WRITE (*,FMT=9010)
139 0 : IER = 2
140 0 : GOTO 90
141 0 : 60 WRITE (*,FMT=9020) IX,XD(IX)
142 0 : IER = 3
143 0 : GOTO 90
144 0 : 70 WRITE (*,FMT=9030) IY,YD(IY)
145 0 : IER = 4
146 0 : GOTO 90
147 0 : 80 WRITE (*,FMT=9040)
148 0 : IER = 5
149 0 : 90 WRITE (*,FMT=9050) NXD,NYD,NIP
150 0 : return
151 : ! Format statements for error messages
152 : 9000 FORMAT (1X,/,'*** RGBI3P Error 1: NXD = 1 or less')
153 : 9010 FORMAT (1X,/,'*** RGBI3P Error 2: NYD = 1 or less')
154 : 9020 FORMAT (1X,/,'*** RGBI3P Error 3: Identical XD values or',
155 : & ' XD values out of sequence',/,' IX =',I6,', XD(IX) =',
156 : & E11.3)
157 : 9030 FORMAT (1X,/,'*** RGBI3P Error 4: Identical YD values or',
158 : & ' YD values out of sequence',/,' IY =',I6,', YD(IY) =',
159 : & E11.3)
160 : 9040 FORMAT (1X,/,'*** RGBI3P Error 5: NIP = 0 or less')
161 : 9050 FORMAT (' NXD =',I5,', NYD =',I5,', NIP =',I5,/)
162 : end subroutine do_RGBI3P_db
163 :
164 :
165 0 : subroutine do_RGSF3P_db(MD,NXD,NYD,XD,YD,ZD,NXI,XI,NYI,YI, ZI,IER, WK)
166 : !
167 : ! Rectangular-grid surface fitting
168 : ! (a master subroutine of the RGBI3P/RGSF3P_db subroutine package)
169 : !
170 : ! Hiroshi Akima
171 : ! U.S. Department of Commerce, NTIA/ITS
172 : ! Version of 1995/08
173 : !
174 : ! This subroutine performs surface fitting by interpolating
175 : ! values of a bivariate function, z(x,y), on a rectangular grid
176 : ! in the x-y plane. It is based on the revised Akima method.
177 : !
178 : ! In this subroutine, the interpolating function is a piecewise
179 : ! function composed of a set of bicubic (bivariate third-degree)
180 : ! polynomials, each applicable to a rectangle of the input grid
181 : ! in the x-y plane. Each polynomial is determined locally.
182 : !
183 : ! This subroutine has the accuracy of a bicubic polynomial, i.e.,
184 : ! it fits the surface accurately when all data points lie on a
185 : ! surface of a bicubic polynomial.
186 : !
187 : ! The grid lines of the input and output data can be unevenly
188 : ! spaced.
189 : !
190 : ! The input arguments are
191 : ! MD = mode of computation
192 : ! = 1 for new XD, YD, or ZD data (default)
193 : ! = 2 for old XD, YD, and ZD data,
194 : ! NXD = number of the input-grid data points in the x
195 : ! coordinate (must be 2 or greater),
196 : ! NYD = number of the input-grid data points in the y
197 : ! coordinate (must be 2 or greater),
198 : ! XD = array of dimension NXD containing the x coordinates
199 : ! of the input-grid data points (must be in a
200 : ! monotonic increasing order),
201 : ! YD = array of dimension NYD containing the y coordinates
202 : ! of the input-grid data points (must be in a
203 : ! monotonic increasing order),
204 : ! ZD = two-dimensional array of dimension NXD*NYD
205 : ! containing the z(x,y) values at the input-grid data
206 : ! points,
207 : ! NXI = number of output grid points in the x coordinate
208 : ! (must be 1 or greater),
209 : ! XI = array of dimension NXI containing the x coordinates
210 : ! of the output grid points,
211 : ! NYI = number of output grid points in the y coordinate
212 : ! (must be 1 or greater),
213 : ! YI = array of dimension NYI containing the y coordinates
214 : ! of the output grid points.
215 : !
216 : ! The output arguments are
217 : ! ZI = two-dimensional array of dimension NXI*NYI, where
218 : ! the interpolated z values at the output grid
219 : ! points are to be stored,
220 : ! IER = error flag
221 : ! = 0 for no error
222 : ! = 1 for NXD = 1 or less
223 : ! = 2 for NYD = 1 or less
224 : ! = 3 for identical XD values or
225 : ! XD values out of sequence
226 : ! = 4 for identical YD values or
227 : ! YD values out of sequence
228 : ! = 5 for NXI = 0 or less
229 : ! = 6 for NYI = 0 or less.
230 : !
231 : ! The other argument is
232 : ! WK = three-dimensional array of dimension 3*NXD*NYD used
233 : ! internally as a work area.
234 : !
235 : ! The very first call to this subroutine and the call with a new
236 : ! XD, YD, or ZD array must be made with MD=1. The call with MD=2
237 : ! must be preceded by another call with the same XD, YD, and ZD
238 : ! arrays. Between the call with MD=2 and its preceding call, the
239 : ! WK array must not be disturbed.
240 : !
241 : ! The constant in the parameter statement below is
242 : ! NIPIMX = maximum number of output points to be processed
243 : ! at a time.
244 : ! The constant value has been selected empirically.
245 : !
246 : ! This subroutine calls the RGPD3P_db, RGLCTN_db, and RGPLNL_db subroutines.
247 : !
248 :
249 : ! Specification statements
250 : ! .. Parameters ..
251 : integer :: NIPIMX
252 : parameter (NIPIMX=51)
253 :
254 : ! .. Scalar Arguments ..
255 : integer :: IER,MD,NXD,NXI,NYD,NYI
256 :
257 : ! .. Array Arguments ..
258 : double precision :: WK(3,NXD,NYD),XD(NXD),XI(NXI),YD(NYD),YI(NYI),ZD(NXD,NYD),ZI(NXI,NYI)
259 :
260 : ! .. Local Scalars ..
261 : integer :: IX,IXI,IY,IYI,NIPI
262 :
263 : ! .. Local Arrays ..
264 0 : double precision :: YII(NIPIMX)
265 : integer :: INXI(NIPIMX),INYI(NIPIMX)
266 :
267 : ! .. External Subroutines ..
268 : EXTERNAL RGLCTN_db,RGPD3P_db,RGPLNL_db
269 :
270 : ! .. Intrinsic Functions ..
271 : INTRINSIC MIN
272 :
273 : ! Preliminary processing
274 : ! Error check
275 0 : if (NXD <= 1) GOTO 60
276 0 : if (NYD <= 1) GOTO 70
277 0 : do IX = 2,NXD
278 0 : if (XD(IX) <= XD(IX-1)) GOTO 80
279 : end do
280 0 : do IY = 2,NYD
281 0 : if (YD(IY) <= YD(IY-1)) GOTO 90
282 : end do
283 0 : if (NXI <= 0) GOTO 100
284 0 : if (NYI <= 0) GOTO 110
285 0 : IER = 0
286 : ! Calculation
287 : ! Estimates partial derivatives at all input-grid data points
288 : ! (for MD=1).
289 0 : if (MD /= 2) then
290 0 : call RGPD3P_db(NXD,NYD,XD,YD,ZD, WK)
291 : ! call RGPD3P_db(NXD,NYD,XD,YD,ZD, PDD)
292 : end if
293 : ! Outermost DO-loop with respect to the y coordinate of the output
294 : ! grid points
295 0 : do IYI = 1,NYI
296 0 : do IXI = 1,NIPIMX
297 0 : YII(IXI) = YI(IYI)
298 : end do
299 : ! Second DO-loop with respect to the x coordinate of the output
300 : ! grid points
301 : ! Processes NIPIMX output-grid points, at most, at a time.
302 0 : do IXI = 1,NXI,NIPIMX
303 0 : NIPI = MIN(NXI- (IXI-1),NIPIMX)
304 : ! Locates the output-grid points.
305 0 : call RGLCTN_db(NXD,NYD,XD,YD,NIPI,XI(IXI),YII, INXI,INYI)
306 : ! call RGLCTN_db(NXD,NYD,XD,YD,NIP,XI,YI, INXI,INYI)
307 : ! Calculates the z values at the output-grid points.
308 0 : call RGPLNL_db(NXD,NYD,XD,YD,ZD,WK,NIPI,XI(IXI),YII,INXI,INYI, ZI(IXI,IYI))
309 : ! call RGPLNL_db(NXD,NYD,XD,YD,ZD,PDD,NIP,XI,YI,INXI,INYI, ZI)
310 : end do
311 : end do
312 0 : return
313 : ! Error exit
314 0 : 60 WRITE (*,FMT=9000)
315 0 : IER = 1
316 0 : GOTO 120
317 0 : 70 WRITE (*,FMT=9010)
318 0 : IER = 2
319 0 : GOTO 120
320 0 : 80 WRITE (*,FMT=9020) IX,XD(IX)
321 0 : IER = 3
322 0 : GOTO 120
323 0 : 90 WRITE (*,FMT=9030) IY,YD(IY)
324 0 : IER = 4
325 0 : GOTO 120
326 0 : 100 WRITE (*,FMT=9040)
327 0 : IER = 5
328 0 : GOTO 120
329 0 : 110 WRITE (*,FMT=9050)
330 0 : IER = 6
331 0 : 120 WRITE (*,FMT=9060) NXD,NYD,NXI,NYI
332 0 : return
333 : ! Format statements for error messages
334 : 9000 FORMAT (1X,/,'*** RGSF3P_db Error 1: NXD = 1 or less')
335 : 9010 FORMAT (1X,/,'*** RGSF3P_db Error 2: NYD = 1 or less')
336 : 9020 FORMAT (1X,/,'*** RGSF3P_db Error 3: Identical XD values or',
337 : & ' XD values out of sequence',/,' IX =',I6,', XD(IX) =',E11.3)
338 : 9030 FORMAT (1X,/,'*** RGSF3P_db Error 4: Identical YD values or',
339 : & ' YD values out of sequence',/,' IY =',I6,', YD(IY) =',E11.3)
340 : 9040 FORMAT (1X,/,'*** RGSF3P_db Error 5: NXI = 0 or less')
341 : 9050 FORMAT (1X,/,'*** RGSF3P_db Error 6: NYI = 0 or less')
342 : 9060 FORMAT (' NXD =',I5,', NYD =',I5,', NXI =',I5,', NYI =',I5,/)
343 : end subroutine do_RGSF3P_db
344 :
345 :
346 0 : subroutine RGPD3P_db(NXD,NYD,XD,YD,ZD, PDD)
347 : !
348 : ! Partial derivatives of a bivariate function on a rectangular grid
349 : ! (a supporting subroutine of the RGBI3P/RGSF3P_db subroutine package)
350 : !
351 : ! Hiroshi Akima
352 : ! U.S. Department of Commerce, NTIA/ITS
353 : ! Version of 1995/08
354 : !
355 : ! This subroutine estimates three partial derivatives, zx, zy, and
356 : ! zxy, of a bivariate function, z(x,y), on a rectangular grid in
357 : ! the x-y plane. It is based on the revised Akima method that has
358 : ! the accuracy of a bicubic polynomial.
359 : !
360 : ! The input arguments are
361 : ! NXD = number of the input-grid data points in the x
362 : ! coordinate (must be 2 or greater),
363 : ! NYD = number of the input-grid data points in the y
364 : ! coordinate (must be 2 or greater),
365 : ! XD = array of dimension NXD containing the x coordinates
366 : ! of the input-grid data points (must be in a
367 : ! monotonic increasing order),
368 : ! YD = array of dimension NYD containing the y coordinates
369 : ! of the input-grid data points (must be in a
370 : ! monotonic increasing order),
371 : ! ZD = two-dimensional array of dimension NXD*NYD
372 : ! containing the z(x,y) values at the input-grid data
373 : ! points.
374 : !
375 : ! The output argument is
376 : ! PDD = three-dimensional array of dimension 3*NXD*NYD,
377 : ! where the estimated zx, zy, and zxy values at the
378 : ! input-grid data points are to be stored.
379 : !
380 :
381 : ! Specification statements
382 : ! .. Scalar Arguments ..
383 : integer :: NXD,NYD
384 :
385 : ! .. Array Arguments ..
386 : double precision :: PDD(3,NXD,NYD),XD(NXD),YD(NYD),ZD(NXD,NYD)
387 :
388 : ! .. Local Scalars ..
389 0 : double precision :: B00,B00X,B00Y,B01,B10,B11,CX1,CX2,CX3,CY1,CY2,
390 0 : & CY3,DISF,DNM,DZ00,DZ01,DZ02,DZ03,DZ10,DZ11,DZ12,
391 0 : & DZ13,DZ20,DZ21,DZ22,DZ23,DZ30,DZ31,DZ32,DZ33,
392 0 : & DZX10,DZX20,DZX30,DZXY11,DZXY12,DZXY13,DZXY21,
393 0 : & DZXY22,DZXY23,DZXY31,DZXY32,DZXY33,DZY01,DZY02,
394 0 : & DZY03,EPSLN,PEZX,PEZXY,PEZY,SMPEF,SMPEI,SMWTF,
395 0 : & SMWTI,SX,SXX,SXXY,SXXYY,SXY,SXYY,SXYZ,SXZ,SY,SYY,
396 0 : & SYZ,SZ,VOLF,WT,X0,X1,X2,X3,XX1,XX2,XX3,Y0,Y1,Y2,
397 0 : & Y3,Z00,Z01,Z02,Z03,Z10,Z11,Z12,Z13,Z20,Z21,Z22,
398 0 : & Z23,Z30,Z31,Z32,Z33,ZXDI,ZXYDI,ZYDI,ZZ0,ZZ1,ZZ2
399 : integer :: IPEX,IPEY,IX0,IX1,IX2,IX3,IY0,IY1,IY2,IY3,JPEXY,JXY,NX0,NY0
400 :
401 : ! .. Local Arrays ..
402 0 : double precision :: B00XA(4),B00YA(4),B01A(4),B10A(4),CXA(3,4),
403 0 : & CYA(3,4),SXA(4),SXXA(4),SYA(4),SYYA(4),XA(3,4),
404 0 : & YA(3,4),Z0IA(3,4),ZI0A(3,4)
405 : integer :: IDLT(3,4)
406 :
407 : ! .. Intrinsic Functions ..
408 : INTRINSIC MAX
409 :
410 : ! .. Statement Functions ..
411 : double precision :: Z2F,Z3F
412 :
413 : ! Data statements
414 : DATA ((IDLT(JXY,JPEXY),JPEXY=1,4),JXY=1,3)/-3,-2,-1,1,-2,-1,1,2,-1,1,2,3/
415 :
416 : ! Statement Function definitions
417 : Z2F(XX1,XX2,ZZ0,ZZ1) = (ZZ1-ZZ0)*XX2/XX1 + ZZ0
418 : Z3F(XX1,XX2,XX3,ZZ0,ZZ1,ZZ2) = ((ZZ2-ZZ0)* (XX3-XX1)/XX2-
419 : & (ZZ1-ZZ0)* (XX3-XX2)/XX1)*
420 : & (XX3/ (XX2-XX1)) + ZZ0
421 :
422 : ! Calculation
423 : ! Initial setting of some local variables
424 0 : NX0 = MAX(4,NXD)
425 0 : NY0 = MAX(4,NYD)
426 : ! Double DO-loop with respect to the input grid points
427 0 : do 60 IY0 = 1,NYD
428 0 : do 50 IX0 = 1,NXD
429 0 : X0 = XD(IX0)
430 0 : Y0 = YD(IY0)
431 0 : Z00 = ZD(IX0,IY0)
432 : ! Part 1. Estimation of ZXDI
433 : ! Initial setting
434 0 : SMPEF = 0.0
435 0 : SMWTF = 0.0
436 0 : SMPEI = 0.0
437 0 : SMWTI = 0.0
438 : ! DO-loop with respect to the primary estimate
439 0 : do 10 IPEX = 1,4
440 : ! Selects necessary grid points in the x direction.
441 0 : IX1 = IX0 + IDLT(1,IPEX)
442 0 : IX2 = IX0 + IDLT(2,IPEX)
443 0 : IX3 = IX0 + IDLT(3,IPEX)
444 : if ((IX1 < 1) .or. (IX2 < 1) .or. (IX3 < 1) .or.
445 0 : & (IX1 > NX0) .or. (IX2 > NX0) .or.
446 : & (IX3 > NX0)) GOTO 10
447 : ! Selects and/or supplements the x and z values.
448 0 : X1 = XD(IX1) - X0
449 0 : Z10 = ZD(IX1,IY0)
450 0 : if (NXD >= 4) then
451 0 : X2 = XD(IX2) - X0
452 0 : X3 = XD(IX3) - X0
453 0 : Z20 = ZD(IX2,IY0)
454 0 : Z30 = ZD(IX3,IY0)
455 0 : else if (NXD == 3) then
456 0 : X2 = XD(IX2) - X0
457 0 : Z20 = ZD(IX2,IY0)
458 0 : X3 = 2*XD(3) - XD(2) - X0
459 0 : Z30 = Z3F(X1,X2,X3,Z00,Z10,Z20)
460 0 : else if (NXD == 2) then
461 0 : X2 = 2*XD(2) - XD(1) - X0
462 0 : Z20 = Z2F(X1,X2,Z00,Z10)
463 0 : X3 = 2*XD(1) - XD(2) - X0
464 0 : Z30 = Z2F(X1,X3,Z00,Z10)
465 : end if
466 0 : DZX10 = (Z10-Z00)/X1
467 0 : DZX20 = (Z20-Z00)/X2
468 0 : DZX30 = (Z30-Z00)/X3
469 : ! Calculates the primary estimate of partial derivative zx as
470 : ! the coefficient of the bicubic polynomial.
471 0 : CX1 = X2*X3/ ((X1-X2)* (X1-X3))
472 0 : CX2 = X3*X1/ ((X2-X3)* (X2-X1))
473 0 : CX3 = X1*X2/ ((X3-X1)* (X3-X2))
474 0 : PEZX = CX1*DZX10 + CX2*DZX20 + CX3*DZX30
475 : ! Calculates the volatility factor and distance factor in the x
476 : ! direction for the primary estimate of zx.
477 0 : SX = X1 + X2 + X3
478 0 : SZ = Z00 + Z10 + Z20 + Z30
479 0 : SXX = X1*X1 + X2*X2 + X3*X3
480 0 : SXZ = X1*Z10 + X2*Z20 + X3*Z30
481 0 : DNM = 4.0*SXX - SX*SX
482 0 : B00 = (SXX*SZ-SX*SXZ)/DNM
483 0 : B10 = (4.0*SXZ-SX*SZ)/DNM
484 0 : DZ00 = Z00 - B00
485 0 : DZ10 = Z10 - (B00+B10*X1)
486 0 : DZ20 = Z20 - (B00+B10*X2)
487 0 : DZ30 = Z30 - (B00+B10*X3)
488 0 : VOLF = DZ00**2 + DZ10**2 + DZ20**2 + DZ30**2
489 0 : DISF = SXX
490 : ! Calculates the EPSLN value, which is used to decide whether or
491 : ! not the volatility factor is essentially zero.
492 0 : EPSLN = (Z00**2+Z10**2+Z20**2+Z30**2)*1.0E-12
493 : ! Accumulates the weighted primary estimates of zx and their
494 : ! weights.
495 0 : if (VOLF > EPSLN) then
496 : ! - For a finite weight.
497 0 : WT = 1.0/ (VOLF*DISF)
498 0 : SMPEF = SMPEF + WT*PEZX
499 0 : SMWTF = SMWTF + WT
500 : else
501 : ! - For an infinite weight.
502 0 : SMPEI = SMPEI + PEZX
503 0 : SMWTI = SMWTI + 1.0
504 : end if
505 : ! Saves the necessary values for estimating zxy
506 0 : XA(1,IPEX) = X1
507 0 : XA(2,IPEX) = X2
508 0 : XA(3,IPEX) = X3
509 0 : ZI0A(1,IPEX) = Z10
510 0 : ZI0A(2,IPEX) = Z20
511 0 : ZI0A(3,IPEX) = Z30
512 0 : CXA(1,IPEX) = CX1
513 0 : CXA(2,IPEX) = CX2
514 0 : CXA(3,IPEX) = CX3
515 0 : SXA(IPEX) = SX
516 0 : SXXA(IPEX) = SXX
517 0 : B00XA(IPEX) = B00
518 0 : B10A(IPEX) = B10
519 0 : 10 continue
520 : ! Calculates the final estimate of zx.
521 0 : if (SMWTI < 0.5) then
522 : ! - When no infinite weights exist.
523 0 : ZXDI = SMPEF/SMWTF
524 : else
525 : ! - When infinite weights exist.
526 0 : ZXDI = SMPEI/SMWTI
527 : end if
528 : ! End of Part 1.
529 : ! Part 2. Estimation of ZYDI
530 : ! Initial setting
531 0 : SMPEF = 0.0
532 0 : SMWTF = 0.0
533 0 : SMPEI = 0.0
534 0 : SMWTI = 0.0
535 : ! DO-loop with respect to the primary estimate
536 0 : do 20 IPEY = 1,4
537 : ! Selects necessary grid points in the y direction.
538 0 : IY1 = IY0 + IDLT(1,IPEY)
539 0 : IY2 = IY0 + IDLT(2,IPEY)
540 0 : IY3 = IY0 + IDLT(3,IPEY)
541 : if ((IY1 < 1) .or. (IY2 < 1) .or. (IY3 < 1) .or.
542 0 : & (IY1 > NY0) .or. (IY2 > NY0) .or.
543 : & (IY3 > NY0)) GOTO 20
544 : ! Selects and/or supplements the y and z values.
545 0 : Y1 = YD(IY1) - Y0
546 0 : Z01 = ZD(IX0,IY1)
547 0 : if (NYD >= 4) then
548 0 : Y2 = YD(IY2) - Y0
549 0 : Y3 = YD(IY3) - Y0
550 0 : Z02 = ZD(IX0,IY2)
551 0 : Z03 = ZD(IX0,IY3)
552 0 : else if (NYD == 3) then
553 0 : Y2 = YD(IY2) - Y0
554 0 : Z02 = ZD(IX0,IY2)
555 0 : Y3 = 2*YD(3) - YD(2) - Y0
556 0 : Z03 = Z3F(Y1,Y2,Y3,Z00,Z01,Z02)
557 0 : else if (NYD == 2) then
558 0 : Y2 = 2*YD(2) - YD(1) - Y0
559 0 : Z02 = Z2F(Y1,Y2,Z00,Z01)
560 0 : Y3 = 2*YD(1) - YD(2) - Y0
561 0 : Z03 = Z2F(Y1,Y3,Z00,Z01)
562 : end if
563 0 : DZY01 = (Z01-Z00)/Y1
564 0 : DZY02 = (Z02-Z00)/Y2
565 0 : DZY03 = (Z03-Z00)/Y3
566 : ! Calculates the primary estimate of partial derivative zy as
567 : ! the coefficient of the bicubic polynomial.
568 0 : CY1 = Y2*Y3/ ((Y1-Y2)* (Y1-Y3))
569 0 : CY2 = Y3*Y1/ ((Y2-Y3)* (Y2-Y1))
570 0 : CY3 = Y1*Y2/ ((Y3-Y1)* (Y3-Y2))
571 0 : PEZY = CY1*DZY01 + CY2*DZY02 + CY3*DZY03
572 : ! Calculates the volatility factor and distance factor in the y
573 : ! direction for the primary estimate of zy.
574 0 : SY = Y1 + Y2 + Y3
575 0 : SZ = Z00 + Z01 + Z02 + Z03
576 0 : SYY = Y1*Y1 + Y2*Y2 + Y3*Y3
577 0 : SYZ = Y1*Z01 + Y2*Z02 + Y3*Z03
578 0 : DNM = 4.0*SYY - SY*SY
579 0 : B00 = (SYY*SZ-SY*SYZ)/DNM
580 0 : B01 = (4.0*SYZ-SY*SZ)/DNM
581 0 : DZ00 = Z00 - B00
582 0 : DZ01 = Z01 - (B00+B01*Y1)
583 0 : DZ02 = Z02 - (B00+B01*Y2)
584 0 : DZ03 = Z03 - (B00+B01*Y3)
585 0 : VOLF = DZ00**2 + DZ01**2 + DZ02**2 + DZ03**2
586 0 : DISF = SYY
587 : ! Calculates the EPSLN value, which is used to decide whether or
588 : ! not the volatility factor is essentially zero.
589 0 : EPSLN = (Z00**2+Z01**2+Z02**2+Z03**2)*1.0E-12
590 : ! Accumulates the weighted primary estimates of zy and their
591 : ! weights.
592 0 : if (VOLF > EPSLN) then
593 : ! - For a finite weight.
594 0 : WT = 1.0/ (VOLF*DISF)
595 0 : SMPEF = SMPEF + WT*PEZY
596 0 : SMWTF = SMWTF + WT
597 : else
598 : ! - For an infinite weight.
599 0 : SMPEI = SMPEI + PEZY
600 0 : SMWTI = SMWTI + 1.0
601 : end if
602 : ! Saves the necessary values for estimating zxy
603 0 : YA(1,IPEY) = Y1
604 0 : YA(2,IPEY) = Y2
605 0 : YA(3,IPEY) = Y3
606 0 : Z0IA(1,IPEY) = Z01
607 0 : Z0IA(2,IPEY) = Z02
608 0 : Z0IA(3,IPEY) = Z03
609 0 : CYA(1,IPEY) = CY1
610 0 : CYA(2,IPEY) = CY2
611 0 : CYA(3,IPEY) = CY3
612 0 : SYA(IPEY) = SY
613 0 : SYYA(IPEY) = SYY
614 0 : B00YA(IPEY) = B00
615 0 : B01A(IPEY) = B01
616 0 : 20 continue
617 : ! Calculates the final estimate of zy.
618 0 : if (SMWTI < 0.5) then
619 : ! - When no infinite weights exist.
620 0 : ZYDI = SMPEF/SMWTF
621 : else
622 : ! - When infinite weights exist.
623 0 : ZYDI = SMPEI/SMWTI
624 : end if
625 : ! End of Part 2.
626 : ! Part 3. Estimation of ZXYDI
627 : ! Initial setting
628 0 : SMPEF = 0.0
629 0 : SMWTF = 0.0
630 0 : SMPEI = 0.0
631 0 : SMWTI = 0.0
632 : ! Outer DO-loops with respect to the primary estimates in the x
633 : ! direction
634 0 : do 40 IPEX = 1,4
635 0 : IX1 = IX0 + IDLT(1,IPEX)
636 0 : IX2 = IX0 + IDLT(2,IPEX)
637 0 : IX3 = IX0 + IDLT(3,IPEX)
638 : if ((IX1 < 1) .or. (IX2 < 1) .or. (IX3 < 1) .or.
639 0 : & (IX1 > NX0) .or. (IX2 > NX0) .or.
640 : & (IX3 > NX0)) GOTO 40
641 : ! Retrieves the necessary values for estimating zxy in the x
642 : ! direction.
643 0 : X1 = XA(1,IPEX)
644 0 : X2 = XA(2,IPEX)
645 0 : X3 = XA(3,IPEX)
646 0 : Z10 = ZI0A(1,IPEX)
647 0 : Z20 = ZI0A(2,IPEX)
648 0 : Z30 = ZI0A(3,IPEX)
649 0 : CX1 = CXA(1,IPEX)
650 0 : CX2 = CXA(2,IPEX)
651 0 : CX3 = CXA(3,IPEX)
652 0 : SX = SXA(IPEX)
653 0 : SXX = SXXA(IPEX)
654 0 : B00X = B00XA(IPEX)
655 0 : B10 = B10A(IPEX)
656 : ! Inner DO-loops with respect to the primary estimates in the y
657 : ! direction
658 0 : do 30 IPEY = 1,4
659 0 : IY1 = IY0 + IDLT(1,IPEY)
660 0 : IY2 = IY0 + IDLT(2,IPEY)
661 0 : IY3 = IY0 + IDLT(3,IPEY)
662 : if ((IY1 < 1) .or. (IY2 < 1) .or.
663 : & (IY3 < 1) .or. (IY1 > NY0) .or.
664 0 : & (IY2 > NY0) .or. (IY3 > NY0)) GOTO 30
665 : ! Retrieves the necessary values for estimating zxy in the y
666 : ! direction.
667 0 : Y1 = YA(1,IPEY)
668 0 : Y2 = YA(2,IPEY)
669 0 : Y3 = YA(3,IPEY)
670 0 : Z01 = Z0IA(1,IPEY)
671 0 : Z02 = Z0IA(2,IPEY)
672 0 : Z03 = Z0IA(3,IPEY)
673 0 : CY1 = CYA(1,IPEY)
674 0 : CY2 = CYA(2,IPEY)
675 0 : CY3 = CYA(3,IPEY)
676 0 : SY = SYA(IPEY)
677 0 : SYY = SYYA(IPEY)
678 0 : B00Y = B00YA(IPEY)
679 0 : B01 = B01A(IPEY)
680 : ! Selects and/or supplements the z values.
681 0 : if (NYD >= 4) then
682 0 : Z11 = ZD(IX1,IY1)
683 0 : Z12 = ZD(IX1,IY2)
684 0 : Z13 = ZD(IX1,IY3)
685 0 : if (NXD >= 4) then
686 0 : Z21 = ZD(IX2,IY1)
687 0 : Z22 = ZD(IX2,IY2)
688 0 : Z23 = ZD(IX2,IY3)
689 0 : Z31 = ZD(IX3,IY1)
690 0 : Z32 = ZD(IX3,IY2)
691 0 : Z33 = ZD(IX3,IY3)
692 0 : else if (NXD == 3) then
693 0 : Z21 = ZD(IX2,IY1)
694 0 : Z22 = ZD(IX2,IY2)
695 0 : Z23 = ZD(IX2,IY3)
696 0 : Z31 = Z3F(X1,X2,X3,Z01,Z11,Z21)
697 0 : Z32 = Z3F(X1,X2,X3,Z02,Z12,Z22)
698 0 : Z33 = Z3F(X1,X2,X3,Z03,Z13,Z23)
699 0 : else if (NXD == 2) then
700 0 : Z21 = Z2F(X1,X2,Z01,Z11)
701 0 : Z22 = Z2F(X1,X2,Z02,Z12)
702 0 : Z23 = Z2F(X1,X2,Z03,Z13)
703 0 : Z31 = Z2F(X1,X3,Z01,Z11)
704 0 : Z32 = Z2F(X1,X3,Z02,Z12)
705 0 : Z33 = Z2F(X1,X3,Z03,Z13)
706 : end if
707 0 : else if (NYD == 3) then
708 0 : Z11 = ZD(IX1,IY1)
709 0 : Z12 = ZD(IX1,IY2)
710 0 : Z13 = Z3F(Y1,Y2,Y3,Z10,Z11,Z12)
711 0 : if (NXD >= 4) then
712 0 : Z21 = ZD(IX2,IY1)
713 0 : Z22 = ZD(IX2,IY2)
714 0 : Z31 = ZD(IX3,IY1)
715 0 : Z32 = ZD(IX3,IY2)
716 0 : else if (NXD == 3) then
717 0 : Z21 = ZD(IX2,IY1)
718 0 : Z22 = ZD(IX2,IY2)
719 0 : Z31 = Z3F(X1,X2,X3,Z01,Z11,Z21)
720 0 : Z32 = Z3F(X1,X2,X3,Z02,Z12,Z22)
721 0 : else if (NXD == 2) then
722 0 : Z21 = Z2F(X1,X2,Z01,Z11)
723 0 : Z22 = Z2F(X1,X2,Z02,Z12)
724 0 : Z31 = Z2F(X1,X3,Z01,Z11)
725 0 : Z32 = Z2F(X1,X3,Z02,Z12)
726 : end if
727 0 : Z23 = Z3F(Y1,Y2,Y3,Z20,Z21,Z22)
728 0 : Z33 = Z3F(Y1,Y2,Y3,Z30,Z31,Z32)
729 0 : else if (NYD == 2) then
730 0 : Z11 = ZD(IX1,IY1)
731 0 : Z12 = Z2F(Y1,Y2,Z10,Z11)
732 0 : Z13 = Z2F(Y1,Y3,Z10,Z11)
733 0 : if (NXD >= 4) then
734 0 : Z21 = ZD(IX2,IY1)
735 0 : Z31 = ZD(IX3,IY1)
736 0 : else if (NXD == 3) then
737 0 : Z21 = ZD(IX2,IY1)
738 0 : Z31 = Z3F(X1,X2,X3,Z01,Z11,Z21)
739 0 : else if (NXD == 2) then
740 0 : Z21 = Z2F(X1,X2,Z01,Z11)
741 0 : Z31 = Z2F(X1,X3,Z01,Z11)
742 : end if
743 0 : Z22 = Z2F(Y1,Y2,Z20,Z21)
744 0 : Z23 = Z2F(Y1,Y3,Z20,Z21)
745 0 : Z32 = Z2F(Y1,Y2,Z30,Z31)
746 0 : Z33 = Z2F(Y1,Y3,Z30,Z31)
747 : end if
748 : ! Calculates the primary estimate of partial derivative zxy as
749 : ! the coefficient of the bicubic polynomial.
750 0 : DZXY11 = (Z11-Z10-Z01+Z00)/ (X1*Y1)
751 0 : DZXY12 = (Z12-Z10-Z02+Z00)/ (X1*Y2)
752 0 : DZXY13 = (Z13-Z10-Z03+Z00)/ (X1*Y3)
753 0 : DZXY21 = (Z21-Z20-Z01+Z00)/ (X2*Y1)
754 0 : DZXY22 = (Z22-Z20-Z02+Z00)/ (X2*Y2)
755 0 : DZXY23 = (Z23-Z20-Z03+Z00)/ (X2*Y3)
756 0 : DZXY31 = (Z31-Z30-Z01+Z00)/ (X3*Y1)
757 0 : DZXY32 = (Z32-Z30-Z02+Z00)/ (X3*Y2)
758 0 : DZXY33 = (Z33-Z30-Z03+Z00)/ (X3*Y3)
759 : PEZXY = CX1* (CY1*DZXY11+CY2*DZXY12+CY3*DZXY13) +
760 : & CX2* (CY1*DZXY21+CY2*DZXY22+CY3*DZXY23) +
761 0 : & CX3* (CY1*DZXY31+CY2*DZXY32+CY3*DZXY33)
762 : ! Calculates the volatility factor and distance factor in the x
763 : ! and y directions for the primary estimate of zxy.
764 0 : B00 = (B00X+B00Y)/2.0
765 0 : SXY = SX*SY
766 0 : SXXY = SXX*SY
767 0 : SXYY = SX*SYY
768 0 : SXXYY = SXX*SYY
769 : SXYZ = X1* (Y1*Z11+Y2*Z12+Y3*Z13) +
770 : & X2* (Y1*Z21+Y2*Z22+Y3*Z23) +
771 0 : & X3* (Y1*Z31+Y2*Z32+Y3*Z33)
772 0 : B11 = (SXYZ-B00*SXY-B10*SXXY-B01*SXYY)/SXXYY
773 0 : DZ00 = Z00 - B00
774 0 : DZ01 = Z01 - (B00+B01*Y1)
775 0 : DZ02 = Z02 - (B00+B01*Y2)
776 0 : DZ03 = Z03 - (B00+B01*Y3)
777 0 : DZ10 = Z10 - (B00+B10*X1)
778 0 : DZ11 = Z11 - (B00+B01*Y1+X1* (B10+B11*Y1))
779 0 : DZ12 = Z12 - (B00+B01*Y2+X1* (B10+B11*Y2))
780 0 : DZ13 = Z13 - (B00+B01*Y3+X1* (B10+B11*Y3))
781 0 : DZ20 = Z20 - (B00+B10*X2)
782 0 : DZ21 = Z21 - (B00+B01*Y1+X2* (B10+B11*Y1))
783 0 : DZ22 = Z22 - (B00+B01*Y2+X2* (B10+B11*Y2))
784 0 : DZ23 = Z23 - (B00+B01*Y3+X2* (B10+B11*Y3))
785 0 : DZ30 = Z30 - (B00+B10*X3)
786 0 : DZ31 = Z31 - (B00+B01*Y1+X3* (B10+B11*Y1))
787 0 : DZ32 = Z32 - (B00+B01*Y2+X3* (B10+B11*Y2))
788 0 : DZ33 = Z33 - (B00+B01*Y3+X3* (B10+B11*Y3))
789 : VOLF = DZ00**2 + DZ01**2 + DZ02**2 + DZ03**2 +
790 : & DZ10**2 + DZ11**2 + DZ12**2 + DZ13**2 +
791 : & DZ20**2 + DZ21**2 + DZ22**2 + DZ23**2 +
792 0 : & DZ30**2 + DZ31**2 + DZ32**2 + DZ33**2
793 0 : DISF = SXX*SYY
794 : ! Calculates EPSLN.
795 : EPSLN = (Z00**2+Z01**2+Z02**2+Z03**2+Z10**2+
796 : & Z11**2+Z12**2+Z13**2+Z20**2+Z21**2+Z22**2+
797 : & Z23**2+Z30**2+Z31**2+Z32**2+Z33**2)*
798 0 : & 1.0E-12
799 : ! Accumulates the weighted primary estimates of zxy and their
800 : ! weights.
801 0 : if (VOLF > EPSLN) then
802 : ! - For a finite weight.
803 0 : WT = 1.0/ (VOLF*DISF)
804 0 : SMPEF = SMPEF + WT*PEZXY
805 0 : SMWTF = SMWTF + WT
806 : else
807 : ! - For an infinite weight.
808 0 : SMPEI = SMPEI + PEZXY
809 0 : SMWTI = SMWTI + 1.0
810 : end if
811 0 : 30 continue
812 0 : 40 continue
813 : ! Calculates the final estimate of zxy.
814 0 : if (SMWTI < 0.5) then
815 : ! - When no infinite weights exist.
816 0 : ZXYDI = SMPEF/SMWTF
817 : else
818 : ! - When infinite weights exist.
819 0 : ZXYDI = SMPEI/SMWTI
820 : end if
821 : ! End of Part 3
822 0 : PDD(1,IX0,IY0) = ZXDI
823 0 : PDD(2,IX0,IY0) = ZYDI
824 0 : PDD(3,IX0,IY0) = ZXYDI
825 0 : 50 continue
826 0 : 60 continue
827 0 : return
828 : end subroutine RGPD3P_db
829 :
830 :
831 0 : subroutine RGLCTN_db(NXD,NYD,XD,YD,NIP,XI,YI, INXI,INYI)
832 : !
833 : ! Location of the desired points in a rectangular grid
834 : ! (a supporting subroutine of the RGBI3P/RGSF3P_db subroutine package)
835 : !
836 : ! Hiroshi Akima
837 : ! U.S. Department of Commerce, NTIA/ITS
838 : ! Version of 1995/08
839 : !
840 : ! This subroutine locates the desired points in a rectangular grid
841 : ! in the x-y plane.
842 : !
843 : ! The grid lines can be unevenly spaced.
844 : !
845 : ! The input arguments are
846 : ! NXD = number of the input-grid data points in the x
847 : ! coordinate (must be 2 or greater),
848 : ! NYD = number of the input-grid data points in the y
849 : ! coordinate (must be 2 or greater),
850 : ! XD = array of dimension NXD containing the x coordinates
851 : ! of the input-grid data points (must be in a
852 : ! monotonic increasing order),
853 : ! YD = array of dimension NYD containing the y coordinates
854 : ! of the input-grid data points (must be in a
855 : ! monotonic increasing order),
856 : ! NIP = number of the output points to be located (must be
857 : ! 1 or greater),
858 : ! XI = array of dimension NIP containing the x coordinates
859 : ! of the output points to be located,
860 : ! YI = array of dimension NIP containing the y coordinates
861 : ! of the output points to be located.
862 : !
863 : ! The output arguments are
864 : ! INXI = integer array of dimension NIP where the interval
865 : ! numbers of the XI array elements are to be stored,
866 : ! INYI = integer array of dimension NIP where the interval
867 : ! numbers of the YI array elements are to be stored.
868 : ! The interval numbers are between 0 and NXD and between 0 and NYD,
869 : ! respectively.
870 : !
871 :
872 : ! Specification statements
873 : ! .. Scalar Arguments ..
874 : integer :: NIP,NXD,NYD
875 :
876 : ! .. Array Arguments ..
877 : double precision :: XD(NXD),XI(NIP),YD(NYD),YI(NIP)
878 : integer :: INXI(NIP),INYI(NIP)
879 :
880 : ! .. Local Scalars ..
881 0 : double precision :: XII,YII
882 : integer :: IIP,IMD,IMN,IMX,IXD,IYD,NINTX,NINTY
883 :
884 : ! DO-loop with respect to IIP, which is the point number of the
885 : ! output point
886 0 : do IIP = 1,NIP
887 0 : XII = XI(IIP)
888 0 : YII = YI(IIP)
889 : ! Checks if the x coordinate of the IIPth output point, XII, is
890 : ! in a new interval. (NINTX is the new-interval flag.)
891 0 : if (IIP == 1) then
892 : NINTX = 1
893 : else
894 0 : NINTX = 0
895 0 : if (IXD == 0) then
896 0 : if (XII > XD(1)) NINTX = 1
897 0 : else if (IXD < NXD) then
898 0 : if ((XII < XD(IXD)) .or. (XII > XD(IXD+1))) NINTX = 1
899 : else
900 0 : if (XII < XD(NXD)) NINTX = 1
901 : end if
902 : end if
903 : ! Locates the output point by binary search if XII is in a new
904 : ! interval. Determines IXD for which XII lies between XD(IXD)
905 : ! and XD(IXD+1).
906 : if (NINTX == 1) then
907 0 : if (XII <= XD(1)) then
908 : IXD = 0
909 0 : else if (XII < XD(NXD)) then
910 0 : IMN = 1
911 0 : IMX = NXD
912 0 : IMD = (IMN+IMX)/2
913 0 : 10 if (XII >= XD(IMD)) then
914 : IMN = IMD
915 : else
916 0 : IMX = IMD
917 : end if
918 0 : IMD = (IMN+IMX)/2
919 0 : if (IMD > IMN) GOTO 10
920 : IXD = IMD
921 : else
922 : IXD = NXD
923 : end if
924 : end if
925 0 : INXI(IIP) = IXD
926 : ! Checks if the y coordinate of the IIPth output point, YII, is
927 : ! in a new interval. (NINTY is the new-interval flag.)
928 0 : if (IIP == 1) then
929 : NINTY = 1
930 : else
931 0 : NINTY = 0
932 0 : if (IYD == 0) then
933 0 : if (YII > YD(1)) NINTY = 1
934 0 : else if (IYD < NYD) then
935 0 : if ((YII < YD(IYD)) .or. (YII > YD(IYD+1))) NINTY = 1
936 : else
937 0 : if (YII < YD(NYD)) NINTY = 1
938 : end if
939 : end if
940 : ! Locates the output point by binary search if YII is in a new
941 : ! interval. Determines IYD for which YII lies between YD(IYD)
942 : ! and YD(IYD+1).
943 : if (NINTY == 1) then
944 0 : if (YII <= YD(1)) then
945 : IYD = 0
946 0 : else if (YII < YD(NYD)) then
947 0 : IMN = 1
948 0 : IMX = NYD
949 0 : IMD = (IMN+IMX)/2
950 0 : 20 if (YII >= YD(IMD)) then
951 : IMN = IMD
952 : else
953 0 : IMX = IMD
954 : end if
955 0 : IMD = (IMN+IMX)/2
956 0 : if (IMD > IMN) GOTO 20
957 : IYD = IMD
958 : else
959 : IYD = NYD
960 : end if
961 : end if
962 0 : INYI(IIP) = IYD
963 : end do
964 0 : return
965 : end subroutine RGLCTN_db
966 :
967 :
968 0 : subroutine RGPLNL_db(NXD,NYD,XD,YD,ZD,PDD,NIP,XI,YI,INXI,INYI, ZI)
969 : !
970 : ! Polynomials for rectangular-grid bivariate interpolation and
971 : ! surface fitting
972 : ! (a supporting subroutine of the RGBI3P/RGSF3P_db subroutine package)
973 : !
974 : ! Hiroshi Akima
975 : ! U.S. Department of Commerce, NTIA/ITS
976 : ! Version of 1995/08
977 : !
978 : ! This subroutine determines a polynomial in x and y for a rectangle
979 : ! of the input grid in the x-y plane and calculates the z value for
980 : ! the desired points by evaluating the polynomial for rectangular-
981 : ! grid bivariate interpolation and surface fitting.
982 : !
983 : ! The input arguments are
984 : ! NXD = number of the input-grid data points in the x
985 : ! coordinate (must be 2 or greater),
986 : ! NYD = number of the input-grid data points in the y
987 : ! coordinate (must be 2 or greater),
988 : ! XD = array of dimension NXD containing the x coordinates
989 : ! of the input-grid data points (must be in a
990 : ! monotonic increasing order),
991 : ! YD = array of dimension NYD containing the y coordinates
992 : ! of the input-grid data points (must be in a
993 : ! monotonic increasing order),
994 : ! ZD = two-dimensional array of dimension NXD*NYD
995 : ! containing the z(x,y) values at the input-grid data
996 : ! points,
997 : ! PDD = three-dimensional array of dimension 3*NXD*NYD
998 : ! containing the estimated zx, zy, and zxy values
999 : ! at the input-grid data points,
1000 : ! NIP = number of the output points at which interpolation
1001 : ! is to be performed,
1002 : ! XI = array of dimension NIP containing the x coordinates
1003 : ! of the output points,
1004 : ! YI = array of dimension NIP containing the y coordinates
1005 : ! of the output points,
1006 : ! INXI = integer array of dimension NIP containing the
1007 : ! interval numbers of the input grid intervals in the
1008 : ! x direction where the x coordinates of the output
1009 : ! points lie,
1010 : ! INYI = integer array of dimension NIP containing the
1011 : ! interval numbers of the input grid intervals in the
1012 : ! y direction where the y coordinates of the output
1013 : ! points lie.
1014 : !
1015 : ! The output argument is
1016 : ! ZI = array of dimension NIP, where the interpolated z
1017 : ! values at the output points are to be stored.
1018 : !
1019 :
1020 : ! Specification statements
1021 : ! .. Scalar Arguments ..
1022 : integer :: NIP,NXD,NYD
1023 :
1024 : ! .. Array Arguments ..
1025 : double precision :: PDD(3,NXD,NYD),XD(NXD),XI(NIP),YD(NYD),YI(NIP),ZD(NXD,NYD),ZI(NIP)
1026 : integer :: INXI(NIP),INYI(NIP)
1027 :
1028 : ! .. Local Scalars ..
1029 0 : double precision :: A,B,C,D,DX,DXSQ,DY,DYSQ,P00,P01,P02,P03,P10,P11,
1030 0 : & P12,P13,P20,P21,P22,P23,P30,P31,P32,P33,Q0,Q1,Q2,
1031 0 : & Q3,U,V,X0,XII,Y0,YII,Z00,Z01,Z0DX,Z0DY,Z10,Z11,
1032 0 : & Z1DX,Z1DY,ZDXDY,ZII,ZX00,ZX01,ZX0DY,ZX10,ZX11,
1033 0 : & ZX1DY,ZXY00,ZXY01,ZXY10,ZXY11,ZY00,ZY01,ZY0DX,
1034 0 : & ZY10,ZY11,ZY1DX
1035 : integer :: IIP,IXD0,IXD1,IXDI,IXDIPV,IYD0,IYD1,IYDI,IYDIPV
1036 :
1037 : ! .. Intrinsic Functions ..
1038 : INTRINSIC MAX
1039 :
1040 : ! Calculation
1041 : ! Outermost DO-loop with respect to the output point
1042 0 : do IIP = 1,NIP
1043 0 : XII = XI(IIP)
1044 0 : YII = YI(IIP)
1045 0 : if (IIP == 1) then
1046 : IXDIPV = -1
1047 : IYDIPV = -1
1048 : else
1049 0 : IXDIPV = IXDI
1050 0 : IYDIPV = IYDI
1051 : end if
1052 0 : IXDI = INXI(IIP)
1053 0 : IYDI = INYI(IIP)
1054 : ! Retrieves the z and partial derivative values at the origin of
1055 : ! the coordinate for the rectangle.
1056 0 : if (IXDI /= IXDIPV .or. IYDI /= IYDIPV) then
1057 0 : IXD0 = MAX(1,IXDI)
1058 0 : IYD0 = MAX(1,IYDI)
1059 0 : X0 = XD(IXD0)
1060 0 : Y0 = YD(IYD0)
1061 0 : Z00 = ZD(IXD0,IYD0)
1062 0 : ZX00 = PDD(1,IXD0,IYD0)
1063 0 : ZY00 = PDD(2,IXD0,IYD0)
1064 0 : ZXY00 = PDD(3,IXD0,IYD0)
1065 : end if
1066 : ! Case 1. When the rectangle is inside the data area in both the
1067 : ! x and y directions.
1068 0 : if ((IXDI > 0.and.IXDI < NXD) .and. (IYDI > 0.and.IYDI < NYD)) then
1069 : ! Retrieves the z and partial derivative values at the other three
1070 : ! vertices of the rectangle.
1071 0 : if (IXDI /= IXDIPV .or. IYDI /= IYDIPV) then
1072 0 : IXD1 = IXD0 + 1
1073 0 : DX = XD(IXD1) - X0
1074 0 : DXSQ = DX*DX
1075 0 : IYD1 = IYD0 + 1
1076 0 : DY = YD(IYD1) - Y0
1077 0 : DYSQ = DY*DY
1078 0 : Z10 = ZD(IXD1,IYD0)
1079 0 : Z01 = ZD(IXD0,IYD1)
1080 0 : Z11 = ZD(IXD1,IYD1)
1081 0 : ZX10 = PDD(1,IXD1,IYD0)
1082 0 : ZX01 = PDD(1,IXD0,IYD1)
1083 0 : ZX11 = PDD(1,IXD1,IYD1)
1084 0 : ZY10 = PDD(2,IXD1,IYD0)
1085 0 : ZY01 = PDD(2,IXD0,IYD1)
1086 0 : ZY11 = PDD(2,IXD1,IYD1)
1087 0 : ZXY10 = PDD(3,IXD1,IYD0)
1088 0 : ZXY01 = PDD(3,IXD0,IYD1)
1089 0 : ZXY11 = PDD(3,IXD1,IYD1)
1090 : ! Calculates the polynomial coefficients.
1091 0 : Z0DX = (Z10-Z00)/DX
1092 0 : Z1DX = (Z11-Z01)/DX
1093 0 : Z0DY = (Z01-Z00)/DY
1094 0 : Z1DY = (Z11-Z10)/DY
1095 0 : ZX0DY = (ZX01-ZX00)/DY
1096 0 : ZX1DY = (ZX11-ZX10)/DY
1097 0 : ZY0DX = (ZY10-ZY00)/DX
1098 0 : ZY1DX = (ZY11-ZY01)/DX
1099 0 : ZDXDY = (Z1DY-Z0DY)/DX
1100 0 : A = ZDXDY - ZX0DY - ZY0DX + ZXY00
1101 0 : B = ZX1DY - ZX0DY - ZXY10 + ZXY00
1102 0 : C = ZY1DX - ZY0DX - ZXY01 + ZXY00
1103 0 : D = ZXY11 - ZXY10 - ZXY01 + ZXY00
1104 0 : P00 = Z00
1105 0 : P01 = ZY00
1106 0 : P02 = (2.0* (Z0DY-ZY00)+Z0DY-ZY01)/DY
1107 0 : P03 = (-2.0*Z0DY+ZY01+ZY00)/DYSQ
1108 0 : P10 = ZX00
1109 0 : P11 = ZXY00
1110 0 : P12 = (2.0* (ZX0DY-ZXY00)+ZX0DY-ZXY01)/DY
1111 0 : P13 = (-2.0*ZX0DY+ZXY01+ZXY00)/DYSQ
1112 0 : P20 = (2.0* (Z0DX-ZX00)+Z0DX-ZX10)/DX
1113 0 : P21 = (2.0* (ZY0DX-ZXY00)+ZY0DX-ZXY10)/DX
1114 0 : P22 = (3.0* (3.0*A-B-C)+D)/ (DX*DY)
1115 0 : P23 = (-6.0*A+2.0*B+3.0*C-D)/ (DX*DYSQ)
1116 0 : P30 = (-2.0*Z0DX+ZX10+ZX00)/DXSQ
1117 0 : P31 = (-2.0*ZY0DX+ZXY10+ZXY00)/DXSQ
1118 0 : P32 = (-6.0*A+3.0*B+2.0*C-D)/ (DXSQ*DY)
1119 0 : P33 = (2.0* (2.0*A-B-C)+D)/ (DXSQ*DYSQ)
1120 : end if
1121 : ! Evaluates the polynomial.
1122 0 : U = XII - X0
1123 0 : V = YII - Y0
1124 0 : Q0 = P00 + V* (P01+V* (P02+V*P03))
1125 0 : Q1 = P10 + V* (P11+V* (P12+V*P13))
1126 0 : Q2 = P20 + V* (P21+V* (P22+V*P23))
1127 0 : Q3 = P30 + V* (P31+V* (P32+V*P33))
1128 0 : ZII = Q0 + U* (Q1+U* (Q2+U*Q3))
1129 : ! End of Case 1
1130 : ! Case 2. When the rectangle is inside the data area in the x
1131 : ! direction but outside in the y direction.
1132 0 : else if ((IXDI > 0.and.IXDI < NXD) .and. (IYDI <= 0.or.IYDI >= NYD)) then
1133 : ! Retrieves the z and partial derivative values at the other
1134 : ! vertex of the semi-infinite rectangle.
1135 0 : if (IXDI /= IXDIPV .or. IYDI /= IYDIPV) then
1136 0 : IXD1 = IXD0 + 1
1137 0 : DX = XD(IXD1) - X0
1138 0 : DXSQ = DX*DX
1139 0 : Z10 = ZD(IXD1,IYD0)
1140 0 : ZX10 = PDD(1,IXD1,IYD0)
1141 0 : ZY10 = PDD(2,IXD1,IYD0)
1142 0 : ZXY10 = PDD(3,IXD1,IYD0)
1143 : ! Calculates the polynomial coefficients.
1144 0 : Z0DX = (Z10-Z00)/DX
1145 0 : ZY0DX = (ZY10-ZY00)/DX
1146 0 : P00 = Z00
1147 0 : P01 = ZY00
1148 0 : P10 = ZX00
1149 0 : P11 = ZXY00
1150 0 : P20 = (2.0* (Z0DX-ZX00)+Z0DX-ZX10)/DX
1151 0 : P21 = (2.0* (ZY0DX-ZXY00)+ZY0DX-ZXY10)/DX
1152 0 : P30 = (-2.0*Z0DX+ZX10+ZX00)/DXSQ
1153 0 : P31 = (-2.0*ZY0DX+ZXY10+ZXY00)/DXSQ
1154 : end if
1155 : ! Evaluates the polynomial.
1156 0 : U = XII - X0
1157 0 : V = YII - Y0
1158 0 : Q0 = P00 + V*P01
1159 0 : Q1 = P10 + V*P11
1160 0 : Q2 = P20 + V*P21
1161 0 : Q3 = P30 + V*P31
1162 0 : ZII = Q0 + U* (Q1+U* (Q2+U*Q3))
1163 : ! End of Case 2
1164 : ! Case 3. When the rectangle is outside the data area in the x
1165 : ! direction but inside in the y direction.
1166 0 : else if ((IXDI <= 0.or.IXDI >= NXD) .and. (IYDI > 0.and.IYDI < NYD)) then
1167 : ! Retrieves the z and partial derivative values at the other
1168 : ! vertex of the semi-infinite rectangle.
1169 0 : if (IXDI /= IXDIPV .or. IYDI /= IYDIPV) then
1170 0 : IYD1 = IYD0 + 1
1171 0 : DY = YD(IYD1) - Y0
1172 0 : DYSQ = DY*DY
1173 0 : Z01 = ZD(IXD0,IYD1)
1174 0 : ZX01 = PDD(1,IXD0,IYD1)
1175 0 : ZY01 = PDD(2,IXD0,IYD1)
1176 0 : ZXY01 = PDD(3,IXD0,IYD1)
1177 : ! Calculates the polynomial coefficients.
1178 0 : Z0DY = (Z01-Z00)/DY
1179 0 : ZX0DY = (ZX01-ZX00)/DY
1180 0 : P00 = Z00
1181 0 : P01 = ZY00
1182 0 : P02 = (2.0* (Z0DY-ZY00)+Z0DY-ZY01)/DY
1183 0 : P03 = (-2.0*Z0DY+ZY01+ZY00)/DYSQ
1184 0 : P10 = ZX00
1185 0 : P11 = ZXY00
1186 0 : P12 = (2.0* (ZX0DY-ZXY00)+ZX0DY-ZXY01)/DY
1187 0 : P13 = (-2.0*ZX0DY+ZXY01+ZXY00)/DYSQ
1188 : end if
1189 : ! Evaluates the polynomial.
1190 0 : U = XII - X0
1191 0 : V = YII - Y0
1192 0 : Q0 = P00 + V* (P01+V* (P02+V*P03))
1193 0 : Q1 = P10 + V* (P11+V* (P12+V*P13))
1194 0 : ZII = Q0 + U*Q1
1195 : ! End of Case 3
1196 : ! Case 4. When the rectangle is outside the data area in both the
1197 : ! x and y direction.
1198 0 : else if ((IXDI <= 0.or.IXDI >= NXD) .and. (IYDI <= 0.or.IYDI >= NYD)) then
1199 : ! Calculates the polynomial coefficients.
1200 0 : if (IXDI /= IXDIPV .or. IYDI /= IYDIPV) then
1201 0 : P00 = Z00
1202 0 : P01 = ZY00
1203 0 : P10 = ZX00
1204 0 : P11 = ZXY00
1205 : end if
1206 : ! Evaluates the polynomial.
1207 0 : U = XII - X0
1208 0 : V = YII - Y0
1209 0 : Q0 = P00 + V*P01
1210 0 : Q1 = P10 + V*P11
1211 0 : ZII = Q0 + U*Q1
1212 : end if
1213 : ! End of Case 4
1214 0 : ZI(IIP) = ZII
1215 : end do
1216 0 : return
1217 : end subroutine RGPLNL_db
|
