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
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 : ! Accepts cached lookup table data instead of loading from file
38 : !****************************
39 0 : subroutine calculate_bolometric(teff, log_g, metallicity, R, d, bolometric_magnitude, &
40 : bolometric_flux, wavelengths, fluxes, sed_filepath, interpolation_radius, &
41 0 : lu_file_names, lu_teff, lu_logg, lu_meta)
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 : ! Cached lookup table data (passed in from colors_settings)
48 : character(len=100), intent(in) :: lu_file_names(:)
49 : real(dp), intent(in) :: lu_teff(:), lu_logg(:), lu_meta(:)
50 :
51 : character(len=32) :: interpolation_method
52 :
53 0 : interpolation_method = 'Hermite' ! or 'Linear' / 'KNN' later
54 :
55 : ! Quantify how far (teff, log_g, metallicity) is from the grid points
56 : interpolation_radius = compute_interp_radius(teff, log_g, metallicity, &
57 0 : lu_teff, lu_logg, lu_meta)
58 :
59 0 : select case (interpolation_method)
60 : case ('Hermite', 'hermite', 'HERMITE')
61 : call construct_sed_hermite(teff, log_g, metallicity, R, d, lu_file_names, &
62 : lu_teff, lu_logg, lu_meta, sed_filepath, &
63 0 : wavelengths, fluxes)
64 :
65 : case ('Linear', 'linear', 'LINEAR')
66 : call construct_sed_linear(teff, log_g, metallicity, R, d, lu_file_names, &
67 : lu_teff, lu_logg, lu_meta, sed_filepath, &
68 0 : wavelengths, fluxes)
69 :
70 : case ('KNN', 'knn', 'Knn')
71 : call construct_sed_knn(teff, log_g, metallicity, R, d, lu_file_names, &
72 : lu_teff, lu_logg, lu_meta, sed_filepath, &
73 0 : wavelengths, fluxes)
74 :
75 : case default
76 : ! Fallback: Hermite
77 : call construct_sed_hermite(teff, log_g, metallicity, R, d, lu_file_names, &
78 : lu_teff, lu_logg, lu_meta, sed_filepath, &
79 0 : wavelengths, fluxes)
80 : end select
81 :
82 : ! Calculate bolometric flux and magnitude
83 0 : call calculate_bolometric_phot(wavelengths, fluxes, bolometric_magnitude, bolometric_flux)
84 0 : end subroutine calculate_bolometric
85 :
86 : !****************************
87 : ! Calculate Bolometric Magnitude and Flux
88 : !****************************
89 0 : subroutine calculate_bolometric_phot(wavelengths, fluxes, bolometric_magnitude, bolometric_flux)
90 : real(dp), dimension(:), intent(inout) :: wavelengths, fluxes
91 : real(dp), intent(out) :: bolometric_magnitude, bolometric_flux
92 : integer :: i
93 :
94 : ! Validate inputs and replace invalid values with 0
95 0 : do i = 1, size(wavelengths) - 1
96 0 : if (wavelengths(i) <= 0.0d0 .or. fluxes(i) < 0.0d0) then
97 0 : fluxes(i) = 0.0d0
98 : end if
99 : end do
100 :
101 : ! Integrate to get bolometric flux
102 0 : call romberg_integration(wavelengths, fluxes, bolometric_flux)
103 :
104 : ! Validate and calculate magnitude
105 0 : if (bolometric_flux <= 0.0d0) then
106 0 : print *, "Error: Flux integration resulted in non-positive value."
107 0 : bolometric_magnitude = 99.0d0
108 0 : return
109 0 : else if (bolometric_flux < 1.0d-10) then
110 0 : print *, "Warning: Flux value is very small, precision might be affected."
111 : end if
112 :
113 0 : bolometric_magnitude = flux_to_magnitude(bolometric_flux)
114 : end subroutine calculate_bolometric_phot
115 :
116 : !****************************
117 : ! Convert Flux to Magnitude
118 : !****************************
119 0 : real(dp) function flux_to_magnitude(flux)
120 : real(dp), intent(in) :: flux
121 0 : if (flux <= 0.0d0) then
122 0 : print *, "Error: Flux must be positive to calculate magnitude."
123 0 : flux_to_magnitude = 99.0d0
124 : else
125 0 : flux_to_magnitude = -2.5d0 * log10(flux)
126 : end if
127 0 : end function flux_to_magnitude
128 :
129 : !--------------------------------------------------------------------
130 : ! Scalar metric: distance to nearest grid point in normalized space
131 : !--------------------------------------------------------------------
132 0 : real(dp) function compute_interp_radius(teff, log_g, metallicity, &
133 0 : lu_teff, lu_logg, lu_meta)
134 :
135 : real(dp), intent(in) :: teff, log_g, metallicity
136 : real(dp), intent(in) :: lu_teff(:), lu_logg(:), lu_meta(:)
137 :
138 : real(dp) :: teff_min, teff_max, logg_min, logg_max, meta_min, meta_max
139 : real(dp) :: teff_range, logg_range, meta_range
140 : real(dp) :: norm_teff, norm_logg, norm_meta
141 : real(dp) :: grid_teff, grid_logg, grid_meta
142 : real(dp) :: d, d_min
143 : integer :: i, n
144 :
145 : logical :: use_teff, use_logg, use_meta
146 : real(dp), parameter :: eps = 1.0d-12
147 :
148 : ! Detect dummy columns (entire axis is 0 or ±999)
149 : use_teff = .not. (all(lu_teff == 0.0d0) .or. &
150 : all(lu_teff == 999.0d0) .or. &
151 0 : all(lu_teff == -999.0d0))
152 :
153 : use_logg = .not. (all(lu_logg == 0.0d0) .or. &
154 : all(lu_logg == 999.0d0) .or. &
155 0 : all(lu_logg == -999.0d0))
156 :
157 : use_meta = .not. (all(lu_meta == 0.0d0) .or. &
158 : all(lu_meta == 999.0d0) .or. &
159 0 : all(lu_meta == -999.0d0))
160 :
161 : ! Compute min/max for valid axes
162 0 : if (use_teff) then
163 0 : teff_min = minval(lu_teff)
164 0 : teff_max = maxval(lu_teff)
165 0 : teff_range = max(teff_max - teff_min, eps)
166 0 : norm_teff = (teff - teff_min) / teff_range
167 : end if
168 :
169 0 : if (use_logg) then
170 0 : logg_min = minval(lu_logg)
171 0 : logg_max = maxval(lu_logg)
172 0 : logg_range = max(logg_max - logg_min, eps)
173 0 : norm_logg = (log_g - logg_min) / logg_range
174 : end if
175 :
176 0 : if (use_meta) then
177 0 : meta_min = minval(lu_meta)
178 0 : meta_max = maxval(lu_meta)
179 0 : meta_range = max(meta_max - meta_min, eps)
180 0 : norm_meta = (metallicity - meta_min) / meta_range
181 : end if
182 :
183 : ! Find minimum distance to any grid point
184 0 : d_min = huge(1.0d0)
185 0 : n = size(lu_teff)
186 :
187 0 : do i = 1, n
188 0 : d = 0.0d0
189 :
190 0 : if (use_teff) then
191 0 : grid_teff = (lu_teff(i) - teff_min) / teff_range
192 0 : d = d + (norm_teff - grid_teff)**2
193 : end if
194 :
195 0 : if (use_logg) then
196 0 : grid_logg = (lu_logg(i) - logg_min) / logg_range
197 0 : d = d + (norm_logg - grid_logg)**2
198 : end if
199 :
200 0 : if (use_meta) then
201 0 : grid_meta = (lu_meta(i) - meta_min) / meta_range
202 0 : d = d + (norm_meta - grid_meta)**2
203 : end if
204 :
205 0 : d = sqrt(d)
206 0 : if (d < d_min) d_min = d
207 : end do
208 :
209 0 : compute_interp_radius = d_min
210 :
211 0 : end function compute_interp_radius
212 :
213 : end module bolometric
|
