Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2026 Niall Miller & 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 : ! linear interpolation for SEDs
21 : !
22 : ! data-loading strategy selected by rq%cube_loaded:
23 : ! .true. -> use the preloaded 4-D flux cube on the handle
24 : ! .false. -> load individual SED files via the lookup table (fallback)
25 :
26 : module linear_interp
27 : use const_def, only: dp
28 : use colors_def, only: Colors_General_Info
29 : use colors_utils, only: dilute_flux, find_containing_cell, find_interval, &
30 : find_nearest_point, find_bracket_index, &
31 : load_sed_cached, load_stencil
32 : use utils_lib, only: mesa_error
33 : implicit none
34 :
35 : private
36 : public :: construct_sed_linear, trilinear_interp
37 :
38 : contains
39 :
40 : ! main entry point -- construct a SED using trilinear interpolation
41 : ! strategy controlled by rq%cube_loaded (set at init)
42 0 : subroutine construct_sed_linear(rq, teff, log_g, metallicity, R, d, &
43 : stellar_model_dir, wavelengths, fluxes)
44 : type(Colors_General_Info), intent(inout) :: rq
45 : real(dp), intent(in) :: teff, log_g, metallicity, R, d
46 : character(len=*), intent(in) :: stellar_model_dir
47 : real(dp), dimension(:), allocatable, intent(out) :: wavelengths, fluxes
48 :
49 : integer :: n_lambda
50 0 : real(dp), dimension(:), allocatable :: interp_flux, diluted_flux
51 :
52 0 : if (rq%cube_loaded) then
53 : ! fast path: use preloaded cube from handle
54 0 : n_lambda = size(rq%cube_wavelengths)
55 :
56 0 : allocate (wavelengths(n_lambda))
57 0 : wavelengths = rq%cube_wavelengths
58 :
59 : ! vectorised interpolation -- cell location computed once, reused across n_lambda
60 0 : allocate (interp_flux(n_lambda))
61 : call trilinear_interp_vector(teff, log_g, metallicity, &
62 : rq%cube_teff_grid, rq%cube_logg_grid, &
63 : rq%cube_meta_grid, &
64 0 : rq%cube_flux, n_lambda, interp_flux)
65 : else
66 : ! fallback path: load individual SED files from lookup table
67 : call construct_sed_from_files(rq, teff, log_g, metallicity, &
68 0 : stellar_model_dir, interp_flux, wavelengths)
69 0 : n_lambda = size(wavelengths)
70 : end if
71 :
72 0 : allocate (diluted_flux(n_lambda))
73 0 : call dilute_flux(interp_flux, R, d, diluted_flux)
74 0 : fluxes = diluted_flux
75 :
76 0 : end subroutine construct_sed_linear
77 :
78 : ! fallback: build a 2x2x2 sub-cube from SED files, then trilinear-interpolate
79 : ! unlike hermite, no derivative context needed -- stencil is exactly the 2x2x2 cell corners
80 0 : subroutine construct_sed_from_files(rq, teff, log_g, metallicity, &
81 : stellar_model_dir, interp_flux, wavelengths)
82 : use colors_utils, only: resolve_path, build_grid_to_lu_map
83 : type(Colors_General_Info), intent(inout) :: rq
84 : real(dp), intent(in) :: teff, log_g, metallicity
85 : character(len=*), intent(in) :: stellar_model_dir
86 : real(dp), dimension(:), allocatable, intent(out) :: interp_flux, wavelengths
87 :
88 : integer :: i_t, i_g, i_m ! bracketing indices in unique grids
89 : integer :: lo_t, hi_t, lo_g, hi_g, lo_m, hi_m ! stencil bounds
90 : integer :: nt, ng, nm, n_lambda
91 : character(len=512) :: resolved_dir
92 : logical :: need_reload
93 :
94 0 : resolved_dir = trim(resolve_path(stellar_model_dir))
95 :
96 : ! ensure the grid-to-lu mapping exists (built once, then reused)
97 0 : if (.not. rq%grid_map_built) call build_grid_to_lu_map(rq)
98 :
99 : ! find bracketing cell in the unique grids
100 0 : call find_bracket_index(rq%u_teff, teff, i_t)
101 0 : call find_bracket_index(rq%u_logg, log_g, i_g)
102 0 : call find_bracket_index(rq%u_meta, metallicity, i_m)
103 :
104 : ! check if the stencil cache is still valid for this cell
105 0 : need_reload = .true.
106 : if (rq%stencil_valid .and. &
107 : i_t == rq%stencil_i_t .and. &
108 0 : i_g == rq%stencil_i_g .and. &
109 : i_m == rq%stencil_i_m) then
110 : need_reload = .false.
111 : end if
112 :
113 : if (need_reload) then
114 : ! trilinear needs exactly the 2x2x2 cell corners -- no extension
115 0 : nt = size(rq%u_teff)
116 0 : ng = size(rq%u_logg)
117 0 : nm = size(rq%u_meta)
118 :
119 0 : if (nt < 2) then
120 0 : lo_t = 1; hi_t = 1
121 : else
122 0 : lo_t = i_t
123 0 : hi_t = min(nt, i_t + 1)
124 : end if
125 :
126 0 : if (ng < 2) then
127 0 : lo_g = 1; hi_g = 1
128 : else
129 0 : lo_g = i_g
130 0 : hi_g = min(ng, i_g + 1)
131 : end if
132 :
133 0 : if (nm < 2) then
134 0 : lo_m = 1; hi_m = 1
135 : else
136 0 : lo_m = i_m
137 0 : hi_m = min(nm, i_m + 1)
138 : end if
139 :
140 : ! load SEDs for every stencil point (using memory cache)
141 0 : call load_stencil(rq, resolved_dir, lo_t, hi_t, lo_g, hi_g, lo_m, hi_m)
142 :
143 : ! store subgrid arrays on the handle
144 0 : if (allocated(rq%stencil_teff)) deallocate (rq%stencil_teff)
145 0 : if (allocated(rq%stencil_logg)) deallocate (rq%stencil_logg)
146 0 : if (allocated(rq%stencil_meta)) deallocate (rq%stencil_meta)
147 :
148 0 : allocate (rq%stencil_teff(hi_t - lo_t + 1))
149 0 : allocate (rq%stencil_logg(hi_g - lo_g + 1))
150 0 : allocate (rq%stencil_meta(hi_m - lo_m + 1))
151 0 : rq%stencil_teff = rq%u_teff(lo_t:hi_t)
152 0 : rq%stencil_logg = rq%u_logg(lo_g:hi_g)
153 0 : rq%stencil_meta = rq%u_meta(lo_m:hi_m)
154 :
155 0 : rq%stencil_i_t = i_t
156 0 : rq%stencil_i_g = i_g
157 0 : rq%stencil_i_m = i_m
158 0 : rq%stencil_valid = .true.
159 : end if
160 :
161 0 : n_lambda = size(rq%stencil_wavelengths)
162 0 : allocate (wavelengths(n_lambda))
163 0 : wavelengths = rq%stencil_wavelengths
164 :
165 0 : allocate (interp_flux(n_lambda))
166 : call trilinear_interp_vector(teff, log_g, metallicity, &
167 : rq%stencil_teff, rq%stencil_logg, rq%stencil_meta, &
168 0 : rq%stencil_fluxes, n_lambda, interp_flux)
169 :
170 0 : end subroutine construct_sed_from_files
171 :
172 : ! vectorised trilinear interpolation over all wavelengths
173 : ! cell location depends only on (teff, logg, meta) -- computed once, reused across n_lambda
174 0 : subroutine trilinear_interp_vector(x_val, y_val, z_val, &
175 0 : x_grid, y_grid, z_grid, &
176 0 : f_values_4d, n_lambda, result_flux)
177 : real(dp), intent(in) :: x_val, y_val, z_val
178 : real(dp), intent(in) :: x_grid(:), y_grid(:), z_grid(:)
179 : real(dp), intent(in) :: f_values_4d(:, :, :, :) ! (nx, ny, nz, n_lambda)
180 : integer, intent(in) :: n_lambda
181 : real(dp), intent(out) :: result_flux(n_lambda)
182 :
183 : integer :: i_x, i_y, i_z, lam
184 : real(dp) :: t_x, t_y, t_z
185 : integer :: nx, ny, nz
186 : real(dp) :: c000, c001, c010, c011, c100, c101, c110, c111
187 : real(dp) :: c00, c01, c10, c11, c0, c1
188 : real(dp) :: lin_result, log_result
189 : real(dp), parameter :: tiny_value = 1.0e-10_dp
190 :
191 0 : nx = size(x_grid)
192 0 : ny = size(y_grid)
193 0 : nz = size(z_grid)
194 :
195 : ! locate the cell once
196 : call find_containing_cell(x_val, y_val, z_val, x_grid, y_grid, z_grid, &
197 0 : i_x, i_y, i_z, t_x, t_y, t_z)
198 :
199 : ! boundary safety check
200 0 : if (i_x < 1) i_x = 1
201 0 : if (i_y < 1) i_y = 1
202 0 : if (i_z < 1) i_z = 1
203 0 : if (i_x >= nx) i_x = max(1, nx - 1)
204 0 : if (i_y >= ny) i_y = max(1, ny - 1)
205 0 : if (i_z >= nz) i_z = max(1, nz - 1)
206 :
207 : ! clamp interpolation parameters to [0,1]
208 0 : t_x = max(0.0_dp, min(1.0_dp, t_x))
209 0 : t_y = max(0.0_dp, min(1.0_dp, t_y))
210 0 : t_z = max(0.0_dp, min(1.0_dp, t_z))
211 :
212 : ! loop over wavelengths with the same cell location
213 0 : do lam = 1, n_lambda
214 : ! get the 8 corners of the cube with safety floor
215 0 : c000 = max(tiny_value, f_values_4d(i_x, i_y, i_z, lam))
216 0 : c001 = max(tiny_value, f_values_4d(i_x, i_y, i_z + 1, lam))
217 0 : c010 = max(tiny_value, f_values_4d(i_x, i_y + 1, i_z, lam))
218 0 : c011 = max(tiny_value, f_values_4d(i_x, i_y + 1, i_z + 1, lam))
219 0 : c100 = max(tiny_value, f_values_4d(i_x + 1, i_y, i_z, lam))
220 0 : c101 = max(tiny_value, f_values_4d(i_x + 1, i_y, i_z + 1, lam))
221 0 : c110 = max(tiny_value, f_values_4d(i_x + 1, i_y + 1, i_z, lam))
222 0 : c111 = max(tiny_value, f_values_4d(i_x + 1, i_y + 1, i_z + 1, lam))
223 :
224 : ! standard linear interpolation first (safer)
225 0 : c00 = c000*(1.0_dp - t_x) + c100*t_x
226 0 : c01 = c001*(1.0_dp - t_x) + c101*t_x
227 0 : c10 = c010*(1.0_dp - t_x) + c110*t_x
228 0 : c11 = c011*(1.0_dp - t_x) + c111*t_x
229 :
230 0 : c0 = c00*(1.0_dp - t_y) + c10*t_y
231 0 : c1 = c01*(1.0_dp - t_y) + c11*t_y
232 :
233 0 : lin_result = c0*(1.0_dp - t_z) + c1*t_z
234 :
235 : ! if valid, try log-space interpolation (smoother for flux)
236 0 : if (lin_result > tiny_value) then
237 0 : c00 = log(c000)*(1.0_dp - t_x) + log(c100)*t_x
238 0 : c01 = log(c001)*(1.0_dp - t_x) + log(c101)*t_x
239 0 : c10 = log(c010)*(1.0_dp - t_x) + log(c110)*t_x
240 0 : c11 = log(c011)*(1.0_dp - t_x) + log(c111)*t_x
241 :
242 0 : c0 = c00*(1.0_dp - t_y) + c10*t_y
243 0 : c1 = c01*(1.0_dp - t_y) + c11*t_y
244 :
245 0 : log_result = c0*(1.0_dp - t_z) + c1*t_z
246 :
247 : ! only use log-space result if valid
248 0 : if (log_result == log_result) then ! NaN check
249 0 : lin_result = exp(log_result)
250 : end if
251 : end if
252 :
253 : ! final sanity check -- fall back to nearest neighbour
254 0 : if (lin_result /= lin_result .or. lin_result <= 0.0_dp) then
255 : call find_nearest_point(x_val, y_val, z_val, x_grid, y_grid, z_grid, &
256 0 : i_x, i_y, i_z)
257 0 : lin_result = max(tiny_value, f_values_4d(i_x, i_y, i_z, lam))
258 : end if
259 :
260 0 : result_flux(lam) = lin_result
261 : end do
262 :
263 0 : end subroutine trilinear_interp_vector
264 :
265 : ! scalar trilinear interpolation (external callers / single-wavelength use)
266 : ! retained for backward compatibility
267 0 : function trilinear_interp(x_val, y_val, z_val, x_grid, y_grid, z_grid, f_values) result(f_interp)
268 : real(dp), intent(in) :: x_val, y_val, z_val
269 : real(dp), intent(in) :: x_grid(:), y_grid(:), z_grid(:)
270 : real(dp), intent(in) :: f_values(:, :, :)
271 : real(dp) :: f_interp
272 : real(dp) :: log_result
273 : integer :: i_x, i_y, i_z
274 : real(dp) :: t_x, t_y, t_z
275 : real(dp) :: c000, c001, c010, c011, c100, c101, c110, c111
276 : real(dp) :: c00, c01, c10, c11, c0, c1
277 : real(dp), parameter :: tiny_value = 1.0e-10_dp
278 :
279 : ! find containing cell and parameter values using binary search
280 : call find_containing_cell(x_val, y_val, z_val, x_grid, y_grid, z_grid, &
281 0 : i_x, i_y, i_z, t_x, t_y, t_z)
282 :
283 : ! boundary safety check
284 0 : if (i_x < lbound(x_grid, 1)) i_x = lbound(x_grid, 1)
285 0 : if (i_y < lbound(y_grid, 1)) i_y = lbound(y_grid, 1)
286 0 : if (i_z < lbound(z_grid, 1)) i_z = lbound(z_grid, 1)
287 0 : if (i_x >= ubound(x_grid, 1)) i_x = ubound(x_grid, 1) - 1
288 0 : if (i_y >= ubound(y_grid, 1)) i_y = ubound(y_grid, 1) - 1
289 0 : if (i_z >= ubound(z_grid, 1)) i_z = ubound(z_grid, 1) - 1
290 :
291 : ! clamp interpolation parameters to [0,1]
292 0 : t_x = max(0.0_dp, MIN(1.0_dp, t_x))
293 0 : t_y = max(0.0_dp, MIN(1.0_dp, t_y))
294 0 : t_z = max(0.0_dp, MIN(1.0_dp, t_z))
295 :
296 : ! get the corners of the cube with safety checks
297 0 : c000 = max(tiny_value, f_values(i_x, i_y, i_z))
298 0 : c001 = max(tiny_value, f_values(i_x, i_y, i_z + 1))
299 0 : c010 = max(tiny_value, f_values(i_x, i_y + 1, i_z))
300 0 : c011 = max(tiny_value, f_values(i_x, i_y + 1, i_z + 1))
301 0 : c100 = max(tiny_value, f_values(i_x + 1, i_y, i_z))
302 0 : c101 = max(tiny_value, f_values(i_x + 1, i_y, i_z + 1))
303 0 : c110 = max(tiny_value, f_values(i_x + 1, i_y + 1, i_z))
304 0 : c111 = max(tiny_value, f_values(i_x + 1, i_y + 1, i_z + 1))
305 :
306 : ! try standard linear interpolation first (safer)
307 0 : c00 = c000*(1.0_dp - t_x) + c100*t_x
308 0 : c01 = c001*(1.0_dp - t_x) + c101*t_x
309 0 : c10 = c010*(1.0_dp - t_x) + c110*t_x
310 0 : c11 = c011*(1.0_dp - t_x) + c111*t_x
311 :
312 0 : c0 = c00*(1.0_dp - t_y) + c10*t_y
313 0 : c1 = c01*(1.0_dp - t_y) + c11*t_y
314 :
315 0 : f_interp = c0*(1.0_dp - t_z) + c1*t_z
316 :
317 : ! if valid, try log-space interpolation (smoother for flux)
318 0 : if (f_interp > tiny_value) then
319 0 : c00 = log(c000)*(1.0_dp - t_x) + log(c100)*t_x
320 0 : c01 = log(c001)*(1.0_dp - t_x) + log(c101)*t_x
321 0 : c10 = log(c010)*(1.0_dp - t_x) + log(c110)*t_x
322 0 : c11 = log(c011)*(1.0_dp - t_x) + log(c111)*t_x
323 :
324 0 : c0 = c00*(1.0_dp - t_y) + c10*t_y
325 0 : c1 = c01*(1.0_dp - t_y) + c11*t_y
326 :
327 0 : log_result = c0*(1.0_dp - t_z) + c1*t_z
328 :
329 : ! only use the log-space result if it's valid
330 0 : if (log_result == log_result) then ! NaN check
331 0 : f_interp = EXP(log_result)
332 : end if
333 : end if
334 :
335 : ! final sanity check
336 0 : if (f_interp /= f_interp .or. f_interp <= 0.0_dp) then
337 0 : call find_nearest_point(x_val, y_val, z_val, x_grid, y_grid, z_grid, i_x, i_y, i_z)
338 0 : f_interp = max(tiny_value, f_values(i_x, i_y, i_z))
339 : end if
340 0 : end function trilinear_interp
341 :
342 : end module linear_interp
|
