Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2025 Niall Miller & The MESA Team
4 : ! Modified to include AB and ST magnitude systems.
5 : !
6 : ! ***********************************************************************
7 :
8 : module synthetic
9 :
10 : ! Added clight for AB conversions
11 : use const_def, only: dp, clight
12 : use utils_lib, only: mkdir, folder_exists
13 : use colors_utils, only: remove_dat, romberg_integration, load_filter, load_vega_sed, load_lookup_table
14 : use knn_interp, only: interpolate_array
15 :
16 : implicit none
17 :
18 : private
19 : public :: calculate_synthetic
20 :
21 : contains
22 :
23 : !****************************
24 : ! Calculate Synthetic Photometry Using SED and Filter
25 : !****************************
26 0 : real(dp) function calculate_synthetic(temperature, gravity, metallicity, ierr, &
27 0 : wavelengths, fluxes, filter_wavelengths, &
28 : filter_trans, &
29 : filter_filepath, vega_filepath, &
30 : filter_name, make_sed, colors_results_directory, &
31 : mag_system)
32 : ! Input arguments
33 : real(dp), intent(in) :: temperature, gravity, metallicity
34 : character(len=*), intent(in) :: filter_filepath, filter_name, vega_filepath, colors_results_directory
35 : character(len=*), intent(in) :: mag_system ! NEW: Arguments for 'Vega', 'AB', or 'ST'
36 : integer, intent(out) :: ierr
37 : character(len=1000) :: line
38 :
39 : real(dp), dimension(:), intent(inout) :: wavelengths, fluxes
40 : real(dp), dimension(:), allocatable, intent(inout) :: filter_wavelengths, filter_trans
41 : logical, intent(in) :: make_sed
42 :
43 : ! Local variables
44 0 : real(dp), dimension(:), allocatable :: convolved_flux
45 : character(len=100) :: csv_file
46 : real(dp) :: synthetic_flux, zero_point_flux
47 : integer :: max_size, i
48 : real(dp) :: wv, fl, cf, fwv, ftr
49 :
50 0 : if (.not. folder_exists(trim(colors_results_directory))) call mkdir(trim(colors_results_directory))
51 0 : csv_file = trim(colors_results_directory)//'/'//trim(remove_dat(filter_name))//'_SED.csv'
52 0 : ierr = 0
53 :
54 : ! Load filter data
55 0 : call load_filter(filter_filepath, filter_wavelengths, filter_trans)
56 :
57 : ! Perform SED convolution (Source Object)
58 0 : allocate (convolved_flux(size(wavelengths)))
59 0 : call convolve_sed(wavelengths, fluxes, filter_wavelengths, filter_trans, convolved_flux)
60 :
61 : ! Write SED to CSV if requested
62 0 : if (make_sed) then
63 : max_size = max(size(wavelengths), size(filter_wavelengths), &
64 0 : size(fluxes), size(convolved_flux), size(filter_trans))
65 :
66 0 : open (unit=10, file=csv_file, status='REPLACE', action='write', iostat=ierr)
67 0 : if (ierr /= 0) then
68 0 : print *, "Error opening file for writing"
69 0 : return
70 : end if
71 :
72 0 : write (10, '(A)') "wavelengths,fluxes,convolved_flux,filter_wavelengths,filter_trans"
73 :
74 0 : do i = 1, max_size
75 0 : wv = 0.0_dp; fl = 0.0_dp; cf = 0.0_dp; fwv = 0.0_dp; ftr = 0.0_dp
76 0 : if (i <= size(wavelengths)) wv = wavelengths(i)
77 0 : if (i <= size(fluxes)) fl = fluxes(i)
78 0 : if (i <= size(convolved_flux)) cf = convolved_flux(i)
79 0 : if (i <= size(filter_wavelengths)) fwv = filter_wavelengths(i)
80 0 : if (i <= size(filter_trans)) ftr = filter_trans(i)
81 :
82 : write (line, '(ES14.6, ",", ES14.6, ",", ES14.6, ",", ES14.6, ",", ES14.6)') &
83 0 : wv, fl, cf, fwv, ftr
84 0 : write (10, '(A)') trim(line)
85 : end do
86 0 : close (10)
87 : end if
88 :
89 : ! ------------------------------------------------------------------
90 : ! Calculate Zero Point Flux based on System Selection
91 : ! ------------------------------------------------------------------
92 0 : select case (trim(mag_system))
93 : case ('VEGA', 'Vega', 'vega')
94 : zero_point_flux = calculate_vega_flux(vega_filepath, filter_wavelengths, filter_trans, &
95 0 : filter_name, make_sed, colors_results_directory)
96 : case ('AB', 'ab')
97 0 : zero_point_flux = calculate_ab_zero_point(filter_wavelengths, filter_trans)
98 : case ('ST', 'st')
99 0 : zero_point_flux = calculate_st_zero_point(filter_wavelengths, filter_trans)
100 : case default
101 0 : print *, "Error: Unknown magnitude system: ", mag_system
102 0 : calculate_synthetic = huge(1.0_dp)
103 0 : return
104 : end select
105 :
106 : ! Calculate synthetic flux (Source Object)
107 : call calculate_synthetic_flux(wavelengths, convolved_flux, synthetic_flux, &
108 0 : filter_wavelengths, filter_trans)
109 :
110 : ! Calculate magnitude using the selected zero point
111 0 : if (zero_point_flux > 0.0_dp) then
112 0 : calculate_synthetic = -2.5d0*log10(synthetic_flux/zero_point_flux)
113 : else
114 0 : print *, "Error: Zero point flux is zero, magnitude calculation is invalid."
115 0 : calculate_synthetic = huge(1.0_dp)
116 : end if
117 :
118 : ! Clean up
119 0 : deallocate (convolved_flux)
120 0 : end function calculate_synthetic
121 :
122 : !-----------------------------------------------------------------------
123 : ! Internal functions for synthetic photometry
124 : !-----------------------------------------------------------------------
125 :
126 : !****************************
127 : ! Convolve SED With Filter
128 : !****************************
129 0 : subroutine convolve_sed(wavelengths, fluxes, filter_wavelengths, filter_trans, convolved_flux)
130 : real(dp), dimension(:), intent(inout) :: wavelengths, fluxes
131 : real(dp), dimension(:), intent(inout) :: filter_wavelengths, filter_trans
132 : real(dp), dimension(:), allocatable, intent(out) :: convolved_flux
133 : real(dp), dimension(:), allocatable :: interpolated_filter
134 : integer :: n
135 :
136 0 : n = size(wavelengths)
137 0 : allocate (interpolated_filter(n))
138 0 : call interpolate_array(filter_wavelengths, filter_trans, wavelengths, interpolated_filter)
139 0 : convolved_flux = fluxes*interpolated_filter
140 0 : deallocate (interpolated_filter)
141 0 : end subroutine convolve_sed
142 :
143 : !****************************
144 : ! Calculate Synthetic Flux (Integration)
145 : !****************************
146 0 : subroutine calculate_synthetic_flux(wavelengths, fluxes, synthetic_flux, &
147 0 : filter_wavelengths, filter_trans)
148 :
149 : real(dp), dimension(:), intent(in) :: wavelengths, fluxes
150 : real(dp), dimension(:), intent(inout) :: filter_wavelengths, filter_trans
151 : real(dp), intent(out) :: synthetic_flux
152 : integer :: i
153 : real(dp) :: integrated_flux, integrated_filter
154 :
155 : ! Validate inputs
156 0 : do i = 1, size(wavelengths) - 1
157 0 : if (wavelengths(i) <= 0.0_dp .or. fluxes(i) < 0.0_dp) then
158 0 : print *, "synthetic Invalid input at index", i, ":", wavelengths(i), fluxes(i)
159 0 : stop
160 : end if
161 : end do
162 :
163 0 : call romberg_integration(wavelengths, fluxes*wavelengths, integrated_flux)
164 : call romberg_integration(filter_wavelengths, &
165 0 : filter_trans*filter_wavelengths, integrated_filter)
166 :
167 0 : if (integrated_filter > 0.0_dp) then
168 0 : synthetic_flux = integrated_flux/integrated_filter
169 : else
170 0 : print *, "Error: Integrated filter transmission is zero."
171 0 : synthetic_flux = -1.0_dp
172 0 : return
173 : end if
174 : end subroutine calculate_synthetic_flux
175 :
176 : !****************************
177 : ! Calculate Vega Flux for Zero Point
178 : !****************************
179 0 : function calculate_vega_flux(vega_filepath, filt_wave, filt_trans, &
180 : filter_name, make_sed, colors_results_directory) result(vega_flux)
181 : character(len=*), intent(in) :: vega_filepath, filter_name, colors_results_directory
182 : character(len=100) :: output_csv
183 : real(dp), dimension(:), intent(inout) :: filt_wave, filt_trans
184 : real(dp) :: vega_flux
185 : real(dp) :: int_flux, int_filter
186 0 : real(dp), allocatable :: vega_wave(:), vega_flux_arr(:), conv_flux(:)
187 : logical, intent(in) :: make_sed
188 : integer :: i, max_size, ierr
189 : real(dp) :: wv, fl, cf, fwv, ftr
190 : character(len=1000) :: line
191 :
192 : ! Load the Vega SED
193 : call load_vega_sed(vega_filepath, vega_wave, vega_flux_arr)
194 :
195 : ! Convolve the Vega SED with the filter transmission
196 0 : allocate (conv_flux(size(vega_wave)))
197 0 : call convolve_sed(vega_wave, vega_flux_arr, filt_wave, filt_trans, conv_flux)
198 :
199 : ! Integrate
200 0 : call romberg_integration(vega_wave, vega_wave*conv_flux, int_flux)
201 0 : call romberg_integration(filt_wave, filt_wave*filt_trans, int_filter)
202 :
203 0 : if (int_filter > 0.0_dp) then
204 0 : vega_flux = int_flux/int_filter
205 : else
206 : vega_flux = -1.0_dp
207 : end if
208 :
209 : ! Write Vega SED to CSV if requested
210 0 : if (make_sed) then
211 : max_size = max(size(vega_wave), size(vega_flux_arr), size(conv_flux), &
212 0 : size(filt_wave), size(filt_trans))
213 :
214 0 : if (.not. folder_exists(trim(colors_results_directory))) call mkdir(trim(colors_results_directory))
215 0 : output_csv = trim(colors_results_directory)//'/VEGA_'//trim(remove_dat(filter_name))//'_SED.csv'
216 :
217 0 : open (unit=10, file=output_csv, status='REPLACE', action='write', iostat=ierr)
218 0 : if (ierr /= 0) then
219 0 : print *, "Error opening file for writing"
220 0 : return
221 : end if
222 :
223 0 : write (10, '(A)') "wavelengths,fluxes,convolved_flux,filter_wavelengths,filter_trans"
224 :
225 0 : do i = 1, max_size
226 0 : wv = 0.0_dp; fl = 0.0_dp; cf = 0.0_dp; fwv = 0.0_dp; ftr = 0.0_dp
227 0 : if (i <= size(vega_wave)) wv = vega_wave(i)
228 0 : if (i <= size(vega_flux_arr)) fl = vega_flux_arr(i)
229 0 : if (i <= size(conv_flux)) cf = conv_flux(i)
230 0 : if (i <= size(filt_wave)) fwv = filt_wave(i)
231 0 : if (i <= size(filt_trans)) ftr = filt_trans(i)
232 :
233 : write (line, '(ES14.6, ",", ES14.6, ",", ES14.6, ",", ES14.6, ",", ES14.6)') &
234 0 : wv, fl, cf, fwv, ftr
235 0 : write (10, '(A)') trim(line)
236 : end do
237 0 : close (10)
238 : end if
239 0 : deallocate (conv_flux, vega_wave, vega_flux_arr)
240 0 : end function calculate_vega_flux
241 :
242 : !****************************
243 : ! Calculate AB Zero Point Flux
244 : ! f_nu = 3631 Jy = 3.631e-20 erg/s/cm^2/Hz
245 : ! f_lambda = f_nu * c / lambda^2
246 : !****************************
247 0 : function calculate_ab_zero_point(filt_wave, filt_trans) result(ab_flux)
248 : real(dp), dimension(:), intent(inout) :: filt_wave, filt_trans
249 : real(dp) :: ab_flux
250 : real(dp) :: int_flux, int_filter
251 0 : real(dp), allocatable :: ab_sed_flux(:)
252 : integer :: i
253 :
254 0 : allocate(ab_sed_flux(size(filt_wave)))
255 :
256 : ! Construct AB Spectrum (f_lambda) on the filter wavelength grid
257 : ! Assumes wavelengths are in Angstroms and clight is in Angstroms/sec
258 : ! 3631 Jy = 3.631E-20 erg/s/cm^2/Hz
259 0 : do i = 1, size(filt_wave)
260 0 : if (filt_wave(i) > 0.0_dp) then
261 0 : ab_sed_flux(i) = 3.631d-20 * (clight / (filt_wave(i)**2))
262 : else
263 0 : ab_sed_flux(i) = 0.0_dp
264 : endif
265 : end do
266 :
267 : ! Integrate using same method as source (f_lambda * T * lambda)
268 : ! Note: We multiply by filt_wave inside the integration because the
269 : ! romberg helper expects (flux * lambda)
270 0 : call romberg_integration(filt_wave, ab_sed_flux * filt_trans * filt_wave, int_flux)
271 0 : call romberg_integration(filt_wave, filt_wave * filt_trans, int_filter)
272 :
273 0 : if (int_filter > 0.0_dp) then
274 0 : ab_flux = int_flux / int_filter
275 : else
276 : ab_flux = -1.0_dp
277 : end if
278 :
279 0 : deallocate(ab_sed_flux)
280 0 : end function calculate_ab_zero_point
281 :
282 : !****************************
283 : ! Calculate ST Zero Point Flux
284 : ! f_lambda = 3.63e-9 erg/s/cm^2/A (Constant)
285 : !****************************
286 0 : function calculate_st_zero_point(filt_wave, filt_trans) result(st_flux)
287 : real(dp), dimension(:), intent(inout) :: filt_wave, filt_trans
288 : real(dp) :: st_flux
289 : real(dp) :: int_flux, int_filter
290 0 : real(dp), allocatable :: st_sed_flux(:)
291 :
292 : ! For ST system, flux is constant in wavelength
293 : ! However, to maintain exact consistency with how the source is integrated
294 : ! (numerical integration over the filter grid), we integrate the constant array.
295 :
296 0 : allocate(st_sed_flux(size(filt_wave)))
297 0 : st_sed_flux = 3.63d-9
298 :
299 0 : call romberg_integration(filt_wave, st_sed_flux * filt_trans * filt_wave, int_flux)
300 0 : call romberg_integration(filt_wave, filt_wave * filt_trans, int_filter)
301 :
302 0 : if (int_filter > 0.0_dp) then
303 0 : st_flux = int_flux / int_filter
304 : else
305 : st_flux = -1.0_dp
306 : end if
307 :
308 0 : deallocate(st_sed_flux)
309 0 : end function calculate_st_zero_point
310 :
311 : end module synthetic
|
