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_lib
21 :
22 : use const_def, only: dp, strlen, mesa_dir
23 : use bolometric, only: calculate_bolometric
24 : use synthetic, only: calculate_synthetic
25 : use colors_utils, only: read_strings_from_file, load_lookup_table, load_filter, load_vega_sed, &
26 : resolve_path, load_flux_cube, build_unique_grids, build_grid_to_lu_map
27 : use colors_history, only: how_many_colors_history_columns, data_for_colors_history_columns
28 :
29 : implicit none
30 :
31 : private
32 :
33 : public :: colors_init, colors_shutdown
34 : public :: alloc_colors_handle, alloc_colors_handle_using_inlist, free_colors_handle
35 : public :: colors_ptr
36 : public :: colors_setup_tables, colors_setup_hooks
37 :
38 : public :: calculate_bolometric, calculate_synthetic
39 : public :: how_many_colors_history_columns, data_for_colors_history_columns
40 : ! old bolometric correction stubs that MESA expects (remove later):
41 : public :: get_bc_id_by_name, get_lum_band_by_id, get_abs_mag_by_id
42 : public :: get_bc_by_id, get_bc_name_by_id, get_bc_by_name
43 : public :: get_abs_bolometric_mag, get_abs_mag_by_name, get_bcs_all
44 : public :: get_lum_band_by_name
45 :
46 : contains
47 :
48 : ! call this routine to initialize the colors module.
49 : ! only needs to be done once at start of run.
50 : ! reads data from the 'colors' directory in the data_dir.
51 : ! if use_cache is true and there is a 'colors/cache' directory, it will try that first.
52 : ! if it doesn't find what it needs in the cache,
53 : ! it reads the data and writes the cache for next time.
54 2 : subroutine colors_init(use_cache, colors_cache_dir, ierr)
55 : use colors_def, only: colors_def_init, colors_use_cache, colors_is_initialized
56 : logical, intent(in) :: use_cache
57 : character(len=*), intent(in) :: colors_cache_dir ! blank means use default
58 : integer, intent(out) :: ierr ! 0 means AOK.
59 2 : ierr = 0
60 2 : if (colors_is_initialized) return
61 2 : call colors_def_init(colors_cache_dir)
62 2 : colors_use_cache = use_cache
63 2 : colors_is_initialized = .true.
64 : end subroutine colors_init
65 :
66 2 : subroutine colors_shutdown
67 : use colors_def, only: do_free_colors_tables, colors_is_initialized
68 2 : call do_free_colors_tables()
69 2 : colors_is_initialized = .false.
70 2 : end subroutine colors_shutdown
71 :
72 : ! after colors_init has finished, you can allocate a "handle".
73 0 : integer function alloc_colors_handle(ierr) result(handle)
74 : integer, intent(out) :: ierr ! 0 means AOK.
75 : character(len=0) :: inlist
76 0 : handle = alloc_colors_handle_using_inlist(inlist, ierr)
77 0 : end function alloc_colors_handle
78 :
79 6 : integer function alloc_colors_handle_using_inlist(inlist, ierr) result(handle)
80 : use colors_def, only: Colors_General_Info, do_alloc_colors, colors_is_initialized, get_colors_ptr
81 : use colors_ctrls_io, only: read_namelist
82 : character(len=*), intent(in) :: inlist ! empty means just use defaults.
83 : integer, intent(out) :: ierr ! 0 means AOK.
84 : type(Colors_General_Info), pointer :: rq
85 : ierr = 0
86 : handle = -1
87 2 : if (.not. colors_is_initialized) then
88 0 : ierr = -1
89 0 : return
90 : end if
91 2 : handle = do_alloc_colors(ierr)
92 2 : if (ierr /= 0) return
93 2 : call read_namelist(handle, inlist, ierr)
94 2 : if (ierr /= 0) return
95 2 : call get_colors_ptr(handle, rq, ierr)
96 2 : if (ierr /= 0) return
97 : ! skip colors table setup and hooks if use_colors = .false.
98 2 : if (.not. rq%use_colors) return
99 0 : call colors_setup_tables(handle, ierr)
100 0 : call colors_setup_hooks(handle, ierr)
101 0 : end function alloc_colors_handle_using_inlist
102 :
103 2 : subroutine free_colors_handle(handle)
104 : ! frees the handle and all associated data
105 : use colors_def, only: colors_General_Info, do_free_colors
106 : integer, intent(in) :: handle
107 2 : call do_free_colors(handle)
108 2 : end subroutine free_colors_handle
109 :
110 4 : subroutine colors_ptr(handle, rq, ierr)
111 :
112 : use colors_def, only: Colors_General_Info, get_colors_ptr, colors_is_initialized
113 :
114 : type(colors_General_Info), pointer, intent(out) :: rq
115 : integer, intent(in) :: handle
116 : integer, intent(out):: ierr
117 :
118 2 : if (.not. colors_is_initialized) then
119 0 : ierr = -1
120 0 : return
121 : end if
122 :
123 2 : call get_colors_ptr(handle, rq, ierr)
124 :
125 : end subroutine colors_ptr
126 :
127 1 : subroutine colors_setup_tables(handle, ierr)
128 : use colors_def, only: Colors_General_Info, get_colors_ptr, color_filter_names, num_color_filters
129 : use synthetic, only: compute_vega_zero_point, compute_ab_zero_point, compute_st_zero_point
130 : integer, intent(in) :: handle
131 : integer, intent(out):: ierr
132 :
133 : type(Colors_General_Info), pointer :: rq
134 : character(len=256) :: lookup_file, filter_dir, filter_filepath, vega_filepath
135 1 : REAL, allocatable :: lookup_table(:, :) ! unused but required by load_lookup_table
136 : integer :: i
137 :
138 : ierr = 0
139 1 : call get_colors_ptr(handle, rq, ierr)
140 1 : if (ierr /= 0) return
141 :
142 1 : call read_strings_from_file(rq, color_filter_names, num_color_filters, ierr)
143 1 : if (ierr /= 0) return
144 :
145 : ! load lookup table (stellar atmosphere grid)
146 1 : if (.not. rq%lookup_loaded) then
147 1 : lookup_file = trim(resolve_path(rq%stellar_atm))//'/lookup_table.csv'
148 : call load_lookup_table(lookup_file, lookup_table, &
149 1 : rq%lu_file_names, rq%lu_logg, rq%lu_meta, rq%lu_teff)
150 1 : rq%lookup_loaded = .true.
151 1 : if (allocated(lookup_table)) deallocate (lookup_table)
152 :
153 : ! build unique sorted grids once for fallback interpolation
154 1 : call build_unique_grids(rq)
155 :
156 : ! build the grid-to-lu mapping for O(1) stencil lookups
157 1 : call build_grid_to_lu_map(rq)
158 : end if
159 :
160 : ! load vega SED
161 1 : if (.not. rq%vega_loaded) then
162 1 : vega_filepath = trim(resolve_path(rq%vega_sed))
163 1 : call load_vega_sed(vega_filepath, rq%vega_wavelengths, rq%vega_fluxes)
164 1 : rq%vega_loaded = .true.
165 : end if
166 :
167 : ! load filter transmission curves and precompute zero-points
168 1 : if (.not. rq%filters_loaded) then
169 1 : filter_dir = trim(resolve_path(rq%instrument))
170 :
171 10 : allocate (rq%filters(num_color_filters))
172 :
173 8 : do i = 1, num_color_filters
174 7 : rq%filters(i)%name = color_filter_names(i)
175 7 : filter_filepath = trim(filter_dir)//'/'//trim(color_filter_names(i))
176 7 : call load_filter(filter_filepath, rq%filters(i)%wavelengths, rq%filters(i)%transmission)
177 :
178 : ! precompute zero-points for all magnitude systems
179 : rq%filters(i)%vega_zero_point = compute_vega_zero_point( &
180 : rq%vega_wavelengths, rq%vega_fluxes, &
181 7 : rq%filters(i)%wavelengths, rq%filters(i)%transmission)
182 :
183 : rq%filters(i)%ab_zero_point = compute_ab_zero_point( &
184 7 : rq%filters(i)%wavelengths, rq%filters(i)%transmission)
185 :
186 : rq%filters(i)%st_zero_point = compute_st_zero_point( &
187 8 : rq%filters(i)%wavelengths, rq%filters(i)%transmission)
188 : end do
189 :
190 1 : rq%filters_loaded = .true.
191 : end if
192 :
193 : ! try to load the full flux cube -- if allocation fails, cube_loaded stays
194 : ! .false. and we fall back to loading individual SED files
195 1 : if (.not. rq%cube_loaded) then
196 1 : call load_flux_cube(rq, rq%stellar_atm)
197 : end if
198 :
199 1 : end subroutine colors_setup_tables
200 :
201 1 : subroutine colors_setup_hooks(handle, ierr)
202 : use colors_def, only: colors_General_Info, get_colors_ptr
203 : integer, intent(in) :: handle
204 : integer, intent(out):: ierr
205 :
206 : type(colors_General_Info), pointer :: rq
207 :
208 : ierr = 0
209 1 : call get_colors_ptr(handle, rq, ierr)
210 :
211 : ! TODO: currently does nothing. See kap if this feature is needed
212 :
213 1 : end subroutine colors_setup_hooks
214 :
215 : ! bolometric correction stubs (legacy mesa interface)
216 :
217 0 : real(dp) function get_bc_by_name(name, log_Teff, log_g, M_div_h, ierr)
218 : character(len=*), intent(in) :: name
219 : real(dp), intent(in) :: log_Teff ! log10 of surface temp
220 : real(dp), intent(in) :: log_g ! log_10 of surface gravity
221 : real(dp), intent(in) :: M_div_h ! [M/H]
222 : integer, intent(inout) :: ierr
223 :
224 0 : get_bc_by_name = -99.9d0
225 0 : ierr = 0
226 0 : end function get_bc_by_name
227 :
228 0 : real(dp) function get_bc_by_id(id, log_Teff, log_g, M_div_h, ierr)
229 : integer, intent(in) :: id
230 : real(dp), intent(in) :: log_Teff ! log10 of surface temp
231 : real(dp), intent(in) :: log_g ! log_10 of surface gravity
232 : real(dp), intent(in) :: M_div_h ! [M/H]
233 : integer, intent(inout) :: ierr
234 :
235 0 : get_bc_by_id = -99.9d0
236 0 : ierr = 0
237 0 : end function get_bc_by_id
238 :
239 0 : integer function get_bc_id_by_name(name, ierr)
240 : character(len=*), intent(in) :: name
241 : integer, intent(inout) :: ierr
242 :
243 0 : get_bc_id_by_name = -1
244 0 : ierr = 0
245 0 : end function get_bc_id_by_name
246 :
247 0 : character(len=strlen) function get_bc_name_by_id(id, ierr)
248 : integer, intent(in) :: id
249 : integer, intent(inout) :: ierr
250 :
251 0 : get_bc_name_by_id = ''
252 0 : ierr = 0
253 0 : end function get_bc_name_by_id
254 :
255 0 : real(dp) function get_abs_bolometric_mag(lum)
256 : use const_def, only: dp
257 : real(dp), intent(in) :: lum ! luminosity in lsun
258 :
259 0 : get_abs_bolometric_mag = -99.9d0
260 0 : end function get_abs_bolometric_mag
261 :
262 0 : real(dp) function get_abs_mag_by_name(name, log_Teff, log_g, M_div_h, lum, ierr)
263 : character(len=*), intent(in) :: name
264 : real(dp), intent(in) :: log_Teff ! log10 of surface temp
265 : real(dp), intent(in) :: M_div_h ! [M/H]
266 : real(dp), intent(in) :: log_g ! log_10 of surface gravity
267 : real(dp), intent(in) :: lum ! luminosity in lsun
268 : integer, intent(inout) :: ierr
269 :
270 0 : ierr = 0
271 0 : get_abs_mag_by_name = -99.9d0
272 0 : end function get_abs_mag_by_name
273 :
274 0 : real(dp) function get_abs_mag_by_id(id, log_Teff, log_g, M_div_h, lum, ierr)
275 : integer, intent(in) :: id
276 : real(dp), intent(in) :: log_Teff ! log10 of surface temp
277 : real(dp), intent(in) :: log_g ! log_10 of surface gravity
278 : real(dp), intent(in) :: M_div_h ! [M/H]
279 : real(dp), intent(in) :: lum ! luminosity in lsun
280 : integer, intent(inout) :: ierr
281 :
282 0 : ierr = 0
283 0 : get_abs_mag_by_id = -99.9d0
284 0 : end function get_abs_mag_by_id
285 :
286 0 : subroutine get_bcs_all(log_Teff, log_g, M_div_h, results, ierr)
287 : real(dp), intent(in) :: log_Teff ! log10 of surface temp
288 : real(dp), intent(in) :: M_div_h ! [M/H]
289 : real(dp), dimension(:), intent(out) :: results
290 : real(dp), intent(in) :: log_g
291 : integer, intent(inout) :: ierr
292 :
293 0 : ierr = 0
294 0 : results(:) = -99.d0
295 0 : end subroutine get_bcs_all
296 :
297 0 : real(dp) function get_lum_band_by_name(name, log_Teff, log_g, M_div_h, lum, ierr)
298 : character(len=*), intent(in) :: name
299 : real(dp), intent(in) :: log_Teff ! log10 of surface temp
300 : real(dp), intent(in) :: M_div_h ! [M/H]
301 : real(dp), intent(in) :: log_g ! log_10 of surface gravity
302 : real(dp), intent(in) :: lum ! luminosity in lsun
303 : integer, intent(inout) :: ierr
304 :
305 0 : ierr = 0
306 0 : get_lum_band_by_name = -99.d0
307 0 : end function get_lum_band_by_name
308 :
309 0 : real(dp) function get_lum_band_by_id(id, log_Teff, log_g, M_div_h, lum, ierr)
310 : integer, intent(in) :: id
311 : real(dp), intent(in) :: log_Teff ! log10 of surface temp
312 : real(dp), intent(in) :: log_g ! log_10 of surface gravity
313 : real(dp), intent(in) :: M_div_h ! [M/H]
314 : real(dp), intent(in) :: lum ! luminosity in lsun
315 : integer, intent(inout) :: ierr
316 :
317 0 : ierr = 0
318 0 : get_lum_band_by_id = -99.d0
319 0 : end function get_lum_band_by_id
320 :
321 : end module colors_lib
|
