Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2026 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 bolometric
21 :
22 : use const_def, only: dp, f_bol_zero_point
23 : use colors_def, only: Colors_General_Info
24 : use colors_utils, only: simpson_integration
25 : use hermite_interp, only: construct_sed_hermite
26 : use hermite_interp_bounded, only: construct_sed_hermite_bounded
27 : use linear_interp, only: construct_sed_linear
28 : use knn_interp, only: construct_sed_knn
29 : use utils_lib, only: is_inf, is_nan, mesa_error
30 :
31 : implicit none
32 :
33 : private
34 : public :: calculate_bolometric, calculate_bolometric_phot
35 :
36 : contains
37 :
38 : ! rq carries cached lookup table data and the preloaded flux cube (if available)
39 8 : subroutine calculate_bolometric(rq, teff, log_g, metallicity, R, d, bolometric_magnitude, &
40 : bolometric_flux, wavelengths, fluxes, sed_filepath, interpolation_radius)
41 : type(Colors_General_Info), intent(inout) :: rq
42 : real(dp), intent(in) :: teff, log_g, metallicity, R, d
43 : character(len=*), intent(in) :: sed_filepath
44 : real(dp), intent(out) :: bolometric_magnitude, bolometric_flux, interpolation_radius
45 : real(dp), dimension(:), allocatable, intent(out) :: wavelengths, fluxes
46 :
47 : character(len=32) :: interpolation_method
48 :
49 4 : interpolation_method = 'Hermite_bounded' ! or 'Linear' / 'KNN' / 'Hermite' / 'Hermite_bounded'
50 :
51 : ! how far (teff, log_g, metallicity) is from the nearest grid point
52 : interpolation_radius = compute_interp_radius(teff, log_g, metallicity, &
53 4 : rq%lu_teff, rq%lu_logg, rq%lu_meta)
54 :
55 0 : select case (interpolation_method)
56 : case ('Hermite', 'hermite', 'HERMITE')
57 : call construct_sed_hermite(rq, teff, log_g, metallicity, R, d, &
58 0 : sed_filepath, wavelengths, fluxes)
59 :
60 : case ('Linear', 'linear', 'LINEAR')
61 : call construct_sed_linear(rq, teff, log_g, metallicity, R, d, &
62 0 : sed_filepath, wavelengths, fluxes)
63 :
64 : case ('KNN', 'knn', 'Knn')
65 : call construct_sed_knn(rq, teff, log_g, metallicity, R, d, &
66 0 : sed_filepath, wavelengths, fluxes)
67 :
68 : case ('Hermite_bounded', 'hermite_bounded', 'HERMITE_BOUNDED')
69 : call construct_sed_hermite_bounded(rq, teff, log_g, metallicity, R, d, &
70 4 : sed_filepath, wavelengths, fluxes)
71 :
72 : case default
73 : ! fallback: bounded hermite (switches to linear in regions where hermite overshoots.)
74 : call construct_sed_hermite_bounded(rq, teff, log_g, metallicity, R, d, &
75 4 : sed_filepath, wavelengths, fluxes)
76 : end select
77 :
78 4 : call calculate_bolometric_phot(wavelengths, fluxes, bolometric_magnitude, bolometric_flux)
79 4 : end subroutine calculate_bolometric
80 :
81 5 : subroutine calculate_bolometric_phot(wavelengths, fluxes, bolometric_magnitude, bolometric_flux)
82 : real(dp), dimension(:), intent(in) :: wavelengths, fluxes
83 : real(dp), intent(out) :: bolometric_magnitude, bolometric_flux
84 5 : real(dp), allocatable :: clean_fluxes(:)
85 : integer :: i
86 :
87 5 : if (size(wavelengths) /= size(fluxes) .or. size(wavelengths) < 2) then
88 0 : write (*, *) 'colors: invalid array sizes in calculate_bolometric_phot'
89 0 : call mesa_error(__FILE__, __LINE__)
90 : end if
91 :
92 15 : allocate (clean_fluxes(size(fluxes)))
93 4809 : clean_fluxes = fluxes
94 :
95 4804 : do i = 1, size(wavelengths)
96 : if (is_nan(wavelengths(i)) .or. is_inf(wavelengths(i)) .or. wavelengths(i) <= 0.0d0 .or. &
97 4804 : is_nan(clean_fluxes(i)) .or. is_inf(clean_fluxes(i)) .or. clean_fluxes(i) < 0.0d0) then
98 0 : clean_fluxes(i) = 0.0d0
99 : end if
100 : end do
101 :
102 5 : call simpson_integration(wavelengths, clean_fluxes, bolometric_flux)
103 :
104 5 : if (bolometric_flux <= 0.0d0) then
105 0 : print *, "Error: Flux integration resulted in non-positive value."
106 0 : bolometric_magnitude = 99.0d0
107 : return
108 5 : else if (bolometric_flux < 1.0d-10) then
109 0 : print *, "Warning: Flux value is very small, precision might be affected."
110 : end if
111 :
112 5 : bolometric_magnitude = flux_to_magnitude(bolometric_flux)
113 5 : end subroutine calculate_bolometric_phot
114 :
115 5 : real(dp) function flux_to_magnitude(flux)
116 : real(dp), intent(in) :: flux ! must be positive; checked by the caller
117 5 : flux_to_magnitude = -2.5d0*log10(flux/f_bol_zero_point)
118 5 : end function flux_to_magnitude
119 :
120 : ! scalar metric: distance to nearest grid point in normalized space
121 4 : real(dp) function compute_interp_radius(teff, log_g, metallicity, &
122 4 : lu_teff, lu_logg, lu_meta)
123 :
124 : real(dp), intent(in) :: teff, log_g, metallicity
125 : real(dp), intent(in) :: lu_teff(:), lu_logg(:), lu_meta(:)
126 :
127 : real(dp) :: teff_min, teff_max, logg_min, logg_max, meta_min, meta_max
128 : real(dp) :: teff_range, logg_range, meta_range
129 : real(dp) :: norm_teff, norm_logg, norm_meta
130 : real(dp) :: grid_teff, grid_logg, grid_meta
131 : real(dp) :: d, d_min
132 : integer :: i, n
133 :
134 : logical :: use_teff, use_logg, use_meta
135 : real(dp), parameter :: eps = 1.0d-12
136 :
137 : ! detect dummy columns (entire axis is 0 or ±999)
138 : use_teff = .not. (all(lu_teff == 0.0d0) .or. &
139 : all(lu_teff == 999.0d0) .or. &
140 12 : all(lu_teff == -999.0d0))
141 :
142 : use_logg = .not. (all(lu_logg == 0.0d0) .or. &
143 : all(lu_logg == 999.0d0) .or. &
144 12 : all(lu_logg == -999.0d0))
145 :
146 : use_meta = .not. (all(lu_meta == 0.0d0) .or. &
147 : all(lu_meta == 999.0d0) .or. &
148 12 : all(lu_meta == -999.0d0))
149 :
150 4 : if (use_teff) then
151 15236 : teff_min = minval(lu_teff)
152 15236 : teff_max = maxval(lu_teff)
153 4 : teff_range = max(teff_max - teff_min, eps)
154 4 : norm_teff = (teff - teff_min)/teff_range
155 : end if
156 :
157 4 : if (use_logg) then
158 15236 : logg_min = minval(lu_logg)
159 15236 : logg_max = maxval(lu_logg)
160 4 : logg_range = max(logg_max - logg_min, eps)
161 4 : norm_logg = (log_g - logg_min)/logg_range
162 : end if
163 :
164 4 : if (use_meta) then
165 15236 : meta_min = minval(lu_meta)
166 15236 : meta_max = maxval(lu_meta)
167 4 : meta_range = max(meta_max - meta_min, eps)
168 4 : norm_meta = (metallicity - meta_min)/meta_range
169 : end if
170 :
171 : ! find minimum distance to any grid point
172 4 : d_min = huge(1.0d0)
173 4 : n = size(lu_teff)
174 :
175 15236 : do i = 1, n
176 15232 : d = 0.0d0
177 :
178 15232 : if (use_teff) then
179 15232 : grid_teff = (lu_teff(i) - teff_min)/teff_range
180 15232 : d = d + (norm_teff - grid_teff)**2
181 : end if
182 :
183 15232 : if (use_logg) then
184 15232 : grid_logg = (lu_logg(i) - logg_min)/logg_range
185 15232 : d = d + (norm_logg - grid_logg)**2
186 : end if
187 :
188 15232 : if (use_meta) then
189 15232 : grid_meta = (lu_meta(i) - meta_min)/meta_range
190 15232 : d = d + (norm_meta - grid_meta)**2
191 : end if
192 :
193 15232 : d = sqrt(d)
194 4 : if (d < d_min) d_min = d
195 : end do
196 :
197 4 : compute_interp_radius = d_min
198 :
199 4 : end function compute_interp_radius
200 :
201 : end module bolometric
|
