Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module interp_2d_lib_sg
21 : ! single precision library for 2D interpolation
22 :
23 : implicit none
24 :
25 : contains ! the procedure interface for the library
26 : ! client programs should only call these routines.
27 :
28 : ! Contents
29 :
30 : ! Rectangular-grid of data points
31 :
32 : ! interp_RGBI3P_sg -- point interpolation (Akima)
33 : ! interp_RGSF3P_sg -- surface interpolation (Akima)
34 :
35 : ! interp_mkbicub_sg -- bicubic splines
36 :
37 : ! Scattered set of data points
38 :
39 : ! interp_CS2VAL_sg -- point interpolation (Renka)
40 : ! interp_CS2GRD_sg -- point interpolation with gradients (Renka)
41 :
42 : ! see the documents in refs directory for more info.
43 :
44 : ! ***********************************************************************
45 : ! ***********************************************************************
46 : ! ***********************************************************************
47 : ! ***********************************************************************
48 :
49 : ! Rectangular-grid bivariate interpolation and surface fitting
50 : ! from ACM Algorithm 760., ACM Trans. Math. Software (22) 1996, 357-361.
51 : ! Hiroshi Akima
52 : ! U.S. Department of Commerce, NTIA/ITS
53 : ! Version of 1995/08
54 :
55 : ! NOTE: versions for scattered data follow.
56 :
57 : ! see interp_2d_lib_db for double precision version.
58 :
59 0 : subroutine interp_RGBI3P_sg(MD, NXD, NYD, XD, YD, ZD, NIP, XI, YI, ZI, ierr, WK)
60 : integer, intent(in) :: MD, NXD, NYD, NIP
61 : real, intent(in) :: XD(NXD), YD(NYD), ZD(NXD, NYD), XI(NIP), YI(NIP)
62 : real, intent(inout) :: ZI(NIP), WK(3, NXD, NYD)
63 : integer, intent(out) :: ierr
64 : !
65 : ! This subroutine performs interpolation of a bivariate function,
66 : ! z(x,y), on a rectangular grid in the x-y plane. It is based on
67 : ! the revised Akima method.
68 : !
69 : ! In this subroutine, the interpolating function is a piecewise
70 : ! function composed of a set of bicubic (bivariate third-degree)
71 : ! polynomials, each applicable to a rectangle of the input grid
72 : ! in the x-y plane. Each polynomial is determined locally.
73 : !
74 : ! This subroutine has the accuracy of a bicubic polynomial, i.e.,
75 : ! it interpolates accurately when all data points lie on a
76 : ! surface of a bicubic polynomial.
77 : !
78 : ! The grid lines can be unevenly spaced.
79 : !
80 : ! The input arguments are
81 : ! MD = mode of computation
82 : ! = 1 for new XD, YD, or ZD data (default)
83 : ! = 2 for old XD, YD, and ZD data,
84 : ! NXD = number of the input-grid data points in the x
85 : ! coordinate (must be 2 or greater),
86 : ! NYD = number of the input-grid data points in the y
87 : ! coordinate (must be 2 or greater),
88 : ! XD = array of dimension NXD containing the x coordinates
89 : ! of the input-grid data points (must be in a
90 : ! monotonic increasing order),
91 : ! YD = array of dimension NYD containing the y coordinates
92 : ! of the input-grid data points (must be in a
93 : ! monotonic increasing order),
94 : ! ZD = two-dimensional array of dimension NXD*NYD
95 : ! containing the z(x,y) values at the input-grid data
96 : ! points,
97 : ! NIP = number of the output points at which interpolation
98 : ! of the z value is desired (must be 1 or greater),
99 : ! XI = array of dimension NIP containing the x coordinates
100 : ! of the output points,
101 : ! YI = array of dimension NIP containing the y coordinates
102 : ! of the output points.
103 : !
104 : ! The output arguments are
105 : ! ZI = array of dimension NIP where the interpolated z
106 : ! values at the output points are to be stored,
107 : ! ierr = error flag
108 : ! = 0 for no errors
109 : ! = 1 for NXD = 1 or less
110 : ! = 2 for NYD = 1 or less
111 : ! = 3 for identical XD values or
112 : ! XD values out of sequence
113 : ! = 4 for identical YD values or
114 : ! YD values out of sequence
115 : ! = 5 for NIP = 0 or less.
116 : !
117 : ! The other argument is
118 : ! WK = three dimensional array of dimension 3*NXD*NYD used
119 : ! internally as a work area.
120 : !
121 : ! The very first call to this subroutine and the call with a new
122 : ! XD, YD, and ZD array must be made with MD=1. The call with MD=2
123 : ! must be preceded by another call with the same XD, YD, and ZD
124 : ! arrays. Between the call with MD=2 and its preceding call, the
125 : ! WK array must not be disturbed.
126 : !
127 0 : call do_RGBI3P_sg(MD, NXD, NYD, XD, YD, ZD, NIP, XI, YI, ZI, ierr, WK)
128 0 : end subroutine interp_RGBI3P_sg
129 :
130 0 : subroutine interp_RGSF3P_sg(MD, NXD, NYD, XD, YD, ZD, NXI, XI, NYI, YI, ZI, ierr, WK)
131 : integer, intent(in) :: MD, NXD, NYD, NXI, NYI
132 : real, intent(in) :: XD(NXD), YD(NYD), ZD(NXD, NYD), XI(NXI), YI(NYI)
133 : real, intent(inout) :: ZI(NXI, NYI), WK(3, NXD, NYD)
134 : integer, intent(out) :: ierr
135 : !
136 : ! Rectangular-grid surface fitting
137 : ! (a master subroutine of the RGBI3P/RGSF3P_sg subroutine package)
138 : !
139 : ! Hiroshi Akima
140 : ! U.S. Department of Commerce, NTIA/ITS
141 : ! Version of 1995/08
142 : !
143 : ! This subroutine performs surface fitting by interpolating
144 : ! values of a bivariate function, z(x,y), on a rectangular grid
145 : ! in the x-y plane. It is based on the revised Akima method.
146 : !
147 : ! In this subroutine, the interpolating function is a piecewise
148 : ! function composed of a set of bicubic (bivariate third-degree)
149 : ! polynomials, each applicable to a rectangle of the input grid
150 : ! in the x-y plane. Each polynomial is determined locally.
151 : !
152 : ! This subroutine has the accuracy of a bicubic polynomial, i.e.,
153 : ! it fits the surface accurately when all data points lie on a
154 : ! surface of a bicubic polynomial.
155 : !
156 : ! The grid lines of the input and output data can be unevenly
157 : ! spaced.
158 : !
159 : ! The input arguments are
160 : ! MD = mode of computation
161 : ! = 1 for new XD, YD, or ZD data (default)
162 : ! = 2 for old XD, YD, and ZD data,
163 : ! NXD = number of the input-grid data points in the x
164 : ! coordinate (must be 2 or greater),
165 : ! NYD = number of the input-grid data points in the y
166 : ! coordinate (must be 2 or greater),
167 : ! XD = array of dimension NXD containing the x coordinates
168 : ! of the input-grid data points (must be in a
169 : ! monotonic increasing order),
170 : ! YD = array of dimension NYD containing the y coordinates
171 : ! of the input-grid data points (must be in a
172 : ! monotonic increasing order),
173 : ! ZD = two-dimensional array of dimension NXD*NYD
174 : ! containing the z(x,y) values at the input-grid data
175 : ! points,
176 : ! NXI = number of output grid points in the x coordinate
177 : ! (must be 1 or greater),
178 : ! XI = array of dimension NXI containing the x coordinates
179 : ! of the output grid points,
180 : ! NYI = number of output grid points in the y coordinate
181 : ! (must be 1 or greater),
182 : ! YI = array of dimension NYI containing the y coordinates
183 : ! of the output grid points.
184 : !
185 : ! The output arguments are
186 : ! ZI = two-dimensional array of dimension NXI*NYI, where
187 : ! the interpolated z values at the output grid
188 : ! points are to be stored,
189 : ! ierr = error flag
190 : ! = 0 for no error
191 : ! = 1 for NXD = 1 or less
192 : ! = 2 for NYD = 1 or less
193 : ! = 3 for identical XD values or
194 : ! XD values out of sequence
195 : ! = 4 for identical YD values or
196 : ! YD values out of sequence
197 : ! = 5 for NXI = 0 or less
198 : ! = 6 for NYI = 0 or less.
199 : !
200 : ! The other argument is
201 : ! WK = three-dimensional array of dimension 3*NXD*NYD used
202 : ! internally as a work area.
203 : !
204 : ! The very first call to this subroutine and the call with a new
205 : ! XD, YD, or ZD array must be made with MD=1. The call with MD=2
206 : ! must be preceded by another call with the same XD, YD, and ZD
207 : ! arrays. Between the call with MD=2 and its preceding call, the
208 : ! WK array must not be disturbed.
209 :
210 0 : call do_RGSF3P_sg(MD, NXD, NYD, XD, YD, ZD, NXI, XI, NYI, YI, ZI, ierr, WK)
211 0 : end subroutine interp_RGSF3P_sg
212 :
213 : ! ***********************************************************************
214 : ! ***********************************************************************
215 : ! ***********************************************************************
216 : ! ***********************************************************************
217 :
218 : ! cubic Shepard method for bivariate interpolation of scattered data.
219 : ! from ACM Algorithm 790., ACM Trans. Math. Software (25) 1999, 70-73.
220 : ! Robert J. Renka
221 : ! Dept. of Computer Science
222 : ! Univ. of North Texas
223 : ! renka@cs.unt.edu
224 :
225 : ! use CSHEP2_sg to set up the interpolation information.
226 : ! use CS2VAL_sg to evaluate it.
227 : ! use CS2GRD_sg to get value and derivatives.
228 : ! see interp_2d_lib_db for double precision versions.
229 :
230 0 : subroutine interp_CSHEP2_sg(N, X, Y, F, NC, NW, NR, LCELL, LNEXT, XMIN, YMIN, DX, DY, RMAX, RW, A, ierr)
231 : integer, intent(in) :: N, NC, NW, NR
232 : integer, intent(out) :: LCELL(NR, NR), LNEXT(N), ierr
233 : real, intent(in) :: X(N), Y(N), F(N)
234 : real, intent(inout) :: XMIN, YMIN, DX, DY, RMAX, RW(N), A(9, N)
235 : !
236 : ! Here's the [slightly edited] documentation.
237 : !
238 : ! This subroutine computes a set of parameters defining a
239 : ! C2 (twice continuously differentiable) bivariate function
240 : ! C(X,Y) which interpolates data values F at a set of N
241 : ! arbitrarily distributed points (X,Y) in the plane (nodes).
242 : ! The interpolant C may be evaluated at an arbitrary point
243 : ! by function CS2VAL_sg, and its first partial derivatives are
244 : ! computed by Subroutine CS2GRD_sg.
245 : !
246 : ! The interpolation scheme is a modified Cubic Shepard
247 : ! method:
248 : !
249 : ! C = [W(1)*C(1)+W(2)*C(2)+..+W(N)*C(N)]/[W(1)+W(2)+..+W(N)]
250 : !
251 : ! for bivariate functions W(k) and C(k). The nodal func-
252 : ! tions are given by
253 : !
254 : ! C(k)(x,y) = A(1,k)*(x-X(k))**3 +
255 : ! A(2,k)*(x-X(k))**2*(y-Y(k)) +
256 : ! A(3,k)*(x-X(k))*(y-Y(k))**2 +
257 : ! A(4,k)*(y-Y(k))**3 + A(5,k)*(x-X(k))**2 +
258 : ! A(6,k)*(x-X(k))*(y-Y(k)) + A(7,k)*(y-Y(k))**2
259 : ! + A(8,k)*(x-X(k)) + A(9,k)*(y-Y(k)) + F(k) .
260 : !
261 : ! Thus, C(k) is a cubic function which interpolates the data
262 : ! value at node k. Its coefficients A(,k) are obtained by a
263 : ! weighted least squares fit to the closest NC data points
264 : ! with weights similar to W(k). Note that the radius of
265 : ! influence for the least squares fit is fixed for each k,
266 : ! but varies with k.
267 : !
268 : ! The weights are taken to be
269 : !
270 : ! W(k)(x,y) = ( (R(k)-D(k))+ / R(k)*D(k) )**3 ,
271 : !
272 : ! where (R(k)-D(k))+ = 0 if R(k) < D(k), and D(k)(x,y) is
273 : ! the Euclidean distance between (x,y) and (X(k),Y(k)). The
274 : ! radius of influence R(k) varies with k and is chosen so
275 : ! that NW nodes are within the radius. Note that W(k) is
276 : ! not defined at node (X(k),Y(k)), but C(x,y) has limit F(k)
277 : ! as (x,y) approaches (X(k),Y(k)).
278 : ! On input:
279 : !
280 : ! N = Number of nodes and data values. N .GE. 10.
281 : !
282 : ! X,Y = Arrays of length N containing the Cartesian
283 : ! coordinates of the nodes.
284 : !
285 : ! F = Array of length N containing the data values
286 : ! in one-to-one correspondence with the nodes.
287 : !
288 : ! NC = Number of data points to be used in the least
289 : ! squares fit for coefficients defining the nodal
290 : ! functions C(k). Values found to be optimal for
291 : ! test data sets ranged from 11 to 25. A recom-
292 : ! mended value for general data sets is NC = 17.
293 : ! For nodes lying on (or close to) a rectangular
294 : ! grid, the recommended value is NC = 11. In any
295 : ! case, NC must be in the range 9 to Min(40,N-1).
296 : !
297 : ! NW = Number of nodes within (and defining) the radii
298 : ! of influence R(k) which enter into the weights
299 : ! W(k). For N sufficiently large, a recommended
300 : ! value is NW = 30. In general, NW should be
301 : ! about 1.5*NC. 1 .LE. NW .LE. Min(40,N-1).
302 : !
303 : ! NR = Number of rows and columns in the cell grid de-
304 : ! fined in Subroutine STORE2. A rectangle con-
305 : ! taining the nodes is partitioned into cells in
306 : ! order to increase search efficiency. NR =
307 : ! Sqrt(N/3) is recommended. NR .GE. 1.
308 : !
309 : ! The above parameters are not altered by this routine.
310 : !
311 : ! LCELL = Array of length .GE. NR**2.
312 : !
313 : ! LNEXT = Array of length .GE. N.
314 : !
315 : ! RW = Array of length .GE. N.
316 : !
317 : ! A = Array of length .GE. 9N.
318 : !
319 : ! On output:
320 : !
321 : ! LCELL = NR by NR array of nodal indexes associated
322 : ! with cells. Refer to Subroutine STORE2.
323 : !
324 : ! LNEXT = Array of length N containing next-node
325 : ! indexes. Refer to Subroutine STORE2.
326 : !
327 : ! XMIN,YMIN,DX,DY = Minimum nodal coordinates and cell
328 : ! dimensions. Refer to Subroutine
329 : ! STORE2.
330 : !
331 : ! RMAX = Largest element in RW -- maximum radius R(k).
332 : !
333 : ! RW = Array containing the the radii R(k) which enter
334 : ! into the weights W(k).
335 : !
336 : ! A = 9 by N array containing the coefficients for
337 : ! cubic nodal function C(k) in column k.
338 : !
339 : ! Note that the output parameters described above are not
340 : ! defined unless ierr = 0.
341 : !
342 : ! ierr = Error indicator:
343 : ! ierr = 0 if no errors were encountered.
344 : ! ierr = 1 if N, NC, NW, or NR is outside its
345 : ! valid range.
346 : ! ierr = 2 if duplicate nodes were encountered.
347 : ! ierr = 3 if all nodes are collinear.
348 :
349 0 : call do_CSHEP2_sg(N, X, Y, F, NC, NW, NR, LCELL, LNEXT, XMIN, YMIN, DX, DY, RMAX, RW, A, ierr)
350 :
351 0 : end subroutine interp_CSHEP2_sg
352 :
353 0 : real function interp_CS2VAL_sg(PX, PY, N, X, Y, F, NR, LCELL, LNEXT, XMIN, YMIN, DX, DY, RMAX, RW, A, ierr)
354 : integer, intent(in) :: N, NR, LCELL(NR, NR), LNEXT(N)
355 : real, intent(in) :: PX, PY, X(N), Y(N), F(N), XMIN, YMIN, DX, DY, RMAX, RW(N), A(9, N)
356 : integer, intent(out) :: ierr
357 : !
358 : ! This function returns the value C(PX,PY), where C is the
359 : ! weighted sum of cubic nodal functions defined in Subrou-
360 : ! tine CSHEP2_sg. CS2GRD_sg may be called to compute a gradient
361 : ! of C along with the value, and/or to test for errors.
362 : ! CS2HES may be called to compute a value, first partial
363 : ! derivatives, and second partial derivatives at a point.
364 : !
365 : ! On input:
366 : !
367 : ! PX,PY = Cartesian coordinates of the point P at
368 : ! which C is to be evaluated.
369 : !
370 : ! N = Number of nodes and data values defining C.
371 : ! N .GE. 10.
372 : !
373 : ! X,Y,F = Arrays of length N containing the nodes and
374 : ! data values interpolated by C.
375 : !
376 : ! NR = Number of rows and columns in the cell grid.
377 : ! Refer to Subroutine STORE2. NR .GE. 1.
378 : !
379 : ! LCELL = NR by NR array of nodal indexes associated
380 : ! with cells. Refer to Subroutine STORE2.
381 : !
382 : ! LNEXT = Array of length N containing next-node
383 : ! indexes. Refer to Subroutine STORE2.
384 : !
385 : ! XMIN,YMIN,DX,DY = Minimum nodal coordinates and cell
386 : ! dimensions. DX and DY must be
387 : ! positive. Refer to Subroutine
388 : ! STORE2.
389 : !
390 : ! RMAX = Largest element in RW -- maximum radius R(k).
391 : !
392 : ! RW = Array containing the the radii R(k) which enter
393 : ! into the weights W(k) defining C.
394 : !
395 : ! A = 9 by N array containing the coefficients for
396 : ! cubic nodal function C(k) in column k.
397 : !
398 : ! Input parameters are not altered by this function. The
399 : ! parameters other than PX and PY should be input unaltered
400 : ! from their values on output from CSHEP2. This function
401 : ! should not be called if a nonzero error flag was returned
402 : ! by CSHEP2_sg.
403 : !
404 : ! On output:
405 : !
406 : ! CS2VAL_sg = Function value C(PX,PY) unless N, NR, DX,
407 : ! DY, or RMAX is invalid, in which case ierr /= 0.
408 : !
409 : ! ierr = Error indicator:
410 : ! ierr = 0 if no errors were encountered.
411 : !
412 : real :: do_CS2VAL_sg
413 0 : interp_CS2VAL_sg = do_CS2VAL_sg(PX, PY, N, X, Y, F, NR, LCELL, LNEXT, XMIN, YMIN, DX, DY, RMAX, RW, A, ierr)
414 0 : end function interp_CS2VAL_sg
415 :
416 0 : subroutine interp_CS2GRD_sg(PX, PY, N, X, Y, F, NR, LCELL, LNEXT, XMIN, YMIN, DX, DY, RMAX, RW, A, C, CX, CY, ierr)
417 : integer, intent(in) :: N, NR, LCELL(NR, NR), LNEXT(N)
418 : real, intent(in) :: PX, PY, X(N), Y(N), F(N), XMIN, YMIN, DX, DY, RMAX, RW(N), A(9, N)
419 : real, intent(out) :: C, CX, CY
420 : integer, intent(out) :: ierr
421 : !
422 : ! This subroutine computes the value and gradient at P =
423 : ! (PX,PY) of the interpolatory function C defined in Sub-
424 : ! routine CSHEP2_sg. C is a weighted sum of cubic nodal
425 : ! functions.
426 : !
427 : ! On input:
428 : !
429 : ! PX,PY = Cartesian coordinates of the point P at
430 : ! which C and its partial derivatives are
431 : ! to be evaluated.
432 : !
433 : ! N = Number of nodes and data values defining C.
434 : ! N .GE. 10.
435 : !
436 : ! X,Y,F = Arrays of length N containing the nodes and
437 : ! data values interpolated by C.
438 : !
439 : ! NR = Number of rows and columns in the cell grid.
440 : ! Refer to Subroutine STORE2. NR .GE. 1.
441 : !
442 : ! LCELL = NR by NR array of nodal indexes associated
443 : ! with cells. Refer to Subroutine STORE2.
444 : !
445 : ! LNEXT = Array of length N containing next-node
446 : ! indexes. Refer to Subroutine STORE2.
447 : !
448 : ! XMIN,YMIN,DX,DY = Minimum nodal coordinates and cell
449 : ! dimensions. DX and DY must be
450 : ! positive. Refer to Subroutine
451 : ! STORE2.
452 : !
453 : ! RMAX = Largest element in RW -- maximum radius R(k).
454 : !
455 : ! RW = Array of length N containing the the radii R(k)
456 : ! which enter into the weights W(k) defining C.
457 : !
458 : ! A = 9 by N array containing the coefficients for
459 : ! cubic nodal function C(k) in column k.
460 : !
461 : ! Input parameters are not altered by this subroutine.
462 : ! The parameters other than PX and PY should be input
463 : ! unaltered from their values on output from CSHEP2_sg. This
464 : ! subroutine should not be called if a nonzero error flag
465 : ! was returned by CSHEP2_sg.
466 : !
467 : ! On output:
468 : !
469 : ! C = Value of C at (PX,PY) unless ierr .EQ. 1, in
470 : ! which case no values are returned.
471 : !
472 : ! CX,CY = First partial derivatives of C at (PX,PY)
473 : ! unless ierr .EQ. 1.
474 : !
475 : ! ierr = Error indicator:
476 : ! ierr = 0 if no errors were encountered.
477 : ! ierr = 1 if N, NR, DX, DY or RMAX is invalid.
478 : ! ierr = 2 if no errors were encountered but
479 : ! (PX,PY) is not within the radius R(k)
480 : ! for any node k (and thus C=CX=CY=0).
481 :
482 0 : call do_CS2GRD_sg(PX, PY, N, X, Y, F, NR, LCELL, LNEXT, XMIN, YMIN, DX, DY, RMAX, RW, A, C, CX, CY, ierr)
483 0 : end subroutine interp_CS2GRD_sg
484 :
485 : ! ***********************************************************************
486 : ! ***********************************************************************
487 : ! ***********************************************************************
488 : ! ***********************************************************************
489 :
490 : ! bicubic splines and hermite
491 :
492 : ! from PSPLINE by Doug McCune (version as of February, 2004)
493 : ! PSPLINE Home Page:
494 : ! http://w3.pppl.gov/NTCC/PSPLINE
495 : ! Doug McCune, Princeton University
496 : ! dmccune@pppl.gov
497 :
498 : ! use interp_mkbicub_sg to set up the interpolation information
499 : ! use interp_evbicub_sg to evaluate it
500 : ! see interp_2d_lib_db for double precision versions.
501 :
502 1 : subroutine interp_mkbicub_sg(x, nx, y, ny, f1, nf2, &
503 1 : ibcxmin, bcxmin, ibcxmax, bcxmax, ibcymin, bcymin, ibcymax, bcymax, ilinx, iliny, ierr)
504 :
505 : use bicub_sg
506 :
507 : integer, intent(in) :: nx ! length of x vector
508 : integer, intent(in) :: ny ! length of y vector
509 : real, intent(in) :: x(:) ! (nx) ! x vector, strict ascending
510 : real, intent(in) :: y(:) ! (ny) ! y vector, strict ascending
511 : integer, intent(in) :: nf2 ! 2nd dimension of f, nf2.ge.nx
512 : real, intent(inout), pointer :: f1(:) ! =(4,nf2,ny) ! data & spline coefficients
513 : !
514 : ! on input: f(1,i,j) = f(x(i),y(j))
515 : ! on output: f(1,i,j) unchanged
516 : ! f(2,i,j) = d2f/dx2(x(i),y(j))
517 : ! f(3,i,j) = d2f/dy2(x(i),y(j))
518 : ! f(4,i,j) = d4f/dx2dy2(x(i),y(j))
519 : !
520 : ! and the interpolation formula for (x,y) in (x(i),x(i+1))^(y(j),y(j+1))
521 : ! is:
522 : ! hx = x(i+1)-x(i) hy = y(j+1)-y(j)
523 : ! dxp= (x-x(i))/hx dxm= 1-dxp dxp,dxm in (0,1)
524 : ! dyp= (x-x(i))/hx dym= 1-dyp dyp,dym in (0,1)
525 : ! dx3p = dxp**3-dxp dx3m = dxm**3-dxm dxp3,dxm3 in (0,1)
526 : !
527 : ! finterp = dxm*(dym*f(1,i,j)+dyp*f(1,i,j+1))
528 : ! +dxp*(dym*f(1,i+1,j)+dyp*f(1,i+1,j+1))
529 : ! +1/6*hx**2*
530 : ! dx3m*(dym*f(2,i,j)+dyp*f(2,i,j+1))
531 : ! +dx3p*(dym*f(2,i+1,j)+dyp*f(2,i+1,j+1))
532 : ! +1/6*hy**2*
533 : ! dxm*(dy3m*f(3,i,j)+dy3p*f(3,i,j+1))
534 : ! +dxp*(dy3m*f(3,i+1,j)+dy3p*f(3,i+1,j+1))
535 : ! +1/36*hx**2*hy**2*
536 : ! dx3m*(dym*f(4,i,j)+dyp*f(4,i,j+1))
537 : ! +dx3p*(dym*f(4,i+1,j)+dyp*f(4,i+1,j+1))
538 : !
539 : ! where the f(2:4,*,*) are cleverly evaluated to assure
540 : ! (a) finterp is continuous and twice differentiable across all
541 : ! grid cell boundaries, and
542 : ! (b) all boundary conditions are satisfied.
543 : !
544 : ! input bdy condition data:
545 : integer, intent(in) :: ibcxmin ! bc flag for x=xmin
546 : real, intent(in) :: bcxmin(:) ! (ny) ! bc data vs. y at x=xmin
547 : integer, intent(in) :: ibcxmax ! bc flag for x=xmax
548 : real, intent(in) :: bcxmax(:) ! (ny) ! bc data vs. y at x=xmax
549 : integer, intent(in) :: ibcymin ! bc flag for y=ymin
550 : real, intent(in) :: bcymin(:) ! (nx) ! bc data vs. x at y=ymin
551 : integer, intent(in) :: ibcymax ! bc flag for y=ymax
552 : real, intent(in) :: bcymax(:) ! (nx) ! bc data vs. x at y=ymax
553 : !
554 : ! with interpretation:
555 : ! ibcxmin -- indicator for boundary condition at x(1):
556 : ! bcxmin(...) -- boundary condition data
557 : ! =-1 -- periodic boundary condition
558 : ! =0 -- use "not a knot"
559 : ! =1 -- match slope, specified at x(1),th(ith) by bcxmin(ith)
560 : ! =2 -- match 2nd derivative, specified at x(1),th(ith) by bcxmin(ith)
561 : ! =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all th(j)
562 : ! =4 -- boundary condition is d2f/dx2=0 at x(1), all th(j)
563 : ! =5 -- match 1st derivative to 1st divided difference
564 : ! =6 -- match 2nd derivative to 2nd divided difference
565 : ! =7 -- match 3rd derivative to 3rd divided difference
566 : ! (for more detailed definition of B!s 5-7, see the
567 : ! comments of subroutine mkspline)
568 : ! NOTE bcxmin(...) referenced ONLY if ibcxmin=1 or ibcxmin=2
569 : !
570 : ! ibcxmax -- indicator for boundary condition at x(nx):
571 : ! bcxmax(...) -- boundary condition data
572 : ! (interpretation as with ibcxmin, bcxmin)
573 : ! NOTE: if ibcxmin=-1, ibcxmax is ignored! ...and the B! is periodic.
574 : !
575 : ! and analogous interpretation for ibcymin,bcymin,ibcymax,bcymax
576 : ! (df/dy or d2f/dy2 boundary conditions at y=ymin and y=ymax).
577 : !
578 : ! output linear grid flags and completion code (ierr=0 is normal):
579 : !
580 : integer, intent(out) :: ilinx ! =1: x grid is "nearly" equally spaced
581 : integer, intent(out) :: iliny ! =1: y grid is "nearly" equally spaced
582 : ! ilinx and iliny are set to zero if corresponding grids are not equally
583 : ! spaced
584 : !
585 : integer, intent(out) :: ierr ! =0 on exit if there is no error.
586 : !
587 : ! if there is an error, ierr is set and a message is output on unit 6.
588 : ! these are considered programming errors in the calling routine.
589 : !
590 : ! possible errors:
591 : ! x(...) not strict ascending
592 : ! y(...) not strict ascending
593 : ! nx .lt. 4
594 : ! ny .lt. 4
595 : ! invalid boundary condition flag
596 : !
597 : !-----------------------
598 :
599 1 : call do_mkbicub(x, nx, y, ny, f1, nf2, ibcxmin, bcxmin, ibcxmax, bcxmax, ibcymin, bcymin, ibcymax, bcymax, ilinx, iliny, ierr)
600 :
601 1 : end subroutine interp_mkbicub_sg
602 :
603 6 : subroutine interp_evbicub_sg(xget, yget, x, nx, y, ny, ilinx, iliny, f1, inf2, ict, fval, ierr)
604 :
605 : use bicub_sg
606 :
607 : ! evaluate bicubic spline interpolant (single precision version)
608 : !
609 : ! use mkbicub to set up spline coefficients
610 : !
611 : ! evaluate a 2d cubic Spline interpolant on a rectilinear
612 : ! grid -- this is C2 in both directions.
613 : !
614 : ! this subroutine calls two subroutines:
615 : ! herm2xy -- find cell containing (xget,yget)
616 : ! fvbicub -- evaluate interpolant function and (optionally) derivatives
617 :
618 : ! input arguments:
619 : ! ================
620 : !
621 : integer, intent(in) :: nx, ny ! grid sizes
622 : real, intent(in) :: xget, yget ! target of this interpolation
623 : real, intent(in) :: x(:) ! (nx) ! ordered x grid
624 : real, intent(in) :: y(:) ! (ny) ! ordered y grid
625 : integer, intent(in) :: ilinx ! ilinx=1 => assume x evenly spaced
626 : integer, intent(in) :: iliny ! iliny=1 => assume y evenly spaced
627 : !
628 : integer, intent(in) :: inf2
629 : real, intent(inout), pointer :: f1(:) ! function data
630 : !
631 : ! f 2nd dimension inf2 must be .ge. nx
632 : ! contents of f:
633 : !
634 : ! f(0,i,j) = f @ x(i),y(j)
635 : ! f(1,i,j) = d2f/dx2 @ x(i),y(j)
636 : ! f(2,i,j) = d2f/dy2 @ x(i),y(j)
637 : ! f(3,i,j) = d4f/dx2dy2 @ x(i),y(j)
638 : !
639 : ! (these are spline coefficients selected for continuous 2-
640 : ! diffentiability, see mkbicub[w].for)
641 : !
642 : integer, intent(in) :: ict(6) ! code specifying output desired
643 : !
644 : ! ict(1)=1 -- return f (0, don't)
645 : ! ict(2)=1 -- return df/dx (0, don't)
646 : ! ict(3)=1 -- return df/dy (0, don't)
647 : ! ict(4)=1 -- return d2f/dx2 (0, don't)
648 : ! ict(5)=1 -- return d2f/dy2 (0, don't)
649 : ! ict(6)=1 -- return d2f/dxdy (0, don't)
650 : ! the number of non zero values ict(1:6)
651 : ! determines the number of outputs...
652 : !
653 : ! new dmc December 2005 -- access to higher derivatives (even if not
654 : ! continuous-- but can only go up to 3rd derivatives on any one coordinate.
655 : ! if ict(1)=3 -- want 3rd derivatives
656 : ! ict(2)=1 for d3f/dx3
657 : ! ict(3)=1 for d3f/dx2dy
658 : ! ict(4)=1 for d3f/dxdy2
659 : ! ict(5)=1 for d3f/dy3
660 : ! number of non-zero values ict(2:5) gives no. of outputs
661 : ! if ict(1)=4 -- want 4th derivatives
662 : ! ict(2)=1 for d4f/dx3dy
663 : ! ict(3)=1 for d4f/dx2dy2
664 : ! ict(4)=1 for d4f/dxdy3
665 : ! number of non-zero values ict(2:4) gives no. of outputs
666 : ! if ict(1)=5 -- want 5th derivatives
667 : ! ict(2)=1 for d5f/dx3dy2
668 : ! ict(3)=1 for d5f/dx2dy3
669 : ! number of non-zero values ict(2:3) gives no. of outputs
670 : ! if ict(1)=6 -- want 6th derivatives
671 : ! d6f/dx3dy3 -- one value is returned.
672 : !
673 : ! output arguments:
674 : ! =================
675 : !
676 : real, intent(inout) :: fval(6) ! output data
677 : integer, intent(out) :: ierr ! error code =0 ==> no error
678 : !
679 : ! fval(1) receives the first output (depends on ict(...) spec)
680 : ! fval(2) receives the second output (depends on ict(...) spec)
681 : ! fval(3) receives the third output (depends on ict(...) spec)
682 : ! fval(4) receives the fourth output (depends on ict(...) spec)
683 : ! fval(5) receives the fourth output (depends on ict(...) spec)
684 : ! fval(6) receives the fourth output (depends on ict(...) spec)
685 : !
686 : ! examples:
687 : ! on input ict = [1,1,1,0,0,1]
688 : ! on output fval= [f,df/dx,df/dy,d2f/dxdy], elements 5 & 6 not referenced.
689 : !
690 : ! on input ict = [1,0,0,0,0,0]
691 : ! on output fval= [f] ... elements 2 -- 6 never referenced.
692 : !
693 : ! on input ict = [0,0,0,1,1,0]
694 : ! on output fval= [d2f/dx2,d2f/dy2] ... elements 3 -- 6 never referenced.
695 : !
696 : ! on input ict = [0,0,1,0,0,0]
697 : ! on output fval= [df/dy] ... elements 2 -- 6 never referenced.
698 : !
699 : ! ierr -- completion code: 0 means OK
700 : !-------------------
701 :
702 : integer :: i, j
703 : integer :: ii(1), jj(1)
704 42 : real :: xparam(1), yparam(1), hx(1), hxi(1), hy(1), hyi(1)
705 6 : call herm2xy(xget, yget, x, nx, y, ny, ilinx, iliny, i, j, xparam(1), yparam(1), hx(1), hxi(1), hy(1), hyi(1), ierr)
706 6 : if (ierr /= 0) return
707 6 : ii(1) = i
708 6 : jj(1) = j
709 6 : call fvbicub(ict, 1, 1, fval, ii, jj, xparam, yparam, hx, hxi, hy, hyi, f1, inf2, ny)
710 :
711 : end subroutine interp_evbicub_sg
712 :
713 : ! ***********************************************************************
714 : ! ***********************************************************************
715 :
716 : ! 2d piecewise monotonic interpolation -- values only, no slopes.
717 : ! does 4 1d interpolations in x followed by 1 1d interpolation in y.
718 : ! use interp_mkbipm_db to set up the interpolation information
719 : ! use interp_evbipm_db to evaluate it
720 :
721 1 : subroutine interp_mkbipm_sg(x, nx, y, ny, f1, nf2, ierr)
722 : use bipm_sg, only: do_mkbipm_sg
723 : integer, intent(in) :: nx ! length of x vector
724 : integer, intent(in) :: ny ! length of y vector
725 : real, intent(in), pointer :: x(:) ! (nx) ! x vector, strict ascending
726 : real, intent(in), pointer :: y(:) ! (ny) ! y vector, strict ascending
727 : integer, intent(in) :: nf2 ! 2nd dimension of f, nf2.ge.nx
728 : real, intent(inout), pointer :: f1(:) ! =(4,nf2,ny) ! data & interpolant coefficients
729 : integer, intent(out) :: ierr ! =0 on exit if there is no error.
730 1 : call do_mkbipm_sg(x, nx, y, ny, f1, nf2, ierr)
731 1 : end subroutine interp_mkbipm_sg
732 :
733 6 : subroutine interp_evbipm_sg(xget, yget, x, nx, y, ny, f1, nf2, z, ierr)
734 : use bipm_sg, only: do_evbipm_sg
735 : integer, intent(in) :: nx, ny
736 : real, intent(in) :: xget, yget ! target of this interpolation
737 : real, intent(in), pointer :: x(:) ! (nx) ! ordered x grid
738 : real, intent(in), pointer :: y(:) ! (ny) ! ordered y grid
739 : integer, intent(in) :: nf2
740 : real, intent(in), pointer :: f1(:) ! =(4,nf2,ny) ! function data
741 : real, intent(out) :: z
742 : integer, intent(out) :: ierr ! error code =0 ==> no error
743 6 : call do_evbipm_sg(xget, yget, x, nx, y, ny, f1, nf2, z, ierr)
744 6 : end subroutine interp_evbipm_sg
745 :
746 : end module interp_2d_lib_sg
|
