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_1d_lib
21 : use const_def, only: dp
22 : use auto_diff
23 :
24 : implicit none
25 :
26 : contains
27 :
28 : ! this routine is a simply wrapper for making an interpolant and then using it.
29 46 : subroutine interpolate_vector( &
30 46 : n_old, x_old, n_new, x_new, v_old, v_new, interp_vec, nwork, work1, str, ierr)
31 : integer, intent(in) :: n_old, n_new
32 : real(dp), intent(in) :: x_old(:) !(n_old)
33 : real(dp), intent(in) :: v_old(:) !(n_old)
34 : real(dp), intent(in) :: x_new(:) !(n_new)
35 : real(dp), intent(inout) :: v_new(:) ! (n_new)
36 : interface
37 : subroutine interp_vec(x, nx, f1, nwork, work1, str, ierr) ! make cubic interpolant
38 : ! e.g., interp_pm, interp_m3a, interp_m3b, or interp_m3q
39 : use const_def, only: dp
40 : implicit none
41 : integer, intent(in) :: nx ! length of x vector
42 : real(dp), intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
43 : real(dp), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
44 : integer, intent(in) :: nwork
45 : real(dp), intent(inout), pointer :: work1(:) ! =(n_old, nwork)
46 : character (len=*) :: str ! for debugging
47 : integer, intent(out) :: ierr
48 : end subroutine interp_vec
49 : end interface
50 : integer, intent(in) :: nwork
51 : real(dp), intent(inout), pointer :: work1(:) ! =(n_old, nwork)
52 : character (len=*) :: str ! for debugging
53 : integer, intent(out) :: ierr
54 46 : real(dp), pointer :: f1(:), f(:,:)
55 : integer :: k
56 46 : ierr = 0
57 46 : allocate(f1(4*n_old), stat=ierr)
58 46 : if (ierr /= 0) return
59 46 : f(1:4,1:n_old) => f1(1:4*n_old)
60 48484 : do k=1,n_old
61 48484 : f(1,k) = v_old(k)
62 : end do
63 46 : call interp_vec(x_old, n_old, f1, nwork, work1, str, ierr) ! make interpolant
64 46 : if (ierr /= 0) then
65 0 : deallocate(f1)
66 0 : return
67 : end if
68 46 : call interp_values(x_old, n_old, f1, n_new, x_new, v_new, ierr)
69 46 : deallocate(f1)
70 46 : end subroutine interpolate_vector
71 :
72 : ! autodiff version of above
73 0 : subroutine interpolate_vector_autodiff( &
74 0 : n_old, x_old, n_new, x_new, v_old, v_new, interp_vec_autodiff, nwork, work1, str, ierr)
75 : integer, intent(in) :: n_old, n_new
76 : type(auto_diff_real_2var_order1), intent(in) :: x_old(:) !(n_old)
77 : type(auto_diff_real_2var_order1), intent(in) :: v_old(:) !(n_old)
78 : type(auto_diff_real_2var_order1), intent(in) :: x_new(:) !(n_new)
79 : type(auto_diff_real_2var_order1), intent(inout) :: v_new(:) ! (n_new)
80 : interface
81 : subroutine interp_vec_autodiff(x, nx, f1, nwork, work1, str, ierr) ! make cubic interpolant
82 : ! e.g., interp_pm, interp_m3a, interp_m3b, or interp_m3q
83 : use const_def, only: dp
84 : use auto_diff
85 : implicit none
86 : integer, intent(in) :: nx ! length of x vector
87 : type(auto_diff_real_2var_order1), intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
88 : type(auto_diff_real_2var_order1), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
89 : integer, intent(in) :: nwork
90 : type(auto_diff_real_2var_order1), intent(inout), pointer :: work1(:) ! =(n_old, nwork)
91 : character (len=*) :: str ! for debugging
92 : integer, intent(out) :: ierr
93 : end subroutine interp_vec_autodiff
94 : end interface
95 : integer, intent(in) :: nwork
96 : type(auto_diff_real_2var_order1), intent(inout), pointer :: work1(:) ! =(n_old, nwork)
97 : character (len=*) :: str ! for debugging
98 : integer, intent(out) :: ierr
99 0 : type(auto_diff_real_2var_order1), pointer :: f1(:), f(:,:)
100 : integer :: k
101 0 : ierr = 0
102 0 : allocate(f1(4*n_old), stat=ierr)
103 0 : if (ierr /= 0) return
104 0 : f(1:4,1:n_old) => f1(1:4*n_old)
105 0 : do k=1,n_old
106 0 : f(1,k) = v_old(k)
107 : end do
108 0 : call interp_vec_autodiff(x_old, n_old, f1, nwork, work1, str, ierr) ! make interpolant
109 0 : if (ierr /= 0) then
110 0 : deallocate(f1)
111 0 : return
112 : end if
113 0 : call interp_values_autodiff(x_old, n_old, f1, n_new, x_new, v_new, ierr)
114 0 : deallocate(f1)
115 0 : end subroutine interpolate_vector_autodiff
116 :
117 : ! this routine is a simply wrapper for making an interpolant with interp_pm and then using it.
118 1 : subroutine interpolate_vector_pm( &
119 1 : n_old, x_old, n_new, x_new, v_old, v_new, work1, str, ierr)
120 : use interp_1d_def, only: pm_work_size
121 : integer, intent(in) :: n_old, n_new
122 : real(dp), intent(in) :: x_old(:) !(n_old)
123 : real(dp), intent(in) :: v_old(:) !(n_old)
124 : real(dp), intent(in) :: x_new(:) !(n_new)
125 : real(dp), intent(inout) :: v_new(:) ! (n_new)
126 : real(dp), intent(inout), pointer :: work1(:) ! =(n_old, pm_work_size)
127 : character (len=*) :: str ! for debugging
128 : integer, intent(out) :: ierr
129 1 : real(dp), pointer :: f1(:), f(:,:)
130 : integer :: k
131 1 : ierr = 0
132 1 : allocate(f1(4*n_old), stat=ierr)
133 1 : if (ierr /= 0) return
134 1 : f(1:4,1:n_old) => f1(1:4*n_old)
135 4 : do k=1,n_old
136 4 : f(1,k) = v_old(k)
137 : end do
138 1 : call interp_pm(x_old, n_old, f1, pm_work_size, work1, str, ierr) ! make interpolant
139 1 : if (ierr /= 0) then
140 0 : deallocate(f1)
141 0 : return
142 : end if
143 1 : call interp_values(x_old, n_old, f1, n_new, x_new, v_new, ierr)
144 1 : deallocate(f1)
145 1 : end subroutine interpolate_vector_pm
146 :
147 :
148 0 : subroutine interp_4_to_1(pdqm1, pdq00, pdqp1, ndq00, pfm1, pf00, pfp1, pfp2, nf00, str, ierr)
149 : ! 4 points in, 1 point out
150 : ! piecewise monotonic cubic interpolation
151 : use interp_1d_def, only: pm_work_size
152 : real(dp), intent(in) :: pdqm1, pdq00, pdqp1 ! spacing between input points
153 : real(dp), intent(in) :: ndq00
154 : real(dp), intent(in) :: pfm1, pf00, pfp1, pfp2 ! previous values
155 : real(dp), intent(out) :: nf00 ! new value at nk
156 : character (len=*) :: str ! for debugging
157 : integer, intent(out) :: ierr
158 : integer, parameter :: n_old=4, n_new=1
159 0 : real(dp) :: x_old(n_old), v_old(n_old), x_new(n_new), v_new(n_new)
160 0 : real(dp), target :: work1_ary(n_old*pm_work_size)
161 : real(dp), pointer :: work1(:)
162 0 : work1 => work1_ary
163 : ierr = 0
164 0 : x_old(1) = 0d0
165 0 : x_old(2) = pdqm1
166 0 : x_old(3) = pdqm1+pdq00
167 0 : x_old(4) = pdqm1+pdq00+pdqp1
168 0 : v_old(1) = pfm1
169 0 : v_old(2) = pf00
170 0 : v_old(3) = pfp1
171 0 : v_old(4) = pfp2
172 0 : x_new(1) = ndq00
173 : call interpolate_vector_pm( &
174 0 : n_old, x_old, n_new, x_new, v_old, v_new, work1, str, ierr)
175 0 : nf00 = v_new(1)
176 0 : end subroutine interp_4_to_1
177 :
178 :
179 0 : subroutine interp_3_to_1(pdqm1, pdq00, ndqm1, pfm1, pf00, pfp1, nf00, str, ierr)
180 : ! 3 points in, 1 point out
181 : ! piecewise monotonic quadratic interpolation
182 : use interp_1d_def, only: pm_work_size
183 : real(dp), intent(in) :: pdqm1, pdq00 ! spacing between input points
184 : real(dp), intent(in) :: ndqm1 ! new spacing to nk
185 : real(dp), intent(in) :: pfm1, pf00, pfp1 ! previous values at pkm1, pk, pkp1
186 : real(dp), intent(out) :: nf00 ! new value at nk
187 : character (len=*) :: str ! for debugging
188 : integer, intent(out) :: ierr
189 : integer, parameter :: n_old=3, n_new=1
190 0 : real(dp) :: x_old(n_old), v_old(n_old), x_new(n_new), v_new(n_new)
191 0 : real(dp), target :: work1_ary(n_old*pm_work_size)
192 : real(dp), pointer :: work1(:)
193 0 : work1 => work1_ary
194 : ierr = 0
195 0 : x_old(1) = 0d0
196 0 : x_old(2) = pdqm1
197 0 : x_old(3) = pdqm1+pdq00
198 0 : v_old(1) = pfm1
199 0 : v_old(2) = pf00
200 0 : v_old(3) = pfp1
201 0 : x_new(1) = ndqm1
202 : call interpolate_vector_pm( &
203 0 : n_old, x_old, n_new, x_new, v_old, v_new, work1, str, ierr)
204 0 : nf00 = v_new(1)
205 0 : end subroutine interp_3_to_1
206 :
207 :
208 1 : subroutine interp_3_to_2(pdqm1, pdq00, ndqm1, ndq00, pfm1, pf00, pfp1, nf00, nfp1, str, ierr)
209 : ! 3 points in, 2 points out
210 : ! piecewise monotonic quadratic interpolation
211 : use interp_1d_def, only: pm_work_size
212 : real(dp), intent(in) :: pdqm1, pdq00 ! previous spacing to pk and pkp1
213 : real(dp), intent(in) :: ndqm1, ndq00 ! new spacing to nk and nk+1
214 : real(dp), intent(in) :: pfm1, pf00, pfp1 ! previous values at pkm1, pk, pkp1
215 : real(dp), intent(out) :: nf00, nfp1 ! new values at nk, nk+1
216 : character (len=*) :: str ! for debugging
217 : integer, intent(out) :: ierr
218 : integer, parameter :: n_old=3, n_new=2
219 14 : real(dp) :: x_old(n_old), v_old(n_old), x_new(n_new), v_new(n_new)
220 10 : real(dp), target :: work1_ary(n_old*pm_work_size)
221 : real(dp), pointer :: work1(:)
222 1 : work1 => work1_ary
223 : ierr = 0
224 1 : x_old(1) = 0d0
225 1 : x_old(2) = pdqm1
226 1 : x_old(3) = pdqm1+pdq00
227 1 : v_old(1) = pfm1
228 1 : v_old(2) = pf00
229 1 : v_old(3) = pfp1
230 1 : x_new(1) = ndqm1
231 1 : x_new(2) = ndqm1+ndq00
232 : call interpolate_vector_pm( &
233 1 : n_old, x_old, n_new, x_new, v_old, v_new, work1, str, ierr)
234 1 : nf00 = v_new(1)
235 1 : nfp1 = v_new(2)
236 1 : end subroutine interp_3_to_2
237 :
238 :
239 : ! general routines
240 :
241 : ! these routines use previously created interpolant information (f)
242 : ! the interpolant can come from either the piecewise monotonic routines, or
243 : ! from the monotonicity preserving routines -- they use the same format for f.
244 :
245 105 : subroutine interp_values(init_x, nx, f1, nv, x, vals, ierr)
246 : use interp_1d_def
247 : use interp_1d_misc
248 : real(dp), intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
249 : integer, intent(in) :: nx ! length of init_x vector
250 : real(dp), intent(in), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
251 : integer, intent(in) :: nv ! length of new x vector and vals vector
252 : real(dp), intent(in) :: x(:) ! (nv) ! locations where want interpolated values
253 : ! strictly monotonic in same way as init_x
254 : ! values out of range of init_x's are clipped to boundaries of init_x's
255 : real(dp), intent(inout) :: vals(:) ! (nv)
256 : integer, intent(out) :: ierr ! 0 means AOK
257 105 : call do_interp_values(init_x, nx, f1, nv, x, vals, ierr)
258 105 : end subroutine interp_values
259 :
260 0 : subroutine interp_values_autodiff(init_x, nx, f1, nv, x, vals, ierr)
261 105 : use interp_1d_def
262 : use interp_1d_misc
263 : type(auto_diff_real_2var_order1), intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
264 : integer, intent(in) :: nx ! length of init_x vector
265 : type(auto_diff_real_2var_order1), intent(in), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
266 : integer, intent(in) :: nv ! length of new x vector and vals vector
267 : type(auto_diff_real_2var_order1), intent(in) :: x(:) ! (nv) ! locations where want interpolated values
268 : ! strictly monotonic in same way as init_x
269 : ! values out of range of init_x's are clipped to boundaries of init_x's
270 : type(auto_diff_real_2var_order1), intent(inout) :: vals(:) ! (nv)
271 : integer, intent(out) :: ierr ! 0 means AOK
272 0 : call do_interp_values_autodiff(init_x, nx, f1, nv, x, vals, ierr)
273 0 : end subroutine interp_values_autodiff
274 :
275 :
276 6232 : subroutine interp_value(init_x, nx, f1, xval, val, ierr)
277 0 : use interp_1d_def
278 : use interp_1d_misc
279 : real(dp), intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
280 : integer, intent(in) :: nx ! length of init_x vector
281 : real(dp), intent(in), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
282 : real(dp), intent(in) :: xval ! location where want interpolated value
283 : real(dp), intent(out) :: val
284 : integer, intent(out) :: ierr ! 0 means AOK
285 : integer, parameter :: nv = 1
286 12464 : real(dp) :: x(nv), vals(nv)
287 6232 : x(1) = xval
288 6232 : call do_interp_values(init_x, nx, f1, nv, x, vals, ierr)
289 6232 : val = vals(1)
290 6232 : end subroutine interp_value
291 :
292 :
293 0 : subroutine interp_values_and_slopes(init_x, nx, f1, nv, x, vals, slopes, ierr)
294 6232 : use interp_1d_def
295 : use interp_1d_misc
296 : real(dp), intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
297 : integer, intent(in) :: nx ! length of init_x vector
298 : real(dp), intent(in), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
299 : integer, intent(in) :: nv ! length of new x vector and vals and slopes vectors
300 : real(dp), intent(in) :: x(:) ! (nv) ! locations where want interpolated values
301 : ! strictly monotonic in same way as init_x
302 : ! values out of range of init_x's are clipped to boundaries of init_x's
303 : real(dp), intent(inout) :: vals(:) ! (nv)
304 : real(dp), intent(inout) :: slopes(:) ! (nv)
305 : integer, intent(out) :: ierr ! 0 means AOK
306 0 : call do_interp_values_and_slopes(init_x, nx, f1, nv, x, vals, slopes, ierr)
307 0 : end subroutine interp_values_and_slopes
308 :
309 :
310 16 : subroutine interp_value_and_slope(init_x, nx, f1, xval, val, slope, ierr)
311 0 : use interp_1d_def
312 : use interp_1d_misc
313 : real(dp), intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
314 : integer, intent(in) :: nx ! length of init_x vector
315 : real(dp), intent(in), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
316 : real(dp), intent(in) :: xval ! location where want interpolated values
317 : real(dp), intent(out) :: val, slope
318 : integer, intent(out) :: ierr ! 0 means AOK
319 : integer, parameter :: nv = 1
320 64 : real(dp) :: x(nv), vals(nv), slopes(nv)
321 16 : x(1) = xval
322 16 : call do_interp_values_and_slopes(init_x, nx, f1, nv, x, vals, slopes, ierr)
323 16 : val = vals(1)
324 16 : slope = slopes(1)
325 16 : end subroutine interp_value_and_slope
326 :
327 :
328 0 : subroutine interp2_values_and_slopes( &
329 0 : init_x, nx, f1_1, f1_2, nv, x, vals_1, slopes_1, vals_2, slopes_2, ierr)
330 16 : use interp_1d_def
331 : use interp_1d_misc
332 : real(dp), intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
333 : integer, intent(in) :: nx ! length of init_x vector
334 : real(dp), intent(in), pointer :: f1_1(:), f1_2(:) ! =(4,nx) ! data & interpolation coefficients
335 : integer, intent(in) :: nv ! length of new x vector and vals and slopes vectors
336 : real(dp), intent(in) :: x(:) ! (nv) ! locations where want interpolated values
337 : ! strictly monotonic in same way as init_x
338 : ! values out of range of init_x's are clipped to boundaries of init_x's
339 : real(dp), intent(inout) :: vals_1(:), vals_2(:) ! (nv)
340 : real(dp), intent(inout) :: slopes_1(:), slopes_2(:) ! (nv)
341 : integer, intent(out) :: ierr ! 0 means AOK
342 : call do_interp2_values_and_slopes( &
343 0 : init_x, nx, f1_1, f1_2, nv, x, vals_1, slopes_1, vals_2, slopes_2, ierr)
344 0 : end subroutine interp2_values_and_slopes
345 :
346 :
347 0 : subroutine interp2_value_and_slope(init_x, nx, f1_1, f1_2, xval, val_1, slope_1, val_2, slope_2, ierr)
348 0 : use interp_1d_def
349 : use interp_1d_misc
350 : real(dp), intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
351 : integer, intent(in) :: nx ! length of init_x vector
352 : real(dp), intent(in), pointer :: f1_1(:), f1_2(:) ! =(4,nx) ! data & interpolation coefficients
353 : real(dp), intent(in) :: xval ! location where want interpolated values
354 : real(dp), intent(out) :: val_1, slope_1, val_2, slope_2
355 : integer, intent(out) :: ierr ! 0 means AOK
356 : integer, parameter :: nv = 1
357 0 : real(dp) :: x(nv), vals_1(nv), slopes_1(nv), vals_2(nv), slopes_2(nv)
358 0 : x(1) = xval
359 : call do_interp2_values_and_slopes( &
360 0 : init_x, nx, f1_1, f1_2, nv, x, vals_1, slopes_1, vals_2, slopes_2, ierr)
361 0 : val_1 = vals_1(1)
362 0 : slope_1 = slopes_1(1)
363 0 : val_2 = vals_2(1)
364 0 : slope_2 = slopes_2(1)
365 0 : end subroutine interp2_value_and_slope
366 :
367 :
368 0 : subroutine interp3_values_and_slopes( &
369 0 : init_x, nx, f1_1, f1_2, f1_3, nv, x, &
370 0 : vals_1, slopes_1, vals_2, slopes_2, vals_3, slopes_3, ierr)
371 0 : use interp_1d_def
372 : use interp_1d_misc
373 : real(dp), intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
374 : integer, intent(in) :: nx ! length of init_x vector
375 : real(dp), intent(in), pointer :: f1_1(:), f1_2(:), f1_3(:) ! =(4,nx) ! data & interpolation coefficients
376 : integer, intent(in) :: nv ! length of new x vector and vals and slopes vectors
377 : real(dp), intent(in) :: x(:) ! (nv) ! locations where want interpolated values
378 : ! strictly monotonic in same way as init_x
379 : ! values out of range of init_x's are clipped to boundaries of init_x's
380 : real(dp), intent(inout) :: vals_1(:), vals_2(:), vals_3(:) ! (nv)
381 : real(dp), intent(inout) :: slopes_1(:), slopes_2(:), slopes_3(:) ! (nv)
382 : integer, intent(out) :: ierr ! 0 means AOK
383 : call do_interp3_values_and_slopes( &
384 : init_x, nx, f1_1, f1_2, f1_3, nv, x, &
385 0 : vals_1, slopes_1, vals_2, slopes_2, vals_3, slopes_3, ierr)
386 0 : end subroutine interp3_values_and_slopes
387 :
388 :
389 0 : subroutine interp3_value_and_slope( &
390 0 : init_x, nx, f1_1, f1_2, f1_3, xval, &
391 : val_1, slope_1, val_2, slope_2, val_3, slope_3, ierr)
392 0 : use interp_1d_def
393 : use interp_1d_misc
394 : real(dp), intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
395 : integer, intent(in) :: nx ! length of init_x vector
396 : real(dp), intent(in), pointer :: f1_1(:), f1_2(:), f1_3(:) ! =(4,nx) ! data & interpolation coefficients
397 : real(dp), intent(in) :: xval ! location where want interpolated values
398 : real(dp), intent(out) :: val_1, slope_1, val_2, slope_2, val_3, slope_3
399 : integer, intent(out) :: ierr ! 0 means AOK
400 : integer, parameter :: nv = 1
401 0 : real(dp) :: x(nv), vals_1(nv), slopes_1(nv), vals_2(nv), slopes_2(nv), vals_3(nv), slopes_3(nv)
402 0 : x(1) = xval
403 : call do_interp3_values_and_slopes( &
404 : init_x, nx, f1_1, f1_2, f1_3, nv, x, &
405 0 : vals_1, slopes_1, vals_2, slopes_2, vals_3, slopes_3, ierr)
406 0 : val_1 = vals_1(1)
407 0 : slope_1 = slopes_1(1)
408 0 : val_2 = vals_2(1)
409 0 : slope_2 = slopes_2(1)
410 0 : val_3 = vals_3(1)
411 0 : slope_3 = slopes_3(1)
412 0 : end subroutine interp3_value_and_slope
413 :
414 :
415 0 : subroutine interp6_values_and_slopes( &
416 0 : init_x, nx, f1_1, f1_2, f1_3, f1_4, f1_5, f1_6, nv, x, &
417 0 : vals_1, slopes_1, vals_2, slopes_2, vals_3, slopes_3, &
418 0 : vals_4, slopes_4, vals_5, slopes_5, vals_6, slopes_6, &
419 : ierr)
420 0 : use interp_1d_misc
421 : real(dp), intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
422 : integer, intent(in) :: nx ! length of init_x vector
423 : real(dp), intent(in), pointer, dimension(:) :: f1_1, f1_2, f1_3, f1_4, f1_5, f1_6
424 : integer, intent(in) :: nv ! length of new x vector and vals vector
425 : real(dp), intent(in) :: x(:) ! (nv) ! locations where want interpolated values
426 : ! strictly monotonic in same way as init_x
427 : ! values out of range of init_x's are clipped to boundaries of init_x's
428 : real(dp), intent(inout), dimension(:) :: &
429 : vals_1, slopes_1, vals_2, slopes_2, vals_3, slopes_3, &
430 : vals_4, slopes_4, vals_5, slopes_5, vals_6, slopes_6
431 : integer, intent(out) :: ierr ! 0 means aok
432 : ierr = 0
433 : call do_interp6_values_and_slopes( &
434 : init_x, nx, f1_1, f1_2, f1_3, f1_4, f1_5, f1_6, nv, x, &
435 : vals_1, slopes_1, vals_2, slopes_2, vals_3, slopes_3, &
436 : vals_4, slopes_4, vals_5, slopes_5, vals_6, slopes_6, &
437 0 : ierr)
438 0 : end subroutine interp6_values_and_slopes
439 :
440 :
441 0 : subroutine interp6_value_and_slope( &
442 0 : init_x, nx, f1_1, f1_2, f1_3, f1_4, f1_5, f1_6, xval, &
443 : val_1, slope_1, val_2, slope_2, val_3, slope_3, &
444 : val_4, slope_4, val_5, slope_5, val_6, slope_6, &
445 : ierr)
446 0 : use interp_1d_misc
447 : real(dp), intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
448 : integer, intent(in) :: nx ! length of init_x vector
449 : real(dp), intent(in), pointer, dimension(:) :: f1_1, f1_2, f1_3, f1_4, f1_5, f1_6
450 : real(dp), intent(in) :: xval ! location where want interpolated values
451 : real(dp), intent(out) :: val_1, slope_1, val_2, slope_2, val_3, slope_3, &
452 : val_4, slope_4, val_5, slope_5, val_6, slope_6
453 : integer, intent(out) :: ierr ! 0 means AOK
454 :
455 : integer, parameter :: nv = 1
456 0 : real(dp), dimension(nv) :: x, &
457 0 : vals_1, slopes_1, vals_2, slopes_2, vals_3, slopes_3, &
458 0 : vals_4, slopes_4, vals_5, slopes_5, vals_6, slopes_6
459 : ierr = 0
460 0 : x(1) = xval
461 : call do_interp6_values_and_slopes( &
462 : init_x, nx, f1_1, f1_2, f1_3, f1_4, f1_5, f1_6, nv, x, &
463 : vals_1, slopes_1, vals_2, slopes_2, vals_3, slopes_3, &
464 : vals_4, slopes_4, vals_5, slopes_5, vals_6, slopes_6, &
465 0 : ierr)
466 0 : val_1 = vals_1(1); slope_1 = slopes_1(1)
467 0 : val_2 = vals_2(1); slope_2 = slopes_2(1)
468 0 : val_3 = vals_3(1); slope_3 = slopes_3(1)
469 0 : val_4 = vals_4(1); slope_4 = slopes_4(1)
470 0 : val_5 = vals_5(1); slope_5 = slopes_5(1)
471 0 : val_6 = vals_6(1); slope_6 = slopes_6(1)
472 0 : end subroutine interp6_value_and_slope
473 :
474 :
475 2 : subroutine integrate_values(init_x, nx, f1, nv, x, vals, ierr)
476 0 : use interp_1d_def
477 : use interp_1d_misc
478 : real(dp), intent(in) :: init_x(:) ! (nx) ! junction points, strictly increasing
479 : integer, intent(in) :: nx ! length of init_x vector
480 : real(dp), intent(in), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
481 : integer, intent(in) :: nv ! length of new x vector and vals vector
482 : real(dp), intent(in) :: x(:) ! (nv)
483 : ! strictly monotonic in same way as init_x
484 : ! NOTE: no extrapolation allowed -- x's must be within range of init_x's
485 : real(dp), intent(inout) :: vals(:) ! (nv)
486 : ! for i > 1, vals(i) = integral of interpolating poly from x(i-1) to x(i)
487 : ! vals(1) = 0
488 : integer, intent(out) :: ierr ! 0 means AOK
489 :
490 2 : call do_integrate_values(init_x, nx, f1, nv, x, vals, ierr)
491 :
492 2 : end subroutine integrate_values
493 :
494 :
495 : ! piecewise monotonic routines
496 :
497 : ! the following produce piecewise monotonic interpolants rather than monotonicity preserving
498 : ! this stricter limit never introduces interpolated values exceeding the given values,
499 : ! even in places where the given values are not monotonic.
500 : ! the downside is reduced accuracy on smooth data compared to the mp routines.
501 :
502 :
503 : ! Steffen, M., "A simple method for monotonic interpolation in one dimension",
504 : ! Astron. Astrophys., (239) 1990, 443-450.
505 :
506 :
507 76990 : subroutine interp_pm(x, nx, f1, nwork, work1, str, ierr) ! make piecewise monotonic cubic interpolant
508 2 : use interp_1d_def
509 : use interp_1d_pm
510 : integer, intent(in) :: nx ! length of x vector (>= 2)
511 : real(dp), intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
512 : real(dp), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
513 : integer, intent(in) :: nwork ! nwork must be >= pm_work_size (see interp_1d_def)
514 : real(dp), intent(inout), pointer :: work1(:) ! =(nx, nwork)
515 : character (len=*) :: str ! for debugging
516 : integer, intent(out) :: ierr
517 76990 : call mk_pmcub(x, nx, f1, .false., nwork, work1, str, ierr)
518 76990 : end subroutine interp_pm
519 :
520 0 : subroutine interp_pm_autodiff(x, nx, f1, nwork, work1, str, ierr) ! make piecewise monotonic cubic interpolant
521 : use interp_1d_def
522 : use interp_1d_pm_autodiff
523 : integer, intent(in) :: nx ! length of x vector (>= 2)
524 : type(auto_diff_real_2var_order1), intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
525 : type(auto_diff_real_2var_order1), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
526 : integer, intent(in) :: nwork ! nwork must be >= pm_work_size (see interp_1d_def)
527 : type(auto_diff_real_2var_order1), intent(inout), pointer :: work1(:) ! =(nx, nwork)
528 : character (len=*) :: str ! for debugging
529 : integer, intent(out) :: ierr
530 0 : call mk_pmcub_autodiff(x, nx, f1, .false., nwork, work1, str, ierr)
531 0 : end subroutine interp_pm_autodiff
532 :
533 0 : subroutine interp_pm_slopes_only(x, nx, f1, nwork, work1, str, ierr)
534 : ! identical to interp_pm, but only calculates slopes and stores them in f(2,:)
535 : ! this is a little faster for the special case in which you just want the slopes at x
536 0 : use interp_1d_def
537 : use interp_1d_pm
538 : integer, intent(in) :: nx ! length of x vector (>= 2)
539 : real(dp), intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
540 : real(dp), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
541 : integer, intent(in) :: nwork ! nwork must be >= pm_work_size (see interp_1d_def)
542 : real(dp), intent(inout), pointer :: work1(:) ! =(nx, nwork)
543 : character (len=*) :: str ! for debugging
544 : integer, intent(out) :: ierr
545 0 : call mk_pmcub(x, nx, f1, .true., nwork, work1, str, ierr)
546 0 : end subroutine interp_pm_slopes_only
547 :
548 :
549 1 : subroutine interp_4pt_pm(x, y, a)
550 : ! returns coefficients for monotonic cubic interpolation from x(2) to x(3)
551 : real(dp), intent(in) :: x(4) ! junction points, strictly monotonic
552 : real(dp), intent(in) :: y(4) ! data values at x's
553 : real(dp), intent(inout) :: a(3) ! coefficients
554 1 : real(dp) :: h1, h2, h3, s1, s2, s3, p2, p3, as2, ss2, yp2, yp3
555 : ! for x(2) <= x <= x(3) and dx = x-x(2),
556 : ! y(x) = y(2) + dx*(a(1) + dx*(a(2) + dx*a(3)))
557 1 : h1 = x(2)-x(1)
558 1 : h2 = x(3)-x(2)
559 1 : h3 = x(4)-x(3)
560 1 : s1 = (y(2)-y(1))/h1
561 1 : s2 = (y(3)-y(2))/h2
562 1 : s3 = (y(4)-y(3))/h3
563 1 : p2 = (s1*h2+s2*h1)/(h1+h2)
564 1 : p3 = (s2*h3+s3*h2)/(h2+h3)
565 1 : as2 = abs(s2)
566 1 : ss2 = sign(1d0, s2)
567 1 : yp2 = (sign(1d0, s1)+ss2)*min(abs(s1), as2, 0.5d0*abs(p2))
568 1 : yp3 = (ss2+sign(1d0, s3))*min(as2, abs(s3), 0.5d0*abs(p3))
569 1 : a(1) = yp2
570 1 : a(2) = (3*s2-2*yp2-yp3)/h2
571 1 : a(3) = (yp2+yp3-2*s2)/(h2*h2)
572 1 : end subroutine interp_4pt_pm
573 :
574 :
575 0 : subroutine interp_pm_on_uniform_grid(dx, nx, f1, nwork, work1, str, ierr)
576 : ! make piecewise monotonic cubic interpolant on uniformly spaced mesh
577 : use interp_1d_def
578 : use interp_1d_pm
579 : real(dp), intent(in) :: dx ! grid spacing
580 : integer, intent(in) :: nx ! length of vector (>= 2)
581 : real(dp), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
582 : integer, intent(in) :: nwork ! nwork must be >= pm_work_size (see interp_1d_def)
583 : real(dp), intent(inout), pointer :: work1(:) ! =(nx, nwork)
584 : character (len=*) :: str ! for debugging
585 : integer, intent(out) :: ierr
586 0 : call mk_pmcub_uniform(dx, nx, f1, .false., nwork, work1, str, ierr)
587 0 : end subroutine interp_pm_on_uniform_grid
588 :
589 :
590 :
591 : ! monotonicity preserving routines
592 :
593 : ! Huynh, H.T., "Accurate Monotone Cubic Interpolation", SIAM J Numer. Anal. (30) 1993, 57-100.
594 :
595 : ! Suresh, A, and H.T. Huynh, "Accurate Monotonicity-Preserving Schemes with Runge-Kutta
596 : ! Time Stepping", JCP (136) 1997, 83-99.
597 :
598 :
599 0 : subroutine interp_m3(x, nx, f1, which, nwork, work1, str, ierr)
600 : ! make monotonicity preserving cubic interpolant on arbitrarily spaced grid
601 : use interp_1d_def
602 : use interp_1d_mp
603 : integer, intent(in) :: nx ! length of x vector (>= 4)
604 : real(dp), intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
605 : real(dp), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
606 : integer, intent(in) :: which ! average, quartic, or super_bee
607 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
608 : real(dp), intent(inout), pointer :: work1(:) ! =(nx, nwork)
609 : character (len=*) :: str ! for debugging
610 : integer, intent(out) :: ierr
611 0 : call m3(x, nx, f1, which, .false., nwork, work1, str, ierr)
612 0 : end subroutine interp_m3
613 :
614 :
615 0 : subroutine interp_m3a(x, nx, f1, nwork, work1, str, ierr)
616 : ! make monotonicity preserving cubic interpolant on arbitrarily spaced grid
617 : use interp_1d_def
618 : use interp_1d_mp
619 : integer, intent(in) :: nx ! length of x vector (>= 4)
620 : real(dp), intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
621 : real(dp), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
622 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
623 : real(dp), intent(inout), pointer :: work1(:) ! =(nx, nwork)
624 : character (len=*) :: str ! for debugging
625 : integer, intent(out) :: ierr
626 0 : call m3(x, nx, f1, average, .false., nwork, work1, str, ierr)
627 0 : end subroutine interp_m3a
628 :
629 :
630 1315 : subroutine interp_m3q(x, nx, f1, nwork, work1, str, ierr)
631 : ! make monotonicity preserving cubic interpolant on arbitrarily spaced grid
632 : use interp_1d_def
633 : use interp_1d_mp
634 : integer, intent(in) :: nx ! length of x vector (>= 4)
635 : real(dp), intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
636 : real(dp), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
637 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
638 : real(dp), intent(inout), pointer :: work1(:) ! =(nx, nwork)
639 : character (len=*) :: str ! for debugging
640 : integer, intent(out) :: ierr
641 1315 : call m3(x, nx, f1, quartic, .false., nwork, work1, str, ierr)
642 1315 : end subroutine interp_m3q
643 :
644 :
645 0 : subroutine interp_m3b(x, nx, f1, nwork, work1, str, ierr)
646 : ! make monotonicity preserving cubic interpolant on arbitrarily spaced grid
647 : use interp_1d_def
648 : use interp_1d_mp
649 : integer, intent(in) :: nx ! length of x vector (>= 4)
650 : real(dp), intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
651 : real(dp), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
652 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
653 : real(dp), intent(inout), pointer :: work1(:) ! =(nx, nwork)
654 : character (len=*) :: str ! for debugging
655 : integer, intent(out) :: ierr
656 0 : call m3(x, nx, f1, super_bee, .false., nwork, work1, str, ierr)
657 0 : end subroutine interp_m3b
658 :
659 :
660 0 : subroutine interp_m3_on_uniform_grid(dx, nx, f1, which, nwork, work1, str, ierr)
661 : ! make monotonicity preserving cubic interpolant on uniformly spaced grid
662 : use interp_1d_def
663 : use interp_1d_mp
664 : real(dp), intent(in) :: dx ! the grid spacing
665 : integer, intent(in) :: nx ! length of x vector (>= 4)
666 : real(dp), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
667 : integer, intent(in) :: which ! average, quartic, or super_bee
668 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
669 : real(dp), intent(inout), pointer :: work1(:) ! =(nx, nwork)
670 : character (len=*) :: str ! for debugging
671 : integer, intent(out) :: ierr
672 0 : call m3_on_uniform_grid(dx, nx, f1, which, .false., nwork, work1, str, ierr)
673 0 : end subroutine interp_m3_on_uniform_grid
674 :
675 :
676 0 : subroutine interp_m3a_on_uniform_grid(dx, nx, f1, nwork, work1, str, ierr)
677 : ! make monotonicity preserving cubic interpolant on uniformly spaced grid
678 : use interp_1d_def
679 : use interp_1d_mp
680 : real(dp), intent(in) :: dx ! the grid spacing
681 : integer, intent(in) :: nx ! length of x vector (>= 4)
682 : real(dp), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
683 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
684 : real(dp), intent(inout), pointer :: work1(:) ! =(nx, nwork)
685 : character (len=*) :: str ! for debugging
686 : integer, intent(out) :: ierr
687 0 : call m3_on_uniform_grid(dx, nx, f1, average, .false., nwork, work1, str, ierr)
688 0 : end subroutine interp_m3a_on_uniform_grid
689 :
690 :
691 0 : subroutine interp_m3b_on_uniform_grid(dx, nx, f1, nwork, work1, str, ierr)
692 : ! make monotonicity preserving cubic interpolant on uniformly spaced grid
693 : use interp_1d_def
694 : use interp_1d_mp
695 : real(dp), intent(in) :: dx ! the grid spacing
696 : integer, intent(in) :: nx ! length of x vector (>= 4)
697 : real(dp), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
698 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
699 : real(dp), intent(inout), pointer :: work1(:) ! =(nx, nwork)
700 : character (len=*) :: str ! for debugging
701 : integer, intent(out) :: ierr
702 0 : call m3_on_uniform_grid(dx, nx, f1, super_bee, .false., nwork, work1, str, ierr)
703 0 : end subroutine interp_m3b_on_uniform_grid
704 :
705 :
706 0 : subroutine interp_m3q_on_uniform_grid(dx, nx, f1, nwork, work1, str, ierr)
707 : ! make monotonicity preserving cubic interpolant on uniformly spaced grid
708 : use interp_1d_def
709 : use interp_1d_mp
710 : real(dp), intent(in) :: dx ! the grid spacing
711 : integer, intent(in) :: nx ! length of x vector (>= 4)
712 : real(dp), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
713 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
714 : real(dp), intent(inout), pointer :: work1(:) ! =(nx, nwork)
715 : character (len=*) :: str ! for debugging
716 : integer, intent(out) :: ierr
717 0 : call m3_on_uniform_grid(dx, nx, f1, quartic, .false., nwork, work1, str, ierr)
718 0 : end subroutine interp_m3q_on_uniform_grid
719 :
720 :
721 0 : subroutine interp_m3_autodiff(x, nx, f1, which, nwork, work1, str, ierr)
722 : ! make monotonicity preserving cubic interpolant on arbitrarily spaced grid
723 : use interp_1d_def
724 : use interp_1d_mp_autodiff
725 : integer, intent(in) :: nx ! length of x vector (>= 4)
726 : type(auto_diff_real_2var_order1), intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
727 : type(auto_diff_real_2var_order1), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
728 : integer, intent(in) :: which ! average, quartic, or super_bee
729 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
730 : type(auto_diff_real_2var_order1), intent(inout), pointer :: work1(:) ! =(nx, nwork)
731 : character (len=*) :: str ! for debugging
732 : integer, intent(out) :: ierr
733 0 : call m3_autodiff(x, nx, f1, which, .false., nwork, work1, str, ierr)
734 0 : end subroutine interp_m3_autodiff
735 :
736 :
737 0 : subroutine interp_m3a_autodiff(x, nx, f1, nwork, work1, str, ierr)
738 : ! make monotonicity preserving cubic interpolant on arbitrarily spaced grid
739 0 : use interp_1d_def
740 : use interp_1d_mp_autodiff
741 : integer, intent(in) :: nx ! length of x vector (>= 4)
742 : type(auto_diff_real_2var_order1), intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
743 : type(auto_diff_real_2var_order1), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
744 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
745 : type(auto_diff_real_2var_order1), intent(inout), pointer :: work1(:) ! =(nx, nwork)
746 : character (len=*) :: str ! for debugging
747 : integer, intent(out) :: ierr
748 0 : call m3_autodiff(x, nx, f1, average, .false., nwork, work1, str, ierr)
749 0 : end subroutine interp_m3a_autodiff
750 :
751 :
752 0 : subroutine interp_m3q_autodiff(x, nx, f1, nwork, work1, str, ierr)
753 : ! make monotonicity preserving cubic interpolant on arbitrarily spaced grid
754 0 : use interp_1d_def
755 : use interp_1d_mp_autodiff
756 : integer, intent(in) :: nx ! length of x vector (>= 4)
757 : type(auto_diff_real_2var_order1), intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
758 : type(auto_diff_real_2var_order1), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
759 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
760 : type(auto_diff_real_2var_order1), intent(inout), pointer :: work1(:) ! =(nx, nwork)
761 : character (len=*) :: str ! for debugging
762 : integer, intent(out) :: ierr
763 0 : call m3_autodiff(x, nx, f1, quartic, .false., nwork, work1, str, ierr)
764 0 : end subroutine interp_m3q_autodiff
765 :
766 :
767 0 : subroutine interp_m3b_autodiff(x, nx, f1, nwork, work1, str, ierr)
768 : ! make monotonicity preserving cubic interpolant on arbitrarily spaced grid
769 0 : use interp_1d_def
770 : use interp_1d_mp_autodiff
771 : integer, intent(in) :: nx ! length of x vector (>= 4)
772 : type(auto_diff_real_2var_order1), intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
773 : type(auto_diff_real_2var_order1), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
774 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
775 : type(auto_diff_real_2var_order1), intent(inout), pointer :: work1(:) ! =(nx, nwork)
776 : character (len=*) :: str ! for debugging
777 : integer, intent(out) :: ierr
778 0 : call m3_autodiff(x, nx, f1, super_bee, .false., nwork, work1, str, ierr)
779 0 : end subroutine interp_m3b_autodiff
780 :
781 :
782 0 : subroutine interp_m3_on_uniform_grid_autodiff(dx, nx, f1, which, nwork, work1, str, ierr)
783 : ! make monotonicity preserving cubic interpolant on uniformly spaced grid
784 0 : use interp_1d_def
785 : use interp_1d_mp_autodiff
786 : type(auto_diff_real_2var_order1), intent(in) :: dx ! the grid spacing
787 : integer, intent(in) :: nx ! length of x vector (>= 4)
788 : type(auto_diff_real_2var_order1), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
789 : integer, intent(in) :: which ! average, quartic, or super_bee
790 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
791 : type(auto_diff_real_2var_order1), intent(inout), pointer :: work1(:) ! =(nx, nwork)
792 : character (len=*) :: str ! for debugging
793 : integer, intent(out) :: ierr
794 0 : call m3_on_uniform_grid_autodiff(dx, nx, f1, which, .false., nwork, work1, str, ierr)
795 0 : end subroutine interp_m3_on_uniform_grid_autodiff
796 :
797 :
798 0 : subroutine interp_m3a_on_uniform_grid_autodiff(dx, nx, f1, nwork, work1, str, ierr)
799 : ! make monotonicity preserving cubic interpolant on uniformly spaced grid
800 0 : use interp_1d_def
801 : use interp_1d_mp_autodiff
802 : type(auto_diff_real_2var_order1), intent(in) :: dx ! the grid spacing
803 : integer, intent(in) :: nx ! length of x vector (>= 4)
804 : type(auto_diff_real_2var_order1), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
805 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
806 : type(auto_diff_real_2var_order1), intent(inout), pointer :: work1(:) ! =(nx, nwork)
807 : character (len=*) :: str ! for debugging
808 : integer, intent(out) :: ierr
809 0 : call m3_on_uniform_grid_autodiff(dx, nx, f1, average, .false., nwork, work1, str, ierr)
810 0 : end subroutine interp_m3a_on_uniform_grid_autodiff
811 :
812 :
813 0 : subroutine interp_m3b_on_uniform_grid_autodiff(dx, nx, f1, nwork, work1, str, ierr)
814 : ! make monotonicity preserving cubic interpolant on uniformly spaced grid
815 0 : use interp_1d_def
816 : use interp_1d_mp_autodiff
817 : type(auto_diff_real_2var_order1), intent(in) :: dx ! the grid spacing
818 : integer, intent(in) :: nx ! length of x vector (>= 4)
819 : type(auto_diff_real_2var_order1), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
820 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
821 : type(auto_diff_real_2var_order1), intent(inout), pointer :: work1(:) ! =(nx, nwork)
822 : character (len=*) :: str ! for debugging
823 : integer, intent(out) :: ierr
824 0 : call m3_on_uniform_grid_autodiff(dx, nx, f1, super_bee, .false., nwork, work1, str, ierr)
825 0 : end subroutine interp_m3b_on_uniform_grid_autodiff
826 :
827 :
828 0 : subroutine interp_m3q_on_uniform_grid_autodiff(dx, nx, f1, nwork, work1, str, ierr)
829 : ! make monotonicity preserving cubic interpolant on uniformly spaced grid
830 0 : use interp_1d_def
831 : use interp_1d_mp_autodiff
832 : type(auto_diff_real_2var_order1), intent(in) :: dx ! the grid spacing
833 : integer, intent(in) :: nx ! length of x vector (>= 4)
834 : type(auto_diff_real_2var_order1), intent(inout), pointer :: f1(:) ! =(4,nx) ! data & interpolation coefficients
835 : integer, intent(in) :: nwork ! nwork must be >= mp_work_size (see interp_1d_def)
836 : type(auto_diff_real_2var_order1), intent(inout), pointer :: work1(:) ! =(nx, nwork)
837 : character (len=*) :: str ! for debugging
838 : integer, intent(out) :: ierr
839 0 : call m3_on_uniform_grid_autodiff(dx, nx, f1, quartic, .false., nwork, work1, str, ierr)
840 0 : end subroutine interp_m3q_on_uniform_grid_autodiff
841 :
842 :
843 :
844 : ! simple Lagrange polynomial weights
845 :
846 0 : subroutine interp_weights_lp3(x, x1, x2, x3, c1, c2, c3)
847 :
848 : real(dp), intent(in) :: x, x1, x2, x3
849 : real(dp), intent(out) :: c1, c2, c3
850 :
851 0 : real(dp) :: dx1, dx2, dx3, dx12i, dx13i, dx23i
852 :
853 0 : dx1 = x-x1
854 0 : dx2 = x-x2
855 0 : dx3 = x-x3
856 :
857 0 : dx12i = 1d0/(x1-x2)
858 0 : dx13i = 1d0/(x1-x3)
859 0 : dx23i = 1d0/(x2-x3)
860 :
861 0 : c1 = dx2*dx3*dx12i*dx13i
862 0 : c2 = -dx1*dx3*dx12i*dx23i
863 0 : c3 = dx1*dx2*dx13i*dx23i
864 :
865 0 : end subroutine interp_weights_lp3
866 :
867 :
868 0 : subroutine interp_weights_lp4(x, x1, x2, x3, x4, c1, c2, c3, c4)
869 :
870 : real(dp), intent(in) :: x, x1, x2, x3, x4
871 : real(dp), intent(out) :: c1, c2, c3, c4
872 :
873 0 : real(dp) :: dx1, dx2, dx3, dx4, dx12i, dx13i, dx14i, dx23i, dx24i, dx34i
874 :
875 0 : dx1 = x-x1
876 0 : dx2 = x-x2
877 0 : dx3 = x-x3
878 0 : dx4 = x-x4
879 :
880 0 : dx12i = 1d0/(x1-x2)
881 0 : dx13i = 1d0/(x1-x3)
882 0 : dx14i = 1d0/(x1-x4)
883 0 : dx23i = 1d0/(x2-x3)
884 0 : dx24i = 1d0/(x2-x4)
885 0 : dx34i = 1d0/(x3-x4)
886 :
887 0 : c1 = dx2*dx3*dx4*dx12i*dx13i*dx14i
888 0 : c2 = -dx1*dx3*dx4*dx12i*dx23i*dx24i
889 0 : c3 = dx1*dx2*dx4*dx13i*dx23i*dx34i
890 0 : c4 = -dx1*dx2*dx3*dx14i*dx24i*dx34i
891 :
892 0 : end subroutine interp_weights_lp4
893 :
894 : end module interp_1d_lib
|
