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_pm_sg ! piecewise monotonic algorithms
21 :
22 : implicit none
23 :
24 : contains
25 :
26 : ! the following produce piecewise monotonic interpolants
27 : ! rather than monotonicity preserving
28 : ! this stricter limit never introduces interpolated values exceeding the given values,
29 : ! even in places where the given values are not monotonic.
30 : ! the downside is reduced accuracy on smooth data compared to the mp routines.
31 :
32 :
33 46 : subroutine mk_pmcub_sg(x, nx, f1, slope_only, nwork, work1, str, ierr)
34 : ! make piecewise monotonic cubic interpolant
35 : use interp_1d_def
36 : integer, intent(in) :: nx ! length of x vector (nx >= 2)
37 : real, intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
38 : real, intent(inout), pointer :: f1(:) ! =(4, nx) ! data & interpolation coefficients
39 : logical, intent(in) :: slope_only
40 : integer, intent(in) :: nwork ! nwork must be >= pm_work_size (see interp_1d_def)
41 : real, intent(inout), pointer :: work1(:) ! =(nx, nwork)
42 : character (len=*) :: str
43 : integer, intent(out) :: ierr
44 :
45 46 : real, dimension(:), pointer :: h, s, p
46 : integer :: i
47 : logical, parameter :: dbg = .true.
48 46 : real, pointer :: f(:,:) ! (4, nx) ! data & interpolation coefficients
49 46 : f(1:4,1:nx) => f1(1:4*nx)
50 :
51 : include 'formats'
52 :
53 46 : ierr = 0
54 :
55 46 : if (nx < 2) then
56 7 : return
57 : end if
58 :
59 46 : if (nx == 2) then
60 0 : call mk_pmlinear_sg(x, f1, slope_only, nwork, work1, str, ierr)
61 0 : return
62 : end if
63 :
64 46 : if (nx == 3) then
65 1 : call mk_pmquad_sg(x, f1, slope_only, nwork, work1, str, ierr)
66 1 : return
67 : end if
68 :
69 45 : if (nx == 4) then
70 6 : call mk_pmcub4_sg(x, f1, slope_only, nwork, work1, str, ierr)
71 6 : return
72 : end if
73 :
74 39 : if (nwork < pm_work_size) then
75 0 : ierr = -1
76 0 : return
77 : end if
78 :
79 39 : h(1:nx) => work1(1:nx)
80 39 : s(1:nx) => work1(1+nx:2*nx)
81 39 : p(1:nx) => work1(1+2*nx:3*nx)
82 :
83 : if (dbg) then
84 4494 : h(:) = 0; s(:) = 0; p(:) = 0
85 : end if
86 :
87 1485 : do i=1,nx-1
88 1485 : h(i) = x(i+1) - x(i) ! width of interval
89 : end do
90 1485 : do i = 1, nx-1
91 1485 : if (h(i) == 0) then
92 : write(*, '(a,1x,2i5,1x,a)') &
93 0 : 'same interpolation x values at', i, i+1, 'for ' // trim(str)
94 0 : ierr = -1
95 0 : return
96 : end if
97 : end do
98 :
99 1485 : do i=1,nx-1
100 1485 : s(i) = (f(1,i+1) - f(1,i)) / h(i) ! slope across interval
101 : end do
102 1446 : do i=2,nx-1
103 1446 : p(i) = (s(i-1)*h(i) + s(i)*h(i-1))/(h(i-1)+h(i))
104 : ! slope at i of parabola through i-1, i, and i+1
105 : end do
106 1446 : do i=2,nx-1
107 : f(2,i) = (sign(1.0, s(i-1))+sign(1.0, s(i)))* &
108 1446 : min(abs(s(i-1)), abs(s(i)), 0.5*abs(p(i)))
109 : ! "safe" slope at i to ensure monotonic -- see Steffen's paper for explanation.
110 : end do
111 :
112 39 : p(1) = s(1)*(1 + h(1) / (h(1) + h(2))) - s(2) * h(1) / (h(1) + h(2))
113 : ! slope at 1 of parabola through 1st 3 points
114 39 : if (p(1)*s(1) <= 0) then
115 0 : f(2, 1) = 0
116 39 : else if (abs(p(1)) > 2*abs(s(1))) then
117 0 : f(2, 1) = 2*s(1)
118 : else
119 39 : f(2, 1) = p(1)
120 : end if
121 :
122 : p(nx) = s(nx-1)*(1 + h(nx-1) / (h(nx-1) + h(nx-2))) &
123 39 : - s(nx-2)*h(nx-1) / (h(nx-1) + h(nx-2))
124 : ! slope at nx of parabola through last 3 points
125 39 : if (p(nx)*s(nx-1) <= 0) then
126 0 : f(2,nx) = 0
127 39 : else if (abs(p(nx)) > 2*abs(s(nx-1))) then
128 0 : f(2,nx) = 2*s(nx-1)
129 : else
130 39 : f(2,nx) = p(nx)
131 : end if
132 :
133 39 : if (slope_only) return
134 :
135 1485 : do i=1,nx-1
136 1446 : f(3,i) = (3*s(i) - 2*f(2,i) - f(2,i+1)) / h(i)
137 1485 : f(4,i) = (f(2,i) + f(2,i+1) - 2*s(i)) / (h(i)*h(i))
138 : end do
139 39 : f(3,nx) = 0
140 39 : f(4,nx) = 0
141 :
142 46 : end subroutine mk_pmcub_sg
143 :
144 :
145 : ! optimize special case for nx = 4
146 6 : subroutine mk_pmcub4_sg(x, f1, slope_only, nwork, work1, str, ierr)
147 : use interp_1d_def
148 : integer, parameter :: nx = 4 ! length of x vector (nx >= 2)
149 : real, intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
150 : real, intent(inout), pointer :: f1(:) ! =(4, nx) ! data & interpolation coefficients
151 : logical, intent(in) :: slope_only
152 : integer, intent(in) :: nwork ! nwork must be >= pm_work_size (see interp_1d_def)
153 : real, intent(inout), pointer :: work1(:) ! =(nx, nwork)
154 : character (len=*) :: str
155 : integer, intent(out) :: ierr
156 :
157 6 : real, dimension(:), pointer :: h, s, p
158 : integer :: i
159 6 : real, pointer :: f(:,:) ! (4, nx) ! data & interpolation coefficients
160 6 : f(1:4,1:nx) => f1(1:4*nx)
161 :
162 6 : ierr = 0
163 :
164 6 : if (nwork < pm_work_size) then
165 0 : ierr = -1
166 0 : return
167 : end if
168 :
169 6 : h(1:nx) => work1(1:nx)
170 6 : s(1:nx) => work1(1+nx:2*nx)
171 6 : p(1:nx) => work1(1+2*nx:3*nx)
172 :
173 24 : do i=1,nx-1
174 24 : h(i) = x(i+1) - x(i) ! width of interval
175 : end do
176 24 : do i = 1, nx-1
177 24 : if (h(i) == 0) then
178 : write(*, '(a,1x,2i5,1x,a)') &
179 0 : 'same interpolation x values at', i, i+1, 'for ' // trim(str)
180 0 : ierr = -1
181 0 : return
182 : end if
183 : end do
184 :
185 24 : do i=1,nx-1
186 24 : s(i) = (f(1,i+1) - f(1,i)) / h(i) ! slope across interval
187 : end do
188 18 : do i=2,nx-1
189 18 : p(i) = (s(i-1)*h(i) + s(i)*h(i-1))/(h(i-1)+h(i))
190 : ! slope at i of parabola through i-1, i, and i+1
191 : end do
192 18 : do i=2,nx-1
193 : f(2,i) = (sign(1.0, s(i-1))+sign(1.0, s(i)))* &
194 18 : min(abs(s(i-1)), abs(s(i)), 0.5*abs(p(i)))
195 : ! "safe" slope at i to ensure monotonic -- see Steffen's paper for explanation.
196 : end do
197 :
198 6 : p(1) = s(1)*(1 + h(1) / (h(1) + h(2))) - s(2) * h(1) / (h(1) + h(2))
199 : ! slope at 1 of parabola through 1st 3 points
200 6 : if (p(1)*s(1) <= 0) then
201 0 : f(2, 1) = 0
202 6 : else if (abs(p(1)) > 2*abs(s(1))) then
203 0 : f(2, 1) = 2*s(1)
204 : else
205 6 : f(2, 1) = p(1)
206 : end if
207 :
208 : p(nx) = s(nx-1)*(1 + h(nx-1) / (h(nx-1) + h(nx-2))) &
209 6 : - s(nx-2)*h(nx-1) / (h(nx-1) + h(nx-2))
210 : ! slope at nx of parabola through last 3 points
211 6 : if (p(nx)*s(nx-1) <= 0) then
212 0 : f(2, nx) = 0
213 6 : else if (abs(p(nx)) > 2*abs(s(nx-1))) then
214 0 : f(2, nx) = 2*s(nx-1)
215 : else
216 6 : f(2, nx) = p(nx)
217 : end if
218 :
219 6 : if (slope_only) return
220 :
221 24 : do i=1,nx-1
222 18 : f(3,i) = (3*s(i) - 2*f(2,i) - f(2,i+1)) / h(i)
223 24 : f(4,i) = (f(2,i) + f(2,i+1) - 2*s(i)) / (h(i)*h(i))
224 : end do
225 6 : f(3,nx) = 0
226 6 : f(4,nx) = 0
227 :
228 6 : end subroutine mk_pmcub4_sg
229 :
230 :
231 : ! optimize special case for nx = 3
232 1 : subroutine mk_pmquad_sg(x, f1, slope_only, nwork, work1, str, ierr)
233 : ! make piecewise monotonic quadratic interpolant
234 : use interp_1d_def
235 : integer, parameter :: nx = 3
236 : real, intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
237 : real, intent(inout), pointer :: f1(:) ! =(4, nx) ! data & interpolation coefficients
238 : logical, intent(in) :: slope_only
239 : integer, intent(in) :: nwork ! nwork must be >= pm_work_size (see interp_1d_def)
240 : real, intent(inout), pointer :: work1(:) ! =(nx, nwork)
241 : character (len=*) :: str
242 : integer, intent(out) :: ierr
243 :
244 1 : real, dimension(:), pointer :: h, s, p
245 : integer :: i
246 1 : real, pointer :: f(:,:) ! (4, nx) ! data & interpolation coefficients
247 1 : f(1:4,1:nx) => f1(1:4*nx)
248 :
249 1 : if (nwork < pm_work_size) then
250 0 : ierr = -1
251 0 : return
252 : end if
253 1 : ierr = 0
254 :
255 1 : h(1:nx) => work1(1:nx)
256 1 : s(1:nx) => work1(1+nx:2*nx)
257 1 : p(1:nx) => work1(1+2*nx:3*nx)
258 :
259 3 : do i=1,nx-1
260 2 : h(i) = x(i+1) - x(i) ! width of interval
261 3 : if (h(i) == 0) then
262 : write(*, '(a,1x,2i5,1x,a)') &
263 0 : 'same interpolation x values at', i, i+1, 'for ' // trim(str)
264 0 : ierr = -1
265 0 : return
266 : end if
267 : end do
268 :
269 3 : do i=1,nx-1
270 3 : s(i) = (f(1,i+1) - f(1,i)) / h(i) ! slope across interval
271 : end do
272 2 : do i=2,nx-1
273 2 : p(i) = (s(i-1)*h(i) + s(i)*h(i-1))/(h(i-1)+h(i))
274 : ! slope at i of parabola through i-1, i, and i+1
275 : end do
276 2 : do i=2,nx-1
277 : f(2,i) = (sign(1.0, s(i-1))+sign(1.0, s(i)))* &
278 2 : min(abs(s(i-1)), abs(s(i)), 0.5*abs(p(i)))
279 : ! "safe" slope at i to ensure monotonic -- see Steffen's paper for explanation.
280 : end do
281 :
282 1 : p(1) = s(1)*(1 + h(1) / (h(1) + h(2))) - s(2) * h(1) / (h(1) + h(2))
283 : ! slope at 1 of parabola through 1st 3 points
284 1 : if (sign(1.0,p(1)) /= sign(1.0,s(1))) then
285 0 : f(2, 1) = 0
286 1 : else if (abs(p(1)) > 2*abs(s(1))) then
287 0 : f(2, 1) = 2*s(1)
288 : else
289 1 : f(2, 1) = p(1)
290 : end if
291 :
292 : p(nx) = s(nx-1)*(1 + h(nx-1) / (h(nx-1) + h(nx-2))) &
293 1 : - s(nx-2)*h(nx-1) / (h(nx-1) + h(nx-2))
294 : ! slope at nx of parabola through last 3 points
295 1 : if (sign(1.0,p(nx)) /= sign(1.0,s(nx-1))) then
296 0 : f(2, nx) = 0
297 1 : else if (abs(p(nx)) > 2*abs(s(nx-1))) then
298 0 : f(2, nx) = 2*s(nx-1)
299 : else
300 1 : f(2, nx) = p(nx)
301 : end if
302 :
303 1 : if (slope_only) return
304 :
305 3 : do i=1,nx-1
306 2 : f(3,i) = (3*s(i) - 2*f(2,i) - f(2,i+1)) / h(i)
307 3 : f(4,i) = (f(2,i) + f(2,i+1) - 2*s(i)) / (h(i)*h(i))
308 : end do
309 1 : f(3,nx) = 0
310 1 : f(4,nx) = 0
311 :
312 1 : end subroutine mk_pmquad_sg
313 :
314 :
315 : ! optimize special case for nx = 2
316 0 : subroutine mk_pmlinear_sg(x, f1, slope_only, nwork, work1, str, ierr)
317 : use interp_1d_def
318 : integer, parameter :: nx = 2
319 : real, intent(in) :: x(:) ! (nx) ! junction points, strictly monotonic
320 : real, intent(inout), pointer :: f1(:) ! =(4, nx) ! data & interpolation coefficients
321 : logical, intent(in) :: slope_only
322 : integer, intent(in) :: nwork ! nwork must be >= pm_work_size (see interp_1d_def)
323 : real, intent(inout), pointer :: work1(:) ! =(nx, nwork)
324 : character (len=*) :: str
325 : integer, intent(out) :: ierr
326 :
327 0 : real :: h, s
328 0 : real, pointer :: f(:,:) ! (4, nx) ! data & interpolation coefficients
329 0 : f(1:4,1:nx) => f1(1:4*nx)
330 :
331 0 : ierr = 0
332 :
333 0 : if (nwork < pm_work_size) then
334 0 : ierr = -1
335 0 : return
336 : end if
337 :
338 0 : h = x(2) - x(1) ! width of interval
339 0 : if (h == 0) then
340 : write(*, '(a,1x,2i5,1x,a)') &
341 0 : 'same interpolation x values at', 1, 2, 'for ' // trim(str)
342 0 : ierr = -1
343 0 : return
344 : end if
345 :
346 0 : s = (f(1, 2) - f(1, 1)) / h ! slope across interval
347 0 : f(2, 1) = s
348 0 : f(2, 2) = 0
349 :
350 0 : if (slope_only) return
351 :
352 0 : f(3, 1:2) = 0
353 0 : f(4, 1:2) = 0
354 :
355 0 : end subroutine mk_pmlinear_sg
356 :
357 :
358 0 : subroutine mk_pmcub_uniform_sg(dx, n, f1, slope_only, nwork, work1, str, ierr)
359 : ! make piecewise monotonic cubic interpolant on unit spaced mesh
360 : use interp_1d_def
361 : integer, intent(in) :: n ! length of vector
362 : real, intent(in) :: dx
363 : real, intent(inout), pointer :: f1(:) ! =(4, n) ! data & interpolation coefficients
364 : logical, intent(in) :: slope_only
365 : integer, intent(in) :: nwork ! nwork must be >= pm_work_size (see interp_1d_def)
366 : real, intent(inout), pointer :: work1(:) ! =(nx, nwork)
367 : character (len=*) :: str
368 : integer, intent(out) :: ierr
369 :
370 0 : real, dimension(:), pointer :: s, p
371 0 : real :: x(2)
372 : integer :: i
373 0 : real, pointer :: f(:,:) ! (4, n) ! data & interpolation coefficients
374 0 : f(1:4,1:n) => f1(1:4*n)
375 :
376 0 : ierr = 0
377 :
378 0 : if (n < 2) then
379 0 : return
380 : end if
381 :
382 0 : if (n == 2) then
383 0 : x(1) = 0
384 0 : x(2) = dx
385 0 : call mk_pmlinear_sg(x, f1, slope_only, nwork, work1, str, ierr)
386 0 : return
387 : end if
388 :
389 0 : if (dx == 0) then
390 0 : ierr = -1
391 0 : return
392 : end if
393 :
394 0 : if (nwork < pm_work_size) then
395 0 : ierr = -1
396 0 : return
397 : end if
398 : ierr = 0
399 :
400 0 : s(1:n) => work1(1:n)
401 0 : p(1:n) => work1(1+n:2*n)
402 :
403 : ierr = 0
404 :
405 0 : do i=1,n-1
406 0 : s(i) = f(1,i+1) - f(1,i) ! slope across interval
407 : end do
408 0 : do i=2,n-1
409 0 : p(i) = 0.5*(s(i-1) + s(i))
410 : ! slope at i of parabola through i-1, i, and i+1
411 : end do
412 0 : do i=2,n-1
413 : f(2,i) = (sign(1.0, s(i-1))+sign(1.0, s(i)))* &
414 0 : min(abs(s(i-1)), abs(s(i)), 0.5*abs(p(i)))
415 : ! "safe" slope at i to ensure monotonic -- see Steffen's paper for explanation.
416 : end do
417 :
418 0 : p(1) = 1.5 * s(1) - 0.5 * s(2)
419 : ! slope at 1 of parabola through 1st 3 points
420 0 : if (p(1)*s(1) <= 0) then
421 0 : f(2, 1) = 0
422 0 : else if (abs(p(1)) > 2*abs(s(1))) then
423 0 : f(2, 1) = 2*s(1)
424 : else
425 0 : f(2, 1) = p(1)
426 : end if
427 :
428 0 : p(n) = 1.5 * s(n-1) - 0.5 * s(n-2)
429 : ! slope at n of parabola through last 3 points
430 0 : if (p(n)*s(n-1) <= 0) then
431 0 : f(2, n) = 0
432 0 : else if (abs(p(n)) > 2*abs(s(n-1))) then
433 0 : f(2, n) = 2*s(n-1)
434 : else
435 0 : f(2, n) = p(n)
436 : end if
437 :
438 0 : f(2, 1:n) = f(2, 1:n) / dx
439 :
440 0 : if (slope_only) return
441 :
442 : ! 2nd and 3rd derivatives
443 0 : do i=1,n-1
444 0 : f(3,i) = (3*s(i) - 2*f(2,i) - f(2,i+1)) / dx
445 0 : f(4,i) = (f(2,i) + f(2,i+1) - 2*s(i)) / (dx*dx)
446 : end do
447 0 : f(3,n) = (3*f(1, n-1) - 3*f(1, n) + (f(2, n-1) + 2*f(2, n)) * dx) / (dx*dx)
448 0 : f(4,n) = (-2*f(1, n-1) + 2*f(1, n) - (f(2, n-1) + f(2, n))*dx) / (dx*dx*dx)
449 :
450 0 : end subroutine mk_pmcub_uniform_sg
451 :
452 :
453 : end module interp_1d_pm_sg
|
