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