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 : module colors_utils
21 : use const_def, only: dp, strlen, mesa_dir
22 : use colors_def, only: Colors_General_Info, sed_mem_cache_cap
23 : use utils_lib, only: mesa_error
24 :
25 : implicit none
26 :
27 : public :: dilute_flux, trapezoidal_integration, &
28 : simpson_integration, load_sed, load_filter, load_vega_sed, &
29 : load_lookup_table, remove_dat, load_flux_cube, build_unique_grids, &
30 : build_grid_to_lu_map, &
31 : find_containing_cell, find_interval, find_nearest_point, &
32 : find_bracket_index, load_sed_cached, load_stencil
33 : contains
34 :
35 : ! apply dilution factor (R/d)^2 to convert surface flux to observed flux
36 20 : subroutine dilute_flux(surface_flux, R, d, calibrated_flux)
37 : real(dp), intent(in) :: surface_flux(:)
38 : real(dp), intent(in) :: R, d ! R = stellar radius, d = distance (both in the same units, e.g., cm)
39 : real(dp), intent(out) :: calibrated_flux(:)
40 :
41 20 : if (size(calibrated_flux) /= size(surface_flux)) then
42 0 : print *, "Error in dilute_flux: Output array must have the same size as input array."
43 0 : call mesa_error(__FILE__, __LINE__)
44 : end if
45 :
46 24000 : calibrated_flux = surface_flux*((R/d)**2)
47 20 : end subroutine dilute_flux
48 :
49 0 : subroutine trapezoidal_integration(x, y, result)
50 : real(dp), dimension(:), intent(in) :: x, y
51 : real(dp), intent(out) :: result
52 :
53 : integer :: i, n
54 : real(dp) :: sum
55 :
56 0 : n = size(x)
57 0 : sum = 0.0_dp
58 :
59 0 : if (size(x) /= size(y)) then
60 0 : print *, "Error: x and y arrays must have the same size."
61 0 : call mesa_error(__FILE__, __LINE__)
62 : end if
63 :
64 0 : if (size(x) < 2) then
65 0 : print *, "Error: x and y arrays must have at least 2 elements."
66 0 : call mesa_error(__FILE__, __LINE__)
67 : end if
68 :
69 0 : do i = 1, n - 1
70 0 : sum = sum + 0.5_dp*(x(i + 1) - x(i))*(y(i + 1) + y(i))
71 : end do
72 :
73 0 : result = sum
74 0 : end subroutine trapezoidal_integration
75 :
76 272 : subroutine simpson_integration(x, y, result)
77 : real(dp), dimension(:), intent(in) :: x, y
78 : real(dp), intent(out) :: result
79 :
80 : integer :: i, n
81 : real(dp) :: sum, h1, h2, f1, f2, f0
82 :
83 272 : n = size(x)
84 272 : sum = 0.0_dp
85 :
86 272 : if (size(x) /= size(y)) then
87 0 : print *, "Error: x and y arrays must have the same size."
88 0 : call mesa_error(__FILE__, __LINE__)
89 : end if
90 :
91 272 : if (size(x) < 2) then
92 0 : print *, "Error: x and y arrays must have at least 2 elements."
93 0 : call mesa_error(__FILE__, __LINE__)
94 : end if
95 :
96 : ! adaptive Simpson's rule
97 272 : do i = 1, n - 2, 2
98 194682 : h1 = x(i + 1) - x(i)
99 194682 : h2 = x(i + 2) - x(i + 1)
100 :
101 194682 : f0 = y(i)
102 194682 : f1 = y(i + 1)
103 194682 : f2 = y(i + 2)
104 :
105 : ! Simpson's rule: (h/3) * (f0 + 4f1 + f2)
106 : sum = sum + (h1 + h2)/6.0_dp * ( &
107 : f0*(2.0_dp*h1 - h2)/h1 + &
108 : f1*(h1 + h2)**2/(h1*h2) + &
109 194682 : f2*(2.0_dp*h2 - h1)/h2)
110 :
111 : end do
112 :
113 : ! handle the last interval if n is even (odd number of points)
114 272 : if (MOD(n, 2) == 0) then
115 18 : sum = sum + 0.5_dp*(x(n) - x(n - 1))*(y(n) + y(n - 1))
116 : end if
117 :
118 272 : result = sum
119 272 : end subroutine simpson_integration
120 :
121 1 : subroutine load_vega_sed(filepath, wavelengths, flux)
122 : character(len=*), intent(in) :: filepath
123 : real(dp), dimension(:), allocatable, intent(out) :: wavelengths, flux
124 : character(len=512) :: line
125 : integer :: unit, n_rows, status, i
126 : real(dp) :: temp_wave, temp_flux
127 :
128 1 : open (newunit=unit, file=trim(filepath), status='OLD', action='READ', iostat=status)
129 1 : if (status /= 0) then
130 0 : print *, "Error: Could not open Vega SED file ", trim(filepath)
131 0 : call mesa_error(__FILE__, __LINE__)
132 : end if
133 :
134 1 : read (unit, '(A)', iostat=status) line
135 1 : if (status /= 0) then
136 0 : print *, "Error: Could not read header from Vega SED file ", trim(filepath)
137 0 : call mesa_error(__FILE__, __LINE__)
138 : end if
139 :
140 : n_rows = 0
141 8094 : do
142 8095 : read (unit, '(A)', iostat=status) line
143 8095 : if (status /= 0) exit
144 8094 : n_rows = n_rows + 1
145 : end do
146 :
147 1 : rewind (unit)
148 1 : read (unit, '(A)', iostat=status) line ! skip header again
149 :
150 3 : allocate (wavelengths(n_rows))
151 2 : allocate (flux(n_rows))
152 :
153 1 : i = 0
154 8094 : do
155 8095 : read (unit, *, iostat=status) temp_wave, temp_flux ! ignore any extra columns
156 8095 : if (status /= 0) exit
157 8094 : i = i + 1
158 8094 : wavelengths(i) = temp_wave
159 8094 : flux(i) = temp_flux
160 : end do
161 :
162 1 : close (unit)
163 1 : end subroutine load_vega_sed
164 :
165 7 : subroutine load_filter(directory, filter_wavelengths, filter_trans)
166 : character(len=*), intent(in) :: directory
167 : real(dp), dimension(:), allocatable, intent(out) :: filter_wavelengths, filter_trans
168 :
169 : character(len=512) :: line
170 : integer :: unit, n_rows, status, i
171 : real(dp) :: temp_wavelength, temp_trans
172 :
173 7 : open (newunit=unit, file=trim(directory), status='OLD', action='READ', iostat=status)
174 7 : if (status /= 0) then
175 0 : print *, "Error: Could not open file ", trim(directory)
176 0 : call mesa_error(__FILE__, __LINE__)
177 : end if
178 :
179 7 : read (unit, '(A)', iostat=status) line
180 7 : if (status /= 0) then
181 0 : print *, "Error: Could not read the file", trim(directory)
182 0 : call mesa_error(__FILE__, __LINE__)
183 : end if
184 :
185 : n_rows = 0
186 142 : do
187 149 : read (unit, '(A)', iostat=status) line
188 149 : if (status /= 0) exit
189 142 : n_rows = n_rows + 1
190 : end do
191 :
192 21 : allocate (filter_wavelengths(n_rows))
193 14 : allocate (filter_trans(n_rows))
194 :
195 : ! rewind to the first non-comment line
196 7 : rewind (unit)
197 : do
198 7 : read (unit, '(A)', iostat=status) line
199 7 : if (status /= 0) then
200 0 : print *, "Error: Could not rewind file", trim(directory)
201 0 : call mesa_error(__FILE__, __LINE__)
202 : end if
203 7 : if (line(1:1) /= "#") exit
204 : end do
205 :
206 : i = 0
207 142 : do
208 149 : read (unit, *, iostat=status) temp_wavelength, temp_trans
209 149 : if (status /= 0) exit
210 142 : i = i + 1
211 :
212 142 : filter_wavelengths(i) = temp_wavelength
213 142 : filter_trans(i) = temp_trans
214 : end do
215 :
216 7 : close (unit)
217 7 : end subroutine load_filter
218 :
219 : ! parses a csv lookup table mapping atmosphere grid parameters to SED filenames
220 1 : subroutine load_lookup_table(lookup_file, lookup_table, out_file_names, &
221 : out_logg, out_meta, out_teff)
222 :
223 : character(len=*), intent(in) :: lookup_file
224 : REAL, dimension(:, :), allocatable, intent(out) :: lookup_table
225 : character(len=100), allocatable, intent(inout) :: out_file_names(:)
226 : real(dp), allocatable, intent(inout) :: out_logg(:), out_meta(:), out_teff(:)
227 :
228 : integer :: i, n_rows, status, unit, ios
229 : character(len=512) :: line
230 : character(len=*), parameter :: delimiter = ","
231 1 : character(len=100), allocatable :: columns(:), headers(:)
232 : character(len=256) :: token
233 : integer :: logg_col, meta_col, teff_col
234 :
235 1 : open (newunit=unit, file=lookup_file, status='old', action='read', iostat=status)
236 1 : if (status /= 0) then
237 0 : print *, "Error: Could not open file", lookup_file
238 0 : call mesa_error(__FILE__, __LINE__)
239 : end if
240 :
241 1 : read (unit, '(A)', iostat=status) line
242 1 : if (status /= 0) then
243 0 : print *, "Error: Could not read header line"
244 0 : call mesa_error(__FILE__, __LINE__)
245 : end if
246 :
247 1 : call split_line(line, delimiter, headers)
248 :
249 : ! determine column indices -- try all plausible header name variants
250 1 : logg_col = get_column_index(headers, "logg")
251 1 : if (logg_col < 0) logg_col = get_column_index(headers, "log_g")
252 0 : if (logg_col < 0) logg_col = get_column_index(headers, "log(g)")
253 1 : if (logg_col < 0) logg_col = get_column_index(headers, "log10g")
254 0 : if (logg_col < 0) logg_col = get_column_index(headers, "log10_g")
255 :
256 1 : teff_col = get_column_index(headers, "teff")
257 1 : if (teff_col < 0) teff_col = get_column_index(headers, "t_eff")
258 0 : if (teff_col < 0) teff_col = get_column_index(headers, "t(eff)")
259 0 : if (teff_col < 0) teff_col = get_column_index(headers, "temperature")
260 0 : if (teff_col < 0) teff_col = get_column_index(headers, "temp")
261 :
262 1 : meta_col = get_column_index(headers, "meta")
263 1 : if (meta_col < 0) meta_col = get_column_index(headers, "feh")
264 1 : if (meta_col < 0) meta_col = get_column_index(headers, "fe_h")
265 1 : if (meta_col < 0) meta_col = get_column_index(headers, "[fe/h]")
266 1 : if (meta_col < 0) meta_col = get_column_index(headers, "mh")
267 1 : if (meta_col < 0) meta_col = get_column_index(headers, "[m/h]")
268 1 : if (meta_col < 0) meta_col = get_column_index(headers, "m_h")
269 1 : if (meta_col < 0) meta_col = get_column_index(headers, "z")
270 1 : if (meta_col < 0) meta_col = get_column_index(headers, "logz")
271 1 : if (meta_col < 0) meta_col = get_column_index(headers, "metallicity")
272 1 : if (meta_col < 0) meta_col = get_column_index(headers, "metals")
273 :
274 1 : n_rows = 0
275 3808 : do
276 3809 : read (unit, '(A)', iostat=status) line
277 3809 : if (status /= 0) exit
278 3808 : n_rows = n_rows + 1
279 : end do
280 1 : rewind (unit)
281 :
282 1 : read (unit, '(A)', iostat=status) line ! skip header
283 :
284 3 : allocate (out_file_names(n_rows))
285 5 : allocate (out_logg(n_rows), out_meta(n_rows), out_teff(n_rows))
286 :
287 1 : i = 0
288 : do
289 3809 : read (unit, '(A)', iostat=status) line
290 3809 : if (status /= 0) exit
291 3808 : i = i + 1
292 :
293 3808 : call split_line(line, delimiter, columns)
294 :
295 3808 : out_file_names(i) = columns(1)
296 :
297 : ! robust numeric parsing: never crash on bad/missing values
298 3808 : if (logg_col > 0) then
299 3808 : token = trim(adjustl(columns(logg_col)))
300 3808 : if (len_trim(token) == 0 .or. token == '""') then
301 0 : out_logg(i) = -999.0_dp
302 : else
303 3808 : read(token, *, iostat=ios) out_logg(i)
304 3808 : if (ios /= 0) out_logg(i) = -999.0_dp
305 : end if
306 : else
307 0 : out_logg(i) = -999.0
308 : end if
309 :
310 3808 : if (meta_col > 0) then
311 3808 : token = trim(adjustl(columns(meta_col)))
312 3808 : if (len_trim(token) == 0 .or. token == '""') then
313 0 : out_meta(i) = 0.0_dp
314 : else
315 3808 : read(token, *, iostat=ios) out_meta(i)
316 3808 : if (ios /= 0) out_meta(i) = 0.0_dp
317 : end if
318 : else
319 0 : out_meta(i) = 0.0
320 : end if
321 :
322 3809 : if (teff_col > 0) then
323 3808 : token = trim(adjustl(columns(teff_col)))
324 3808 : if (len_trim(token) == 0 .or. token == '""') then
325 0 : out_teff(i) = 0.0_dp
326 : else
327 3808 : read(token, *, iostat=ios) out_teff(i)
328 3808 : if (ios /= 0) out_teff(i) = 0.0_dp
329 : end if
330 : else
331 0 : out_teff(i) = 0.0
332 : end if
333 : end do
334 :
335 1 : close (unit)
336 :
337 : contains
338 :
339 24 : function get_column_index(headers, target) result(index)
340 : character(len=100), intent(in) :: headers(:)
341 : character(len=*), intent(in) :: target
342 : integer :: index, i
343 : character(len=100) :: clean_header, clean_target
344 :
345 12 : index = -1
346 12 : clean_target = trim(adjustl(target))
347 :
348 159 : do i = 1, size(headers)
349 150 : clean_header = trim(adjustl(headers(i)))
350 159 : if (clean_header == clean_target) then
351 : index = i
352 : exit
353 : end if
354 : end do
355 12 : end function get_column_index
356 :
357 3809 : subroutine split_line(line, delimiter, tokens)
358 : character(len=*), intent(in) :: line, delimiter
359 : character(len=100), allocatable, intent(out) :: tokens(:)
360 : integer :: num_tokens, pos, start, len_delim
361 :
362 3809 : len_delim = len_trim(delimiter)
363 3809 : start = 1
364 3809 : num_tokens = 0
365 3809 : if (allocated(tokens)) deallocate (tokens)
366 :
367 49517 : do
368 53326 : pos = index(line(start:), delimiter)
369 :
370 53326 : if (pos == 0) exit
371 49517 : num_tokens = num_tokens + 1
372 49517 : call append_token(tokens, line(start:start + pos - 2))
373 49517 : start = start + pos + len_delim - 1
374 : end do
375 :
376 3809 : num_tokens = num_tokens + 1
377 3809 : call append_token(tokens, line(start:))
378 3809 : end subroutine split_line
379 :
380 53326 : subroutine append_token(tokens, token)
381 : character(len=*), intent(in) :: token
382 : character(len=100), allocatable, intent(inout) :: tokens(:)
383 53326 : character(len=100), allocatable :: temp(:)
384 : integer :: n
385 :
386 53326 : if (.not. allocated(tokens)) then
387 3809 : allocate (tokens(1))
388 3809 : tokens(1) = token
389 : else
390 49517 : n = size(tokens)
391 148551 : allocate (temp(n))
392 445653 : temp = tokens
393 49517 : deallocate (tokens)
394 148551 : allocate (tokens(n + 1))
395 396136 : tokens(1:n) = temp
396 49517 : tokens(n + 1) = token
397 49517 : deallocate (temp) ! unsure if this is till needed.
398 : end if
399 53326 : end subroutine append_token
400 :
401 : end subroutine load_lookup_table
402 :
403 :
404 :
405 :
406 0 : subroutine load_sed(directory, index, wavelengths, flux)
407 : character(len=*), intent(in) :: directory
408 : integer, intent(in) :: index
409 : real(dp), dimension(:), allocatable, intent(out) :: wavelengths, flux
410 :
411 : character(len=512) :: line
412 : integer :: unit, n_rows, status, i, header_lines
413 : real(dp) :: temp_wavelength, temp_flux
414 :
415 :
416 0 : header_lines = 0
417 0 : open (newunit=unit, file=trim(directory), status='OLD', action='READ', iostat=status)
418 0 : if (status /= 0) then
419 0 : print *, "Error: Could not open file ", trim(directory)
420 0 : call mesa_error(__FILE__, __LINE__)
421 : end if
422 :
423 0 : do
424 0 : read (unit, '(A)', iostat=status) line
425 0 : if (status /= 0) exit
426 0 : if (line(1:1) /= "#") exit
427 0 : header_lines = header_lines + 1
428 : end do
429 :
430 : ! count data rows (we've already read one; count it)
431 : n_rows = 1
432 0 : do
433 0 : read (unit, '(A)', iostat=status) line
434 0 : if (status /= 0) exit
435 0 : n_rows = n_rows + 1
436 : end do
437 :
438 0 : allocate (wavelengths(n_rows))
439 0 : allocate (flux(n_rows))
440 :
441 0 : rewind (unit)
442 : ! skip exactly header_lines lines
443 0 : do i = 1, header_lines
444 0 : read (unit, '(A)', iostat=status) line
445 : end do
446 :
447 : i = 0
448 0 : do
449 0 : read (unit, *, iostat=status) temp_wavelength, temp_flux
450 0 : if (status /= 0) exit
451 0 : i = i + 1
452 0 : wavelengths(i) = temp_wavelength
453 0 : flux(i) = temp_flux
454 : end do
455 :
456 0 : close (unit)
457 :
458 0 : end subroutine load_sed
459 :
460 105 : function remove_dat(path) result(base)
461 : ! returns the portion of the string before the first dot
462 : character(len=*), intent(in) :: path
463 : character(len=strlen) :: base
464 : integer :: first_dot
465 :
466 105 : first_dot = 0
467 210 : do while (first_dot < len_trim(path) .and. path(first_dot + 1:first_dot + 1) /= '.')
468 105 : first_dot = first_dot + 1
469 : end do
470 :
471 105 : if (first_dot < len_trim(path)) then
472 105 : base = path(:first_dot)
473 : else
474 0 : base = path
475 : end if
476 105 : end function remove_dat
477 :
478 1 : function basename(path) result(name)
479 : character(len=*), intent(in) :: path
480 : character(len=strlen) :: name
481 : integer :: i
482 1 : if (len_trim(path) == 0) then
483 0 : name = ''
484 0 : return
485 : end if
486 1 : i = index(path, '/', back=.true.)
487 1 : name = path(i + 1:)
488 : end function basename
489 :
490 : ! resolve a path relative to mesa_dir unless it is already explicit.
491 : ! legacy colors defaults used '/data/...'; if that absolute-looking path
492 : ! does not exist on disk, warn and retry relative to mesa_dir.
493 21 : function resolve_path(path) result(full_path)
494 : use const_def, only: mesa_dir
495 : character(len=*), intent(in) :: path
496 : character(len=512) :: full_path
497 21 : character(len=:), allocatable :: p
498 : logical :: exists
499 : integer :: n
500 :
501 21 : exists = .false.
502 21 : p = trim(adjustl(path))
503 21 : n = len_trim(p)
504 :
505 21 : if (n >= 2 .and. p(1:2) == './') then
506 0 : full_path = p
507 21 : else if (n >= 3 .and. p(1:3) == '../') then
508 0 : full_path = p
509 :
510 21 : else if (n >= 1 .and. p(1:1) == '/') then
511 0 : inquire (file=p, exist=exists)
512 0 : if (.not. exists) inquire (file=p//'/.', exist=exists)
513 :
514 0 : if (exists) then
515 0 : full_path = p
516 : else
517 0 : write (*, '(5a)') 'colors: path ', trim(p), &
518 0 : ' not found; trying ', trim(mesa_dir)//trim(p), ' instead'
519 0 : full_path = trim(mesa_dir)//trim(p)
520 : end if
521 :
522 : else
523 21 : full_path = trim(mesa_dir)//'/'//trim(p)
524 : end if
525 42 : end function resolve_path
526 :
527 1 : subroutine read_strings_from_file(colors_settings, strings, n, ierr)
528 : character(len=512) :: filename
529 : character(len=100), allocatable, intent(out) :: strings(:)
530 : integer, intent(out) :: n, ierr
531 : integer :: unit, i, status
532 : character(len=100) :: line
533 : type(Colors_General_Info), pointer :: colors_settings
534 :
535 1 : ierr = 0
536 :
537 : filename = trim(resolve_path(colors_settings%instrument))//"/"// &
538 1 : trim(basename(colors_settings%instrument))
539 :
540 1 : n = 0
541 1 : open (newunit=unit, file=filename, status='old', action='read', iostat=status)
542 1 : if (status /= 0) then
543 0 : ierr = -1
544 0 : print *, "Error: Could not open file", filename
545 0 : call mesa_error(__FILE__, __LINE__)
546 : end if
547 :
548 7 : do
549 8 : read (unit, '(A)', iostat=status) line
550 8 : if (status /= 0) exit
551 7 : n = n + 1
552 : end do
553 1 : rewind (unit)
554 1 : if (allocated(strings)) deallocate (strings)
555 3 : allocate (strings(n))
556 8 : do i = 1, n
557 8 : read (unit, '(A)') strings(i)
558 : end do
559 1 : close (unit)
560 1 : end subroutine read_strings_from_file
561 :
562 : ! load flux cube from binary file into handle at initialization.
563 : ! if the cube is missing, malformed, or too large to allocate, keep
564 : ! cube_loaded = .false. and fall back to loading individual SED files.
565 : ! validate the header and on-disk size before reading the large 4-D payload
566 : ! so a stale or truncated cube reports a clear message instead of failing
567 : ! later inside the Fortran runtime.
568 1 : subroutine load_flux_cube(rq, stellar_model_dir)
569 : use, intrinsic :: iso_fortran_env, only: int64
570 : type(Colors_General_Info), intent(inout) :: rq
571 : character(len=*), intent(in) :: stellar_model_dir
572 :
573 : character(len=512) :: bin_filename
574 : integer :: unit, status, n_teff, n_logg, n_meta, n_lambda
575 : integer(int64) :: file_bytes, expected_bytes, header_bytes, real_bytes
576 : real(dp) :: cube_mb
577 :
578 1 : rq%cube_loaded = .false.
579 :
580 1 : bin_filename = trim(resolve_path(stellar_model_dir))//'/flux_cube.bin'
581 :
582 : open (newunit=unit, file=trim(bin_filename), status='OLD', &
583 1 : access='STREAM', form='UNFORMATTED', iostat=status)
584 1 : if (status /= 0) then
585 : ! no binary cube available -- will use individual SED files
586 0 : write (*, '(a)') 'colors: no flux_cube.bin found; using slower per-file SED loading'
587 1 : return
588 : end if
589 :
590 1 : read (unit, iostat=status) n_teff, n_logg, n_meta, n_lambda
591 1 : if (status /= 0) then
592 0 : write (*, '(a)') 'colors: could not read flux_cube.bin header; falling back to slower per-file SED loading'
593 0 : close (unit)
594 0 : return
595 : end if
596 :
597 1 : if (n_teff <= 0 .or. n_logg <= 0 .or. n_meta <= 0 .or. n_lambda <= 0) then
598 0 : write (*, '(a)') 'colors: invalid flux_cube.bin header; falling back to slower per-file SED loading'
599 0 : close (unit)
600 0 : return
601 : end if
602 :
603 1 : inquire (file=trim(bin_filename), size=file_bytes, iostat=status)
604 1 : if (status == 0) then
605 1 : header_bytes = 4_int64*int(storage_size(n_teff)/8, int64)
606 1 : real_bytes = int(storage_size(1.0_dp)/8, int64)
607 : expected_bytes = header_bytes + real_bytes*( &
608 : int(n_teff, int64) + int(n_logg, int64) + int(n_meta, int64) + int(n_lambda, int64) + &
609 1 : int(n_teff, int64)*int(n_logg, int64)*int(n_meta, int64)*int(n_lambda, int64))
610 1 : if (file_bytes /= expected_bytes) then
611 : write (*, '(a,i0,a,i0,a)') &
612 0 : 'colors: flux_cube.bin size mismatch (', file_bytes, &
613 0 : ' bytes, expected ', expected_bytes, &
614 0 : '); falling back to slower per-file SED loading'
615 0 : close (unit)
616 0 : return
617 : end if
618 : end if
619 :
620 : ! attempt the large allocation first -- this is the one that may fail.
621 : ! doing it before the small grid allocations avoids partial cleanup.
622 6 : allocate (rq%cube_flux(n_teff, n_logg, n_meta, n_lambda), stat=status)
623 1 : if (status /= 0) then
624 0 : if (allocated(rq%cube_flux)) deallocate (rq%cube_flux)
625 0 : cube_mb = real(n_teff, dp)*n_logg*n_meta*n_lambda*8.0_dp/(1024.0_dp**2)
626 : write (*, '(a,f0.1,a)') &
627 0 : 'colors: flux cube allocation failed (', cube_mb, &
628 0 : ' MB); falling back to slower per-file SED loading'
629 0 : close (unit)
630 0 : return
631 : end if
632 :
633 : ! grid arrays are small -- always expected to succeed
634 3 : allocate (rq%cube_teff_grid(n_teff), stat=status)
635 1 : if (status /= 0) goto 900
636 :
637 3 : allocate (rq%cube_logg_grid(n_logg), stat=status)
638 1 : if (status /= 0) goto 900
639 :
640 3 : allocate (rq%cube_meta_grid(n_meta), stat=status)
641 1 : if (status /= 0) goto 900
642 :
643 3 : allocate (rq%cube_wavelengths(n_lambda), stat=status)
644 1 : if (status /= 0) goto 900
645 :
646 1 : read (unit, iostat=status) rq%cube_teff_grid
647 1 : if (status /= 0) goto 900
648 :
649 1 : read (unit, iostat=status) rq%cube_logg_grid
650 1 : if (status /= 0) goto 900
651 :
652 1 : read (unit, iostat=status) rq%cube_meta_grid
653 1 : if (status /= 0) goto 900
654 :
655 1 : read (unit, iostat=status) rq%cube_wavelengths
656 1 : if (status /= 0) goto 900
657 :
658 1 : read (unit, iostat=status) rq%cube_flux
659 1 : if (status /= 0) goto 900
660 :
661 1 : close (unit)
662 1 : rq%cube_loaded = .true.
663 :
664 1 : cube_mb = real(n_teff, dp)*n_logg*n_meta*n_lambda*8.0_dp/(1024.0_dp**2)
665 : write (*, '(a,i0,a,i0,a,i0,a,i0,a,f0.1,a)') &
666 1 : 'colors: flux cube loaded (', &
667 1 : n_teff, ' x ', n_logg, ' x ', n_meta, ' x ', n_lambda, &
668 2 : ', ', cube_mb, ' MB)'
669 1 : return
670 :
671 : ! error cleanup -- deallocate everything that may have been allocated
672 : 900 continue
673 0 : write (*, '(a)') 'colors: error reading flux_cube.bin; falling back to slower per-file SED loading'
674 0 : if (allocated(rq%cube_flux)) deallocate (rq%cube_flux)
675 0 : if (allocated(rq%cube_teff_grid)) deallocate (rq%cube_teff_grid)
676 0 : if (allocated(rq%cube_logg_grid)) deallocate (rq%cube_logg_grid)
677 0 : if (allocated(rq%cube_meta_grid)) deallocate (rq%cube_meta_grid)
678 0 : if (allocated(rq%cube_wavelengths)) deallocate (rq%cube_wavelengths)
679 0 : close (unit)
680 :
681 : end subroutine load_flux_cube
682 :
683 : ! build unique sorted grids from lookup table and store on handle.
684 : ! called once at init so the fallback interpolation path never rebuilds these.
685 1 : subroutine build_unique_grids(rq)
686 : type(Colors_General_Info), intent(inout) :: rq
687 : logical :: found
688 :
689 1 : if (rq%unique_grids_built) return
690 1 : if (.not. rq%lookup_loaded) return
691 :
692 1 : call extract_unique_sorted(rq%lu_teff, rq%u_teff)
693 1 : call extract_unique_sorted(rq%lu_logg, rq%u_logg)
694 1 : call extract_unique_sorted(rq%lu_meta, rq%u_meta)
695 :
696 1 : rq%unique_grids_built = .true.
697 :
698 : contains
699 :
700 3 : subroutine extract_unique_sorted(arr, unique)
701 : real(dp), intent(in) :: arr(:)
702 : real(dp), allocatable, intent(out) :: unique(:)
703 3 : real(dp), allocatable :: buf(:)
704 : integer :: ii, jj, nn, nnu
705 : real(dp) :: sw
706 :
707 3 : nn = size(arr)
708 9 : allocate (buf(nn))
709 11427 : nnu = 0
710 :
711 11427 : do ii = 1, nn
712 11424 : found = .false.
713 194760 : do jj = 1, nnu
714 194760 : if (abs(arr(ii) - buf(jj)) < 1.0e-10_dp) then
715 11329 : found = .true.
716 11329 : exit
717 : end if
718 : end do
719 11427 : if (.not. found) then
720 95 : nnu = nnu + 1
721 95 : buf(nnu) = arr(ii)
722 : end if
723 : end do
724 :
725 : ! insertion sort (grids are small)
726 95 : do ii = 2, nnu
727 92 : sw = buf(ii)
728 92 : jj = ii - 1
729 1376 : do while (jj >= 1 .and. buf(jj) > sw)
730 1284 : buf(jj + 1) = buf(jj)
731 1370 : jj = jj - 1
732 : end do
733 95 : buf(jj + 1) = sw
734 : end do
735 :
736 9 : allocate (unique(nnu))
737 101 : unique = buf(1:nnu)
738 3 : deallocate (buf)
739 3 : end subroutine extract_unique_sorted
740 :
741 : end subroutine build_unique_grids
742 :
743 : ! build a 3-D mapping from unique-grid indices to lookup-table rows.
744 : ! grid_to_lu(i_t, i_g, i_m) = row in lu_* whose parameters match
745 : ! (u_teff(i_t), u_logg(i_g), u_meta(i_m)).
746 : ! called once at init; replaces the O(n_lu) nearest-neighbour scan
747 : ! that previously ran per corner at runtime.
748 1 : subroutine build_grid_to_lu_map(rq)
749 : type(Colors_General_Info), intent(inout) :: rq
750 :
751 : integer :: nt, ng, nm, n_lu
752 : integer :: it, ig, im, idx
753 : integer :: best_idx
754 : real(dp) :: best_dist, dist
755 : real(dp), parameter :: tol = 1.0e-10_dp
756 :
757 1 : if (rq%grid_map_built) return
758 1 : if (.not. rq%unique_grids_built) return
759 1 : if (.not. rq%lookup_loaded) return
760 :
761 1 : nt = size(rq%u_teff)
762 1 : ng = size(rq%u_logg)
763 1 : nm = size(rq%u_meta)
764 1 : n_lu = size(rq%lu_teff)
765 :
766 5 : allocate (rq%grid_to_lu(nt, ng, nm))
767 6785 : rq%grid_to_lu = 0
768 :
769 77 : do it = 1, nt
770 913 : do ig = 1, ng
771 7600 : do im = 1, nm
772 : best_idx = 1
773 : best_dist = huge(1.0_dp)
774 18222256 : do idx = 1, n_lu
775 : dist = abs(rq%lu_teff(idx) - rq%u_teff(it)) + &
776 : abs(rq%lu_logg(idx) - rq%u_logg(ig)) + &
777 18219376 : abs(rq%lu_meta(idx) - rq%u_meta(im))
778 18219376 : if (dist < best_dist) then
779 358488 : best_dist = dist
780 358488 : best_idx = idx
781 : end if
782 18222256 : if (best_dist < tol) exit ! exact match found
783 : end do
784 7524 : rq%grid_to_lu(it, ig, im) = best_idx
785 : end do
786 : end do
787 : end do
788 :
789 1 : rq%grid_map_built = .true.
790 :
791 : end subroutine build_grid_to_lu_map
792 :
793 20 : subroutine find_containing_cell(x_val, y_val, z_val, x_grid, y_grid, z_grid, &
794 : i_x, i_y, i_z, t_x, t_y, t_z)
795 : real(dp), intent(in) :: x_val, y_val, z_val
796 : real(dp), intent(in) :: x_grid(:), y_grid(:), z_grid(:)
797 : integer, intent(out) :: i_x, i_y, i_z
798 : real(dp), intent(out) :: t_x, t_y, t_z
799 :
800 20 : call find_interval(x_grid, x_val, i_x, t_x)
801 20 : call find_interval(y_grid, y_val, i_y, t_y)
802 20 : call find_interval(z_grid, z_val, i_z, t_z)
803 20 : end subroutine find_containing_cell
804 :
805 : ! find the interval in a sorted array containing a value.
806 : ! returns index i such that x(i) <= val <= x(i+1) and fractional position t in [0,1].
807 : ! detects dummy axes (all zeros / 999 / -999) and collapses them to i=1, t=0.
808 60 : subroutine find_interval(x, val, i, t)
809 : real(dp), intent(in) :: x(:), val
810 : integer, intent(out) :: i
811 : real(dp), intent(out) :: t
812 :
813 : integer :: n, lo, hi, mid
814 : logical :: dummy_axis
815 :
816 60 : n = size(x)
817 :
818 : ! detect dummy axis: all values == 0, 999, or -999
819 200 : dummy_axis = all(x == 0.0_dp) .or. all(x == 999.0_dp) .or. all(x == -999.0_dp)
820 :
821 : if (dummy_axis) then
822 0 : i = 1
823 0 : t = 0.0_dp
824 0 : return
825 : end if
826 :
827 60 : if (val <= x(1)) then
828 2 : i = 1
829 2 : t = 0.0_dp
830 2 : return
831 58 : else if (val >= x(n)) then
832 2 : i = n - 1
833 2 : t = 1.0_dp
834 2 : return
835 : end if
836 :
837 : lo = 1
838 : hi = n
839 299 : do while (hi - lo > 1)
840 243 : mid = (lo + hi)/2
841 299 : if (val >= x(mid)) then
842 : lo = mid
843 : else
844 172 : hi = mid
845 : end if
846 : end do
847 :
848 56 : i = lo
849 56 : if (abs(x(i + 1) - x(i)) < 1.0e-30_dp) then
850 0 : t = 0.0_dp ! degenerate interval -- no interpolation needed
851 : else
852 56 : t = (val - x(i))/(x(i + 1) - x(i))
853 : end if
854 : end subroutine find_interval
855 :
856 0 : subroutine find_nearest_point(x_val, y_val, z_val, x_grid, y_grid, z_grid, &
857 : i_x, i_y, i_z)
858 : real(dp), intent(in) :: x_val, y_val, z_val
859 : real(dp), intent(in) :: x_grid(:), y_grid(:), z_grid(:)
860 : integer, intent(out) :: i_x, i_y, i_z
861 :
862 0 : i_x = minloc(abs(x_val - x_grid), 1)
863 0 : i_y = minloc(abs(y_val - y_grid), 1)
864 0 : i_z = minloc(abs(z_val - z_grid), 1)
865 0 : end subroutine find_nearest_point
866 :
867 : ! returns idx such that grid(idx) <= val < grid(idx+1), clamped to bounds
868 0 : subroutine find_bracket_index(grid, val, idx)
869 : real(dp), intent(in) :: grid(:), val
870 : integer, intent(out) :: idx
871 :
872 : integer :: n, lo, hi, mid
873 :
874 0 : n = size(grid)
875 0 : if (n < 2) then
876 0 : idx = 1
877 0 : return
878 : end if
879 :
880 0 : if (val <= grid(1)) then
881 0 : idx = 1
882 0 : return
883 0 : else if (val >= grid(n)) then
884 0 : idx = n - 1
885 0 : return
886 : end if
887 :
888 : lo = 1
889 : hi = n
890 0 : do while (hi - lo > 1)
891 0 : mid = (lo + hi)/2
892 0 : if (val >= grid(mid)) then
893 : lo = mid
894 : else
895 0 : hi = mid
896 : end if
897 : end do
898 0 : idx = lo
899 : end subroutine find_bracket_index
900 :
901 : ! fallback SED cache and stencil loader
902 :
903 : ! load a stencil sub-cube of SED fluxes for the given index ranges.
904 : ! uses load_sed_cached so repeated visits to the same grid point are
905 : ! served from memory rather than disk.
906 0 : subroutine load_stencil(rq, resolved_dir, lo_t, hi_t, lo_g, hi_g, lo_m, hi_m)
907 : type(Colors_General_Info), intent(inout) :: rq
908 : character(len=*), intent(in) :: resolved_dir
909 : integer, intent(in) :: lo_t, hi_t, lo_g, hi_g, lo_m, hi_m
910 :
911 : integer :: st, sg, sm, n_lambda, lu_idx
912 : integer :: it, ig, im
913 0 : real(dp), dimension(:), allocatable :: sed_flux
914 :
915 0 : st = hi_t - lo_t + 1
916 0 : sg = hi_g - lo_g + 1
917 0 : sm = hi_m - lo_m + 1
918 :
919 : ! free previous stencil flux data
920 0 : if (allocated(rq%stencil_fluxes)) deallocate (rq%stencil_fluxes)
921 :
922 0 : n_lambda = 0
923 :
924 0 : do it = lo_t, hi_t
925 0 : do ig = lo_g, hi_g
926 0 : do im = lo_m, hi_m
927 0 : lu_idx = rq%grid_to_lu(it, ig, im)
928 0 : call load_sed_cached(rq, resolved_dir, lu_idx, sed_flux)
929 :
930 0 : if (n_lambda == 0) then
931 0 : n_lambda = size(sed_flux)
932 0 : allocate (rq%stencil_fluxes(st, sg, sm, n_lambda))
933 0 : else if (size(sed_flux) /= n_lambda) then
934 : ! SED files have inconsistent wavelength counts — this is a
935 : ! fatal grid inconsistency; crash with a clear message.
936 : write (*, '(a,i0,a,i0,a,i0)') &
937 0 : 'colors ERROR: SED at lu_idx=', lu_idx, &
938 0 : ' has ', size(sed_flux), ' wavelength points; expected ', n_lambda
939 : write (*, '(a)') &
940 0 : 'colors: bt-settl (or other) grid files have non-uniform wavelength grids.'
941 0 : call mesa_error(__FILE__, __LINE__)
942 : end if
943 :
944 : rq%stencil_fluxes(it - lo_t + 1, ig - lo_g + 1, im - lo_m + 1, :) = &
945 0 : sed_flux(1:n_lambda)
946 0 : if (allocated(sed_flux)) deallocate (sed_flux)
947 : end do
948 : end do
949 : end do
950 :
951 : ! set stencil wavelengths from the canonical copy on the handle
952 0 : if (allocated(rq%stencil_wavelengths)) deallocate (rq%stencil_wavelengths)
953 0 : allocate (rq%stencil_wavelengths(n_lambda))
954 0 : rq%stencil_wavelengths = rq%fallback_wavelengths(1:n_lambda)
955 :
956 0 : end subroutine load_stencil
957 :
958 : ! retrieve an SED flux from the memory cache, or load from disk on miss.
959 : ! uses a bounded circular buffer (sed_mem_cache_cap slots).
960 : ! on the first disk read, the wavelength array is stored once on the handle
961 : ! as rq%fallback_wavelengths -- all SEDs in a given atmosphere grid share
962 : ! the same wavelengths, so only flux is cached and returned.
963 0 : subroutine load_sed_cached(rq, resolved_dir, lu_idx, flux)
964 : type(Colors_General_Info), intent(inout) :: rq
965 : character(len=*), intent(in) :: resolved_dir
966 : integer, intent(in) :: lu_idx
967 : real(dp), dimension(:), allocatable, intent(out) :: flux
968 :
969 : integer :: slot, n_lam, status
970 : character(len=512) :: filepath
971 0 : real(dp), dimension(:), allocatable :: sed_wave
972 0 : real(dp), dimension(:), allocatable :: flux_interp
973 :
974 : ! initialise the cache on first call
975 0 : if (.not. rq%sed_mcache_init) then
976 0 : allocate (rq%sed_mcache_keys(sed_mem_cache_cap))
977 0 : rq%sed_mcache_keys = 0 ! 0 means empty slot
978 0 : rq%sed_mcache_count = 0
979 0 : rq%sed_mcache_next = 1
980 0 : rq%sed_mcache_nlam = 0
981 0 : rq%sed_mcache_init = .true.
982 : end if
983 :
984 : ! search for a cache hit (linear scan over a small array)
985 0 : do slot = 1, rq%sed_mcache_count
986 0 : if (rq%sed_mcache_keys(slot) == lu_idx) then
987 : ! hit -- return cached flux
988 0 : n_lam = rq%sed_mcache_nlam
989 0 : allocate (flux(n_lam))
990 0 : flux = rq%sed_mcache_data(:, slot)
991 0 : return
992 : end if
993 : end do
994 :
995 : ! miss -- load from disk
996 0 : filepath = trim(resolved_dir)//'/'//trim(rq%lu_file_names(lu_idx))
997 0 : call load_sed(filepath, lu_idx, sed_wave, flux)
998 :
999 : ! store the canonical wavelength array on the handle (once only)
1000 0 : if (.not. rq%fallback_wavelengths_set) then
1001 0 : n_lam = size(sed_wave)
1002 0 : allocate (rq%fallback_wavelengths(n_lam))
1003 0 : rq%fallback_wavelengths = sed_wave
1004 0 : rq%fallback_wavelengths_set = .true.
1005 : end if
1006 :
1007 :
1008 :
1009 : ! store flux in the cache
1010 0 : n_lam = size(flux)
1011 0 : if (rq%sed_mcache_nlam == 0) then
1012 0 : rq%sed_mcache_nlam = n_lam
1013 0 : allocate (rq%sed_mcache_data(n_lam, sed_mem_cache_cap), stat=status)
1014 0 : if (status /= 0) then
1015 0 : write (*, '(a,f0.1,a)') 'colors: SED memory cache allocation failed (', &
1016 0 : real(n_lam, dp)*sed_mem_cache_cap*8.0_dp/1024.0_dp**2, ' MB); cache disabled'
1017 0 : rq%sed_mcache_init = .false. ! fall through to disk-only path
1018 0 : if (allocated(sed_wave)) deallocate (sed_wave)
1019 0 : return
1020 : end if
1021 0 : else if (n_lam /= rq%sed_mcache_nlam) then
1022 : ! non-uniform wavelength grids across SED files: interpolate onto the canonical grid
1023 0 : if (.not. rq%fallback_wavelengths_set) then
1024 0 : allocate (rq%fallback_wavelengths(n_lam))
1025 0 : rq%fallback_wavelengths = sed_wave
1026 0 : rq%fallback_wavelengths_set = .true.
1027 : end if
1028 :
1029 0 : allocate (flux_interp(rq%sed_mcache_nlam), stat=status)
1030 0 : if (status /= 0) then
1031 0 : write (*, '(a)') 'colors: WARNING: could not allocate interpolation buffer; SED cache disabled'
1032 0 : rq%sed_mcache_init = .false.
1033 0 : if (allocated(sed_wave)) deallocate (sed_wave)
1034 0 : return
1035 : end if
1036 :
1037 0 : call interp_linear_internal(rq%fallback_wavelengths, sed_wave, flux, flux_interp)
1038 :
1039 0 : deallocate (flux)
1040 0 : allocate (flux(rq%sed_mcache_nlam))
1041 0 : flux = flux_interp
1042 0 : deallocate (flux_interp)
1043 :
1044 0 : n_lam = rq%sed_mcache_nlam
1045 : end if
1046 :
1047 : ! write to the next slot (circular)
1048 0 : slot = rq%sed_mcache_next
1049 0 : rq%sed_mcache_keys(slot) = lu_idx
1050 0 : rq%sed_mcache_data(:, slot) = flux(1:n_lam)
1051 :
1052 : ! advance the circular pointer
1053 0 : if (rq%sed_mcache_count < sed_mem_cache_cap) &
1054 0 : rq%sed_mcache_count = rq%sed_mcache_count + 1
1055 0 : rq%sed_mcache_next = mod(slot, sed_mem_cache_cap) + 1
1056 :
1057 0 : if (allocated(sed_wave)) deallocate (sed_wave)
1058 0 : end subroutine load_sed_cached
1059 :
1060 :
1061 : ! simple 1D linear interpolation (np.interp-like): clamps to endpoints
1062 0 : subroutine interp_linear_internal(x_out, x_in, y_in, y_out)
1063 : real(dp), intent(in) :: x_out(:), x_in(:), y_in(:)
1064 : real(dp), intent(out) :: y_out(:)
1065 : integer :: i, n_in, n_out, lo, hi, mid
1066 : real(dp) :: t, denom
1067 :
1068 0 : n_in = size(x_in)
1069 0 : n_out = size(x_out)
1070 :
1071 0 : if (n_out <= 0) return
1072 :
1073 0 : if (n_in <= 0) then
1074 0 : y_out = 0.0_dp
1075 : return
1076 : end if
1077 :
1078 0 : if (n_in == 1) then
1079 0 : y_out = y_in(1)
1080 : return
1081 : end if
1082 :
1083 0 : do i = 1, n_out
1084 0 : if (x_out(i) <= x_in(1)) then
1085 0 : y_out(i) = y_in(1)
1086 0 : else if (x_out(i) >= x_in(n_in)) then
1087 0 : y_out(i) = y_in(n_in)
1088 : else
1089 : lo = 1
1090 : hi = n_in
1091 0 : do while (hi - lo > 1)
1092 0 : mid = (lo + hi)/2
1093 0 : if (x_out(i) >= x_in(mid)) then
1094 : lo = mid
1095 : else
1096 0 : hi = mid
1097 : end if
1098 : end do
1099 :
1100 0 : denom = x_in(lo + 1) - x_in(lo)
1101 0 : if (abs(denom) <= 0.0_dp) then
1102 0 : y_out(i) = y_in(lo)
1103 : else
1104 0 : t = (x_out(i) - x_in(lo))/denom
1105 0 : y_out(i) = y_in(lo) + t*(y_in(lo + 1) - y_in(lo))
1106 : end if
1107 : end if
1108 : end do
1109 :
1110 : end subroutine interp_linear_internal
1111 : end module colors_utils
|
