Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2025 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 synthetic
21 :
22 : use const_def, only: dp, clight
23 : use utils_lib, only: mkdir, folder_exists
24 : use colors_utils, only: remove_dat, simpson_integration
25 : use knn_interp, only: interpolate_array
26 :
27 : implicit none
28 :
29 : private
30 : public :: calculate_synthetic
31 : public :: compute_vega_zero_point, compute_ab_zero_point, compute_st_zero_point
32 : public :: zero_filter_outside_support
33 :
34 : contains
35 :
36 : ! calculate synthetic magnitude; uses precomputed zero-point from filter_data
37 29 : real(dp) function calculate_synthetic(temperature, gravity, metallicity, ierr, &
38 29 : wavelengths, fluxes, &
39 29 : filter_wavelengths, filter_trans, &
40 : zero_point_flux, &
41 : filter_name, make_sed, sed_per_model, colors_results_directory, model_number)
42 : real(dp), intent(in) :: temperature, gravity, metallicity
43 : character(len=*), intent(in) :: filter_name, colors_results_directory
44 : integer, intent(out) :: ierr
45 :
46 : real(dp), dimension(:), intent(in) :: wavelengths, fluxes
47 : real(dp), dimension(:), intent(in) :: filter_wavelengths, filter_trans
48 : real(dp), intent(in) :: zero_point_flux ! precomputed at initialization
49 : logical, intent(in) :: make_sed, sed_per_model
50 : integer, intent(in) :: model_number
51 :
52 29 : real(dp), dimension(:), allocatable :: convolved_flux, filter_on_sed_grid
53 : character(len=256) :: csv_file
54 : character(len=20) :: model_str
55 : character(len=1000) :: line
56 : real(dp) :: synthetic_flux
57 : integer :: max_size, i
58 : real(dp) :: wv, fl, cf, fwv, ftr
59 :
60 29 : ierr = 0
61 :
62 87 : allocate (convolved_flux(size(wavelengths)))
63 58 : allocate (filter_on_sed_grid(size(wavelengths)))
64 :
65 29 : call interpolate_array(filter_wavelengths, filter_trans, wavelengths, filter_on_sed_grid)
66 29 : call zero_filter_outside_support(wavelengths, filter_wavelengths, filter_on_sed_grid)
67 :
68 33635 : convolved_flux = fluxes*filter_on_sed_grid
69 :
70 29 : if (make_sed) then
71 0 : if (.not. folder_exists(trim(colors_results_directory))) call mkdir(trim(colors_results_directory))
72 :
73 : ! track model number internally when write_sed_per_model is enabled
74 0 : if (sed_per_model) then
75 :
76 0 : write (model_str, '(I8.8)') model_number
77 0 : csv_file = trim(colors_results_directory)//'/'//trim(remove_dat(filter_name))//'_SED_'//trim(model_str)//'.csv'
78 :
79 : else
80 0 : csv_file = trim(colors_results_directory)//'/'//trim(remove_dat(filter_name))//'_SED.csv'
81 : end if
82 :
83 0 : max_size = max(size(wavelengths), size(filter_wavelengths))
84 :
85 0 : open (unit=10, file=csv_file, status='REPLACE', action='write', iostat=ierr)
86 0 : if (ierr /= 0) then
87 0 : print *, "Error opening file for writing"
88 0 : deallocate (convolved_flux, filter_on_sed_grid)
89 0 : calculate_synthetic = huge(1.0_dp)
90 : return
91 : end if
92 :
93 0 : write (10, '(A)') "wavelengths,fluxes,convolved_flux,filter_wavelengths,filter_trans"
94 :
95 0 : do i = 1, max_size
96 0 : wv = 0.0_dp; fl = 0.0_dp; cf = 0.0_dp; fwv = 0.0_dp; ftr = 0.0_dp
97 0 : if (i <= size(wavelengths)) wv = wavelengths(i)
98 0 : if (i <= size(fluxes)) fl = fluxes(i)
99 0 : if (i <= size(convolved_flux)) cf = convolved_flux(i)
100 0 : if (i <= size(filter_wavelengths)) fwv = filter_wavelengths(i)
101 0 : if (i <= size(filter_trans)) ftr = filter_trans(i)
102 :
103 : write (line, '(ES14.6, ",", ES14.6, ",", ES14.6, ",", ES14.6, ",", ES14.6)') &
104 0 : wv, fl, cf, fwv, ftr
105 0 : write (10, '(A)') trim(line)
106 : end do
107 0 : close (10)
108 : end if
109 :
110 29 : call calculate_synthetic_flux(wavelengths, convolved_flux, filter_on_sed_grid, synthetic_flux)
111 :
112 29 : if (zero_point_flux > 0.0_dp .and. synthetic_flux > 0.0_dp) then
113 29 : calculate_synthetic = -2.5d0*log10(synthetic_flux/zero_point_flux)
114 : else
115 0 : if (zero_point_flux <= 0.0_dp) then
116 0 : print *, "Error: Zero point flux is zero or negative for filter ", trim(filter_name)
117 : end if
118 0 : if (synthetic_flux <= 0.0_dp) then
119 0 : print *, "Error: Synthetic flux is zero or negative for filter ", trim(filter_name)
120 : end if
121 : calculate_synthetic = huge(1.0_dp)
122 : end if
123 :
124 29 : deallocate (convolved_flux, filter_on_sed_grid)
125 29 : end function calculate_synthetic
126 :
127 : ! Ensure an interpolated filter response has compact support.
128 : ! Some interpolation backends extrapolate beyond the tabulated filter range.
129 : ! For photometric passbands, transmission outside the supplied passband must be zero.
130 37 : subroutine zero_filter_outside_support(target_wave, filt_wave, filt_on_target)
131 : real(dp), dimension(:), intent(in) :: target_wave, filt_wave
132 : real(dp), dimension(:), intent(inout) :: filt_on_target
133 :
134 : real(dp) :: fmin, fmax
135 : integer :: i
136 :
137 753 : fmin = minval(filt_wave)
138 753 : fmax = maxval(filt_wave)
139 :
140 90277 : do i = 1, size(target_wave)
141 90277 : if (target_wave(i) < fmin .or. target_wave(i) > fmax) then
142 83086 : filt_on_target(i) = 0.0_dp
143 7154 : else if (filt_on_target(i) < 0.0_dp) then
144 0 : filt_on_target(i) = 0.0_dp
145 : end if
146 : end do
147 37 : end subroutine zero_filter_outside_support
148 :
149 : ! photon-counting synthetic flux (numerator and denominator weighted by lambda)
150 29 : subroutine calculate_synthetic_flux(wavelengths, convolved_flux, filter_on_sed_grid, synthetic_flux)
151 : real(dp), dimension(:), intent(in) :: wavelengths, convolved_flux, filter_on_sed_grid
152 : real(dp), intent(out) :: synthetic_flux
153 :
154 : real(dp) :: integrated_flux, integrated_filter
155 :
156 : ! photon-counting: weight by wavelength
157 33606 : call simpson_integration(wavelengths, convolved_flux*wavelengths, integrated_flux)
158 33606 : call simpson_integration(wavelengths, filter_on_sed_grid*wavelengths, integrated_filter)
159 :
160 29 : if (integrated_filter > 0.0_dp) then
161 29 : synthetic_flux = integrated_flux/integrated_filter
162 : else
163 0 : print *, "Error: Integrated filter transmission is zero."
164 0 : synthetic_flux = -1.0_dp
165 : end if
166 29 : end subroutine calculate_synthetic_flux
167 :
168 : ! vega zero-point -- called once at init, result cached in filter_data
169 8 : real(dp) function compute_vega_zero_point(vega_wave, vega_flux, filt_wave, filt_trans)
170 : real(dp), dimension(:), intent(in) :: vega_wave, vega_flux
171 : real(dp), dimension(:), intent(in) :: filt_wave, filt_trans
172 :
173 : real(dp) :: int_flux, int_filter
174 8 : real(dp), allocatable :: filt_on_vega_grid(:), conv_flux(:)
175 :
176 24 : allocate (filt_on_vega_grid(size(vega_wave)))
177 16 : allocate (conv_flux(size(vega_wave)))
178 :
179 8 : call interpolate_array(filt_wave, filt_trans, vega_wave, filt_on_vega_grid)
180 8 : call zero_filter_outside_support(vega_wave, filt_wave, filt_on_vega_grid)
181 :
182 56679 : conv_flux = vega_flux*filt_on_vega_grid
183 :
184 : ! photon-counting integration
185 56671 : call simpson_integration(vega_wave, vega_wave*conv_flux, int_flux)
186 56671 : call simpson_integration(vega_wave, vega_wave*filt_on_vega_grid, int_filter)
187 :
188 8 : if (int_filter > 0.0_dp) then
189 8 : compute_vega_zero_point = int_flux/int_filter
190 : else
191 : compute_vega_zero_point = -1.0_dp
192 : end if
193 :
194 8 : deallocate (filt_on_vega_grid, conv_flux)
195 8 : end function compute_vega_zero_point
196 :
197 : ! ab zero-point -- f_nu = 3631 Jy, converted to f_lambda = f_nu * c / lambda^2
198 : ! called once at init, result cached in filter_data
199 7 : real(dp) function compute_ab_zero_point(filt_wave, filt_trans)
200 : real(dp), dimension(:), intent(in) :: filt_wave, filt_trans
201 :
202 : real(dp) :: int_flux, int_filter
203 7 : real(dp), allocatable :: ab_sed_flux(:)
204 : integer :: i
205 :
206 21 : allocate (ab_sed_flux(size(filt_wave)))
207 :
208 : ! construct ab spectrum (f_lambda) on the filter wavelength grid
209 : ! 3631 Jy = 3.631e-20 erg/s/cm^2/Hz; clight in cm/s, wavelength in angstroms
210 149 : do i = 1, size(filt_wave)
211 149 : if (filt_wave(i) > 0.0_dp) then
212 142 : ab_sed_flux(i) = 3.631d-20*((clight*1.0d8)/(filt_wave(i)**2))
213 : else
214 0 : ab_sed_flux(i) = 0.0_dp
215 : end if
216 : end do
217 :
218 : ! photon-counting integration
219 149 : call simpson_integration(filt_wave, ab_sed_flux*filt_trans*filt_wave, int_flux)
220 149 : call simpson_integration(filt_wave, filt_wave*filt_trans, int_filter)
221 :
222 7 : if (int_filter > 0.0_dp) then
223 7 : compute_ab_zero_point = int_flux/int_filter
224 : else
225 : compute_ab_zero_point = -1.0_dp
226 : end if
227 :
228 7 : deallocate (ab_sed_flux)
229 7 : end function compute_ab_zero_point
230 :
231 : ! st zero-point -- f_lambda = 3.63e-9 erg/s/cm^2/A (flat spectrum)
232 : ! called once at init, result cached in filter_data
233 7 : real(dp) function compute_st_zero_point(filt_wave, filt_trans)
234 : real(dp), dimension(:), intent(in) :: filt_wave, filt_trans
235 :
236 : real(dp) :: int_flux, int_filter
237 7 : real(dp), allocatable :: st_sed_flux(:)
238 :
239 21 : allocate (st_sed_flux(size(filt_wave)))
240 149 : st_sed_flux = 3.63d-9
241 :
242 : ! photon-counting integration
243 149 : call simpson_integration(filt_wave, st_sed_flux*filt_trans*filt_wave, int_flux)
244 149 : call simpson_integration(filt_wave, filt_wave*filt_trans, int_filter)
245 :
246 7 : if (int_filter > 0.0_dp) then
247 7 : compute_st_zero_point = int_flux/int_filter
248 : else
249 : compute_st_zero_point = -1.0_dp
250 : end if
251 :
252 7 : deallocate (st_sed_flux)
253 7 : end function compute_st_zero_point
254 :
255 : end module synthetic
|
