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 :
33 : contains
34 :
35 : ! calculate synthetic magnitude; uses precomputed zero-point from filter_data
36 105 : real(dp) function calculate_synthetic(temperature, gravity, metallicity, ierr, &
37 105 : wavelengths, fluxes, &
38 105 : filter_wavelengths, filter_trans, &
39 : zero_point_flux, &
40 : filter_name, make_sed, sed_per_model, colors_results_directory, model_number)
41 : real(dp), intent(in) :: temperature, gravity, metallicity
42 : character(len=*), intent(in) :: filter_name, colors_results_directory
43 : integer, intent(out) :: ierr
44 :
45 : real(dp), dimension(:), intent(in) :: wavelengths, fluxes
46 : real(dp), dimension(:), intent(in) :: filter_wavelengths, filter_trans
47 : real(dp), intent(in) :: zero_point_flux ! precomputed at initialization
48 : logical, intent(in) :: make_sed, sed_per_model
49 : integer, intent(in) :: model_number
50 :
51 105 : real(dp), dimension(:), allocatable :: convolved_flux, filter_on_sed_grid
52 : character(len=256) :: csv_file
53 : character(len=20) :: model_str
54 : character(len=1000) :: line
55 : real(dp) :: synthetic_flux
56 : integer :: max_size, i
57 : real(dp) :: wv, fl, cf, fwv, ftr
58 :
59 105 : ierr = 0
60 :
61 315 : allocate (convolved_flux(size(wavelengths)))
62 210 : allocate (filter_on_sed_grid(size(wavelengths)))
63 :
64 105 : call interpolate_array(filter_wavelengths, filter_trans, wavelengths, filter_on_sed_grid)
65 :
66 126105 : convolved_flux = fluxes*filter_on_sed_grid
67 :
68 105 : if (make_sed) then
69 0 : if (.not. folder_exists(trim(colors_results_directory))) call mkdir(trim(colors_results_directory))
70 :
71 : ! track model number internally when write_sed_per_model is enabled
72 0 : if (sed_per_model) then
73 :
74 0 : write (model_str, '(I8.8)') model_number
75 0 : csv_file = trim(colors_results_directory)//'/'//trim(remove_dat(filter_name))//'_SED_'//trim(model_str)//'.csv'
76 :
77 : else
78 0 : csv_file = trim(colors_results_directory)//'/'//trim(remove_dat(filter_name))//'_SED.csv'
79 : end if
80 :
81 0 : max_size = max(size(wavelengths), size(filter_wavelengths))
82 :
83 0 : open (unit=10, file=csv_file, status='REPLACE', action='write', iostat=ierr)
84 0 : if (ierr /= 0) then
85 0 : print *, "Error opening file for writing"
86 0 : deallocate (convolved_flux, filter_on_sed_grid)
87 0 : calculate_synthetic = huge(1.0_dp)
88 : return
89 : end if
90 :
91 0 : write (10, '(A)') "wavelengths,fluxes,convolved_flux,filter_wavelengths,filter_trans"
92 :
93 0 : do i = 1, max_size
94 0 : wv = 0.0_dp; fl = 0.0_dp; cf = 0.0_dp; fwv = 0.0_dp; ftr = 0.0_dp
95 0 : if (i <= size(wavelengths)) wv = wavelengths(i)
96 0 : if (i <= size(fluxes)) fl = fluxes(i)
97 0 : if (i <= size(convolved_flux)) cf = convolved_flux(i)
98 0 : if (i <= size(filter_wavelengths)) fwv = filter_wavelengths(i)
99 0 : if (i <= size(filter_trans)) ftr = filter_trans(i)
100 :
101 : write (line, '(ES14.6, ",", ES14.6, ",", ES14.6, ",", ES14.6, ",", ES14.6)') &
102 0 : wv, fl, cf, fwv, ftr
103 0 : write (10, '(A)') trim(line)
104 : end do
105 0 : close (10)
106 : end if
107 :
108 105 : call calculate_synthetic_flux(wavelengths, convolved_flux, filter_on_sed_grid, synthetic_flux)
109 :
110 105 : if (zero_point_flux > 0.0_dp .and. synthetic_flux > 0.0_dp) then
111 105 : 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 105 : deallocate (convolved_flux, filter_on_sed_grid)
123 105 : end function calculate_synthetic
124 :
125 : ! photon-counting synthetic flux (numerator and denominator weighted by lambda)
126 105 : subroutine calculate_synthetic_flux(wavelengths, convolved_flux, filter_on_sed_grid, synthetic_flux)
127 : real(dp), dimension(:), intent(in) :: wavelengths, convolved_flux, filter_on_sed_grid
128 : real(dp), intent(out) :: synthetic_flux
129 :
130 : real(dp) :: integrated_flux, integrated_filter
131 :
132 : ! photon-counting: weight by wavelength
133 126000 : call simpson_integration(wavelengths, convolved_flux*wavelengths, integrated_flux)
134 126000 : call simpson_integration(wavelengths, filter_on_sed_grid*wavelengths, integrated_filter)
135 :
136 105 : if (integrated_filter > 0.0_dp) then
137 105 : synthetic_flux = integrated_flux/integrated_filter
138 : else
139 0 : print *, "Error: Integrated filter transmission is zero."
140 0 : synthetic_flux = -1.0_dp
141 : end if
142 105 : end subroutine calculate_synthetic_flux
143 :
144 : ! vega zero-point -- called once at init, result cached in filter_data
145 7 : real(dp) function compute_vega_zero_point(vega_wave, vega_flux, filt_wave, filt_trans)
146 : real(dp), dimension(:), intent(in) :: vega_wave, vega_flux
147 : real(dp), dimension(:), intent(in) :: filt_wave, filt_trans
148 :
149 : real(dp) :: int_flux, int_filter
150 7 : real(dp), allocatable :: filt_on_vega_grid(:), conv_flux(:)
151 :
152 21 : allocate (filt_on_vega_grid(size(vega_wave)))
153 14 : allocate (conv_flux(size(vega_wave)))
154 :
155 7 : call interpolate_array(filt_wave, filt_trans, vega_wave, filt_on_vega_grid)
156 :
157 56672 : conv_flux = vega_flux*filt_on_vega_grid
158 :
159 : ! photon-counting integration
160 56665 : call simpson_integration(vega_wave, vega_wave*conv_flux, int_flux)
161 56665 : call simpson_integration(vega_wave, vega_wave*filt_on_vega_grid, int_filter)
162 :
163 7 : if (int_filter > 0.0_dp) then
164 7 : compute_vega_zero_point = int_flux/int_filter
165 : else
166 : compute_vega_zero_point = -1.0_dp
167 : end if
168 :
169 7 : deallocate (filt_on_vega_grid, conv_flux)
170 7 : end function compute_vega_zero_point
171 :
172 : ! ab zero-point -- f_nu = 3631 Jy, converted to f_lambda = f_nu * c / lambda^2
173 : ! called once at init, result cached in filter_data
174 7 : real(dp) function compute_ab_zero_point(filt_wave, filt_trans)
175 : real(dp), dimension(:), intent(in) :: filt_wave, filt_trans
176 :
177 : real(dp) :: int_flux, int_filter
178 7 : real(dp), allocatable :: ab_sed_flux(:)
179 : integer :: i
180 :
181 21 : allocate (ab_sed_flux(size(filt_wave)))
182 :
183 : ! construct ab spectrum (f_lambda) on the filter wavelength grid
184 : ! 3631 Jy = 3.631e-20 erg/s/cm^2/Hz; clight in cm/s, wavelength in angstroms
185 149 : do i = 1, size(filt_wave)
186 149 : if (filt_wave(i) > 0.0_dp) then
187 142 : ab_sed_flux(i) = 3.631d-20*((clight*1.0d8)/(filt_wave(i)**2))
188 : else
189 0 : ab_sed_flux(i) = 0.0_dp
190 : end if
191 : end do
192 :
193 : ! photon-counting integration
194 149 : call simpson_integration(filt_wave, ab_sed_flux*filt_trans*filt_wave, int_flux)
195 149 : call simpson_integration(filt_wave, filt_wave*filt_trans, int_filter)
196 :
197 7 : if (int_filter > 0.0_dp) then
198 7 : compute_ab_zero_point = int_flux/int_filter
199 : else
200 : compute_ab_zero_point = -1.0_dp
201 : end if
202 :
203 7 : deallocate (ab_sed_flux)
204 7 : end function compute_ab_zero_point
205 :
206 : ! st zero-point -- f_lambda = 3.63e-9 erg/s/cm^2/A (flat spectrum)
207 : ! called once at init, result cached in filter_data
208 7 : real(dp) function compute_st_zero_point(filt_wave, filt_trans)
209 : real(dp), dimension(:), intent(in) :: filt_wave, filt_trans
210 :
211 : real(dp) :: int_flux, int_filter
212 7 : real(dp), allocatable :: st_sed_flux(:)
213 :
214 21 : allocate (st_sed_flux(size(filt_wave)))
215 149 : st_sed_flux = 3.63d-9
216 :
217 : ! photon-counting integration
218 149 : call simpson_integration(filt_wave, st_sed_flux*filt_trans*filt_wave, int_flux)
219 149 : call simpson_integration(filt_wave, filt_wave*filt_trans, int_filter)
220 :
221 7 : if (int_filter > 0.0_dp) then
222 7 : compute_st_zero_point = int_flux/int_filter
223 : else
224 : compute_st_zero_point = -1.0_dp
225 : end if
226 :
227 7 : deallocate (st_sed_flux)
228 7 : end function compute_st_zero_point
229 :
230 : end module synthetic
|
