Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2025 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 : ! ***********************************************************************
21 : ! Linear interpolation module for spectral energy distributions (SEDs)
22 : ! ***********************************************************************
23 :
24 : module linear_interp
25 : use const_def, only: dp
26 : use colors_utils, only: dilute_flux
27 : use utils_lib, only: mesa_error
28 : implicit none
29 :
30 : private
31 : public :: construct_sed_linear, trilinear_interp
32 :
33 : contains
34 :
35 : !---------------------------------------------------------------------------
36 : ! Main entry point: Construct a SED using linear interpolation
37 : !---------------------------------------------------------------------------
38 :
39 0 : subroutine construct_sed_linear(teff, log_g, metallicity, R, d, file_names, &
40 : lu_teff, lu_logg, lu_meta, stellar_model_dir, &
41 : wavelengths, fluxes)
42 :
43 : real(dp), intent(in) :: teff, log_g, metallicity, R, d
44 : real(dp), intent(in) :: lu_teff(:), lu_logg(:), lu_meta(:)
45 : character(len=*), intent(in) :: stellar_model_dir
46 : character(len=100), intent(in) :: file_names(:)
47 : real(dp), dimension(:), allocatable, intent(out) :: wavelengths, fluxes
48 :
49 : integer :: i, n_lambda, status, n_teff, n_logg, n_meta
50 0 : real(dp), dimension(:), allocatable :: interp_flux, diluted_flux
51 0 : real(dp), dimension(:, :, :, :), allocatable :: precomputed_flux_cube
52 0 : real(dp), dimension(:, :, :), allocatable :: flux_cube_lambda
53 : real(dp) :: min_flux, max_flux, mean_flux, progress_pct
54 :
55 : ! Parameter grids
56 0 : real(dp), allocatable :: teff_grid(:), logg_grid(:), meta_grid(:)
57 : character(len=256) :: bin_filename, clean_path
58 : logical :: file_exists
59 :
60 : ! Clean up any double slashes in the path
61 0 : clean_path = trim(stellar_model_dir)
62 0 : if (clean_path(len_trim(clean_path):len_trim(clean_path)) == '/') then
63 0 : bin_filename = trim(clean_path)//'flux_cube.bin'
64 : else
65 0 : bin_filename = trim(clean_path)//'/flux_cube.bin'
66 : end if
67 :
68 : ! Check if file exists first
69 0 : INQUIRE (file=bin_filename, EXIST=file_exists)
70 :
71 0 : if (.not. file_exists) then
72 0 : print *, 'Missing required binary file for interpolation'
73 0 : call mesa_error(__FILE__, __LINE__)
74 : end if
75 :
76 : ! Load the data from binary file
77 : call load_binary_data(bin_filename, teff_grid, logg_grid, meta_grid, &
78 0 : wavelengths, precomputed_flux_cube, status)
79 :
80 0 : if (status /= 0) then
81 0 : print *, 'Binary data loading error'
82 0 : call mesa_error(__FILE__, __LINE__)
83 : end if
84 :
85 0 : n_teff = size(teff_grid)
86 0 : n_logg = size(logg_grid)
87 0 : n_meta = size(meta_grid)
88 0 : n_lambda = size(wavelengths)
89 :
90 : ! Allocate space for interpolated flux
91 0 : allocate (interp_flux(n_lambda))
92 0 : allocate (flux_cube_lambda(n_teff, n_logg, n_meta))
93 :
94 : ! Perform trilinear interpolation for each wavelength
95 0 : do i = 1, n_lambda
96 :
97 : ! Extract the 3D grid for this wavelength
98 0 : flux_cube_lambda = precomputed_flux_cube(:, :, :, i)
99 :
100 : ! Simple trilinear interpolation at the target parameters
101 : interp_flux(i) = trilinear_interp(teff, log_g, metallicity, &
102 0 : teff_grid, logg_grid, meta_grid, flux_cube_lambda)
103 : end do
104 :
105 :
106 0 : deallocate(flux_cube_lambda)
107 :
108 : ! Calculate statistics for validation
109 : min_flux = minval(interp_flux)
110 : max_flux = maxval(interp_flux)
111 0 : mean_flux = sum(interp_flux)/n_lambda
112 :
113 : ! Apply distance dilution to get observed flux
114 0 : allocate (diluted_flux(n_lambda))
115 0 : call dilute_flux(interp_flux, R, d, diluted_flux)
116 0 : fluxes = diluted_flux
117 :
118 : ! Calculate statistics after dilution
119 : min_flux = minval(diluted_flux)
120 : max_flux = maxval(diluted_flux)
121 : mean_flux = sum(diluted_flux)/n_lambda
122 :
123 0 : end subroutine construct_sed_linear
124 :
125 : !---------------------------------------------------------------------------
126 : ! Load data from binary file
127 : !---------------------------------------------------------------------------
128 0 : subroutine load_binary_data(filename, teff_grid, logg_grid, meta_grid, &
129 : wavelengths, flux_cube, status)
130 : character(len=*), intent(in) :: filename
131 : real(dp), allocatable, intent(out) :: teff_grid(:), logg_grid(:), meta_grid(:)
132 : real(dp), allocatable, intent(out) :: wavelengths(:)
133 : real(dp), allocatable, intent(out) :: flux_cube(:, :, :, :)
134 : integer, intent(out) :: status
135 :
136 : integer :: unit, n_teff, n_logg, n_meta, n_lambda
137 :
138 0 : unit = 99
139 0 : status = 0
140 :
141 : ! Open the binary file
142 0 : open (unit=unit, file=filename, status='OLD', ACCESS='STREAM', FORM='UNFORMATTED', iostat=status)
143 0 : if (status /= 0) then
144 : !print *, 'Error opening binary file:', trim(filename)
145 : return
146 : end if
147 :
148 : ! Read dimensions
149 0 : read (unit, iostat=status) n_teff, n_logg, n_meta, n_lambda
150 0 : if (status /= 0) then
151 : !print *, 'Error reading dimensions from binary file'
152 0 : close (unit)
153 0 : return
154 : end if
155 :
156 : ! Allocate arrays based on dimensions
157 0 : allocate (teff_grid(n_teff), STAT=status)
158 0 : if (status /= 0) then
159 : !print *, 'Error allocating teff_grid array'
160 0 : close (unit)
161 0 : return
162 : end if
163 :
164 0 : allocate (logg_grid(n_logg), STAT=status)
165 0 : if (status /= 0) then
166 : !print *, 'Error allocating logg_grid array'
167 0 : close (unit)
168 0 : return
169 : end if
170 :
171 0 : allocate (meta_grid(n_meta), STAT=status)
172 0 : if (status /= 0) then
173 : !print *, 'Error allocating meta_grid array'
174 0 : close (unit)
175 0 : return
176 : end if
177 :
178 0 : allocate (wavelengths(n_lambda), STAT=status)
179 0 : if (status /= 0) then
180 : !print *, 'Error allocating wavelengths array'
181 0 : close (unit)
182 0 : return
183 : end if
184 :
185 0 : allocate (flux_cube(n_teff, n_logg, n_meta, n_lambda), STAT=status)
186 0 : if (status /= 0) then
187 : !print *, 'Error allocating flux_cube array'
188 0 : close (unit)
189 0 : return
190 : end if
191 :
192 : ! Read grid arrays
193 0 : read (unit, iostat=status) teff_grid
194 0 : if (status /= 0) then
195 : !print *, 'Error reading teff_grid'
196 : GOTO 999 ! Cleanup and return
197 : end if
198 :
199 0 : read (unit, iostat=status) logg_grid
200 0 : if (status /= 0) then
201 : !print *, 'Error reading logg_grid'
202 : GOTO 999 ! Cleanup and return
203 : end if
204 :
205 0 : read (unit, iostat=status) meta_grid
206 0 : if (status /= 0) then
207 : !print *, 'Error reading meta_grid'
208 : GOTO 999 ! Cleanup and return
209 : end if
210 :
211 0 : read (unit, iostat=status) wavelengths
212 0 : if (status /= 0) then
213 : !print *, 'Error reading wavelengths'
214 : GOTO 999 ! Cleanup and return
215 : end if
216 :
217 : ! Read flux cube
218 0 : read (unit, iostat=status) flux_cube
219 0 : if (status /= 0) then
220 : !print *, 'Error reading flux_cube'
221 : GOTO 999 ! Cleanup and return
222 : end if
223 :
224 : ! Close file and return success
225 0 : close (unit)
226 0 : return
227 :
228 : 999 CONTINUE
229 : ! Cleanup on error
230 0 : close (unit)
231 0 : return
232 :
233 : ! After reading the grid arrays
234 : !print *, 'Teff grid min/max:', minval(teff_grid), maxval(teff_grid)
235 : !print *, 'logg grid min/max:', minval(logg_grid), maxval(logg_grid)
236 : !print *, 'meta grid min/max:', minval(meta_grid), maxval(meta_grid)
237 :
238 : end subroutine load_binary_data
239 :
240 : !---------------------------------------------------------------------------
241 : ! Simple trilinear interpolation function
242 : !---------------------------------------------------------------------------
243 : !---------------------------------------------------------------------------
244 : ! Log-space trilinear interpolation function with normalization
245 : !---------------------------------------------------------------------------
246 0 : function trilinear_interp(x_val, y_val, z_val, x_grid, y_grid, z_grid, f_values) result(f_interp)
247 : real(dp), intent(in) :: x_val, y_val, z_val
248 : real(dp), intent(in) :: x_grid(:), y_grid(:), z_grid(:)
249 : real(dp), intent(in) :: f_values(:, :, :)
250 : real(dp) :: f_interp
251 : ! Compute log-space result
252 : real(dp) :: log_result
253 : integer :: i_x, i_y, i_z
254 : real(dp) :: t_x, t_y, t_z
255 : real(dp) :: c000, c001, c010, c011, c100, c101, c110, c111
256 : real(dp) :: c00, c01, c10, c11, c0, c1
257 : real(dp), parameter :: tiny_value = 1.0e-10_dp
258 :
259 : ! Find containing cell and parameter values using binary search
260 : call find_containing_cell(x_val, y_val, z_val, x_grid, y_grid, z_grid, &
261 0 : i_x, i_y, i_z, t_x, t_y, t_z)
262 :
263 : ! Boundary safety check
264 0 : if (i_x < lbound(x_grid,1)) i_x = lbound(x_grid,1)
265 0 : if (i_y < lbound(y_grid,1)) i_y = lbound(y_grid,1)
266 0 : if (i_z < lbound(z_grid,1)) i_z = lbound(z_grid,1)
267 0 : if (i_x >= ubound(x_grid,1)) i_x = ubound(x_grid,1) - 1
268 0 : if (i_y >= ubound(y_grid,1)) i_y = ubound(y_grid,1) - 1
269 0 : if (i_z >= ubound(z_grid,1)) i_z = ubound(z_grid,1) - 1
270 :
271 : ! Force interpolation parameters to be in [0,1]
272 0 : t_x = max(0.0_dp, MIN(1.0_dp, t_x))
273 0 : t_y = max(0.0_dp, MIN(1.0_dp, t_y))
274 0 : t_z = max(0.0_dp, MIN(1.0_dp, t_z))
275 :
276 : ! Get the corners of the cube with safety checks
277 0 : c000 = max(tiny_value, f_values(i_x, i_y, i_z))
278 0 : c001 = max(tiny_value, f_values(i_x, i_y, i_z + 1))
279 0 : c010 = max(tiny_value, f_values(i_x, i_y + 1, i_z))
280 0 : c011 = max(tiny_value, f_values(i_x, i_y + 1, i_z + 1))
281 0 : c100 = max(tiny_value, f_values(i_x + 1, i_y, i_z))
282 0 : c101 = max(tiny_value, f_values(i_x + 1, i_y, i_z + 1))
283 0 : c110 = max(tiny_value, f_values(i_x + 1, i_y + 1, i_z))
284 0 : c111 = max(tiny_value, f_values(i_x + 1, i_y + 1, i_z + 1))
285 :
286 : ! Try standard linear interpolation first (safer)
287 0 : c00 = c000*(1.0_dp - t_x) + c100*t_x
288 0 : c01 = c001*(1.0_dp - t_x) + c101*t_x
289 0 : c10 = c010*(1.0_dp - t_x) + c110*t_x
290 0 : c11 = c011*(1.0_dp - t_x) + c111*t_x
291 :
292 0 : c0 = c00*(1.0_dp - t_y) + c10*t_y
293 0 : c1 = c01*(1.0_dp - t_y) + c11*t_y
294 :
295 0 : f_interp = c0*(1.0_dp - t_z) + c1*t_z
296 :
297 : ! If the linear result is valid and non-zero, try log space
298 0 : if (f_interp > tiny_value) then
299 : ! Perform log-space interpolation
300 0 : c00 = log(c000)*(1.0_dp - t_x) + log(c100)*t_x
301 0 : c01 = log(c001)*(1.0_dp - t_x) + log(c101)*t_x
302 0 : c10 = log(c010)*(1.0_dp - t_x) + log(c110)*t_x
303 0 : c11 = log(c011)*(1.0_dp - t_x) + log(c111)*t_x
304 :
305 0 : c0 = c00*(1.0_dp - t_y) + c10*t_y
306 0 : c1 = c01*(1.0_dp - t_y) + c11*t_y
307 :
308 0 : log_result = c0*(1.0_dp - t_z) + c1*t_z
309 :
310 : ! Only use the log-space result if it's valid
311 0 : if (log_result == log_result) then ! NaN check
312 0 : f_interp = EXP(log_result)
313 : end if
314 : end if
315 :
316 : ! Final sanity check
317 0 : if (f_interp /= f_interp .or. f_interp <= 0.0_dp) then
318 : ! If we somehow still got an invalid result, use nearest neighbor
319 0 : call find_nearest_point(x_val, y_val, z_val, x_grid, y_grid, z_grid, i_x, i_y, i_z)
320 0 : f_interp = max(tiny_value, f_values(i_x, i_y, i_z))
321 : end if
322 0 : end function trilinear_interp
323 :
324 : !---------------------------------------------------------------------------
325 : ! Find the cell containing the interpolation point
326 : !---------------------------------------------------------------------------
327 0 : subroutine find_containing_cell(x_val, y_val, z_val, x_grid, y_grid, z_grid, &
328 : i_x, i_y, i_z, t_x, t_y, t_z)
329 : real(dp), intent(in) :: x_val, y_val, z_val
330 : real(dp), intent(in) :: x_grid(:), y_grid(:), z_grid(:)
331 : integer, intent(out) :: i_x, i_y, i_z
332 : real(dp), intent(out) :: t_x, t_y, t_z
333 :
334 : ! Find x interval
335 0 : call find_interval(x_grid, x_val, i_x, t_x)
336 :
337 : ! Find y interval
338 0 : call find_interval(y_grid, y_val, i_y, t_y)
339 :
340 : ! Find z interval
341 0 : call find_interval(z_grid, z_val, i_z, t_z)
342 0 : end subroutine find_containing_cell
343 :
344 : !---------------------------------------------------------------------------
345 : ! Find the interval in a sorted array containing a value
346 : !---------------------------------------------------------------------------
347 0 : subroutine find_interval(x, val, i, t)
348 : real(dp), intent(in) :: x(:), val
349 : integer, intent(out) :: i
350 : real(dp), intent(out) :: t
351 :
352 : integer :: n, lo, hi, mid
353 : logical :: dummy_axis
354 :
355 0 : n = size(x)
356 :
357 : ! Detect dummy axis
358 0 : dummy_axis = all(x == 0.0_dp) .or. all(x == 999.0_dp) .or. all(x == -999.0_dp)
359 :
360 : if (dummy_axis) then
361 : ! Collapse: use the first element of the axis, no interpolation
362 0 : i = 1
363 0 : t = 0.0_dp
364 0 : return
365 : end if
366 :
367 : ! --- ORIGINAL CODE BELOW ---
368 0 : if (val <= x(1)) then
369 0 : i = 1
370 0 : t = 0.0_dp
371 0 : return
372 0 : else if (val >= x(n)) then
373 0 : i = n - 1
374 0 : t = 1.0_dp
375 0 : return
376 : end if
377 :
378 : lo = 1
379 : hi = n
380 0 : do while (hi - lo > 1)
381 0 : mid = (lo + hi)/2
382 0 : if (val >= x(mid)) then
383 : lo = mid
384 : else
385 0 : hi = mid
386 : end if
387 : end do
388 :
389 0 : i = lo
390 0 : t = (val - x(i))/(x(i + 1) - x(i))
391 : end subroutine find_interval
392 :
393 : !---------------------------------------------------------------------------
394 : ! Find the nearest grid point
395 : !---------------------------------------------------------------------------
396 0 : subroutine find_nearest_point(x_val, y_val, z_val, x_grid, y_grid, z_grid, &
397 : i_x, i_y, i_z)
398 : real(dp), intent(in) :: x_val, y_val, z_val
399 : real(dp), intent(in) :: x_grid(:), y_grid(:), z_grid(:)
400 : integer, intent(out) :: i_x, i_y, i_z
401 :
402 : integer :: i
403 : real(dp) :: min_dist, dist
404 :
405 : ! Find nearest x grid point
406 0 : min_dist = abs(x_val - x_grid(1))
407 0 : i_x = 1
408 0 : do i = 2, size(x_grid)
409 0 : dist = abs(x_val - x_grid(i))
410 0 : if (dist < min_dist) then
411 0 : min_dist = dist
412 0 : i_x = i
413 : end if
414 : end do
415 :
416 : ! Find nearest y grid point
417 0 : min_dist = abs(y_val - y_grid(1))
418 0 : i_y = 1
419 0 : do i = 2, size(y_grid)
420 0 : dist = abs(y_val - y_grid(i))
421 0 : if (dist < min_dist) then
422 0 : min_dist = dist
423 0 : i_y = i
424 : end if
425 : end do
426 :
427 : ! Find nearest z grid point
428 0 : min_dist = abs(z_val - z_grid(1))
429 0 : i_z = 1
430 0 : do i = 2, size(z_grid)
431 0 : dist = abs(z_val - z_grid(i))
432 0 : if (dist < min_dist) then
433 0 : min_dist = dist
434 0 : i_z = i
435 : end if
436 : end do
437 0 : end subroutine find_nearest_point
438 :
439 : end module linear_interp
|
