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 : ! unit test for the MESA colors module.
21 : !
22 : ! compares synthetic magnitudes (Vega system, Kurucz2003, Johnson filters) and
23 : ! sampled SED flux values against a reference test_output file generated from
24 : ! a known-good state of the code. run via ./ck, which invokes ./rn and
25 : ! compares stdout against test_output using diff -b.
26 : !
27 : ! group 1 - representative stellar types:
28 : ! solar Teff = 5778 K, log g = 4.44, [M/H] = 0.0
29 : ! hot_ms Teff = 15000 K, log g = 4.00, [M/H] = 0.0
30 : ! cool_giant Teff = 4000 K, log g = 2.00, [M/H] = 0.0
31 : !
32 : ! group 2 - grid sweeps (exercises each table dimension independently):
33 : ! vary [M/H]: Teff = 5778, log g = 4.44, [M/H] = -2.0, -1.0, 0.0, 0.5
34 : ! vary log g: Teff = 5778, [M/H] = 0.0, log g = 1.0, 2.5, 4.0, 5.0
35 : ! vary Teff: log g = 4.0, [M/H] = 0.0, Teff = 3500, 6000, 10000, 20000
36 :
37 1 : program test_colors
38 :
39 1 : use const_lib, only: const_init
40 : use math_lib, only: math_init
41 : use colors_lib, only: &
42 : colors_init, colors_shutdown, &
43 : alloc_colors_handle_using_inlist, free_colors_handle, colors_ptr, &
44 : colors_setup_tables, colors_setup_hooks, how_many_colors_history_columns, &
45 : data_for_colors_history_columns, calculate_bolometric
46 : use colors_def, only: Colors_General_Info
47 : use const_def, only: dp, rsun, boltz_sigma
48 : use utils_lib, only: mesa_error
49 : use colors_utils, only: resolve_path
50 :
51 : implicit none
52 :
53 : ! -----------------------------------------------------------------------
54 : ! group 1: representative stellar types
55 : ! -----------------------------------------------------------------------
56 :
57 : integer, parameter :: n_cases = 3
58 :
59 : real(dp), parameter :: test_teff(n_cases) = [5778d0, 15000d0, 4000d0 ]
60 : real(dp), parameter :: test_logg(n_cases) = [4.44d0, 4.0d0, 2.0d0 ]
61 : real(dp), parameter :: test_meta(n_cases) = [0.0d0, 0.0d0, 0.0d0 ]
62 : real(dp), parameter :: test_R(n_cases) = [rsun, 5d0*rsun, 20d0*rsun]
63 :
64 : character(len=12), parameter :: labels(n_cases) = &
65 : ['solar ', 'hot_ms ', 'cool_giant ']
66 :
67 : ! -----------------------------------------------------------------------
68 : ! group 2a: fixed Teff=5778, fixed log g=4.44, varying [M/H]
69 : ! -----------------------------------------------------------------------
70 :
71 : integer, parameter :: n_meta = 4
72 :
73 : real(dp), parameter :: sweep_meta(n_meta) = [-2.0d0, -1.0d0, 0.0d0, 0.5d0]
74 :
75 : ! -----------------------------------------------------------------------
76 : ! group 2b: fixed Teff=5778, fixed [M/H]=0.0, varying log g
77 : ! -----------------------------------------------------------------------
78 :
79 : integer, parameter :: n_logg = 4
80 :
81 : real(dp), parameter :: sweep_logg(n_logg) = [1.0d0, 2.5d0, 4.0d0, 5.0d0]
82 :
83 : ! -----------------------------------------------------------------------
84 : ! group 2c: fixed log g=4.0, fixed [M/H]=0.0, varying Teff
85 : ! -----------------------------------------------------------------------
86 :
87 : integer, parameter :: n_teff = 4
88 :
89 : real(dp), parameter :: sweep_teff(n_teff) = [3500d0, 6000d0, 10000d0, 20000d0]
90 :
91 : ! -----------------------------------------------------------------------
92 : ! shared
93 : ! -----------------------------------------------------------------------
94 :
95 : ! 10 parsecs in cm -> absolute magnitudes
96 : real(dp), parameter :: d_10pc = 3.0857d19
97 :
98 : ! number of sampled SED points printed for the SED comparison
99 : integer, parameter :: n_sed_samples = 20
100 :
101 : ! require the integrated SED flux to stay close to sigma*T^4 scaled to
102 : ! the test distance. large deficits usually mean the wavelength coverage
103 : ! is missing too much UV or IR flux.
104 : real(dp), parameter :: bol_flux_rel_tol = 5d-3
105 :
106 : character(len=32) :: my_mesa_dir
107 : integer :: handle, ierr, n_cols, i, j, k
108 : integer :: model_num
109 : type(Colors_General_Info), pointer :: cs
110 1 : character(len=80), allocatable :: col_names(:)
111 1 : real(dp), allocatable :: col_vals(:)
112 1 : real(dp), allocatable :: wavelengths(:), fluxes(:)
113 : real(dp) :: bol_mag, bol_flux, interp_rad
114 : character(len=256) :: sed_filepath
115 : integer :: n_wav, stride
116 :
117 : ! -----------------------------------------------------------------------
118 : ! module initialization
119 : ! -----------------------------------------------------------------------
120 :
121 1 : my_mesa_dir = '../..'
122 1 : call const_init(my_mesa_dir, ierr)
123 1 : if (ierr /= 0) then
124 0 : write(*,*) 'const_init failed'
125 0 : call mesa_error(__FILE__, __LINE__)
126 : end if
127 :
128 1 : call math_init()
129 :
130 1 : call colors_init(.false., '', ierr)
131 1 : if (ierr /= 0) then
132 0 : write(*,*) 'colors_init failed, ierr =', ierr
133 0 : stop 1
134 : end if
135 :
136 : ! -----------------------------------------------------------------------
137 : ! handle setup: empty inlist string -> defaults (Kurucz2003 + Johnson)
138 : ! -----------------------------------------------------------------------
139 :
140 1 : handle = alloc_colors_handle_using_inlist('', ierr)
141 1 : if (ierr /= 0) then
142 0 : write(*,*) 'alloc_colors_handle_using_inlist failed, ierr =', ierr
143 0 : stop 1
144 : end if
145 :
146 1 : call colors_ptr(handle, cs, ierr)
147 1 : if (ierr /= 0) then
148 0 : write(*,*) 'colors_ptr failed, ierr =', ierr
149 0 : stop 1
150 : end if
151 :
152 : ! enable photometry and explicitly load colors data, since the default
153 : ! namelist leaves use_colors = .false. and skips setup at handle.
154 : ! To do : This could be replaced by reading from an inlist instead in the future.
155 1 : cs%use_colors = .true.
156 1 : cs%mag_system = 'Vega'
157 1 : call colors_setup_tables(handle, ierr)
158 1 : if (ierr /= 0) then
159 0 : write(*,*) 'colors_setup_tables failed, ierr =', ierr
160 0 : stop 1
161 : end if
162 :
163 : ! probably unecessary since this is empty right now.
164 1 : call colors_setup_hooks(handle, ierr)
165 1 : if (ierr /= 0) then
166 0 : write(*,*) 'colors_setup_hooks failed, ierr =', ierr
167 0 : stop 1
168 : end if
169 :
170 1 : n_cols = how_many_colors_history_columns(handle)
171 1 : if (n_cols == 0) then
172 0 : write(*,*) 'how_many_colors_history_columns returned 0'
173 0 : stop 1
174 : end if
175 :
176 5 : allocate(col_names(n_cols), col_vals(n_cols))
177 1 : model_num = 0
178 :
179 : ! -----------------------------------------------------------------------
180 : ! group 1: representative stellar types
181 : ! -----------------------------------------------------------------------
182 :
183 1 : call write_section_header('# Group1 system=Vega grid=Kurucz2003 filters=Johnson')
184 4 : do j = 1, n_cases
185 3 : model_num = model_num + 1
186 : call data_for_colors_history_columns( &
187 : test_teff(j), test_logg(j), test_R(j), test_meta(j), model_num, &
188 3 : handle, n_cols, col_names, col_vals, ierr)
189 3 : if (ierr /= 0) then
190 0 : write(*,*) 'data_for_colors_history_columns failed, group1 case', j, ', ierr =', ierr
191 0 : stop 1
192 : end if
193 3 : write(*,'(a, a)') '# case: ', trim(adjustl(labels(j)))
194 37 : do k = 1, n_cols
195 33 : write(*,'(a40, 1pe23.13)') trim(col_names(k)), col_vals(k)
196 : end do
197 : end do
198 :
199 : ! -----------------------------------------------------------------------
200 : ! group 2a: vary [M/H]
201 : ! -----------------------------------------------------------------------
202 :
203 1 : call write_section_header('# Group2a vary_MH Teff=5778 logg=4.44')
204 5 : do j = 1, n_meta
205 4 : model_num = model_num + 1
206 : call data_for_colors_history_columns( &
207 : 5778d0, 4.44d0, rsun, sweep_meta(j), model_num, &
208 4 : handle, n_cols, col_names, col_vals, ierr)
209 4 : if (ierr /= 0) then
210 0 : write(*,*) 'data_for_colors_history_columns failed, group2a case', j, ', ierr =', ierr
211 0 : stop 1
212 : end if
213 4 : write(*,'(a, f6.2)') '# MH= ', sweep_meta(j)
214 45 : do k = 1, n_cols
215 44 : write(*,'(a40, 1pe23.13)') trim(col_names(k)), col_vals(k)
216 : end do
217 : end do
218 :
219 : ! -----------------------------------------------------------------------
220 : ! group 2b: vary log g
221 : ! -----------------------------------------------------------------------
222 :
223 1 : call write_section_header('# Group2b vary_logg Teff=5778 MH=0.0')
224 5 : do j = 1, n_logg
225 4 : model_num = model_num + 1
226 : call data_for_colors_history_columns( &
227 : 5778d0, sweep_logg(j), rsun, 0.0d0, model_num, &
228 4 : handle, n_cols, col_names, col_vals, ierr)
229 4 : if (ierr /= 0) then
230 0 : write(*,*) 'data_for_colors_history_columns failed, group2b case', j, ', ierr =', ierr
231 0 : stop 1
232 : end if
233 4 : write(*,'(a, f6.2)') '# logg= ', sweep_logg(j)
234 45 : do k = 1, n_cols
235 44 : write(*,'(a40, 1pe23.13)') trim(col_names(k)), col_vals(k)
236 : end do
237 : end do
238 :
239 : ! -----------------------------------------------------------------------
240 : ! group 2c: vary Teff
241 : ! -----------------------------------------------------------------------
242 :
243 1 : call write_section_header('# Group2c vary_Teff logg=4.0 MH=0.0')
244 5 : do j = 1, n_teff
245 4 : model_num = model_num + 1
246 : call data_for_colors_history_columns( &
247 : sweep_teff(j), 4.0d0, rsun, 0.0d0, model_num, &
248 4 : handle, n_cols, col_names, col_vals, ierr)
249 4 : if (ierr /= 0) then
250 0 : write(*,*) 'data_for_colors_history_columns failed, group2c case', j, ', ierr =', ierr
251 0 : stop 1
252 : end if
253 4 : write(*,'(a, f10.1)') '# Teff= ', sweep_teff(j)
254 45 : do k = 1, n_cols
255 44 : write(*,'(a40, 1pe23.13)') trim(col_names(k)), col_vals(k)
256 : end do
257 : end do
258 :
259 1 : sed_filepath = trim(resolve_path(cs%stellar_atm))
260 :
261 : ! -----------------------------------------------------------------------
262 : ! group 3: solar SED sample plus wavelength-coverage sanity checks
263 : ! -----------------------------------------------------------------------
264 :
265 1 : call write_section_header('# Group3 SED sample + wavelength_coverage_sanity')
266 1 : write(*,'(a)') '# SED sample case=solar Teff=5778 logg=4.44 FeH=0.0 columns=wavelength_AA flux_erg_s_cm2_AA'
267 : call calculate_bolometric( &
268 : cs, test_teff(1), test_logg(1), test_meta(1), test_R(1), d_10pc, &
269 1 : bol_mag, bol_flux, wavelengths, fluxes, sed_filepath, interp_rad)
270 :
271 1 : n_wav = size(wavelengths)
272 1 : stride = max(1, n_wav / n_sed_samples)
273 :
274 22 : do i = 1, n_wav, stride
275 22 : write(*,'(1pe23.13, 1x, 1pe23.13)') wavelengths(i), fluxes(i)
276 : end do
277 :
278 1 : write(*,'(a)') ''
279 1 : write(*,'(a, 1pe10.2)') '# wavelength_coverage_sanity logg=4.0 FeH=0.0 rel_tol=', bol_flux_rel_tol
280 5 : do j = 1, n_teff
281 : call check_bolometric_coverage( &
282 5 : cs, sed_filepath, 'Teff', sweep_teff(j), 4.0d0, 0.0d0, rsun, bol_flux_rel_tol)
283 : end do
284 :
285 : ! -----------------------------------------------------------------------
286 : ! cleanup
287 : ! -----------------------------------------------------------------------
288 :
289 1 : deallocate(col_names, col_vals)
290 1 : if (allocated(wavelengths)) deallocate(wavelengths)
291 1 : if (allocated(fluxes)) deallocate(fluxes)
292 :
293 1 : call free_colors_handle(handle)
294 1 : call colors_shutdown()
295 :
296 5 : write(*,*) 'test_colors: passed'
297 :
298 : contains
299 :
300 5 : subroutine write_section_header(title)
301 : character(len=*), intent(in) :: title
302 :
303 5 : write(*,'(a)') ''
304 5 : write(*,'(a)') trim(title)
305 5 : write(*,'(a)') ''
306 5 : end subroutine write_section_header
307 :
308 4 : subroutine check_bolometric_coverage(rq, sed_path, label, teff, log_g, metallicity, radius, rel_tol)
309 : type(Colors_General_Info), intent(inout) :: rq
310 : character(len=*), intent(in) :: sed_path, label
311 : real(dp), intent(in) :: teff, log_g, metallicity, radius, rel_tol
312 :
313 4 : real(dp), allocatable :: local_wavelengths(:), local_fluxes(:)
314 : real(dp) :: local_mag, local_flux, local_interp_rad
315 : real(dp) :: expected_flux, abs_err, rel_err
316 :
317 : call calculate_bolometric( &
318 : rq, teff, log_g, metallicity, radius, d_10pc, &
319 4 : local_mag, local_flux, local_wavelengths, local_fluxes, sed_path, local_interp_rad)
320 :
321 4 : expected_flux = boltz_sigma*teff**4*(radius/d_10pc)**2
322 4 : abs_err = abs(local_flux - expected_flux)
323 4 : rel_err = abs(local_flux - expected_flux)/expected_flux
324 :
325 4 : write(*,'(a, a, a, f10.1)') '# case: ', trim(label), '=', teff
326 4800 : write(*,'(a40, 1pe23.13)') 'Wav_min_AA', minval(local_wavelengths)
327 4800 : write(*,'(a40, 1pe23.13)') 'Wav_max_AA', maxval(local_wavelengths)
328 4 : write(*,'(a40, 1pe23.13)') 'Flux_actual', local_flux
329 4 : write(*,'(a40, 1pe23.13)') 'Flux_expected', expected_flux
330 4 : write(*,'(a40, 1pe23.13)') 'Flux_abserr', abs_err
331 4 : write(*,'(a40, 1pe23.13)') 'Flux_relerr', rel_err
332 :
333 4 : if (rel_err > rel_tol) then
334 : write(*,'(a, a, a, 1pe11.3, a, 1pe11.3)') &
335 0 : 'wavelength coverage sanity check failed for ', trim(label), &
336 0 : ': rel_err=', rel_err, ' > rel_tol=', rel_tol
337 0 : stop 1
338 : end if
339 4 : end subroutine check_bolometric_coverage
340 : end program test_colors
|
