Line data Source code
1 : ! ***********************************************************************
2 : ! Copyright (C) 2026 Niall Miller & The MESA Team
3 : ! ***********************************************************************
4 :
5 : ! ***********************************************************************
6 : ! Hermite interpolation module for spectral energy distributions (SEDs)
7 : ! ***********************************************************************
8 :
9 : module hermite_interp
10 : use const_def, only: dp
11 : use colors_def, only: Colors_General_Info
12 : use colors_utils, only: dilute_flux, find_containing_cell, &
13 : find_nearest_point, find_bracket_index, load_stencil
14 : implicit none
15 :
16 : private
17 : public :: construct_sed_hermite, hermite_tensor_interp3d
18 :
19 : contains
20 :
21 : !---------------------------------------------------------------------------
22 : ! Main entry point: Construct a SED using Hermite tensor interpolation.
23 : ! Data loading strategy is determined by rq%cube_loaded (set at init):
24 : ! cube_loaded = .true. -> use the preloaded 4-D cube on the handle
25 : ! cube_loaded = .false. -> load individual SED files via the lookup table
26 : !---------------------------------------------------------------------------
27 7 : subroutine construct_sed_hermite(rq, teff, log_g, metallicity, R, d, &
28 : stellar_model_dir, wavelengths, fluxes)
29 : type(Colors_General_Info), intent(inout) :: rq
30 : real(dp), intent(in) :: teff, log_g, metallicity, R, d
31 : character(len=*), intent(in) :: stellar_model_dir
32 : real(dp), dimension(:), allocatable, intent(out) :: wavelengths, fluxes
33 :
34 : integer :: n_lambda
35 7 : real(dp), dimension(:), allocatable :: interp_flux
36 :
37 7 : if (rq%cube_loaded) then
38 : ! ---- Fast path: use preloaded cube from handle ----
39 7 : n_lambda = size(rq%cube_wavelengths)
40 :
41 : ! Copy wavelengths to output
42 21 : allocate (wavelengths(n_lambda))
43 4819 : wavelengths = rq%cube_wavelengths
44 :
45 : ! Vectorised interpolation over all wavelengths in one pass —
46 : ! cell location is computed once and reused, no per-wavelength
47 : ! allocation or 3-D slice extraction needed.
48 14 : allocate (interp_flux(n_lambda))
49 : call hermite_interp_vector(teff, log_g, metallicity, &
50 : rq%cube_teff_grid, rq%cube_logg_grid, &
51 : rq%cube_meta_grid, &
52 7 : rq%cube_flux, n_lambda, interp_flux)
53 : else
54 : ! ---- Fallback path: load individual SED files from lookup table ----
55 : call construct_sed_from_files(rq, teff, log_g, metallicity, &
56 0 : stellar_model_dir, interp_flux, wavelengths)
57 0 : n_lambda = size(wavelengths)
58 : end if
59 :
60 : ! Apply distance dilution to get observed flux
61 21 : allocate (fluxes(n_lambda))
62 7 : call dilute_flux(interp_flux, R, d, fluxes)
63 :
64 7 : end subroutine construct_sed_hermite
65 :
66 : !---------------------------------------------------------------------------
67 : ! Fallback: Build a local sub-cube from individual SED files with enough
68 : ! context for Hermite derivative computation, then interpolate all
69 : ! wavelengths in a single pass.
70 : !---------------------------------------------------------------------------
71 0 : subroutine construct_sed_from_files(rq, teff, log_g, metallicity, &
72 : stellar_model_dir, interp_flux, wavelengths)
73 : use colors_utils, only: resolve_path, build_grid_to_lu_map
74 : type(Colors_General_Info), intent(inout) :: rq
75 : real(dp), intent(in) :: teff, log_g, metallicity
76 : character(len=*), intent(in) :: stellar_model_dir
77 : real(dp), dimension(:), allocatable, intent(out) :: interp_flux, wavelengths
78 :
79 : integer :: i_t, i_g, i_m ! bracketing indices in unique grids
80 : integer :: lo_t, hi_t, lo_g, hi_g, lo_m, hi_m ! stencil bounds
81 : integer :: nt, ng, nm, n_lambda
82 : character(len=512) :: resolved_dir
83 : logical :: need_reload
84 :
85 0 : resolved_dir = trim(resolve_path(stellar_model_dir))
86 :
87 : ! Ensure the grid-to-lu mapping exists (built once, then reused)
88 0 : if (.not. rq%grid_map_built) call build_grid_to_lu_map(rq)
89 :
90 : ! Find bracketing cell in the unique grids
91 0 : call find_bracket_index(rq%u_teff, teff, i_t)
92 0 : call find_bracket_index(rq%u_logg, log_g, i_g)
93 0 : call find_bracket_index(rq%u_meta, metallicity, i_m)
94 :
95 : ! Check if the stencil cache is still valid for this cell
96 0 : need_reload = .true.
97 : if (rq%stencil_valid .and. &
98 : i_t == rq%stencil_i_t .and. &
99 0 : i_g == rq%stencil_i_g .and. &
100 : i_m == rq%stencil_i_m) then
101 : need_reload = .false.
102 : end if
103 :
104 : if (need_reload) then
105 : ! Determine the extended stencil bounds:
106 : ! For each axis, include one point before and after the cell
107 : ! when available, so that centred differences match the cube.
108 0 : nt = size(rq%u_teff)
109 0 : ng = size(rq%u_logg)
110 0 : nm = size(rq%u_meta)
111 :
112 0 : if (nt < 2) then
113 0 : lo_t = 1; hi_t = 1
114 : else
115 0 : lo_t = max(1, i_t - 1)
116 0 : hi_t = min(nt, i_t + 2)
117 : end if
118 :
119 0 : if (ng < 2) then
120 0 : lo_g = 1; hi_g = 1
121 : else
122 0 : lo_g = max(1, i_g - 1)
123 0 : hi_g = min(ng, i_g + 2)
124 : end if
125 :
126 0 : if (nm < 2) then
127 0 : lo_m = 1; hi_m = 1
128 : else
129 0 : lo_m = max(1, i_m - 1)
130 0 : hi_m = min(nm, i_m + 2)
131 : end if
132 :
133 : ! Load SEDs for every stencil point (using memory cache)
134 0 : call load_stencil(rq, resolved_dir, lo_t, hi_t, lo_g, hi_g, lo_m, hi_m)
135 :
136 : ! Store subgrid arrays on the handle
137 0 : if (allocated(rq%stencil_teff)) deallocate (rq%stencil_teff)
138 0 : if (allocated(rq%stencil_logg)) deallocate (rq%stencil_logg)
139 0 : if (allocated(rq%stencil_meta)) deallocate (rq%stencil_meta)
140 :
141 0 : allocate (rq%stencil_teff(hi_t - lo_t + 1))
142 0 : allocate (rq%stencil_logg(hi_g - lo_g + 1))
143 0 : allocate (rq%stencil_meta(hi_m - lo_m + 1))
144 0 : rq%stencil_teff = rq%u_teff(lo_t:hi_t)
145 0 : rq%stencil_logg = rq%u_logg(lo_g:hi_g)
146 0 : rq%stencil_meta = rq%u_meta(lo_m:hi_m)
147 :
148 0 : rq%stencil_i_t = i_t
149 0 : rq%stencil_i_g = i_g
150 0 : rq%stencil_i_m = i_m
151 0 : rq%stencil_valid = .true.
152 : end if
153 :
154 : ! Copy wavelengths to output
155 0 : n_lambda = size(rq%stencil_wavelengths)
156 0 : allocate (wavelengths(n_lambda))
157 0 : wavelengths = rq%stencil_wavelengths
158 :
159 : ! Interpolate all wavelengths using precomputed stencil
160 0 : allocate (interp_flux(n_lambda))
161 : call hermite_interp_vector(teff, log_g, metallicity, &
162 : rq%stencil_teff, rq%stencil_logg, rq%stencil_meta, &
163 0 : rq%stencil_fluxes, n_lambda, interp_flux)
164 :
165 0 : end subroutine construct_sed_from_files
166 :
167 : !---------------------------------------------------------------------------
168 : ! Vectorised Hermite interpolation over all wavelengths.
169 : !
170 : ! The cell location (i_x, i_y, i_z, t_x, t_y, t_z) depends only on
171 : ! (teff, logg, meta) and the sub-grids — not on wavelength. Computing
172 : ! it once and reusing across all n_lambda samples eliminates redundant
173 : ! binary searches and basis-function evaluations.
174 : !---------------------------------------------------------------------------
175 7 : subroutine hermite_interp_vector(x_val, y_val, z_val, &
176 7 : x_grid, y_grid, z_grid, &
177 7 : f_values_4d, n_lambda, result_flux)
178 : real(dp), intent(in) :: x_val, y_val, z_val
179 : real(dp), intent(in) :: x_grid(:), y_grid(:), z_grid(:)
180 : real(dp), intent(in) :: f_values_4d(:, :, :, :) ! (nx, ny, nz, n_lambda)
181 : integer, intent(in) :: n_lambda
182 : real(dp), intent(out) :: result_flux(n_lambda)
183 :
184 : integer :: i_x, i_y, i_z
185 : integer :: ix_max, iy_max, iz_max
186 : real(dp) :: t_x, t_y, t_z
187 : real(dp) :: dx, dy, dz
188 : real(dp) :: val, df_dx, df_dy, df_dz
189 : integer :: nx, ny, nz
190 : integer :: ix, iy, iz, lam
191 : real(dp) :: h_x(2), h_y(2), h_z(2)
192 : real(dp) :: hx_d(2), hy_d(2), hz_d(2)
193 : real(dp) :: wx, wy, wz, wxd, wyd, wzd
194 :
195 7 : nx = size(x_grid)
196 7 : ny = size(y_grid)
197 7 : nz = size(z_grid)
198 :
199 : ! Find containing cell (done once for all wavelengths). Singleton axes
200 : ! are valid degenerate dimensions, not an out-of-grid condition.
201 : call find_containing_cell(x_val, y_val, z_val, x_grid, y_grid, z_grid, &
202 7 : i_x, i_y, i_z, t_x, t_y, t_z)
203 :
204 7 : ix_max = 1
205 7 : iy_max = 1
206 7 : iz_max = 1
207 :
208 7 : if (nx == 1) then
209 0 : i_x = 1
210 0 : t_x = 0.0_dp
211 0 : ix_max = 0
212 : end if
213 7 : if (ny == 1) then
214 3 : i_y = 1
215 3 : t_y = 0.0_dp
216 3 : iy_max = 0
217 : end if
218 7 : if (nz == 1) then
219 3 : i_z = 1
220 3 : t_z = 0.0_dp
221 3 : iz_max = 0
222 : end if
223 :
224 : ! If outside a non-degenerate grid dimension, use nearest point for all wavelengths.
225 : if ((nx > 1 .and. (i_x < 1 .or. i_x >= nx)) .or. &
226 7 : (ny > 1 .and. (i_y < 1 .or. i_y >= ny)) .or. &
227 : (nz > 1 .and. (i_z < 1 .or. i_z >= nz))) then
228 :
229 : call find_nearest_point(x_val, y_val, z_val, x_grid, y_grid, z_grid, &
230 0 : i_x, i_y, i_z)
231 0 : do lam = 1, n_lambda
232 0 : result_flux(lam) = f_values_4d(i_x, i_y, i_z, lam)
233 : end do
234 0 : return
235 : end if
236 :
237 : ! Grid cell spacing. Degenerate axes have zero derivative contribution.
238 7 : if (nx > 1) then
239 7 : dx = x_grid(i_x + 1) - x_grid(i_x)
240 : else
241 : dx = 0.0_dp
242 : end if
243 7 : if (ny > 1) then
244 4 : dy = y_grid(i_y + 1) - y_grid(i_y)
245 : else
246 : dy = 0.0_dp
247 : end if
248 7 : if (nz > 1) then
249 4 : dz = z_grid(i_z + 1) - z_grid(i_z)
250 : else
251 : dz = 0.0_dp
252 : end if
253 :
254 : ! Precompute Hermite basis functions (same for all wavelengths)
255 21 : h_x = [h00(t_x), h01(t_x)]
256 21 : hx_d = [h10(t_x), h11(t_x)]
257 21 : h_y = [h00(t_y), h01(t_y)]
258 21 : hy_d = [h10(t_y), h11(t_y)]
259 21 : h_z = [h00(t_z), h01(t_z)]
260 21 : hz_d = [h10(t_z), h11(t_z)]
261 :
262 : ! stencil loop -- weights are invariant over lambda, so lambda is innermost
263 4812 : result_flux = 0.0_dp
264 18 : do iz = 0, iz_max
265 11 : wz = h_z(iz + 1)
266 11 : wzd = hz_d(iz + 1)
267 37 : do iy = 0, iy_max
268 19 : wy = h_y(iy + 1)
269 19 : wyd = hy_d(iy + 1)
270 68 : do ix = 0, ix_max
271 38 : wx = h_x(ix + 1)
272 38 : wxd = hx_d(ix + 1)
273 38443 : do lam = 1, n_lambda
274 38386 : val = f_values_4d(i_x + ix, i_y + iy, i_z + iz, lam)
275 :
276 : call compute_derivatives_at_point_4d( &
277 : f_values_4d, i_x + ix, i_y + iy, i_z + iz, lam, &
278 38386 : nx, ny, nz, x_grid, y_grid, z_grid, df_dx, df_dy, df_dz)
279 :
280 : result_flux(lam) = result_flux(lam) &
281 : + wx*wy*wz*val &
282 : + wxd*wy*wz*dx*df_dx &
283 : + wx*wyd*wz*dy*df_dy &
284 38424 : + wx*wy*wzd*dz*df_dz
285 : end do
286 : end do
287 : end do
288 : end do
289 :
290 : end subroutine hermite_interp_vector
291 :
292 :
293 115158 : subroutine compute_derivatives_at_point_4d(f4d, i, j, k, lam, nx, ny, nz, &
294 38386 : x_grid, y_grid, z_grid, df_dx, df_dy, df_dz)
295 : real(dp), intent(in) :: f4d(:, :, :, :)
296 : integer, intent(in) :: i, j, k, lam, nx, ny, nz
297 : real(dp), intent(in) :: x_grid(:), y_grid(:), z_grid(:)
298 : real(dp), intent(out) :: df_dx, df_dy, df_dz
299 : real(dp) :: h_minus, h_plus, d_minus, d_plus
300 :
301 : ! x derivative
302 38386 : if (nx == 1) then
303 0 : df_dx = 0.0_dp
304 38386 : else if (i > 1 .and. i < nx) then
305 38380 : h_minus = x_grid(i) - x_grid(i - 1)
306 38380 : h_plus = x_grid(i + 1) - x_grid(i)
307 38380 : d_minus = (f4d(i, j, k, lam) - f4d(i - 1, j, k, lam))/h_minus
308 38380 : d_plus = (f4d(i + 1, j, k, lam) - f4d(i, j, k, lam))/h_plus
309 38380 : df_dx = (h_plus*d_minus + h_minus*d_plus)/(h_minus + h_plus)
310 6 : else if (i == 1) then
311 : df_dx = (f4d(i + 1, j, k, lam) - f4d(i, j, k, lam))/ &
312 6 : (x_grid(i + 1) - x_grid(i))
313 : else
314 : df_dx = (f4d(i, j, k, lam) - f4d(i - 1, j, k, lam))/ &
315 0 : (x_grid(i) - x_grid(i - 1))
316 : end if
317 :
318 : ! y derivative
319 38386 : if (ny == 1) then
320 18 : df_dy = 0.0_dp
321 38368 : else if (j > 1 .and. j < ny) then
322 38368 : h_minus = y_grid(j) - y_grid(j - 1)
323 38368 : h_plus = y_grid(j + 1) - y_grid(j)
324 38368 : d_minus = (f4d(i, j, k, lam) - f4d(i, j - 1, k, lam))/h_minus
325 38368 : d_plus = (f4d(i, j + 1, k, lam) - f4d(i, j, k, lam))/h_plus
326 38368 : df_dy = (h_plus*d_minus + h_minus*d_plus)/(h_minus + h_plus)
327 0 : else if (j == 1) then
328 : df_dy = (f4d(i, j + 1, k, lam) - f4d(i, j, k, lam))/ &
329 0 : (y_grid(j + 1) - y_grid(j))
330 : else
331 : df_dy = (f4d(i, j, k, lam) - f4d(i, j - 1, k, lam))/ &
332 0 : (y_grid(j) - y_grid(j - 1))
333 : end if
334 :
335 : ! z derivative
336 38386 : if (nz == 1) then
337 18 : df_dz = 0.0_dp
338 38368 : else if (k > 1 .and. k < nz) then
339 38368 : h_minus = z_grid(k) - z_grid(k - 1)
340 38368 : h_plus = z_grid(k + 1) - z_grid(k)
341 38368 : d_minus = (f4d(i, j, k, lam) - f4d(i, j, k - 1, lam))/h_minus
342 38368 : d_plus = (f4d(i, j, k + 1, lam) - f4d(i, j, k, lam))/h_plus
343 38368 : df_dz = (h_plus*d_minus + h_minus*d_plus)/(h_minus + h_plus)
344 0 : else if (k == 1) then
345 : df_dz = (f4d(i, j, k + 1, lam) - f4d(i, j, k, lam))/ &
346 0 : (z_grid(k + 1) - z_grid(k))
347 : else
348 : df_dz = (f4d(i, j, k, lam) - f4d(i, j, k - 1, lam))/ &
349 0 : (z_grid(k) - z_grid(k - 1))
350 : end if
351 :
352 38386 : end subroutine compute_derivatives_at_point_4d
353 :
354 :
355 1 : function hermite_tensor_interp3d(x_val, y_val, z_val, x_grid, y_grid, &
356 1 : z_grid, f_values) result(f_interp)
357 : real(dp), intent(in) :: x_val, y_val, z_val
358 : real(dp), intent(in) :: x_grid(:), y_grid(:), z_grid(:)
359 : real(dp), intent(in) :: f_values(:, :, :)
360 : real(dp) :: f_interp
361 :
362 : integer :: i_x, i_y, i_z, ix, iy, iz
363 : integer :: ix_max, iy_max, iz_max
364 : integer :: nx, ny, nz
365 : real(dp) :: t_x, t_y, t_z, dx, dy, dz, f_sum
366 : real(dp) :: dx_values(2, 2, 2), dy_values(2, 2, 2), dz_values(2, 2, 2)
367 : real(dp) :: values(2, 2, 2)
368 : real(dp) :: h_x(2), h_y(2), h_z(2)
369 : real(dp) :: hx_d(2), hy_d(2), hz_d(2)
370 :
371 1 : nx = size(x_grid)
372 1 : ny = size(y_grid)
373 1 : nz = size(z_grid)
374 :
375 : call find_containing_cell(x_val, y_val, z_val, x_grid, y_grid, z_grid, &
376 1 : i_x, i_y, i_z, t_x, t_y, t_z)
377 :
378 1 : ix_max = 1
379 1 : iy_max = 1
380 1 : iz_max = 1
381 :
382 1 : if (nx == 1) then
383 0 : i_x = 1
384 0 : t_x = 0.0_dp
385 0 : ix_max = 0
386 : end if
387 1 : if (ny == 1) then
388 1 : i_y = 1
389 1 : t_y = 0.0_dp
390 1 : iy_max = 0
391 : end if
392 1 : if (nz == 1) then
393 1 : i_z = 1
394 1 : t_z = 0.0_dp
395 1 : iz_max = 0
396 : end if
397 :
398 : if ((nx > 1 .and. (i_x < 1 .or. i_x >= nx)) .or. &
399 1 : (ny > 1 .and. (i_y < 1 .or. i_y >= ny)) .or. &
400 : (nz > 1 .and. (i_z < 1 .or. i_z >= nz))) then
401 : call find_nearest_point(x_val, y_val, z_val, x_grid, y_grid, z_grid, &
402 0 : i_x, i_y, i_z)
403 0 : f_interp = f_values(i_x, i_y, i_z)
404 0 : return
405 : end if
406 :
407 1 : if (nx > 1) then
408 1 : dx = x_grid(i_x + 1) - x_grid(i_x)
409 : else
410 : dx = 0.0_dp
411 : end if
412 1 : if (ny > 1) then
413 0 : dy = y_grid(i_y + 1) - y_grid(i_y)
414 : else
415 : dy = 0.0_dp
416 : end if
417 1 : if (nz > 1) then
418 0 : dz = z_grid(i_z + 1) - z_grid(i_z)
419 : else
420 : dz = 0.0_dp
421 : end if
422 :
423 2 : do iz = 0, iz_max
424 3 : do iy = 0, iy_max
425 4 : do ix = 0, ix_max
426 2 : values(ix + 1, iy + 1, iz + 1) = f_values(i_x + ix, i_y + iy, i_z + iz)
427 : call compute_derivatives_at_point(f_values, i_x + ix, i_y + iy, i_z + iz, &
428 : nx, ny, nz, x_grid, y_grid, z_grid, &
429 : dx_values(ix + 1, iy + 1, iz + 1), &
430 : dy_values(ix + 1, iy + 1, iz + 1), &
431 3 : dz_values(ix + 1, iy + 1, iz + 1))
432 : end do
433 : end do
434 : end do
435 :
436 3 : h_x = [h00(t_x), h01(t_x)]
437 3 : hx_d = [h10(t_x), h11(t_x)]
438 3 : h_y = [h00(t_y), h01(t_y)]
439 3 : hy_d = [h10(t_y), h11(t_y)]
440 3 : h_z = [h00(t_z), h01(t_z)]
441 3 : hz_d = [h10(t_z), h11(t_z)]
442 :
443 1 : f_sum = 0.0_dp
444 2 : do iz = 1, iz_max + 1
445 3 : do iy = 1, iy_max + 1
446 4 : do ix = 1, ix_max + 1
447 2 : f_sum = f_sum + h_x(ix)*h_y(iy)*h_z(iz)*values(ix, iy, iz)
448 2 : f_sum = f_sum + hx_d(ix)*h_y(iy)*h_z(iz)*dx*dx_values(ix, iy, iz)
449 2 : f_sum = f_sum + h_x(ix)*hy_d(iy)*h_z(iz)*dy*dy_values(ix, iy, iz)
450 3 : f_sum = f_sum + h_x(ix)*h_y(iy)*hz_d(iz)*dz*dz_values(ix, iy, iz)
451 : end do
452 : end do
453 : end do
454 :
455 1 : f_interp = f_sum
456 : end function hermite_tensor_interp3d
457 :
458 6 : subroutine compute_derivatives_at_point(f, i, j, k, nx, ny, nz, &
459 2 : x_grid, y_grid, z_grid, &
460 : df_dx, df_dy, df_dz)
461 : real(dp), intent(in) :: f(:, :, :)
462 : integer, intent(in) :: i, j, k, nx, ny, nz
463 : real(dp), intent(in) :: x_grid(:), y_grid(:), z_grid(:)
464 : real(dp), intent(out) :: df_dx, df_dy, df_dz
465 : real(dp) :: h_minus, h_plus, d_minus, d_plus
466 :
467 2 : if (nx == 1) then
468 0 : df_dx = 0.0_dp
469 2 : else if (i > 1 .and. i < nx) then
470 2 : h_minus = x_grid(i) - x_grid(i - 1)
471 2 : h_plus = x_grid(i + 1) - x_grid(i)
472 2 : d_minus = (f(i, j, k) - f(i - 1, j, k))/h_minus
473 2 : d_plus = (f(i + 1, j, k) - f(i, j, k))/h_plus
474 2 : df_dx = (h_plus*d_minus + h_minus*d_plus)/(h_minus + h_plus)
475 0 : else if (i == 1) then
476 : df_dx = (f(i + 1, j, k) - f(i, j, k))/ &
477 0 : (x_grid(i + 1) - x_grid(i))
478 : else
479 : df_dx = (f(i, j, k) - f(i - 1, j, k))/ &
480 0 : (x_grid(i) - x_grid(i - 1))
481 : end if
482 :
483 2 : if (ny == 1) then
484 2 : df_dy = 0.0_dp
485 0 : else if (j > 1 .and. j < ny) then
486 0 : h_minus = y_grid(j) - y_grid(j - 1)
487 0 : h_plus = y_grid(j + 1) - y_grid(j)
488 0 : d_minus = (f(i, j, k) - f(i, j - 1, k))/h_minus
489 0 : d_plus = (f(i, j + 1, k) - f(i, j, k))/h_plus
490 0 : df_dy = (h_plus*d_minus + h_minus*d_plus)/(h_minus + h_plus)
491 0 : else if (j == 1) then
492 : df_dy = (f(i, j + 1, k) - f(i, j, k))/ &
493 0 : (y_grid(j + 1) - y_grid(j))
494 : else
495 : df_dy = (f(i, j, k) - f(i, j - 1, k))/ &
496 0 : (y_grid(j) - y_grid(j - 1))
497 : end if
498 :
499 2 : if (nz == 1) then
500 2 : df_dz = 0.0_dp
501 0 : else if (k > 1 .and. k < nz) then
502 0 : h_minus = z_grid(k) - z_grid(k - 1)
503 0 : h_plus = z_grid(k + 1) - z_grid(k)
504 0 : d_minus = (f(i, j, k) - f(i, j, k - 1))/h_minus
505 0 : d_plus = (f(i, j, k + 1) - f(i, j, k))/h_plus
506 0 : df_dz = (h_plus*d_minus + h_minus*d_plus)/(h_minus + h_plus)
507 0 : else if (k == 1) then
508 : df_dz = (f(i, j, k + 1) - f(i, j, k))/ &
509 0 : (z_grid(k + 1) - z_grid(k))
510 : else
511 : df_dz = (f(i, j, k) - f(i, j, k - 1))/ &
512 0 : (z_grid(k) - z_grid(k - 1))
513 : end if
514 :
515 2 : end subroutine compute_derivatives_at_point
516 :
517 : !---------------------------------------------------------------------------
518 : ! Hermite basis functions
519 : !---------------------------------------------------------------------------
520 24 : function h00(t) result(h)
521 : real(dp), intent(in) :: t
522 : real(dp) :: h
523 24 : h = 2.0_dp*t**3 - 3.0_dp*t**2 + 1.0_dp
524 24 : end function h00
525 :
526 24 : function h10(t) result(h)
527 : real(dp), intent(in) :: t
528 : real(dp) :: h
529 24 : h = t**3 - 2.0_dp*t**2 + t
530 24 : end function h10
531 :
532 24 : function h01(t) result(h)
533 : real(dp), intent(in) :: t
534 : real(dp) :: h
535 24 : h = -2.0_dp*t**3 + 3.0_dp*t**2
536 24 : end function h01
537 :
538 24 : function h11(t) result(h)
539 : real(dp), intent(in) :: t
540 : real(dp) :: h
541 24 : h = t**3 - t**2
542 24 : end function h11
543 :
544 : end module hermite_interp
|
