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 : ! Regression test for the MESA colors module.
21 : !
22 : ! Like the other module unit tests (see e.g. neu/test, interp_2d/test),
23 : ! this program only prints computed values -- it does not assert against
24 : ! hardcoded expected values itself. Correctness is verified externally by
25 : ! `make check`, which diffs the printed output against test_output.
26 :
27 1 : program test_colors
28 :
29 1 : use const_lib, only: const_init
30 : use math_lib, only: math_init, pow2
31 : use colors_lib, only: &
32 : colors_init, colors_shutdown, &
33 : alloc_colors_handle_using_inlist, free_colors_handle, colors_ptr, &
34 : colors_setup_tables, colors_setup_hooks, how_many_colors_history_columns, &
35 : data_for_colors_history_columns
36 : use colors_def, only: Colors_General_Info
37 : use const_def, only: dp, pi, pc, rsun, Lsun
38 : use utils_lib, only: mesa_error
39 : use synthetic, only: calculate_synthetic, compute_vega_zero_point
40 : use hermite_interp, only: hermite_tensor_interp3d
41 : use hermite_interp_bounded, only: construct_sed_hermite_bounded
42 : use bolometric, only: calculate_bolometric_phot
43 :
44 : implicit none
45 :
46 : integer, parameter :: n_cases = 4
47 :
48 : real(dp), parameter :: test_teff(n_cases) = &
49 : [5778d0, 15000d0, 4000d0, 5778d0]
50 : real(dp), parameter :: test_logg(n_cases) = &
51 : [4.44d0, 4.0d0, 2.0d0, 4.44d0]
52 : real(dp), parameter :: test_meta(n_cases) = &
53 : [0.0d0, 0.0d0, 0.0d0, -2.0d0]
54 : real(dp), parameter :: test_R(n_cases) = &
55 : [rsun, 5d0*rsun, 20d0*rsun, rsun]
56 :
57 : character(len=12), parameter :: labels(n_cases) = &
58 : ['solar ', 'hot_ms ', 'cool_giant ', 'metal_poor ']
59 :
60 : character(len=32) :: my_mesa_dir
61 : integer :: handle, ierr, n_cols, j, k
62 : integer :: model_num
63 : type(Colors_General_Info), pointer :: cs
64 1 : character(len=80), allocatable :: col_names(:)
65 1 : real(dp), allocatable :: col_vals(:)
66 :
67 1 : my_mesa_dir = '../..'
68 1 : call const_init(my_mesa_dir, ierr)
69 1 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
70 :
71 1 : call math_init()
72 :
73 1 : call colors_init(.false., '', ierr)
74 1 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
75 :
76 1 : handle = alloc_colors_handle_using_inlist('', ierr)
77 1 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
78 :
79 1 : call colors_ptr(handle, cs, ierr)
80 1 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
81 :
82 1 : cs%use_colors = .true.
83 1 : cs%mag_system = 'Vega'
84 :
85 1 : call colors_setup_tables(handle, ierr)
86 1 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
87 :
88 1 : call colors_setup_hooks(handle, ierr)
89 1 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
90 :
91 1 : n_cols = how_many_colors_history_columns(handle)
92 1 : if (n_cols <= 0) call mesa_error(__FILE__, __LINE__)
93 :
94 5 : allocate (col_names(n_cols), col_vals(n_cols))
95 1 : model_num = 0
96 :
97 5 : do j = 1, n_cases
98 4 : model_num = model_num + 1
99 : call data_for_colors_history_columns( &
100 : test_teff(j), test_logg(j), test_R(j), test_meta(j), model_num, &
101 4 : handle, n_cols, col_names, col_vals, ierr)
102 4 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
103 :
104 : ! col_names is only populated by the call above, so the header can
105 : ! only be printed once the first case has run.
106 4 : if (j == 1) then
107 1 : write (*, '(a)', advance='no') 'columns: Teff logg MH R_Rsun'
108 11 : do k = 1, n_cols
109 11 : write (*, '(1x,a)', advance='no') trim(col_names(k))
110 : end do
111 1 : write (*, '(a)') ''
112 : end if
113 :
114 : write (*, '(a,":",1x,4(1pe16.8,1x))', advance='no') &
115 4 : trim(labels(j)), test_teff(j), test_logg(j), test_meta(j), &
116 8 : test_R(j)/rsun
117 44 : do k = 1, n_cols
118 44 : write (*, '(1pe16.8,1x)', advance='no') col_vals(k)
119 : end do
120 9 : write (*, '(a)') ''
121 : end do
122 :
123 1 : call check_filter_support()
124 1 : call check_nonuniform_hermite_derivatives()
125 1 : call check_bounded_hermite_small_negative()
126 1 : call check_bounded_hermite_fallback()
127 1 : call check_bolometric_zero_point()
128 :
129 1 : call free_colors_handle(handle)
130 5 : call colors_shutdown()
131 :
132 : contains
133 :
134 1 : subroutine check_filter_support()
135 : real(dp) :: sed_wave(5), filter_wave(3), sed_flux(5), transmission(3)
136 : real(dp) :: magnitude, zero_point
137 : integer :: local_ierr
138 :
139 : ! Filter is tabulated only over the middle three wavelengths; the
140 : ! bright outer two samples must not leak into either result.
141 1 : sed_wave = [1.0_dp, 2.0_dp, 3.0_dp, 4.0_dp, 5.0_dp]
142 1 : filter_wave = [2.0_dp, 3.0_dp, 4.0_dp]
143 1 : sed_flux = [100.0_dp, 2.0_dp, 2.0_dp, 2.0_dp, 100.0_dp]
144 4 : transmission = 1.0_dp
145 :
146 : magnitude = calculate_synthetic( &
147 : 0.0_dp, 0.0_dp, 0.0_dp, local_ierr, sed_wave, sed_flux, &
148 : filter_wave, transmission, 1.0_dp, 'test.dat', .false., &
149 1 : .false., '.', 0)
150 1 : if (local_ierr /= 0) call mesa_error(__FILE__, __LINE__)
151 :
152 : zero_point = compute_vega_zero_point( &
153 1 : sed_wave, sed_flux, filter_wave, transmission)
154 :
155 : write (*, '(a,":",1x,2(1pe23.13,1x))') &
156 1 : 'filter_support', magnitude, zero_point
157 1 : end subroutine check_filter_support
158 :
159 1 : subroutine check_nonuniform_hermite_derivatives()
160 1 : type(Colors_General_Info) :: bounded_cs
161 : real(dp) :: x_grid(4), singleton(1), values(4, 1, 1), actual
162 1 : real(dp), allocatable :: actual_wave(:), actual_flux(:)
163 : integer :: lam
164 :
165 1 : x_grid = [0.0_dp, 1.0_dp, 3.0_dp, 6.0_dp]
166 1 : singleton = [0.0_dp]
167 5 : values(:, 1, 1) = 100.0_dp + x_grid*x_grid
168 :
169 : actual = hermite_tensor_interp3d( &
170 1 : 1.5_dp, 0.0_dp, 0.0_dp, x_grid, singleton, singleton, values)
171 1 : write (*, '(a,":",1x,1pe23.13)') 'nonuniform_hermite_3d', actual
172 :
173 1 : bounded_cs%cube_loaded = .true.
174 5 : bounded_cs%cube_teff_grid = x_grid
175 2 : bounded_cs%cube_logg_grid = singleton
176 2 : bounded_cs%cube_meta_grid = singleton
177 4 : bounded_cs%cube_wavelengths = [1.0_dp, 2.0_dp, 4.0_dp]
178 1 : allocate (bounded_cs%cube_flux(4, 1, 1, 3))
179 4 : do lam = 1, 3
180 16 : bounded_cs%cube_flux(:, 1, 1, lam) = values(:, 1, 1)
181 : end do
182 :
183 : call construct_sed_hermite_bounded( &
184 : bounded_cs, 1.5_dp, 0.0_dp, 0.0_dp, 1.0_dp, 1.0_dp, '', &
185 1 : actual_wave, actual_flux)
186 : write (*, '(a,":",1x,*(1pe23.13,1x))') &
187 1 : 'nonuniform_hermite_4d', actual_flux
188 1 : end subroutine check_nonuniform_hermite_derivatives
189 :
190 1 : subroutine check_bounded_hermite_small_negative()
191 1 : type(Colors_General_Info) :: bounded_cs
192 1 : real(dp), allocatable :: actual_wave(:), actual_flux(:)
193 :
194 1 : bounded_cs%cube_loaded = .true.
195 4 : bounded_cs%cube_teff_grid = [0.0_dp, 1.0_dp, 2.0_dp]
196 2 : bounded_cs%cube_logg_grid = [0.0_dp]
197 2 : bounded_cs%cube_meta_grid = [0.0_dp]
198 4 : bounded_cs%cube_wavelengths = [1.0_dp, 2.0_dp, 4.0_dp]
199 1 : allocate (bounded_cs%cube_flux(3, 1, 1, 3))
200 4 : bounded_cs%cube_flux(:, 1, 1, 1) = 100.0_dp
201 4 : bounded_cs%cube_flux(:, 1, 1, 2) = [1.0_dp, 1.0_dp, 17.000016_dp]
202 4 : bounded_cs%cube_flux(:, 1, 1, 3) = 100.0_dp
203 :
204 : call construct_sed_hermite_bounded( &
205 : bounded_cs, 0.5_dp, 0.0_dp, 0.0_dp, 1.0_dp, 1.0_dp, '', &
206 1 : actual_wave, actual_flux)
207 : write (*, '(a,":",1x,*(1pe23.13,1x))') &
208 1 : 'bounded_hermite_small_negative', actual_flux
209 1 : end subroutine check_bounded_hermite_small_negative
210 :
211 1 : subroutine check_bounded_hermite_fallback()
212 1 : type(Colors_General_Info) :: bounded_cs
213 1 : real(dp), allocatable :: actual_wave(:), actual_flux(:)
214 : integer :: lam
215 :
216 1 : bounded_cs%cube_loaded = .true.
217 4 : bounded_cs%cube_teff_grid = [0.0_dp, 1.0_dp, 2.0_dp]
218 2 : bounded_cs%cube_logg_grid = [0.0_dp]
219 2 : bounded_cs%cube_meta_grid = [0.0_dp]
220 4 : bounded_cs%cube_wavelengths = [1.0_dp, 2.0_dp, 4.0_dp]
221 1 : allocate (bounded_cs%cube_flux(3, 1, 1, 3))
222 4 : do lam = 1, 3
223 13 : bounded_cs%cube_flux(:, 1, 1, lam) = [1.0_dp, 1.0_dp, 100.0_dp]
224 : end do
225 :
226 : call construct_sed_hermite_bounded( &
227 : bounded_cs, 0.5_dp, 0.0_dp, 0.0_dp, 1.0_dp, 1.0_dp, '', &
228 1 : actual_wave, actual_flux)
229 : write (*, '(a,":",1x,*(1pe23.13,1x))') &
230 1 : 'bounded_hermite_fallback', actual_flux
231 1 : end subroutine check_bounded_hermite_fallback
232 :
233 1 : subroutine check_bolometric_zero_point()
234 : real(dp) :: local_wavelengths(3), local_fluxes(3)
235 : real(dp) :: solar_flux, local_mag, local_bol_flux
236 :
237 1 : solar_flux = Lsun/(4.0_dp*pi*pow2(10.0_dp*pc))
238 1 : local_wavelengths = [1.0_dp, 2.0_dp, 3.0_dp]
239 4 : local_fluxes = solar_flux/2.0_dp
240 :
241 : call calculate_bolometric_phot( &
242 1 : local_wavelengths, local_fluxes, local_mag, local_bol_flux)
243 : write (*, '(a,":",1x,1pe23.13)') &
244 1 : 'bolometric_solar_zero_point', local_mag
245 1 : end subroutine check_bolometric_zero_point
246 :
247 : end program test_colors
|
