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_misc
21 :
22 : use const_def, only: dp
23 : use auto_diff
24 :
25 : implicit none
26 :
27 : contains
28 :
29 2 : subroutine do_integrate_values(init_x, nx, f1, nv, x, vals, ierr)
30 : real(dp), intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
31 : integer, intent(in) :: nx ! length of init_x vector
32 : real(dp), intent(in), pointer :: f1(:) ! =(4, nx) ! data & interpolation coefficients
33 : integer, intent(in) :: nv ! length of new x vector and vals vector
34 : real(dp), intent(in) :: x(:) ! (nv)
35 : ! strictly monotonic in same way as init_x
36 : real(dp), intent(inout) :: vals(:) ! (nv)
37 : ! for i > 1, vals(i) = integral of interpolating poly from x(i-1) to x(i)
38 : ! vals(1) = 0
39 : integer, intent(out) :: ierr ! 0 means aok
40 :
41 : integer :: k_old, k_new
42 2 : real(dp) :: xk_old, xkp1_old, xk_new, xk_prev, sum
43 : logical :: increasing
44 : real(dp), pointer :: f(:,:) ! (4, nx) ! data & interpolation coefficients
45 2 : f(1:4,1:nx) => f1(1:4*nx)
46 :
47 2 : increasing = (init_x(1) < init_x(nx))
48 :
49 : if (increasing .and. (x(1) < init_x(1) .or. x(nv) > init_x(nx)) &
50 2 : .or. ((.not. increasing) .and. (x(1) > init_x(1) .or. x(nv) < init_x(nx)))) then
51 0 : ierr = -1
52 0 : return
53 : end if
54 :
55 2 : ierr = 0
56 :
57 2 : k_old = 1; xk_old = init_x(k_old); xkp1_old = init_x(k_old+1)
58 2 : sum = 0; xk_prev = x(1); vals(1) = 0
59 :
60 14 : do k_new = 2, nv
61 :
62 12 : xk_new = x(k_new)
63 74 : do while ((increasing .and. xk_new > xkp1_old) .or. ((.not. increasing) .and. xk_new < xkp1_old))
64 62 : k_old = k_old + 1
65 62 : if (k_old >= nx) then
66 0 : k_old = k_old - 1
67 0 : xk_new = xkp1_old
68 0 : exit
69 : end if
70 62 : call add_to_integral(k_old - 1, xkp1_old)
71 62 : xk_old = xkp1_old
72 68 : xkp1_old = init_x(k_old+1)
73 : end do
74 :
75 12 : call add_to_integral(k_old, xk_new)
76 12 : vals(k_new) = sum
77 14 : sum = 0
78 :
79 : end do
80 :
81 : contains
82 :
83 74 : subroutine add_to_integral(k, x2)
84 : integer, intent(in) :: k
85 : real(dp), intent(in) :: x2
86 :
87 74 : real(dp) :: x0, x1, a1, a2, d1, d2, area
88 :
89 74 : x0 = init_x(k)
90 74 : x1 = xk_prev
91 74 : if (x1 == x2) return
92 74 : d2 = x2 - x0
93 : a2 = d2*(f(1, k) + d2*(f(2, k)/2 &
94 74 : + d2*(f(3, k)/3 + d2*f(4, k)/4)))
95 74 : if (x1 > x0) then
96 5 : d1 = x1 - x0
97 : a1 = d1*(f(1, k) + d1*(f(2, k)/2 &
98 5 : + d1*(f(3, k)/3 + d1*f(4, k)/4)))
99 5 : area = a2 - a1
100 : else
101 74 : d1 = 0; a1 = 0; area = a2
102 : end if
103 74 : sum = sum + area
104 74 : xk_prev = x2
105 :
106 : end subroutine add_to_integral
107 :
108 :
109 : end subroutine do_integrate_values
110 :
111 :
112 6337 : subroutine do_interp_values(init_x, nx, f1, nv, x, vals, ierr)
113 : real(dp), intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
114 : integer, intent(in) :: nx ! length of init_x vector
115 : real(dp), intent(in), pointer :: f1(:) ! =(4, nx) ! data & interpolation coefficients
116 : integer, intent(in) :: nv ! length of new x vector and vals vector
117 : real(dp), intent(in) :: x(:) ! (nv) ! locations where want interpolated values
118 : ! strictly monotonic in same way as init_x
119 : ! values out of range of init_x's are clipped to boundaries of init_x's
120 : real(dp), intent(inout) :: vals(:) ! (nv)
121 : integer, intent(out) :: ierr ! 0 means aok
122 :
123 : integer :: k_old, k_new
124 6337 : real(dp) :: xk_old, xkp1_old, xk_new, delta
125 : logical :: increasing
126 6337 : real(dp), pointer :: f(:,:) ! (4, nx) ! data & interpolation coefficients
127 :
128 6337 : ierr = 0
129 :
130 6337 : if (nx == 1) then
131 0 : vals(1:nv) = f1(1)
132 : return
133 : end if
134 :
135 6337 : f(1:4,1:nx) => f1(1:4*nx)
136 :
137 6337 : increasing = (init_x(1) < init_x(nx))
138 :
139 6337 : k_old = 1; xk_old = init_x(k_old); xkp1_old = init_x(k_old+1)
140 :
141 20557 : do k_new = 1, nv
142 :
143 14220 : xk_new = x(k_new)
144 14220 : if (increasing) then
145 14204 : if (xk_new > init_x(nx)) then
146 14220 : xk_new = init_x(nx)
147 14200 : else if (xk_new < init_x(1)) then
148 14220 : xk_new = init_x(1)
149 : end if
150 : else ! decreasing
151 16 : if (xk_new < init_x(nx)) then
152 : xk_new = init_x(nx)
153 16 : else if (xk_new > init_x(1)) then
154 14220 : xk_new = init_x(1)
155 : end if
156 : end if
157 67151 : do while ((increasing .and. xk_new > xkp1_old) .or. ((.not. increasing) .and. xk_new < xkp1_old))
158 52931 : k_old = k_old + 1
159 52931 : if (k_old >= nx) then
160 : k_old = k_old - 1
161 : xk_new = xkp1_old
162 : exit
163 : end if
164 52931 : xk_old = xkp1_old
165 52947 : xkp1_old = init_x(k_old+1)
166 : end do
167 :
168 14220 : delta = xk_new - xk_old
169 :
170 : vals(k_new) = &
171 : f(1, k_old) + delta*(f(2, k_old) &
172 20557 : + delta*(f(3, k_old) + delta*f(4, k_old)))
173 :
174 : end do
175 :
176 :
177 6337 : end subroutine do_interp_values
178 :
179 0 : subroutine do_interp_values_autodiff(init_x, nx, f1, nv, x, vals, ierr)
180 : type(auto_diff_real_2var_order1), intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
181 : integer, intent(in) :: nx ! length of init_x vector
182 : type(auto_diff_real_2var_order1), intent(in), pointer :: f1(:) ! =(4, nx) ! data & interpolation coefficients
183 : integer, intent(in) :: nv ! length of new x vector and vals vector
184 : type(auto_diff_real_2var_order1), intent(in) :: x(:) ! (nv) ! locations where want interpolated values
185 : ! strictly monotonic in same way as init_x
186 : ! values out of range of init_x's are clipped to boundaries of init_x's
187 : type(auto_diff_real_2var_order1), intent(inout) :: vals(:) ! (nv)
188 : integer, intent(out) :: ierr ! 0 means aok
189 :
190 : integer :: k_old, k_new
191 : type(auto_diff_real_2var_order1) :: xk_old, xkp1_old, xk_new, delta
192 : logical :: increasing
193 0 : type(auto_diff_real_2var_order1), pointer :: f(:,:) ! (4, nx) ! data & interpolation coefficients
194 :
195 0 : ierr = 0
196 :
197 0 : if (nx == 1) then
198 0 : vals(1:nv) = f1(1)
199 0 : return
200 : end if
201 :
202 0 : f(1:4,1:nx) => f1(1:4*nx)
203 :
204 0 : if(init_x(1) < init_x(nx)) then
205 : increasing = .true.
206 : end if
207 :
208 0 : k_old = 1; xk_old = init_x(k_old); xkp1_old = init_x(k_old+1)
209 :
210 0 : do k_new = 1, nv
211 :
212 0 : xk_new = x(k_new)
213 : if (increasing) then
214 0 : if (xk_new > init_x(nx)) then
215 0 : xk_new = init_x(nx)
216 0 : else if (xk_new < init_x(1)) then
217 0 : xk_new = init_x(1)
218 : end if
219 : else ! decreasing
220 : if (xk_new < init_x(nx)) then
221 : xk_new = init_x(nx)
222 : else if (xk_new > init_x(1)) then
223 : xk_new = init_x(1)
224 : end if
225 : end if
226 0 : do while ((increasing .and. xk_new > xkp1_old) .or. ((.not. increasing) .and. xk_new < xkp1_old))
227 0 : k_old = k_old + 1
228 0 : if (k_old >= nx) then
229 0 : k_old = k_old - 1
230 0 : xk_new = xkp1_old
231 0 : exit
232 : end if
233 0 : xk_old = xkp1_old
234 0 : xkp1_old = init_x(k_old+1)
235 : end do
236 :
237 0 : delta = xk_new - xk_old
238 :
239 : vals(k_new) = &
240 : f(1, k_old) + delta*(f(2, k_old) &
241 0 : + delta*(f(3, k_old) + delta*f(4, k_old)))
242 :
243 : end do
244 :
245 :
246 0 : end subroutine do_interp_values_autodiff
247 :
248 16 : subroutine do_interp_values_and_slopes(init_x, nx, f1, nv, x, vals, slopes, ierr)
249 : real(dp), intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
250 : integer, intent(in) :: nx ! length of init_x vector
251 : real(dp), intent(in), pointer :: f1(:) ! =(4, nx) ! data & interpolation coefficients
252 : integer, intent(in) :: nv ! length of new x vector and vals vector
253 : real(dp), intent(in) :: x(:) ! (nv) ! locations where want interpolated values
254 : ! strictly monotonic in same way as init_x
255 : ! values out of range of init_x's are clipped to boundaries of init_x's
256 : real(dp), intent(inout) :: vals(:) ! (nv)
257 : real(dp), intent(inout) :: slopes(:) ! (nv)
258 : integer, intent(out) :: ierr ! 0 means aok
259 :
260 : integer :: k_old, k_new
261 16 : real(dp) :: xk_old, xkp1_old, xk_new, delta
262 : logical :: increasing
263 16 : real(dp), pointer :: f(:,:) ! (4, nx) ! data & interpolation coefficients
264 16 : f(1:4,1:nx) => f1(1:4*nx)
265 :
266 16 : ierr = 0
267 :
268 16 : if (nx == 1) then
269 0 : vals(1:nv) = f(1,1)
270 0 : slopes(1:nv) = 0
271 : return
272 : end if
273 :
274 16 : k_old = 1; xk_old = init_x(k_old); xkp1_old = init_x(k_old+1)
275 :
276 16 : increasing = (init_x(1) < init_x(nx))
277 :
278 32 : do k_new = 1, nv
279 :
280 16 : xk_new = x(k_new)
281 16 : if (increasing) then
282 16 : if (xk_new > init_x(nx)) then
283 16 : xk_new = init_x(nx)
284 16 : else if (xk_new < init_x(1)) then
285 0 : xk_new = init_x(1)
286 : end if
287 : else ! decreasing
288 0 : if (xk_new < init_x(nx)) then
289 : xk_new = init_x(nx)
290 0 : else if (xk_new > init_x(1)) then
291 0 : xk_new = init_x(1)
292 : end if
293 : end if
294 4272 : do while ((increasing .and. xk_new > xkp1_old) .or. ((.not. increasing) .and. xk_new < xkp1_old))
295 4256 : k_old = k_old + 1
296 4256 : if (k_old >= nx) then
297 : k_old = k_old - 1
298 : xk_new = xkp1_old
299 : exit
300 : end if
301 4256 : xk_old = xkp1_old
302 4256 : xkp1_old = init_x(k_old+1)
303 : end do
304 :
305 16 : delta = xk_new - xk_old
306 :
307 : vals(k_new) = &
308 : f(1, k_old) + delta*(f(2, k_old) &
309 16 : + delta*(f(3, k_old) + delta*f(4, k_old)))
310 :
311 : slopes(k_new) = &
312 32 : f(2, k_old) + 2*delta*(f(3, k_old) + 1.5d0*delta*f(4, k_old))
313 :
314 : end do
315 :
316 :
317 16 : end subroutine do_interp_values_and_slopes
318 :
319 :
320 0 : subroutine do_interp2_values_and_slopes( &
321 0 : init_x, nx, f1_1, f1_2, nv, x, vals_1, slopes_1, vals_2, slopes_2, ierr)
322 : real(dp), intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
323 : integer, intent(in) :: nx ! length of init_x vector
324 : real(dp), intent(in), pointer :: f1_1(:), f1_2(:) ! =(4,nx) ! data & interpolation coefficients
325 : integer, intent(in) :: nv ! length of new x vector and vals vector
326 : real(dp), intent(in) :: x(:) ! (nv) ! locations where want interpolated values
327 : ! strictly monotonic in same way as init_x
328 : ! values out of range of init_x's are clipped to boundaries of init_x's
329 : real(dp), intent(inout) :: vals_1(:), vals_2(:) ! (nv)
330 : real(dp), intent(inout) :: slopes_1(:), slopes_2(:) ! (nv)
331 : integer, intent(out) :: ierr ! 0 means aok
332 :
333 : integer :: k_old, k_new
334 0 : real(dp) :: xk_old, xkp1_old, xk_new, delta
335 : logical :: increasing
336 0 : real(dp), pointer :: f_1(:,:), f_2(:,:) ! (4, nx) ! data & interpolation coefficients
337 0 : f_1(1:4,1:nx) => f1_1(1:4*nx)
338 0 : f_2(1:4,1:nx) => f1_2(1:4*nx)
339 :
340 0 : ierr = 0
341 :
342 0 : if (nx == 1) then
343 0 : vals_1(1:nv) = f_1(1,1)
344 0 : slopes_1(1:nv) = 0
345 0 : vals_2(1:nv) = f_2(1,1)
346 0 : slopes_2(1:nv) = 0
347 : return
348 : end if
349 :
350 0 : k_old = 1; xk_old = init_x(k_old); xkp1_old = init_x(k_old+1)
351 :
352 0 : increasing = (init_x(1) < init_x(nx))
353 :
354 0 : do k_new = 1, nv
355 :
356 0 : xk_new = x(k_new)
357 0 : if (increasing) then
358 0 : if (xk_new > init_x(nx)) then
359 0 : xk_new = init_x(nx)
360 0 : else if (xk_new < init_x(1)) then
361 0 : xk_new = init_x(1)
362 : end if
363 : else ! decreasing
364 0 : if (xk_new < init_x(nx)) then
365 : xk_new = init_x(nx)
366 0 : else if (xk_new > init_x(1)) then
367 0 : xk_new = init_x(1)
368 : end if
369 : end if
370 :
371 0 : do while ((increasing .and. xk_new > xkp1_old) .or. ((.not. increasing) .and. xk_new < xkp1_old))
372 0 : k_old = k_old + 1
373 0 : if (k_old >= nx) then
374 : k_old = k_old - 1
375 : xk_new = xkp1_old
376 : exit
377 : end if
378 0 : xk_old = xkp1_old
379 0 : xkp1_old = init_x(k_old+1)
380 : end do
381 :
382 0 : delta = xk_new - xk_old
383 :
384 : vals_1(k_new) = &
385 : f_1(1, k_old) + delta*(f_1(2, k_old) &
386 0 : + delta*(f_1(3, k_old) + delta*f_1(4, k_old)))
387 :
388 : slopes_1(k_new) = &
389 0 : f_1(2, k_old) + 2*delta*(f_1(3, k_old) + 1.5d0*delta*f_1(4, k_old))
390 :
391 : vals_2(k_new) = &
392 : f_2(1, k_old) + delta*(f_2(2, k_old) &
393 0 : + delta*(f_2(3, k_old) + delta*f_2(4, k_old)))
394 :
395 : slopes_2(k_new) = &
396 0 : f_2(2, k_old) + 2*delta*(f_2(3, k_old) + 1.5d0*delta*f_2(4, k_old))
397 :
398 : end do
399 :
400 :
401 0 : end subroutine do_interp2_values_and_slopes
402 :
403 :
404 0 : subroutine do_interp3_values_and_slopes( &
405 0 : init_x, nx, f1_1, f1_2, f1_3, nv, x, &
406 0 : vals_1, slopes_1, vals_2, slopes_2, vals_3, slopes_3, ierr)
407 : real(dp), intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
408 : integer, intent(in) :: nx ! length of init_x vector
409 : real(dp), intent(in), pointer :: f1_1(:), f1_2(:), f1_3(:) ! =(4,nx) ! data & interpolation coefficients
410 : integer, intent(in) :: nv ! length of new x vector and vals vector
411 : real(dp), intent(in) :: x(:) ! (nv) ! locations where want interpolated values
412 : ! strictly monotonic in same way as init_x
413 : ! values out of range of init_x's are clipped to boundaries of init_x's
414 : real(dp), intent(inout) :: vals_1(:), vals_2(:), vals_3(:) ! (nv)
415 : real(dp), intent(inout) :: slopes_1(:), slopes_2(:), slopes_3(:) ! (nv)
416 : integer, intent(out) :: ierr ! 0 means aok
417 :
418 : integer :: k_old, k_new
419 0 : real(dp) :: xk_old, xkp1_old, xk_new, delta
420 : logical :: increasing
421 0 : real(dp), pointer :: f_1(:,:), f_2(:,:), f_3(:,:) ! (4, nx) ! data & interpolation coefficients
422 0 : f_1(1:4,1:nx) => f1_1(1:4*nx)
423 0 : f_2(1:4,1:nx) => f1_2(1:4*nx)
424 0 : f_3(1:4,1:nx) => f1_3(1:4*nx)
425 :
426 0 : ierr = 0
427 :
428 0 : if (nx == 1) then
429 0 : vals_1(1:nv) = f_1(1,1)
430 0 : slopes_1(1:nv) = 0
431 0 : vals_2(1:nv) = f_2(1,1)
432 0 : slopes_2(1:nv) = 0
433 0 : vals_3(1:nv) = f_3(1,1)
434 0 : slopes_3(1:nv) = 0
435 : return
436 : end if
437 :
438 0 : k_old = 1; xk_old = init_x(k_old); xkp1_old = init_x(k_old+1)
439 :
440 0 : increasing = (init_x(1) < init_x(nx))
441 :
442 0 : do k_new = 1, nv
443 :
444 0 : xk_new = x(k_new)
445 0 : if (increasing) then
446 0 : if (xk_new > init_x(nx)) then
447 0 : xk_new = init_x(nx)
448 0 : else if (xk_new < init_x(1)) then
449 0 : xk_new = init_x(1)
450 : end if
451 : else ! decreasing
452 0 : if (xk_new < init_x(nx)) then
453 : xk_new = init_x(nx)
454 0 : else if (xk_new > init_x(1)) then
455 0 : xk_new = init_x(1)
456 : end if
457 : end if
458 :
459 0 : do while ((increasing .and. xk_new > xkp1_old) .or. ((.not. increasing) .and. xk_new < xkp1_old))
460 0 : k_old = k_old + 1
461 0 : if (k_old >= nx) then
462 : k_old = k_old - 1
463 : xk_new = xkp1_old
464 : exit
465 : end if
466 0 : xk_old = xkp1_old
467 0 : xkp1_old = init_x(k_old+1)
468 : end do
469 :
470 0 : delta = xk_new - xk_old
471 :
472 : vals_1(k_new) = &
473 : f_1(1, k_old) + delta*(f_1(2, k_old) &
474 0 : + delta*(f_1(3, k_old) + delta*f_1(4, k_old)))
475 : slopes_1(k_new) = &
476 0 : f_1(2, k_old) + 2*delta*(f_1(3, k_old) + 1.5d0*delta*f_1(4, k_old))
477 :
478 : vals_2(k_new) = &
479 : f_2(1, k_old) + delta*(f_2(2, k_old) &
480 0 : + delta*(f_2(3, k_old) + delta*f_2(4, k_old)))
481 : slopes_2(k_new) = &
482 0 : f_2(2, k_old) + 2*delta*(f_2(3, k_old) + 1.5d0*delta*f_2(4, k_old))
483 :
484 : vals_3(k_new) = &
485 : f_3(1, k_old) + delta*(f_3(2, k_old) &
486 0 : + delta*(f_3(3, k_old) + delta*f_3(4, k_old)))
487 : slopes_3(k_new) = &
488 0 : f_3(2, k_old) + 2*delta*(f_3(3, k_old) + 1.5d0*delta*f_3(4, k_old))
489 :
490 : end do
491 :
492 :
493 0 : end subroutine do_interp3_values_and_slopes
494 :
495 :
496 0 : subroutine do_interp6_values_and_slopes( &
497 0 : init_x, nx, f1_1, f1_2, f1_3, f1_4, f1_5, f1_6, nv, x, &
498 0 : vals_1, slopes_1, vals_2, slopes_2, vals_3, slopes_3, &
499 0 : vals_4, slopes_4, vals_5, slopes_5, vals_6, slopes_6, &
500 : ierr)
501 : real(dp), intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
502 : integer, intent(in) :: nx ! length of init_x vector
503 : real(dp), intent(in), pointer, dimension(:) :: f1_1, f1_2, f1_3, f1_4, f1_5, f1_6
504 : integer, intent(in) :: nv ! length of new x vector and vals vector
505 : real(dp), intent(in) :: x(:) ! (nv) ! locations where want interpolated values
506 : ! strictly monotonic in same way as init_x
507 : ! values out of range of init_x's are clipped to boundaries of init_x's
508 : real(dp), intent(inout), dimension(:) :: &
509 : vals_1, slopes_1, vals_2, slopes_2, vals_3, slopes_3, &
510 : vals_4, slopes_4, vals_5, slopes_5, vals_6, slopes_6
511 : integer, intent(out) :: ierr ! 0 means aok
512 :
513 : integer :: k_old, k_new
514 0 : real(dp) :: xk_old, xkp1_old, xk_new, delta
515 : logical :: increasing
516 0 : real(dp), pointer, dimension(:,:) :: f_1, f_2, f_3, f_4, f_5, f_6
517 0 : f_1(1:4,1:nx) => f1_1(1:4*nx)
518 0 : f_2(1:4,1:nx) => f1_2(1:4*nx)
519 0 : f_3(1:4,1:nx) => f1_3(1:4*nx)
520 0 : f_4(1:4,1:nx) => f1_4(1:4*nx)
521 0 : f_5(1:4,1:nx) => f1_5(1:4*nx)
522 0 : f_6(1:4,1:nx) => f1_6(1:4*nx)
523 :
524 0 : ierr = 0
525 :
526 0 : if (nx == 1) then
527 0 : vals_1(1:nv) = f_1(1,1); slopes_1(1:nv) = 0
528 0 : vals_2(1:nv) = f_2(1,1); slopes_2(1:nv) = 0
529 0 : vals_3(1:nv) = f_3(1,1); slopes_3(1:nv) = 0
530 0 : vals_4(1:nv) = f_4(1,1); slopes_4(1:nv) = 0
531 0 : vals_5(1:nv) = f_5(1,1); slopes_5(1:nv) = 0
532 0 : vals_6(1:nv) = f_6(1,1); slopes_6(1:nv) = 0
533 : return
534 : end if
535 :
536 0 : k_old = 1; xk_old = init_x(k_old); xkp1_old = init_x(k_old+1)
537 :
538 0 : increasing = (init_x(1) < init_x(nx))
539 :
540 0 : do k_new = 1, nv
541 :
542 0 : xk_new = x(k_new)
543 0 : if (increasing) then
544 0 : if (xk_new > init_x(nx)) then
545 0 : xk_new = init_x(nx)
546 0 : else if (xk_new < init_x(1)) then
547 0 : xk_new = init_x(1)
548 : end if
549 : else ! decreasing
550 0 : if (xk_new < init_x(nx)) then
551 : xk_new = init_x(nx)
552 0 : else if (xk_new > init_x(1)) then
553 0 : xk_new = init_x(1)
554 : end if
555 : end if
556 :
557 0 : do while ((increasing .and. xk_new > xkp1_old) .or. ((.not. increasing) .and. xk_new < xkp1_old))
558 0 : k_old = k_old + 1
559 0 : if (k_old >= nx) then
560 : k_old = k_old - 1
561 : xk_new = xkp1_old
562 : exit
563 : end if
564 0 : xk_old = xkp1_old
565 0 : xkp1_old = init_x(k_old+1)
566 : end do
567 :
568 0 : delta = xk_new - xk_old
569 :
570 : vals_1(k_new) = &
571 : f_1(1, k_old) + delta*(f_1(2, k_old) &
572 0 : + delta*(f_1(3, k_old) + delta*f_1(4, k_old)))
573 : slopes_1(k_new) = &
574 0 : f_1(2, k_old) + 2*delta*(f_1(3, k_old) + 1.5d0*delta*f_1(4, k_old))
575 :
576 : vals_2(k_new) = &
577 : f_2(1, k_old) + delta*(f_2(2, k_old) &
578 0 : + delta*(f_2(3, k_old) + delta*f_2(4, k_old)))
579 : slopes_2(k_new) = &
580 0 : f_2(2, k_old) + 2*delta*(f_2(3, k_old) + 1.5d0*delta*f_2(4, k_old))
581 :
582 : vals_3(k_new) = &
583 : f_3(1, k_old) + delta*(f_3(2, k_old) &
584 0 : + delta*(f_3(3, k_old) + delta*f_3(4, k_old)))
585 : slopes_3(k_new) = &
586 0 : f_3(2, k_old) + 2*delta*(f_3(3, k_old) + 1.5d0*delta*f_3(4, k_old))
587 :
588 : vals_4(k_new) = &
589 : f_4(1, k_old) + delta*(f_4(2, k_old) &
590 0 : + delta*(f_4(3, k_old) + delta*f_4(4, k_old)))
591 : slopes_4(k_new) = &
592 0 : f_4(2, k_old) + 2*delta*(f_4(3, k_old) + 1.5d0*delta*f_4(4, k_old))
593 :
594 : vals_5(k_new) = &
595 : f_5(1, k_old) + delta*(f_5(2, k_old) &
596 0 : + delta*(f_5(3, k_old) + delta*f_5(4, k_old)))
597 : slopes_5(k_new) = &
598 0 : f_5(2, k_old) + 2*delta*(f_5(3, k_old) + 1.5d0*delta*f_5(4, k_old))
599 :
600 : vals_6(k_new) = &
601 : f_6(1, k_old) + delta*(f_6(2, k_old) &
602 0 : + delta*(f_6(3, k_old) + delta*f_6(4, k_old)))
603 : slopes_6(k_new) = &
604 0 : f_6(2, k_old) + 2*delta*(f_6(3, k_old) + 1.5d0*delta*f_6(4, k_old))
605 :
606 : end do
607 :
608 0 : end subroutine do_interp6_values_and_slopes
609 :
610 :
611 2630 : real(dp) function minmod1(f1, f2)
612 : real(dp), intent(in) :: f1, f2
613 2630 : minmod1 = 0.5d0 * (sign(1d0, f1) + sign(1d0, f2)) * min(abs(f1), abs(f2))
614 2630 : end function minmod1
615 :
616 :
617 0 : real(dp) function median1(f1, f2, f3)
618 : real(dp), intent(in) :: f1, f2, f3
619 0 : median1 = f1 + minmod1(f2 - f1, f3 - f1)
620 0 : end function median1
621 :
622 :
623 9205 : subroutine minmod(z, n, f1, f2)
624 : real(dp), intent(inout) :: z(:)
625 : integer, intent(in) :: n ! length of vectors
626 : real(dp), intent(in) :: f1(:), f2(:)
627 91775919 : z(1:n) = 0.5d0 * (sign(1d0, f1(1:n)) + sign(1d0, f2(1:n))) * min(abs(f1(1:n)), abs(f2(1:n)))
628 9205 : end subroutine minmod
629 :
630 :
631 1315 : subroutine minmod4(z, n, f1, f2, f3, f4)
632 : real(dp), intent(inout) :: z(:)
633 : integer, intent(in) :: n ! length of vectors
634 : real(dp), intent(in) :: f1(:), f2(:), f3(:), f4(:)
635 1315 : call minmod(z, n, f1, f2)
636 1315 : call minmod(z, n, z, f3)
637 1315 : call minmod(z, n, z, f4)
638 1315 : end subroutine minmod4
639 :
640 :
641 1315 : subroutine median(z, n, f1, f2, f3)
642 : real(dp), intent(inout) :: z(:)
643 : integer, intent(in) :: n ! length of vectors
644 : real(dp), intent(in) :: f1(:), f2(:), f3(:)
645 26220564 : real(dp), target :: tmp1_ary(n), tmp2_ary(n)
646 1315 : real(dp), pointer :: tmp1(:), tmp2(:)
647 1315 : tmp1 => tmp1_ary
648 1315 : tmp2 => tmp2_ary
649 13110282 : tmp1(1:n) = f2(1:n) - f1(1:n)
650 13110282 : tmp2(1:n) = f3(1:n) - f1(1:n)
651 1315 : call minmod(z(1:n), n, tmp1(1:n), tmp2(1:n))
652 13110282 : z(1:n) = z(1:n) + f1(1:n)
653 1315 : end subroutine median
654 :
655 :
656 0 : type(auto_diff_real_2var_order1) function minmod1_autodiff(f1, f2)
657 : use auto_diff
658 : type(auto_diff_real_2var_order1), intent(in) :: f1, f2
659 0 : minmod1_autodiff = 0.5d0 * (sign(f1) + sign(f2)) * min(abs(f1), abs(f2))
660 0 : end function minmod1_autodiff
661 :
662 :
663 0 : type(auto_diff_real_2var_order1) function median1_autodiff(f1, f2, f3)
664 0 : use auto_diff
665 : type(auto_diff_real_2var_order1), intent(in) :: f1, f2, f3
666 0 : median1_autodiff = f1 + minmod1_autodiff(f2 - f1, f3 - f1)
667 0 : end function median1_autodiff
668 :
669 :
670 0 : subroutine minmod_autodiff(z, n, f1, f2)
671 0 : use auto_diff
672 : type(auto_diff_real_2var_order1), intent(inout) :: z(:)
673 : integer, intent(in) :: n ! length of vectors
674 : integer :: k
675 : type(auto_diff_real_2var_order1), intent(in) :: f1(:), f2(:)
676 0 : do k =1, n
677 0 : z(k) = 0.5d0 * (sign(f1(k)) + sign(f2(k))) * min(abs(f1(k)), abs(f2(k)))
678 : end do
679 0 : end subroutine minmod_autodiff
680 :
681 0 : subroutine minmod4_autodiff(z, n, f1, f2, f3, f4)
682 0 : use auto_diff
683 : type(auto_diff_real_2var_order1), intent(inout) :: z(:)
684 : integer, intent(in) :: n ! length of vectors
685 : type(auto_diff_real_2var_order1), intent(in) :: f1(:), f2(:), f3(:), f4(:)
686 0 : call minmod_autodiff(z, n, f1, f2)
687 0 : call minmod_autodiff(z, n, z, f3)
688 0 : call minmod_autodiff(z, n, z, f4)
689 0 : end subroutine minmod4_autodiff
690 :
691 :
692 0 : subroutine median_autodiff(z, n, f1, f2, f3)
693 0 : use auto_diff
694 : type(auto_diff_real_2var_order1), intent(inout) :: z(:)
695 : integer, intent(in) :: n ! length of vectors
696 : integer :: k
697 : type(auto_diff_real_2var_order1), intent(in) :: f1(:), f2(:), f3(:)
698 0 : type(auto_diff_real_2var_order1), target :: tmp1_ary(n), tmp2_ary(n)
699 : type(auto_diff_real_2var_order1), pointer :: tmp1(:), tmp2(:)
700 0 : tmp1 => tmp1_ary
701 0 : tmp2 => tmp2_ary
702 0 : do k =1,n
703 0 : tmp1(k) = f2(k) - f1(k)
704 0 : tmp2(k) = f3(k) - f1(k)
705 : end do
706 0 : call minmod_autodiff(z, n, tmp1, tmp2)
707 0 : do k =1,n
708 0 : z(k) = z(k) + f1(k)
709 : end do
710 0 : end subroutine median_autodiff
711 :
712 : end module interp_1d_misc
|
