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_def
21 :
22 : use const_def, only: dp
23 :
24 : implicit none
25 :
26 : public
27 :
28 : ! max number of SEDs to keep in the memory cache
29 : ! each slot holds one wavelength array (~1200 doubles ~ 10 KB),
30 : ! so 256 slots ~ 2.5 MB -- negligible even when the full cube
31 : ! cannot be allocated
32 : integer, parameter :: sed_mem_cache_cap = 256
33 :
34 : type :: filter_data
35 : character(len=100) :: name
36 : real(dp), allocatable :: wavelengths(:)
37 : real(dp), allocatable :: transmission(:)
38 : ! precomputed zero-point fluxes, computed once at init
39 : real(dp) :: vega_zero_point = -1.0_dp
40 : real(dp) :: ab_zero_point = -1.0_dp
41 : real(dp) :: st_zero_point = -1.0_dp
42 : end type filter_data
43 :
44 : type :: Colors_General_Info
45 : character(len=256) :: instrument
46 : character(len=256) :: vega_sed
47 : character(len=256) :: stellar_atm
48 : character(len=256) :: colors_results_directory
49 : character(len=256) :: mag_system
50 : real(dp) :: metallicity
51 : real(dp) :: distance
52 : real(dp) :: z_over_x_ref
53 : logical :: make_csv
54 : logical :: sed_per_model
55 : logical :: use_colors
56 : integer :: handle
57 : logical :: in_use
58 :
59 : ! cached lookup table data
60 : logical :: lookup_loaded = .false.
61 : character(len=100), allocatable :: lu_file_names(:)
62 : real(dp), allocatable :: lu_logg(:)
63 : real(dp), allocatable :: lu_meta(:)
64 : real(dp), allocatable :: lu_teff(:)
65 :
66 : ! cached vega SED
67 : logical :: vega_loaded = .false.
68 : real(dp), allocatable :: vega_wavelengths(:)
69 : real(dp), allocatable :: vega_fluxes(:)
70 :
71 : ! cached filter data (includes precomputed zero-points)
72 : logical :: filters_loaded = .false.
73 : type(filter_data), allocatable :: filters(:)
74 :
75 : ! cached flux cube
76 : logical :: cube_loaded = .false.
77 : real(dp), allocatable :: cube_flux(:, :, :, :) ! (n_teff, n_logg, n_meta, n_lambda)
78 : real(dp), allocatable :: cube_teff_grid(:)
79 : real(dp), allocatable :: cube_logg_grid(:)
80 : real(dp), allocatable :: cube_meta_grid(:)
81 : real(dp), allocatable :: cube_wavelengths(:)
82 :
83 : ! unique sorted grids (built once from lookup table at init)
84 : logical :: unique_grids_built = .false.
85 : real(dp), allocatable :: u_teff(:), u_logg(:), u_meta(:)
86 :
87 : ! grid_to_lu(i_t, i_g, i_m) gives the lookup-table row index for
88 : ! (u_teff(i_t), u_logg(i_g), u_meta(i_m)) -- avoids O(n_lu)
89 : ! nearest-neighbour searches at runtime
90 : logical :: grid_map_built = .false.
91 : integer, allocatable :: grid_to_lu(:, :, :)
92 :
93 : ! fallback-path caches (used only when cube_loaded == .false.)
94 :
95 : ! stencil cache: the extended neighbourhood around the current
96 : ! interpolation cell, includes derivative-context points
97 : ! (i-1 .. i+2 per axis, clamped to boundaries) so that
98 : ! hermite_tensor_interp3d gives the same result as the cube path
99 : logical :: stencil_valid = .false.
100 : integer :: stencil_i_t = -1, stencil_i_g = -1, stencil_i_m = -1
101 : real(dp), allocatable :: stencil_fluxes(:, :, :, :) ! (st, sg, sm, n_lambda)
102 : real(dp), allocatable :: stencil_wavelengths(:) ! (n_lambda)
103 : real(dp), allocatable :: stencil_teff(:) ! subgrid values (st)
104 : real(dp), allocatable :: stencil_logg(:) ! subgrid values (sg)
105 : real(dp), allocatable :: stencil_meta(:) ! subgrid values (sm)
106 :
107 : ! canonical wavelength grid for fallback SEDs (set once on first disk
108 : ! read -- all SEDs in a given atmosphere grid share the same wavelengths)
109 : logical :: fallback_wavelengths_set = .false.
110 : real(dp), allocatable :: fallback_wavelengths(:) ! (n_lambda)
111 :
112 : ! bounded SED memory cache (circular buffer, keyed by lu index)
113 : ! avoids re-reading text files for SEDs we've already parsed
114 : logical :: sed_mcache_init = .false.
115 : integer :: sed_mcache_count = 0
116 : integer :: sed_mcache_next = 1
117 : integer :: sed_mcache_nlam = 0
118 : integer, allocatable :: sed_mcache_keys(:) ! (sed_mem_cache_cap)
119 : real(dp), allocatable :: sed_mcache_data(:, :) ! (n_lambda, sed_mem_cache_cap)
120 :
121 : end type Colors_General_Info
122 :
123 : ! Global filter name list (shared across handles)
124 : integer :: num_color_filters
125 : character(len=100), allocatable :: color_filter_names(:)
126 :
127 : integer, parameter :: max_colors_handles = 10
128 : type(Colors_General_Info), target :: colors_handles(max_colors_handles)
129 :
130 : logical :: colors_is_initialized = .false.
131 :
132 : character(len=1000) :: colors_dir, colors_cache_dir, colors_temp_cache_dir
133 : logical :: colors_use_cache = .true.
134 :
135 : contains
136 :
137 2 : subroutine colors_def_init(colors_cache_dir_in)
138 : use utils_lib, only: mkdir
139 : use const_def, only: mesa_data_dir, mesa_caches_dir, mesa_temp_caches_dir, use_mesa_temp_cache
140 : character(*), intent(in) :: colors_cache_dir_in
141 : integer :: i
142 :
143 2 : if (len_trim(colors_cache_dir_in) > 0) then
144 0 : colors_cache_dir = colors_cache_dir_in
145 2 : else if (len_trim(mesa_caches_dir) > 0) then
146 0 : colors_cache_dir = trim(mesa_caches_dir)//'/colors_cache'
147 : else
148 2 : colors_cache_dir = trim(mesa_data_dir)//'/colors_data/cache'
149 : end if
150 2 : call mkdir(colors_cache_dir)
151 :
152 22 : do i = 1, max_colors_handles
153 20 : colors_handles(i)%handle = i
154 20 : colors_handles(i)%in_use = .false.
155 20 : colors_handles(i)%lookup_loaded = .false.
156 20 : colors_handles(i)%vega_loaded = .false.
157 20 : colors_handles(i)%filters_loaded = .false.
158 20 : colors_handles(i)%cube_loaded = .false.
159 20 : colors_handles(i)%unique_grids_built = .false.
160 20 : colors_handles(i)%grid_map_built = .false.
161 20 : colors_handles(i)%stencil_valid = .false.
162 20 : colors_handles(i)%stencil_i_t = -1
163 20 : colors_handles(i)%stencil_i_g = -1
164 20 : colors_handles(i)%stencil_i_m = -1
165 20 : colors_handles(i)%sed_mcache_init = .false.
166 20 : colors_handles(i)%sed_mcache_count = 0
167 20 : colors_handles(i)%sed_mcache_next = 1
168 20 : colors_handles(i)%sed_mcache_nlam = 0
169 22 : colors_handles(i)%fallback_wavelengths_set = .false.
170 : end do
171 :
172 2 : colors_temp_cache_dir = trim(mesa_temp_caches_dir)//'/colors_cache'
173 2 : if (use_mesa_temp_cache) call mkdir(colors_temp_cache_dir)
174 :
175 2 : end subroutine colors_def_init
176 :
177 2 : integer function do_alloc_colors(ierr)
178 : integer, intent(out) :: ierr
179 : integer :: i
180 2 : ierr = 0
181 2 : do_alloc_colors = -1
182 4 : !$omp critical (colors_handle)
183 2 : do i = 1, max_colors_handles
184 2 : if (.not. colors_handles(i)%in_use) then
185 2 : colors_handles(i)%in_use = .true.
186 2 : do_alloc_colors = i
187 2 : exit
188 : end if
189 : end do
190 : !$omp end critical (colors_handle)
191 2 : if (do_alloc_colors == -1) then
192 0 : ierr = -1
193 0 : return
194 : end if
195 2 : if (colors_handles(do_alloc_colors)%handle /= do_alloc_colors) then
196 0 : ierr = -1
197 0 : return
198 : end if
199 : end function do_alloc_colors
200 :
201 2 : subroutine do_free_colors(handle)
202 : integer, intent(in) :: handle
203 2 : if (handle >= 1 .and. handle <= max_colors_handles) then
204 2 : colors_handles(handle)%in_use = .false.
205 2 : call free_colors_cache(handle)
206 : end if
207 2 : end subroutine do_free_colors
208 :
209 22 : subroutine free_colors_cache(handle)
210 : integer, intent(in) :: handle
211 : integer :: i
212 :
213 22 : if (handle < 1 .or. handle > max_colors_handles) return
214 :
215 22 : if (allocated(colors_handles(handle)%lu_file_names)) &
216 1 : deallocate (colors_handles(handle)%lu_file_names)
217 22 : if (allocated(colors_handles(handle)%lu_logg)) &
218 1 : deallocate (colors_handles(handle)%lu_logg)
219 22 : if (allocated(colors_handles(handle)%lu_meta)) &
220 1 : deallocate (colors_handles(handle)%lu_meta)
221 22 : if (allocated(colors_handles(handle)%lu_teff)) &
222 1 : deallocate (colors_handles(handle)%lu_teff)
223 22 : colors_handles(handle)%lookup_loaded = .false.
224 :
225 22 : if (allocated(colors_handles(handle)%vega_wavelengths)) &
226 1 : deallocate (colors_handles(handle)%vega_wavelengths)
227 22 : if (allocated(colors_handles(handle)%vega_fluxes)) &
228 1 : deallocate (colors_handles(handle)%vega_fluxes)
229 22 : colors_handles(handle)%vega_loaded = .false.
230 :
231 22 : if (allocated(colors_handles(handle)%filters)) then
232 8 : do i = 1, size(colors_handles(handle)%filters)
233 7 : if (allocated(colors_handles(handle)%filters(i)%wavelengths)) &
234 7 : deallocate (colors_handles(handle)%filters(i)%wavelengths)
235 7 : if (allocated(colors_handles(handle)%filters(i)%transmission)) &
236 8 : deallocate (colors_handles(handle)%filters(i)%transmission)
237 : end do
238 8 : deallocate (colors_handles(handle)%filters)
239 : end if
240 22 : colors_handles(handle)%filters_loaded = .false.
241 :
242 22 : if (allocated(colors_handles(handle)%cube_flux)) &
243 1 : deallocate (colors_handles(handle)%cube_flux)
244 22 : if (allocated(colors_handles(handle)%cube_teff_grid)) &
245 1 : deallocate (colors_handles(handle)%cube_teff_grid)
246 22 : if (allocated(colors_handles(handle)%cube_logg_grid)) &
247 1 : deallocate (colors_handles(handle)%cube_logg_grid)
248 22 : if (allocated(colors_handles(handle)%cube_meta_grid)) &
249 1 : deallocate (colors_handles(handle)%cube_meta_grid)
250 22 : if (allocated(colors_handles(handle)%cube_wavelengths)) &
251 1 : deallocate (colors_handles(handle)%cube_wavelengths)
252 22 : colors_handles(handle)%cube_loaded = .false.
253 :
254 22 : if (allocated(colors_handles(handle)%u_teff)) &
255 1 : deallocate (colors_handles(handle)%u_teff)
256 22 : if (allocated(colors_handles(handle)%u_logg)) &
257 1 : deallocate (colors_handles(handle)%u_logg)
258 22 : if (allocated(colors_handles(handle)%u_meta)) &
259 1 : deallocate (colors_handles(handle)%u_meta)
260 22 : colors_handles(handle)%unique_grids_built = .false.
261 :
262 22 : if (allocated(colors_handles(handle)%grid_to_lu)) &
263 1 : deallocate (colors_handles(handle)%grid_to_lu)
264 22 : colors_handles(handle)%grid_map_built = .false.
265 :
266 22 : if (allocated(colors_handles(handle)%stencil_fluxes)) &
267 0 : deallocate (colors_handles(handle)%stencil_fluxes)
268 22 : if (allocated(colors_handles(handle)%stencil_wavelengths)) &
269 0 : deallocate (colors_handles(handle)%stencil_wavelengths)
270 22 : if (allocated(colors_handles(handle)%stencil_teff)) &
271 0 : deallocate (colors_handles(handle)%stencil_teff)
272 22 : if (allocated(colors_handles(handle)%stencil_logg)) &
273 0 : deallocate (colors_handles(handle)%stencil_logg)
274 22 : if (allocated(colors_handles(handle)%stencil_meta)) &
275 0 : deallocate (colors_handles(handle)%stencil_meta)
276 22 : colors_handles(handle)%stencil_valid = .false.
277 22 : colors_handles(handle)%stencil_i_t = -1
278 22 : colors_handles(handle)%stencil_i_g = -1
279 22 : colors_handles(handle)%stencil_i_m = -1
280 :
281 22 : if (allocated(colors_handles(handle)%sed_mcache_keys)) &
282 0 : deallocate (colors_handles(handle)%sed_mcache_keys)
283 22 : if (allocated(colors_handles(handle)%sed_mcache_data)) &
284 0 : deallocate (colors_handles(handle)%sed_mcache_data)
285 22 : colors_handles(handle)%sed_mcache_init = .false.
286 22 : colors_handles(handle)%sed_mcache_count = 0
287 22 : colors_handles(handle)%sed_mcache_next = 1
288 22 : colors_handles(handle)%sed_mcache_nlam = 0
289 :
290 22 : if (allocated(colors_handles(handle)%fallback_wavelengths)) &
291 0 : deallocate (colors_handles(handle)%fallback_wavelengths)
292 22 : colors_handles(handle)%fallback_wavelengths_set = .false.
293 :
294 : end subroutine free_colors_cache
295 :
296 28 : subroutine get_colors_ptr(handle, rq, ierr)
297 : integer, intent(in) :: handle
298 : type(Colors_General_Info), pointer, intent(out) :: rq
299 : integer, intent(out):: ierr
300 28 : if (handle < 1 .or. handle > max_colors_handles) then
301 0 : ierr = -1
302 0 : return
303 : end if
304 28 : rq => colors_handles(handle)
305 28 : ierr = 0
306 : end subroutine get_colors_ptr
307 :
308 2 : subroutine do_free_colors_tables
309 : integer :: i
310 :
311 2 : if (allocated(color_filter_names)) deallocate (color_filter_names)
312 :
313 22 : do i = 1, max_colors_handles
314 22 : call free_colors_cache(i)
315 : end do
316 :
317 2 : end subroutine do_free_colors_tables
318 :
319 0 : end module colors_def
|
