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_db
21 : ! real(dp) library for 2d interpolation
22 :
23 : use const_def, only: dp
24 :
25 : implicit none
26 :
27 : contains ! the procedure interface for the library
28 : ! client programs should only call these routines.
29 :
30 : ! contents
31 :
32 : ! rectangular-grid of data points
33 : ! interp_rgbi3p_db -- point interpolation (akima)
34 : ! interp_rgsf3p_db -- surface interpolation (akima)
35 : ! interp_mkbicub_db -- bicubic splines
36 : ! interp_mkbipm_db -- 2d piecewise monotonic
37 :
38 : ! scattered set of data points
39 : ! interp_cs2val_db -- point interpolation (renka)
40 : ! interp_cs2grd_db -- point interpolation with gradients (renka)
41 :
42 : ! ***********************************************************************
43 : ! ***********************************************************************
44 : ! ***********************************************************************
45 : ! ***********************************************************************
46 :
47 : ! rectangular-grid bivariate interpolation and surface fitting
48 : ! from acm algorithm 760., acm trans. math. software (22) 1996, 357-361.
49 : ! hiroshi akima
50 : ! u.s. department of commerce, ntia/its
51 : ! version of 1995/08
52 :
53 : ! interpolate at given points
54 0 : subroutine interp_rgbi3p_db(md, nxd, nyd, xd, yd, zd, nip, xi, yi, zi, ierr, wk)
55 : integer, intent(in) :: md, nxd, nyd, nip
56 : real(dp), intent(in) :: xd(nxd), yd(nyd), zd(nxd, nyd), xi(nip), yi(nip)
57 : real(dp), intent(inout) :: zi(nip), wk(3, nxd, nyd)
58 : integer, intent(out) :: ierr
59 0 : call do_rgbi3p_db(md, nxd, nyd, xd, yd, zd, nip, xi, yi, zi, ierr, wk)
60 0 : end subroutine interp_rgbi3p_db
61 :
62 : ! interpolate surface for given grid in x and y
63 0 : subroutine interp_rgsf3p_db(md, nxd, nyd, xd, yd, zd, nxi, xi, nyi, yi, zi, ierr, wk)
64 : integer, intent(in) :: md, nxd, nyd, nxi, nyi
65 : real(dp), intent(in) :: xd(nxd), yd(nyd), zd(nxd, nyd), xi(nxi), yi(nyi)
66 : real(dp), intent(inout) :: zi(nxi, nyi), wk(3, nxd, nyd)
67 : integer, intent(out) :: ierr
68 0 : call do_rgsf3p_db(md, nxd, nyd, xd, yd, zd, nxi, xi, nyi, yi, zi, ierr, wk)
69 0 : end subroutine interp_rgsf3p_db
70 :
71 : ! ***********************************************************************
72 : ! ***********************************************************************
73 : ! ***********************************************************************
74 : ! ***********************************************************************
75 :
76 : ! cubic Shepard method for bivariate interpolation of scattered data.
77 : ! from acm algorithm 790., acm trans. math. software (25) 1999, 70-73.
78 : ! robert j. renka
79 : ! dept. of computer science
80 : ! univ. of north texas
81 : ! renka@cs.unt.edu
82 :
83 : ! use cshep2_db to set up the interpolation information.
84 : ! use cs2val_db to evaluate it.
85 : ! use cs2grd_db to get value and derivatives.
86 :
87 : ! detailed documentation can be found in interp_2d_lib_sg.f
88 : ! these routines are exactly the same except they are real(dp).
89 :
90 : ! see cshep2_sg for documentation details.
91 0 : subroutine interp_cshep2_db(n, x, y, f, nc, nw, nr, lcell, lnext, xmin, ymin, dx, dy, rmax, rw, a, ierr)
92 : integer, intent(in) :: n, nc, nw, nr
93 : integer, intent(out) :: lcell(nr, nr), lnext(n), ierr
94 : real(dp), intent(in) :: x(n), y(n), f(n)
95 : real(dp), intent(inout) :: xmin, ymin, dx, dy, rmax, rw(n), a(9, n)
96 0 : call do_cshep2_db(n, x, y, f, nc, nw, nr, lcell, lnext, xmin, ymin, dx, dy, rmax, rw, a, ierr)
97 0 : end subroutine interp_cshep2_db
98 :
99 : ! see cs2val_sg for documentation details.
100 0 : real(dp) function interp_cs2val_db(px, py, n, x, y, f, nr, lcell, lnext, xmin, ymin, dx, dy, rmax, rw, a, ierr)
101 : integer, intent(in) :: n, nr, lcell(nr, nr), lnext(n)
102 : real(dp), intent(in) :: px, py, x(n), y(n), f(n), xmin, ymin, dx, dy, rmax, rw(n), a(9, n)
103 : integer, intent(out) :: ierr
104 : real(dp) :: do_cs2val_db
105 0 : interp_cs2val_db = do_cs2val_db(px, py, n, x, y, f, nr, lcell, lnext, xmin, ymin, dx, dy, rmax, rw, a, ierr)
106 0 : end function interp_cs2val_db
107 :
108 : ! see cs2grd_sg for documentation details.
109 0 : subroutine interp_cs2grd_db(px, py, n, x, y, f, nr, lcell, lnext, xmin, ymin, dx, dy, rmax, rw, a, c, cx, cy, ierr)
110 : integer, intent(in) :: n, nr, lcell(nr, nr), lnext(n)
111 : real(dp), intent(in) :: px, py, x(n), y(n), f(n), xmin, ymin, dx, dy, rmax, rw(n), a(9, n)
112 : real(dp), intent(out) :: c, cx, cy
113 : integer, intent(out) :: ierr
114 0 : call do_cs2grd_db(px, py, n, x, y, f, nr, lcell, lnext, xmin, ymin, dx, dy, rmax, rw, a, c, cx, cy, ierr)
115 0 : end subroutine interp_cs2grd_db
116 :
117 : ! ***********************************************************************
118 : ! ***********************************************************************
119 :
120 : ! bicubic splines
121 : ! use interp_mkbicub_db to set up the interpolation information
122 : ! use interp_evbicub_db to evaluate it
123 :
124 : ! detailed documentation can be found in interp_2d_lib_sg.f
125 :
126 : ! see interp_mkbicub_sg for documentation details.
127 33168 : subroutine interp_mkbicub_db(x, nx, y, ny, f1, nf2, &
128 33168 : ibcxmin, bcxmin, ibcxmax, bcxmax, ibcymin, bcymin, ibcymax, bcymax, ilinx, iliny, ierr)
129 : use bicub_db, only: do_mkbicub_db
130 : integer, intent(in) :: nx ! length of x vector
131 : integer, intent(in) :: ny ! length of y vector
132 : real(dp), intent(in) :: x(:) ! (nx) ! x vector, strict ascending
133 : real(dp), intent(in) :: y(:) ! (ny) ! y vector, strict ascending
134 : integer, intent(in) :: nf2 ! 2nd dimension of f, nf2 >= nx
135 : real(dp), intent(inout), pointer :: f1(:) ! =(4,nf2,ny) ! data & spline coefficients
136 : integer, intent(in) :: ibcxmin ! bc flag for x=xmin
137 : real(dp), intent(in) :: bcxmin(:) ! (ny) ! bc data vs. y at x=xmin
138 : integer, intent(in) :: ibcxmax ! bc flag for x=xmax
139 : real(dp), intent(in) :: bcxmax(:) ! (ny) ! bc data vs. y at x=xmax
140 : integer, intent(in) :: ibcymin ! bc flag for y=ymin
141 : real(dp), intent(in) :: bcymin(:) ! (nx) ! bc data vs. x at y=ymin
142 : integer, intent(in) :: ibcymax ! bc flag for y=ymax
143 : real(dp), intent(in) :: bcymax(:) ! (nx) ! bc data vs. x at y=ymax
144 : integer, intent(out) :: ilinx ! =1: x grid is "nearly" equally spaced
145 : integer, intent(out) :: iliny ! =1: y grid is "nearly" equally spaced
146 : integer, intent(out) :: ierr ! =0 on exit if there is no error.
147 : call do_mkbicub_db(x, nx, y, ny, f1, nf2, &
148 33168 : ibcxmin, bcxmin, ibcxmax, bcxmax, ibcymin, bcymin, ibcymax, bcymax, ilinx, iliny, ierr)
149 33168 : end subroutine interp_mkbicub_db
150 :
151 : ! see interp_evbicub_sg for documentation details.
152 140 : subroutine interp_evbicub_db(xget, yget, x, nx, y, ny, ilinx, iliny, f1, inf2, ict, fval, ierr)
153 : use bicub_db, only: fvbicub_db, herm2xy_db
154 : integer, intent(in) :: nx, ny ! grid sizes
155 : real(dp), intent(in) :: xget, yget ! target of this interpolation
156 : real(dp), intent(in) :: x(:) ! (nx) ! ordered x grid
157 : real(dp), intent(in) :: y(:) ! (ny) ! ordered y grid
158 : integer, intent(in) :: ilinx ! ilinx=1 => assume x evenly spaced
159 : integer, intent(in) :: iliny ! iliny=1 => assume y evenly spaced
160 : integer, intent(in) :: inf2
161 : real(dp), intent(inout), pointer :: f1(:) ! function data
162 : integer, intent(in) :: ict(6) ! code specifying output desired
163 : real(dp), intent(inout) :: fval(6) ! output data
164 : integer, intent(out) :: ierr ! error code =0 ==> no error
165 : integer :: i, j
166 : integer :: ii(1), jj(1)
167 980 : real(dp) :: xparam(1), yparam(1), hx(1), hxi(1), hy(1), hyi(1)
168 140 : call herm2xy_db(xget, yget, x, nx, y, ny, ilinx, iliny, i, j, xparam(1), yparam(1), hx(1), hxi(1), hy(1), hyi(1), ierr)
169 140 : if (ierr /= 0) return
170 140 : ii(1) = i
171 140 : jj(1) = j
172 140 : call fvbicub_db(ict, 1, 1, fval, ii, jj, xparam, yparam, hx, hxi, hy, hyi, f1, inf2, ny)
173 : end subroutine interp_evbicub_db
174 :
175 : ! this is used by do_mkbicub_db to get 2nd derivatives d_dx2 and d_dy2
176 0 : subroutine interp_mkspline_db(x, nx, fspl, ibcxmin, bcxmin, ibcxmax, bcxmax, ilinx, ierr)
177 : use bicub_db, only: mkspline_db
178 : integer, intent(in) :: nx ! no. of data points
179 : real(dp), intent(in) :: x(nx) ! x axis data, strict ascending order
180 : real(dp), intent(inout) :: fspl(2, nx) ! f(1,*): data in; f(2,*): coeffs out
181 : integer, intent(in) :: ibcxmin ! b.c. flag @ x=xmin=x(1)
182 : real(dp), intent(in) :: bcxmin ! b.c. data @xmin
183 : integer, intent(in) :: ibcxmax ! b.c. flag @ x=xmax=x(nx)
184 : real(dp), intent(in) :: bcxmax ! b.c. data @xmax
185 : integer, intent(in) :: ilinx ! ilinx=1 => assume x evenly spaced
186 : integer, intent(out) :: ierr ! error code =0 ==> no error
187 0 : call mkspline_db(x, nx, fspl, ibcxmin, bcxmin, ibcxmax, bcxmax, ilinx, ierr)
188 0 : end subroutine interp_mkspline_db
189 :
190 : ! ***********************************************************************
191 : ! ***********************************************************************
192 :
193 : ! 2d piecewise monotonic interpolation -- values only, no slopes.
194 : ! does 4 1d interpolations in x followed by 1 1d interpolation in y.
195 : ! use interp_mkbipm_db to set up the interpolation information
196 : ! use interp_evbipm_db to evaluate it
197 :
198 1 : subroutine interp_mkbipm_db(x, nx, y, ny, f1, nf2, ierr)
199 : use bipm_db, only: do_mkbipm_db
200 : integer, intent(in) :: nx ! length of x vector
201 : integer, intent(in) :: ny ! length of y vector
202 : real(dp), intent(in), pointer :: x(:) ! (nx) ! x vector, strict ascending
203 : real(dp), intent(in), pointer :: y(:) ! (ny) ! y vector, strict ascending
204 : integer, intent(in) :: nf2 ! 2nd dimension of f, nf2 >= nx
205 : real(dp), intent(inout), pointer :: f1(:) ! =(4,nf2,ny) ! data & interpolant coefficients
206 : integer, intent(out) :: ierr ! =0 on exit if there is no error.
207 1 : call do_mkbipm_db(x, nx, y, ny, f1, nf2, ierr)
208 1 : end subroutine interp_mkbipm_db
209 :
210 6 : subroutine interp_evbipm_db(xget, yget, x, nx, y, ny, f1, nf2, z, ierr)
211 : use bipm_db, only: do_evbipm_db
212 : integer, intent(in) :: nx, ny
213 : real(dp), intent(in) :: xget, yget ! target of this interpolation
214 : real(dp), intent(in), pointer :: x(:) ! (nx) ! ordered x grid
215 : real(dp), intent(in), pointer :: y(:) ! (ny) ! ordered y grid
216 : integer, intent(in) :: nf2
217 : real(dp), intent(in), pointer :: f1(:) ! =(4,nf2,ny) ! function data
218 : real(dp), intent(out) :: z
219 : integer, intent(out) :: ierr ! error code =0 ==> no error
220 6 : call do_evbipm_db(xget, yget, x, nx, y, ny, f1, nf2, z, ierr)
221 6 : end subroutine interp_evbipm_db
222 :
223 : ! ***********************************************************************
224 : ! ***********************************************************************
225 :
226 : end module interp_2d_lib_db
|
