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