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_sg
21 :
22 : implicit none
23 :
24 : contains
25 :
26 2 : subroutine do_integrate_values_sg(init_x, nx, f1, nv, x, vals, ierr)
27 : real, intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
28 : integer, intent(in) :: nx ! length of init_x vector
29 : real, intent(in), pointer :: f1(:) ! =(4, nx) ! data & interpolation coefficients
30 : integer, intent(in) :: nv ! length of new x vector and vals vector
31 : real, intent(in) :: x(:) ! (nv)
32 : ! strictly monotonic in same way as init_x
33 : real, intent(inout) :: vals(:) ! (nv)
34 : ! for i > 1, vals(i) = integral of interpolating poly from x(i-1) to x(i)
35 : ! vals(1) = 0
36 : integer, intent(out) :: ierr ! 0 means aok
37 :
38 : integer :: k_old, k_new
39 2 : real :: xk_old, xkp1_old, xk_new, xk_prev, sum
40 : logical :: increasing
41 : real, pointer :: f(:,:) ! (4, nx) ! data & interpolation coefficients
42 2 : f(1:4,1:nx) => f1(1:4*nx)
43 :
44 2 : increasing = (init_x(1) < init_x(nx))
45 :
46 : if (increasing .and. (x(1) < init_x(1) .or. x(nv) > init_x(nx)) &
47 2 : .or. ((.not. increasing) .and. (x(1) > init_x(1) .or. x(nv) < init_x(nx)))) then
48 0 : ierr = -1
49 0 : return
50 : end if
51 :
52 2 : ierr = 0
53 :
54 2 : k_old = 1; xk_old = init_x(k_old); xkp1_old = init_x(k_old+1)
55 2 : sum = 0; xk_prev = x(1); vals(1) = 0
56 :
57 14 : do k_new = 2, nv
58 :
59 12 : xk_new = x(k_new)
60 74 : do while ((increasing .and. xk_new > xkp1_old) .or. ((.not. increasing) .and. xk_new < xkp1_old))
61 62 : k_old = k_old + 1
62 62 : if (k_old >= nx) then
63 0 : k_old = k_old - 1
64 0 : xk_new = xkp1_old
65 0 : exit
66 : end if
67 62 : call add_to_integral(k_old - 1, xkp1_old)
68 62 : xk_old = xkp1_old
69 68 : xkp1_old = init_x(k_old+1)
70 : end do
71 :
72 12 : call add_to_integral(k_old, xk_new)
73 12 : vals(k_new) = sum
74 14 : sum = 0
75 :
76 : end do
77 :
78 : contains
79 :
80 74 : subroutine add_to_integral(k, x2)
81 : integer, intent(in) :: k
82 : real, intent(in) :: x2
83 :
84 74 : real :: x0, x1, a1, a2, d1, d2, area
85 :
86 74 : x0 = init_x(k)
87 74 : x1 = xk_prev
88 74 : if (x1 == x2) return
89 74 : d2 = x2 - x0
90 : a2 = d2*(f(1, k) + d2*(f(2, k)/2 &
91 74 : + d2*(f(3, k)/3 + d2*f(4, k)/4)))
92 74 : if (x1 > x0) then
93 5 : d1 = x1 - x0
94 : a1 = d1*(f(1, k) + d1*(f(2, k)/2 &
95 5 : + d1*(f(3, k)/3 + d1*f(4, k)/4)))
96 5 : area = a2 - a1
97 : else
98 74 : d1 = 0; a1 = 0; area = a2
99 : end if
100 74 : sum = sum + area
101 74 : xk_prev = x2
102 :
103 : end subroutine add_to_integral
104 :
105 :
106 : end subroutine do_integrate_values_sg
107 :
108 :
109 15 : subroutine do_interp_values_sg(init_x, nx, f1, nv, x, vals, ierr)
110 : real, intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
111 : integer, intent(in) :: nx ! length of init_x vector
112 : real, intent(in), pointer :: f1(:) ! =(4, nx) ! data & interpolation coefficients
113 : integer, intent(in) :: nv ! length of new x vector and vals vector
114 : real, intent(in) :: x(:) ! (nv) ! locations where want interpolated values
115 : ! strictly monotonic in same way as init_x
116 : ! values out of range of init_x's are clipped to boundaries of init_x's
117 : real, intent(inout) :: vals(:) ! (nv)
118 : integer, intent(out) :: ierr ! 0 means aok
119 :
120 : integer :: k_old, k_new
121 15 : real :: xk_old, xkp1_old, xk_new, delta
122 : logical :: increasing
123 15 : real, pointer :: f(:,:) ! (4, nx) ! data & interpolation coefficients
124 15 : f(1:4,1:nx) => f1(1:4*nx)
125 :
126 15 : ierr = 0
127 :
128 15 : if (nx == 1) then
129 0 : vals(1:nv) = f1(1)
130 : return
131 : end if
132 :
133 15 : f(1:4,1:nx) => f1(1:4*nx)
134 :
135 15 : increasing = (init_x(1) < init_x(nx))
136 :
137 15 : k_old = 1; xk_old = init_x(k_old); xkp1_old = init_x(k_old+1)
138 :
139 55 : do k_new = 1, nv
140 :
141 40 : xk_new = x(k_new)
142 40 : if (increasing) then
143 24 : if (xk_new > init_x(nx)) then
144 40 : xk_new = init_x(nx)
145 24 : else if (xk_new < init_x(1)) then
146 0 : xk_new = init_x(1)
147 : end if
148 : else ! decreasing
149 16 : if (xk_new < init_x(nx)) then
150 : xk_new = init_x(nx)
151 16 : else if (xk_new > init_x(1)) then
152 0 : xk_new = init_x(1)
153 : end if
154 : end if
155 217 : do while ((increasing .and. xk_new > xkp1_old) .or. ((.not. increasing) .and. xk_new < xkp1_old))
156 177 : k_old = k_old + 1
157 177 : if (k_old >= nx) then
158 : k_old = k_old - 1
159 : xk_new = xkp1_old
160 : exit
161 : end if
162 177 : xk_old = xkp1_old
163 193 : xkp1_old = init_x(k_old+1)
164 : end do
165 :
166 40 : delta = xk_new - xk_old
167 :
168 : vals(k_new) = &
169 : f(1, k_old) + delta*(f(2, k_old) &
170 55 : + delta*(f(3, k_old) + delta*f(4, k_old)))
171 :
172 : end do
173 :
174 :
175 15 : end subroutine do_interp_values_sg
176 :
177 :
178 0 : subroutine do_interp_values_and_slopes_sg(init_x, nx, f1, nv, x, vals, slopes, ierr)
179 : real, intent(in) :: init_x(:) ! (nx) ! junction points, strictly monotonic
180 : integer, intent(in) :: nx ! length of init_x vector
181 : real, intent(in), pointer :: f1(:) ! =(4, nx) ! data & interpolation coefficients
182 : integer, intent(in) :: nv ! length of new x vector and vals vector
183 : real, intent(in) :: x(:) ! (nv) ! locations where want interpolated values
184 : ! strictly monotonic in same way as init_x
185 : ! values out of range of init_x's are clipped to boundaries of init_x's
186 : real, intent(inout) :: vals(:) ! (nv)
187 : real, intent(inout) :: slopes(:) ! (nv)
188 : integer, intent(out) :: ierr ! 0 means aok
189 :
190 : integer :: k_old, k_new
191 0 : real :: xk_old, xkp1_old, xk_new, delta
192 : logical :: increasing
193 0 : real, pointer :: f(:,:) ! (4, nx) ! data & interpolation coefficients
194 0 : f(1:4,1:nx) => f1(1:4*nx)
195 :
196 0 : ierr = 0
197 :
198 0 : k_old = 1; xk_old = init_x(k_old); xkp1_old = init_x(k_old+1)
199 :
200 0 : increasing = (init_x(1) < init_x(nx))
201 :
202 0 : do k_new = 1, nv
203 :
204 0 : xk_new = x(k_new)
205 0 : if (increasing) then
206 0 : if (xk_new > init_x(nx)) then
207 0 : xk_new = init_x(nx)
208 0 : else if (xk_new < init_x(1)) then
209 0 : xk_new = init_x(1)
210 : end if
211 : else ! decreasing
212 0 : if (xk_new < init_x(nx)) then
213 : xk_new = init_x(nx)
214 0 : else if (xk_new > init_x(1)) then
215 0 : xk_new = init_x(1)
216 : end if
217 : end if
218 0 : do while ((increasing .and. xk_new > xkp1_old) .or. ((.not. increasing) .and. xk_new < xkp1_old))
219 0 : k_old = k_old + 1
220 0 : if (k_old >= nx) then
221 : k_old = k_old - 1
222 : xk_new = xkp1_old
223 : exit
224 : end if
225 0 : xk_old = xkp1_old
226 0 : xkp1_old = init_x(k_old+1)
227 : end do
228 :
229 0 : delta = xk_new - xk_old
230 :
231 : vals(k_new) = &
232 : f(1, k_old) + delta*(f(2, k_old) &
233 0 : + delta*(f(3, k_old) + delta*f(4, k_old)))
234 :
235 : slopes(k_new) = &
236 0 : f(2, k_old) + 2*delta*(f(3, k_old) + 1.5*delta*f(4, k_old))
237 :
238 : end do
239 :
240 :
241 0 : end subroutine do_interp_values_and_slopes_sg
242 :
243 :
244 8 : real function minmod1_sg(f1, f2)
245 : real, intent(in) :: f1, f2
246 8 : minmod1_sg = 0.5 * (sign(1.0, f1) + sign(1.0, f2)) * min(abs(f1), abs(f2))
247 8 : end function minmod1_sg
248 :
249 :
250 0 : real function median1_sg(f1, f2, f3)
251 : real, intent(in) :: f1, f2, f3
252 0 : median1_sg = f1 + minmod1_sg(f2 - f1, f3 - f1)
253 0 : end function median1_sg
254 :
255 :
256 28 : subroutine minmod_sg(z, n, f1, f2)
257 : real, intent(inout) :: z(:) ! (n)
258 : integer, intent(in) :: n ! length of vectors
259 : real, intent(in) :: f1(:), f2(:) ! (n)
260 1986 : z(1:n) = 0.5 * (sign(1.0, f1(1:n)) + sign(1.0, f2(1:n))) * min(abs(f1(1:n)), abs(f2(1:n)))
261 28 : end subroutine minmod_sg
262 :
263 :
264 4 : subroutine minmod4_sg(z, n, f1, f2, f3, f4)
265 : real, intent(inout) :: z(:) ! (n)
266 : integer, intent(in) :: n ! length of vectors
267 : real, intent(in) :: f1(:), f2(:), f3(:), f4(:)
268 4 : call minmod_sg(z, n, f1, f2)
269 4 : call minmod_sg(z, n, z, f3)
270 4 : call minmod_sg(z, n, z, f4)
271 4 : end subroutine minmod4_sg
272 :
273 :
274 4 : subroutine median_sg(z, n, f1, f2, f3)
275 : real, intent(out) :: z(:)
276 : integer, intent(in) :: n ! length of vectors
277 : real, intent(in) :: f1(:), f2(:), f3(:)
278 564 : real, target :: tmp1_ary(n), tmp2_ary(n)
279 4 : real, pointer :: tmp1(:), tmp2(:)
280 4 : tmp1 => tmp1_ary
281 4 : tmp2 => tmp2_ary
282 282 : tmp1(1:n) = f2(1:n) - f1(1:n)
283 282 : tmp2(1:n) = f3(1:n) - f1(1:n)
284 4 : call minmod_sg(z(1:n), n, tmp1(1:n), tmp2(1:n))
285 282 : z(1:n) = z(1:n) + f1(1:n)
286 4 : end subroutine median_sg
287 :
288 :
289 : end module interp_1d_misc_sg
|
