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, romberg_integration
25 : use knn_interp, only: interpolate_array
26 :
27 : implicit none
28 :
29 : private
30 : public :: calculate_synthetic
31 : ! Export zero-point computation functions for precomputation at initialization
32 : public :: compute_vega_zero_point, compute_ab_zero_point, compute_st_zero_point
33 :
34 : contains
35 :
36 : !****************************
37 : ! Calculate Synthetic Photometry Using SED and Filter
38 : ! Uses precomputed zero-point from filter data
39 : !****************************
40 0 : real(dp) function calculate_synthetic(temperature, gravity, metallicity, ierr, &
41 0 : wavelengths, fluxes, &
42 0 : filter_wavelengths, filter_trans, &
43 : zero_point_flux, &
44 : filter_name, make_sed, colors_results_directory)
45 : ! Input arguments
46 : real(dp), intent(in) :: temperature, gravity, metallicity
47 : character(len=*), intent(in) :: filter_name, colors_results_directory
48 : integer, intent(out) :: ierr
49 :
50 : real(dp), dimension(:), intent(in) :: wavelengths, fluxes
51 : real(dp), dimension(:), intent(in) :: filter_wavelengths, filter_trans
52 : real(dp), intent(in) :: zero_point_flux ! precomputed at initialization
53 : logical, intent(in) :: make_sed
54 :
55 : ! Local variables
56 : real(dp), dimension(:), allocatable :: convolved_flux, filter_on_sed_grid
57 : character(len=256) :: csv_file
58 : character(len=1000) :: line
59 : real(dp) :: synthetic_flux
60 : integer :: max_size, i
61 : real(dp) :: wv, fl, cf, fwv, ftr
62 :
63 0 : ierr = 0
64 :
65 : ! Allocate working arrays
66 0 : allocate(convolved_flux(size(wavelengths)))
67 0 : allocate(filter_on_sed_grid(size(wavelengths)))
68 :
69 : ! Interpolate filter onto SED wavelength grid
70 0 : call interpolate_array(filter_wavelengths, filter_trans, wavelengths, filter_on_sed_grid)
71 :
72 : ! Convolve SED with filter
73 0 : convolved_flux = fluxes * filter_on_sed_grid
74 :
75 : ! Write SED to CSV if requested
76 0 : if (make_sed) then
77 0 : if (.not. folder_exists(trim(colors_results_directory))) call mkdir(trim(colors_results_directory))
78 0 : csv_file = trim(colors_results_directory)//'/'//trim(remove_dat(filter_name))//'_SED.csv'
79 :
80 0 : max_size = max(size(wavelengths), size(filter_wavelengths))
81 :
82 0 : open (unit=10, file=csv_file, status='REPLACE', action='write', iostat=ierr)
83 0 : if (ierr /= 0) then
84 0 : print *, "Error opening file for writing"
85 0 : deallocate(convolved_flux, filter_on_sed_grid)
86 : return
87 : end if
88 :
89 0 : write (10, '(A)') "wavelengths,fluxes,convolved_flux,filter_wavelengths,filter_trans"
90 :
91 0 : do i = 1, max_size
92 0 : wv = 0.0_dp; fl = 0.0_dp; cf = 0.0_dp; fwv = 0.0_dp; ftr = 0.0_dp
93 0 : if (i <= size(wavelengths)) wv = wavelengths(i)
94 0 : if (i <= size(fluxes)) fl = fluxes(i)
95 0 : if (i <= size(convolved_flux)) cf = convolved_flux(i)
96 0 : if (i <= size(filter_wavelengths)) fwv = filter_wavelengths(i)
97 0 : if (i <= size(filter_trans)) ftr = filter_trans(i)
98 :
99 : write (line, '(ES14.6, ",", ES14.6, ",", ES14.6, ",", ES14.6, ",", ES14.6)') &
100 0 : wv, fl, cf, fwv, ftr
101 0 : write (10, '(A)') trim(line)
102 : end do
103 0 : close (10)
104 : end if
105 :
106 : ! Calculate synthetic flux using photon-counting integration
107 0 : call calculate_synthetic_flux(wavelengths, convolved_flux, filter_on_sed_grid, synthetic_flux)
108 :
109 : ! Calculate magnitude
110 0 : if (zero_point_flux > 0.0_dp .and. synthetic_flux > 0.0_dp) then
111 0 : calculate_synthetic = -2.5d0 * log10(synthetic_flux / zero_point_flux)
112 : else
113 0 : if (zero_point_flux <= 0.0_dp) then
114 0 : print *, "Error: Zero point flux is zero or negative for filter ", trim(filter_name)
115 : end if
116 0 : if (synthetic_flux <= 0.0_dp) then
117 0 : print *, "Error: Synthetic flux is zero or negative for filter ", trim(filter_name)
118 : end if
119 : calculate_synthetic = huge(1.0_dp)
120 : end if
121 :
122 0 : deallocate(convolved_flux, filter_on_sed_grid)
123 0 : end function calculate_synthetic
124 :
125 : !****************************
126 : ! Calculate Synthetic Flux (photon-counting integration)
127 : !****************************
128 0 : subroutine calculate_synthetic_flux(wavelengths, convolved_flux, filter_on_sed_grid, synthetic_flux)
129 : real(dp), dimension(:), intent(in) :: wavelengths, convolved_flux, filter_on_sed_grid
130 : real(dp), intent(out) :: synthetic_flux
131 :
132 : real(dp) :: integrated_flux, integrated_filter
133 :
134 : ! Photon-counting: weight by wavelength
135 0 : call romberg_integration(wavelengths, convolved_flux * wavelengths, integrated_flux)
136 0 : call romberg_integration(wavelengths, filter_on_sed_grid * wavelengths, integrated_filter)
137 :
138 0 : if (integrated_filter > 0.0_dp) then
139 0 : synthetic_flux = integrated_flux / integrated_filter
140 : else
141 0 : print *, "Error: Integrated filter transmission is zero."
142 0 : synthetic_flux = -1.0_dp
143 : end if
144 0 : end subroutine calculate_synthetic_flux
145 :
146 : !****************************
147 : ! Compute Vega Zero Point Flux
148 : ! Called once at initialization, result cached in filter_data
149 : !****************************
150 7 : real(dp) function compute_vega_zero_point(vega_wave, vega_flux, filt_wave, filt_trans)
151 : real(dp), dimension(:), intent(in) :: vega_wave, vega_flux
152 : real(dp), dimension(:), intent(in) :: filt_wave, filt_trans
153 :
154 : real(dp) :: int_flux, int_filter
155 7 : real(dp), allocatable :: filt_on_vega_grid(:), conv_flux(:)
156 :
157 7 : allocate(filt_on_vega_grid(size(vega_wave)))
158 7 : allocate(conv_flux(size(vega_wave)))
159 :
160 : ! Interpolate filter onto Vega wavelength grid
161 7 : call interpolate_array(filt_wave, filt_trans, vega_wave, filt_on_vega_grid)
162 :
163 : ! Convolve Vega with filter
164 56672 : conv_flux = vega_flux * filt_on_vega_grid
165 :
166 : ! Photon-counting integration
167 56665 : call romberg_integration(vega_wave, vega_wave * conv_flux, int_flux)
168 56665 : call romberg_integration(vega_wave, vega_wave * filt_on_vega_grid, int_filter)
169 :
170 7 : if (int_filter > 0.0_dp) then
171 7 : compute_vega_zero_point = int_flux / int_filter
172 : else
173 : compute_vega_zero_point = -1.0_dp
174 : end if
175 :
176 7 : deallocate(filt_on_vega_grid, conv_flux)
177 7 : end function compute_vega_zero_point
178 :
179 : !****************************
180 : ! Compute AB Zero Point Flux
181 : ! f_nu = 3631 Jy = 3.631e-20 erg/s/cm^2/Hz
182 : ! f_lambda = f_nu * c / lambda^2
183 : ! Called once at initialization, result cached in filter_data
184 : !****************************
185 7 : real(dp) function compute_ab_zero_point(filt_wave, filt_trans)
186 : real(dp), dimension(:), intent(in) :: filt_wave, filt_trans
187 :
188 : real(dp) :: int_flux, int_filter
189 7 : real(dp), allocatable :: ab_sed_flux(:)
190 : integer :: i
191 :
192 7 : allocate(ab_sed_flux(size(filt_wave)))
193 :
194 : ! Construct AB spectrum (f_lambda) on the filter wavelength grid
195 : ! 3631 Jy = 3.631E-20 erg/s/cm^2/Hz
196 : ! clight in cm/s, wavelength in Angstroms, need to convert
197 149 : do i = 1, size(filt_wave)
198 149 : if (filt_wave(i) > 0.0_dp) then
199 142 : ab_sed_flux(i) = 3.631d-20 * ((clight * 1.0d8) / (filt_wave(i)**2))
200 : else
201 0 : ab_sed_flux(i) = 0.0_dp
202 : end if
203 : end do
204 :
205 : ! Photon-counting integration
206 149 : call romberg_integration(filt_wave, ab_sed_flux * filt_trans * filt_wave, int_flux)
207 149 : call romberg_integration(filt_wave, filt_wave * filt_trans, int_filter)
208 :
209 7 : if (int_filter > 0.0_dp) then
210 7 : compute_ab_zero_point = int_flux / int_filter
211 : else
212 : compute_ab_zero_point = -1.0_dp
213 : end if
214 :
215 7 : deallocate(ab_sed_flux)
216 7 : end function compute_ab_zero_point
217 :
218 : !****************************
219 : ! Compute ST Zero Point Flux
220 : ! f_lambda = 3.63e-9 erg/s/cm^2/A (Constant)
221 : ! Called once at initialization, result cached in filter_data
222 : !****************************
223 7 : real(dp) function compute_st_zero_point(filt_wave, filt_trans)
224 : real(dp), dimension(:), intent(in) :: filt_wave, filt_trans
225 :
226 : real(dp) :: int_flux, int_filter
227 7 : real(dp), allocatable :: st_sed_flux(:)
228 :
229 7 : allocate(st_sed_flux(size(filt_wave)))
230 149 : st_sed_flux = 3.63d-9
231 :
232 : ! Photon-counting integration
233 149 : call romberg_integration(filt_wave, st_sed_flux * filt_trans * filt_wave, int_flux)
234 149 : call romberg_integration(filt_wave, filt_wave * filt_trans, int_filter)
235 :
236 7 : if (int_filter > 0.0_dp) then
237 7 : compute_st_zero_point = int_flux / int_filter
238 : else
239 : compute_st_zero_point = -1.0_dp
240 : end if
241 :
242 7 : deallocate(st_sed_flux)
243 7 : end function compute_st_zero_point
244 :
245 : end module synthetic
|
