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 bolometric
21 :
22 : use const_def, only: dp
23 : use colors_utils, only: romberg_integration, load_lookup_table
24 : use hermite_interp, only: construct_sed_hermite
25 : use linear_interp, only: construct_sed_linear
26 : use knn_interp, only: construct_sed_knn
27 :
28 : implicit none
29 :
30 : private
31 : public :: calculate_bolometric
32 :
33 : contains
34 :
35 : !****************************
36 : ! Calculate Bolometric Photometry Using Multiple SEDs
37 : !****************************
38 0 : subroutine calculate_bolometric(teff, log_g, metallicity, R, d, bolometric_magnitude, &
39 : bolometric_flux, wavelengths, fluxes, sed_filepath, interpolation_radius)
40 : real(dp), intent(in) :: teff, log_g, metallicity, R, d
41 : character(len=*), intent(in) :: sed_filepath
42 : real(dp), intent(out) :: bolometric_magnitude, bolometric_flux, interpolation_radius
43 : real(dp), dimension(:), allocatable, intent(out) :: wavelengths, fluxes
44 :
45 0 : real(dp), allocatable :: lu_logg(:), lu_meta(:), lu_teff(:)
46 :
47 0 : REAL, dimension(:, :), allocatable :: lookup_table
48 0 : character(len=100), allocatable :: file_names(:)
49 : character(len=256) :: lookup_file
50 : character(len=32) :: interpolation_method
51 :
52 :
53 :
54 0 : lookup_file = trim(sed_filepath)//'/lookup_table.csv'
55 :
56 : call load_lookup_table(lookup_file, lookup_table, file_names, lu_logg, lu_meta, lu_teff)
57 :
58 0 : interpolation_method = 'Hermite' ! or 'Linear' / 'KNN' later
59 :
60 : ! Quantify how far (teff, log_g, metallicity) is from the grid points
61 : interpolation_radius = compute_interp_radius(teff, log_g, metallicity, &
62 0 : lu_teff, lu_logg, lu_meta)
63 :
64 0 : select case (interpolation_method)
65 : case ('Hermite','hermite','HERMITE')
66 : call construct_sed_hermite(teff, log_g, metallicity, R, d, file_names, &
67 : lu_teff, lu_logg, lu_meta, sed_filepath, &
68 0 : wavelengths, fluxes)
69 :
70 : case ('Linear','linear','LINEAR')
71 : call construct_sed_linear(teff, log_g, metallicity, R, d, file_names, &
72 : lu_teff, lu_logg, lu_meta, sed_filepath, &
73 0 : wavelengths, fluxes)
74 :
75 : case ('KNN','knn','Knn')
76 : call construct_sed_knn(teff, log_g, metallicity, R, d, file_names, &
77 : lu_teff, lu_logg, lu_meta, sed_filepath, &
78 0 : wavelengths, fluxes)
79 :
80 : case default
81 : ! Fallback: Hermite
82 : call construct_sed_hermite(teff, log_g, metallicity, R, d, file_names, &
83 : lu_teff, lu_logg, lu_meta, sed_filepath, &
84 0 : wavelengths, fluxes)
85 : end select
86 :
87 : ! Calculate bolometric flux and magnitude
88 0 : call calculate_bolometric_phot(wavelengths, fluxes, bolometric_magnitude, bolometric_flux)
89 0 : end subroutine calculate_bolometric
90 :
91 :
92 : !****************************
93 : ! Calculate Bolometric Magnitude and Flux
94 : !****************************
95 0 : subroutine calculate_bolometric_phot(wavelengths, fluxes, bolometric_magnitude, bolometric_flux)
96 : real(dp), dimension(:), intent(inout) :: wavelengths, fluxes
97 : real(dp), intent(out) :: bolometric_magnitude, bolometric_flux
98 : integer :: i
99 :
100 : ! Validate inputs and replace invalid wavelengths with 0
101 0 : do i = 1, size(wavelengths) - 1
102 0 : if (wavelengths(i) <= 0.0d0 .or. fluxes(i) < 0.0d0) then
103 0 : fluxes(i) = 0.0d0 ! Replace invalid wavelength with 0
104 : end if
105 : end do
106 :
107 : ! Call Romberg integration
108 0 : call romberg_integration(wavelengths, fluxes, bolometric_flux)
109 :
110 : ! Validate integration result
111 0 : if (bolometric_flux <= 0.0d0) then
112 0 : print *, "Error: Flux integration resulted in non-positive value."
113 0 : bolometric_magnitude = 99.0d0
114 0 : return
115 : end if
116 :
117 : ! Calculate bolometric magnitude
118 : if (bolometric_flux <= 0.0d0) then
119 : print *, "Error: Flux integration resulted in non-positive value."
120 : bolometric_magnitude = 99.0d0
121 : return
122 0 : else if (bolometric_flux < 1.0d-10) then
123 0 : print *, "Warning: Flux value is very small, precision might be affected."
124 : end if
125 :
126 0 : bolometric_magnitude = flux_to_magnitude(bolometric_flux)
127 : end subroutine calculate_bolometric_phot
128 :
129 : !****************************
130 : ! Convert Flux to Magnitude
131 : !****************************
132 0 : real(dp) function flux_to_magnitude(flux)
133 : real(dp), intent(in) :: flux
134 0 : if (flux <= 0.0d0) then
135 0 : print *, "Error: Flux must be positive to calculate magnitude."
136 0 : flux_to_magnitude = 99.0d0 ! Return an error value
137 : else
138 0 : flux_to_magnitude = -2.5d0*log10(flux)
139 : end if
140 0 : end function flux_to_magnitude
141 :
142 :
143 :
144 :
145 :
146 :
147 :
148 :
149 : !--------------------------------------------------------------------
150 : ! Scalar metric: distance to nearest grid point in normalized space
151 : !
152 : ! - Uses lu_teff, lu_logg, lu_meta as the available atmosphere grid.
153 : ! - Normalize each dimension to [0,1] using min/max of the grid.
154 : ! - Compute Euclidean distance to the nearest grid point in that
155 : ! normalized space.
156 : !
157 : ! interp_radius ~ 0 => sitting very close to an atmosphere point
158 : ! interp_radius ~ O(1) => deep in-between points / extrapolating
159 : !--------------------------------------------------------------------
160 :
161 0 : real(dp) function compute_interp_radius(teff, log_g, metallicity, &
162 0 : lu_teff, lu_logg, lu_meta)
163 :
164 : real(dp), intent(in) :: teff, log_g, metallicity
165 : real(dp), intent(in) :: lu_teff(:), lu_logg(:), lu_meta(:)
166 :
167 : real(dp) :: teff_min, teff_max, logg_min, logg_max, meta_min, meta_max
168 : real(dp) :: teff_range, logg_range, meta_range
169 : real(dp) :: norm_teff, norm_logg, norm_meta
170 : real(dp) :: grid_teff, grid_logg, grid_meta
171 : real(dp) :: d, d_min
172 : integer :: i, n
173 :
174 : logical :: use_teff, use_logg, use_meta
175 : real(dp), parameter :: eps = 1.0d-12
176 :
177 : ! ---------------------------------------------------------
178 : ! Detect dummy columns (entire axis is 0 or 999 or -999)
179 : ! ---------------------------------------------------------
180 : use_teff = .not. ( all(lu_teff == 0.0d0) .or. &
181 : all(lu_teff == 999.0d0) .or. &
182 0 : all(lu_teff == -999.0d0) )
183 :
184 : use_logg = .not. ( all(lu_logg == 0.0d0) .or. &
185 : all(lu_logg == 999.0d0) .or. &
186 0 : all(lu_logg == -999.0d0) )
187 :
188 : use_meta = .not. ( all(lu_meta == 0.0d0) .or. &
189 : all(lu_meta == 999.0d0) .or. &
190 0 : all(lu_meta == -999.0d0) )
191 :
192 : ! ---------------------------------------------------------
193 : ! Compute min/max only for VALID axes
194 : ! ---------------------------------------------------------
195 :
196 0 : if (use_teff) then
197 0 : teff_min = minval(lu_teff)
198 0 : teff_max = maxval(lu_teff)
199 0 : teff_range = max(teff_max - teff_min, eps)
200 0 : norm_teff = (teff - teff_min)/teff_range
201 : endif
202 :
203 0 : if (use_logg) then
204 0 : logg_min = minval(lu_logg)
205 0 : logg_max = maxval(lu_logg)
206 0 : logg_range = max(logg_max - logg_min, eps)
207 0 : norm_logg = (log_g - logg_min)/logg_range
208 : endif
209 :
210 0 : if (use_meta) then
211 0 : meta_min = minval(lu_meta)
212 0 : meta_max = maxval(lu_meta)
213 0 : meta_range = max(meta_max - meta_min, eps)
214 0 : norm_meta = (metallicity - meta_min)/meta_range
215 : endif
216 :
217 : ! ---------------------------------------------------------
218 : ! Compute minimum distance with dimension-dropping
219 : ! ---------------------------------------------------------
220 0 : d_min = huge(1.0d0)
221 0 : n = size(lu_teff)
222 :
223 0 : do i = 1, n
224 :
225 0 : d = 0.0d0
226 :
227 0 : if (use_teff) then
228 0 : grid_teff = (lu_teff(i) - teff_min)/teff_range
229 0 : d = d + (norm_teff - grid_teff)**2
230 : end if
231 :
232 0 : if (use_logg) then
233 0 : grid_logg = (lu_logg(i) - logg_min)/logg_range
234 0 : d = d + (norm_logg - grid_logg)**2
235 : end if
236 :
237 0 : if (use_meta) then
238 0 : grid_meta = (lu_meta(i) - meta_min)/meta_range
239 0 : d = d + (norm_meta - grid_meta)**2
240 : end if
241 :
242 0 : d = sqrt(d)
243 :
244 0 : if (d < d_min) d_min = d
245 : end do
246 :
247 0 : compute_interp_radius = d_min
248 :
249 0 : end function compute_interp_radius
250 :
251 :
252 :
253 :
254 :
255 : end module bolometric
|
