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