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: is_inf, is_nan, 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 1 : 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 1 : real(dp), dimension(:), allocatable :: interp_flux, diluted_flux
51 :
52 1 : if (rq%cube_loaded) then
53 : ! fast path: use preloaded cube from handle
54 1 : n_lambda = size(rq%cube_wavelengths)
55 :
56 3 : allocate (wavelengths(n_lambda))
57 5 : wavelengths = rq%cube_wavelengths
58 :
59 : ! vectorised interpolation -- cell location computed once, reused across n_lambda
60 2 : 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 1 : 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 3 : allocate (diluted_flux(n_lambda))
73 1 : call dilute_flux(interp_flux, R, d, diluted_flux)
74 6 : fluxes = diluted_flux
75 :
76 1 : 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 1 : subroutine trilinear_interp_vector(x_val, y_val, z_val, &
175 1 : x_grid, y_grid, z_grid, &
176 1 : 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, nearest_x, nearest_y, nearest_z, lam
184 : integer :: nx, ny, nz
185 : integer :: ix_max, iy_max, iz_max
186 : integer :: ox, oy, oz
187 : real(dp) :: t_x, t_y, t_z
188 : real(dp) :: w_x(0:1), w_y(0:1), w_z(0:1)
189 : real(dp) :: accum, weight
190 : real(dp), parameter :: tiny_value = 1.0d-300
191 :
192 1 : nx = size(x_grid)
193 1 : ny = size(y_grid)
194 1 : nz = size(z_grid)
195 :
196 : if (nx < 1 .or. ny < 1 .or. nz < 1 .or. n_lambda < 1 .or. &
197 : size(f_values_4d, 1) /= nx .or. size(f_values_4d, 2) /= ny .or. &
198 1 : size(f_values_4d, 3) /= nz .or. size(f_values_4d, 4) /= n_lambda) then
199 0 : write (*, *) 'colors: invalid array sizes in trilinear_interp_vector'
200 0 : call mesa_error(__FILE__, __LINE__)
201 : end if
202 :
203 : call find_containing_cell(x_val, y_val, z_val, x_grid, y_grid, z_grid, &
204 1 : i_x, i_y, i_z, t_x, t_y, t_z)
205 :
206 1 : call setup_axis(nx, x_val, x_grid, i_x, t_x, ix_max, w_x)
207 1 : call setup_axis(ny, y_val, y_grid, i_y, t_y, iy_max, w_y)
208 1 : call setup_axis(nz, z_val, z_grid, i_z, t_z, iz_max, w_z)
209 :
210 4 : do lam = 1, n_lambda
211 3 : accum = 0.0_dp
212 :
213 6 : do oz = 0, iz_max
214 9 : do oy = 0, iy_max
215 12 : do ox = 0, ix_max
216 6 : weight = w_x(ox)*w_y(oy)*w_z(oz)
217 9 : accum = accum + weight*f_values_4d(i_x + ox, i_y + oy, i_z + oz, lam)
218 : end do
219 : end do
220 : end do
221 :
222 3 : if (is_nan(accum) .or. is_inf(accum) .or. accum <= 0.0_dp) then
223 : call find_nearest_point(x_val, y_val, z_val, x_grid, y_grid, z_grid, &
224 0 : nearest_x, nearest_y, nearest_z)
225 0 : accum = f_values_4d(nearest_x, nearest_y, nearest_z, lam)
226 : end if
227 :
228 3 : if (is_nan(accum) .or. is_inf(accum) .or. accum < 0.0_dp) then
229 0 : write (*, *) 'colors: linear and nearest-grid interpolation returned bad flux'
230 0 : call mesa_error(__FILE__, __LINE__)
231 : end if
232 4 : result_flux(lam) = max(tiny_value, accum)
233 : end do
234 :
235 1 : end subroutine trilinear_interp_vector
236 :
237 : ! scalar trilinear interpolation (external callers / single-wavelength use)
238 : ! retained for backward compatibility
239 0 : function trilinear_interp(x_val, y_val, z_val, x_grid, y_grid, z_grid, f_values) result(f_interp)
240 : real(dp), intent(in) :: x_val, y_val, z_val
241 : real(dp), intent(in) :: x_grid(:), y_grid(:), z_grid(:)
242 : real(dp), intent(in) :: f_values(:, :, :)
243 : real(dp) :: f_interp
244 :
245 : integer :: i_x, i_y, i_z
246 : integer :: nx, ny, nz
247 : integer :: ix_max, iy_max, iz_max
248 : integer :: ox, oy, oz
249 : real(dp) :: t_x, t_y, t_z
250 : real(dp) :: w_x(0:1), w_y(0:1), w_z(0:1)
251 : real(dp) :: accum, weight
252 : real(dp), parameter :: tiny_value = 1.0d-300
253 :
254 0 : nx = size(x_grid)
255 0 : ny = size(y_grid)
256 0 : nz = size(z_grid)
257 :
258 : if (nx < 1 .or. ny < 1 .or. nz < 1 .or. &
259 0 : size(f_values, 1) /= nx .or. size(f_values, 2) /= ny .or. &
260 : size(f_values, 3) /= nz) then
261 0 : write (*, *) 'colors: invalid array sizes in trilinear_interp'
262 0 : call mesa_error(__FILE__, __LINE__)
263 : end if
264 :
265 : call find_containing_cell(x_val, y_val, z_val, x_grid, y_grid, z_grid, &
266 0 : i_x, i_y, i_z, t_x, t_y, t_z)
267 :
268 0 : call setup_axis(nx, x_val, x_grid, i_x, t_x, ix_max, w_x)
269 0 : call setup_axis(ny, y_val, y_grid, i_y, t_y, iy_max, w_y)
270 0 : call setup_axis(nz, z_val, z_grid, i_z, t_z, iz_max, w_z)
271 :
272 0 : accum = 0.0_dp
273 :
274 0 : do oz = 0, iz_max
275 0 : do oy = 0, iy_max
276 0 : do ox = 0, ix_max
277 0 : weight = w_x(ox)*w_y(oy)*w_z(oz)
278 0 : accum = accum + weight*f_values(i_x + ox, i_y + oy, i_z + oz)
279 : end do
280 : end do
281 : end do
282 :
283 0 : if (is_nan(accum) .or. is_inf(accum) .or. accum <= 0.0_dp) then
284 : call find_nearest_point(x_val, y_val, z_val, x_grid, y_grid, z_grid, &
285 0 : i_x, i_y, i_z)
286 0 : accum = f_values(i_x, i_y, i_z)
287 : end if
288 :
289 0 : if (is_nan(accum) .or. is_inf(accum) .or. accum < 0.0_dp) then
290 0 : write (*, *) 'colors: linear and nearest-grid interpolation returned bad flux'
291 0 : call mesa_error(__FILE__, __LINE__)
292 : end if
293 0 : f_interp = max(tiny_value, accum)
294 :
295 0 : end function trilinear_interp
296 :
297 3 : subroutine setup_axis(n, x_val, grid, i, t, i_max, w)
298 : integer, intent(in) :: n
299 : real(dp), intent(in) :: x_val
300 : real(dp), intent(in) :: grid(:)
301 : integer, intent(inout) :: i
302 : real(dp), intent(inout) :: t
303 : integer, intent(out) :: i_max
304 : real(dp), intent(out) :: w(0:1)
305 :
306 3 : if (n <= 1) then
307 2 : i = 1
308 2 : t = 0.0_dp
309 2 : i_max = 0
310 2 : w(0) = 1.0_dp
311 2 : w(1) = 0.0_dp
312 2 : return
313 : end if
314 :
315 1 : if (x_val <= grid(1)) then
316 0 : i = 1
317 0 : t = 0.0_dp
318 1 : else if (x_val >= grid(n)) then
319 0 : i = n - 1
320 0 : t = 1.0_dp
321 : else
322 1 : if (i < 1) i = 1
323 1 : if (i >= n) i = n - 1
324 1 : t = max(0.0_dp, min(1.0_dp, t))
325 : end if
326 :
327 1 : i_max = 1
328 1 : w(0) = 1.0_dp - t
329 1 : w(1) = t
330 : end subroutine setup_axis
331 :
332 : end module linear_interp
|
