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