Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2025 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 synthetic
21 :
22 : use const_def, only: dp
23 : use utils_lib, only: mkdir, folder_exists
24 : use colors_utils, only: remove_dat, romberg_integration, load_filter, load_vega_sed, load_lookup_table
25 : use knn_interp, only: interpolate_array
26 :
27 : implicit none
28 :
29 : private
30 : public :: calculate_synthetic
31 :
32 : contains
33 :
34 : !****************************
35 : ! Calculate Synthetic Photometry Using SED and Filter
36 : !****************************
37 0 : real(dp) function calculate_synthetic(temperature, gravity, metallicity, ierr, &
38 0 : wavelengths, fluxes, filter_wavelengths, &
39 : filter_trans, &
40 : filter_filepath, vega_filepath, &
41 : filter_name, make_sed, colors_results_directory)
42 : ! Input arguments
43 : real(dp), intent(in) :: temperature, gravity, metallicity
44 : character(len=*), intent(in) :: filter_filepath, filter_name, vega_filepath, colors_results_directory
45 : integer, intent(out) :: ierr
46 : character(len=1000) :: line
47 :
48 : real(dp), dimension(:), intent(inout) :: wavelengths, fluxes
49 : real(dp), dimension(:), allocatable, intent(inout) :: filter_wavelengths, filter_trans
50 : logical, intent(in) :: make_sed
51 :
52 : ! Local variables
53 0 : real(dp), dimension(:), allocatable :: convolved_flux
54 : character(len=100) :: csv_file
55 0 : real(dp) :: synthetic_flux, vega_flux
56 : integer :: max_size, i
57 0 : real(dp) :: wv, fl, cf, fwv, ftr
58 :
59 0 : if (.not. folder_exists(trim(colors_results_directory))) call mkdir(trim(colors_results_directory))
60 0 : csv_file = trim(colors_results_directory)//'/'//trim(remove_dat(filter_name))//'_SED.csv'
61 0 : ierr = 0
62 :
63 : ! Load filter data
64 0 : call load_filter(filter_filepath, filter_wavelengths, filter_trans)
65 :
66 :
67 : ! Perform SED convolution
68 0 : allocate (convolved_flux(size(wavelengths)))
69 0 : call convolve_sed(wavelengths, fluxes, filter_wavelengths, filter_trans, convolved_flux)
70 :
71 : ! Write SED to CSV if requested
72 0 : if (make_sed) then
73 : ! Determine the maximum size among all arrays
74 : max_size = max(size(wavelengths), size(filter_wavelengths), &
75 0 : size(fluxes), size(convolved_flux), size(filter_trans))
76 :
77 : ! Open the CSV file for writing
78 0 : open (unit=10, file=csv_file, status='REPLACE', action='write', iostat=ierr)
79 0 : if (ierr /= 0) then
80 0 : print *, "Error opening file for writing"
81 : return
82 : end if
83 :
84 : ! Write headers to the CSV file
85 0 : write (10, '(A)') "wavelengths,fluxes,convolved_flux,filter_wavelengths,filter_trans"
86 :
87 : ! Loop through data and safely write values, ensuring no out-of-bounds errors
88 0 : do i = 1, max_size
89 : ! Initialize values to zero in case they are out of bounds
90 0 : wv = 0.0_dp
91 0 : fl = 0.0_dp
92 0 : cf = 0.0_dp
93 0 : fwv = 0.0_dp
94 0 : ftr = 0.0_dp
95 :
96 : ! Assign actual values only if within valid indices
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 the formatted output
104 : write (line, '(ES14.6, ",", ES14.6, ",", ES14.6, ",", ES14.6, ",", ES14.6)') &
105 0 : wv, fl, cf, fwv, ftr
106 0 : write (10, '(A)') trim(line)
107 : end do
108 :
109 : ! Close the file
110 0 : close (10)
111 : end if
112 :
113 : ! Calculate Vega flux for zero point calibration
114 : vega_flux = calculate_vega_flux(vega_filepath, filter_wavelengths, filter_trans, &
115 0 : filter_name, make_sed, colors_results_directory)
116 :
117 : ! Calculate synthetic flux
118 : call calculate_synthetic_flux(wavelengths, convolved_flux, synthetic_flux, &
119 0 : filter_wavelengths, filter_trans)
120 :
121 : ! Calculate magnitude using Vega zero point
122 0 : if (vega_flux > 0.0_dp) then
123 0 : calculate_synthetic = -2.5d0*log10(synthetic_flux/vega_flux)
124 : else
125 0 : print *, "Error: Vega flux is zero, magnitude calculation is invalid."
126 0 : calculate_synthetic = huge(1.0_dp)
127 : end if
128 :
129 : ! Clean up
130 0 : deallocate (convolved_flux)
131 0 : end function calculate_synthetic
132 :
133 : !-----------------------------------------------------------------------
134 : ! Internal functions for synthetic photometry
135 : !-----------------------------------------------------------------------
136 :
137 : !****************************
138 : ! Convolve SED With Filter
139 : !****************************
140 0 : subroutine convolve_sed(wavelengths, fluxes, filter_wavelengths, filter_trans, convolved_flux)
141 : real(dp), dimension(:), intent(inout) :: wavelengths, fluxes
142 : real(dp), dimension(:), intent(inout) :: filter_wavelengths, filter_trans
143 : real(dp), dimension(:), allocatable, intent(out) :: convolved_flux
144 : real(dp), dimension(:), allocatable :: interpolated_filter
145 : integer :: n
146 :
147 0 : n = size(wavelengths)
148 :
149 : ! Allocate arrays
150 0 : allocate (interpolated_filter(n))
151 :
152 : ! Interpolate the filter transmission onto the wavelengths array
153 0 : call interpolate_array(filter_wavelengths, filter_trans, wavelengths, interpolated_filter)
154 :
155 : ! Perform convolution (element-wise multiplication)
156 0 : convolved_flux = fluxes*interpolated_filter
157 :
158 : ! Deallocate temporary arrays
159 0 : deallocate (interpolated_filter)
160 0 : end subroutine convolve_sed
161 :
162 : !****************************
163 : ! Calculate Synthetic Flux
164 : !****************************
165 0 : subroutine calculate_synthetic_flux(wavelengths, fluxes, synthetic_flux, &
166 0 : filter_wavelengths, filter_trans)
167 :
168 : real(dp), dimension(:), intent(in) :: wavelengths, fluxes
169 : real(dp), dimension(:), intent(inout) :: filter_wavelengths, filter_trans
170 : real(dp), intent(out) :: synthetic_flux
171 : integer :: i
172 0 : real(dp) :: integrated_flux, integrated_filter
173 :
174 : ! Validate inputs
175 0 : do i = 1, size(wavelengths) - 1
176 0 : if (wavelengths(i) <= 0.0_dp .or. fluxes(i) < 0.0_dp) then
177 0 : print *, "synthetic Invalid input at index", i, ":", wavelengths(i), fluxes(i)
178 0 : stop
179 : end if
180 : end do
181 :
182 0 : call romberg_integration(wavelengths, fluxes*wavelengths, integrated_flux)
183 : call romberg_integration(filter_wavelengths, &
184 0 : filter_trans*filter_wavelengths, integrated_filter)
185 : ! Store the total flux
186 0 : if (integrated_filter > 0.0_dp) then
187 0 : synthetic_flux = integrated_flux/integrated_filter
188 : else
189 0 : print *, "Error: Integrated filter transmission is zero."
190 0 : synthetic_flux = -1.0_dp
191 0 : return
192 : end if
193 : end subroutine calculate_synthetic_flux
194 :
195 : !****************************
196 : ! Calculate Vega Flux for Zero Point
197 : !****************************
198 0 : function calculate_vega_flux(vega_filepath, filt_wave, filt_trans, &
199 0 : filter_name, make_sed, colors_results_directory) result(vega_flux)
200 : character(len=*), intent(in) :: vega_filepath, filter_name, colors_results_directory
201 : character(len=100) :: output_csv
202 : real(dp), dimension(:), intent(inout) :: filt_wave, filt_trans
203 : real(dp) :: vega_flux
204 0 : real(dp) :: int_flux, int_filter
205 0 : real(dp), allocatable :: vega_wave(:), vega_flux_arr(:), conv_flux(:)
206 : logical, intent(in) :: make_sed
207 : integer :: i, unit, max_size
208 0 : real(dp) :: wv, fl, cf, fwv, ftr
209 : integer:: ierr
210 : character(len=1000) :: line
211 :
212 : ! Load the Vega SED
213 : call load_vega_sed(vega_filepath, vega_wave, vega_flux_arr)
214 :
215 : ! Convolve the Vega SED with the filter transmission
216 0 : allocate (conv_flux(size(vega_wave)))
217 0 : call convolve_sed(vega_wave, vega_flux_arr, filt_wave, filt_trans, conv_flux)
218 :
219 : ! Integrate the convolved Vega SED and the filter transmission
220 0 : call romberg_integration(vega_wave, vega_wave*conv_flux, int_flux)
221 0 : call romberg_integration(filt_wave, filt_wave*filt_trans, int_filter)
222 :
223 0 : if (int_filter > 0.0_dp) then
224 0 : vega_flux = int_flux/int_filter
225 : else
226 : vega_flux = -1.0_dp
227 : end if
228 :
229 : ! Write Vega SED to CSV if requested
230 0 : if (make_sed) then
231 : ! Determine the maximum size among all arrays
232 : max_size = max(size(vega_wave), size(vega_flux_arr), size(conv_flux), &
233 0 : size(filt_wave), size(filt_trans))
234 :
235 0 : if (.not. folder_exists(trim(colors_results_directory))) call mkdir(trim(colors_results_directory))
236 0 : output_csv = trim(colors_results_directory)//'/VEGA_'//trim(remove_dat(filter_name))//'_SED.csv'
237 :
238 : ! Open the CSV file for writing
239 0 : open (unit=10, file=output_csv, status='REPLACE', action='write', iostat=ierr)
240 0 : if (ierr /= 0) then
241 0 : print *, "Error opening file for writing"
242 0 : return
243 : end if
244 :
245 0 : write (10, '(A)') "wavelengths,fluxes,convolved_flux,filter_wavelengths,filter_trans"
246 :
247 : ! Loop through data and safely write values, ensuring no out-of-bounds errors
248 0 : do i = 1, max_size
249 : ! Initialize values to zero in case they are out of bounds
250 0 : wv = 0.0_dp
251 0 : fl = 0.0_dp
252 0 : cf = 0.0_dp
253 0 : fwv = 0.0_dp
254 0 : ftr = 0.0_dp
255 :
256 : ! Assign actual values only if within valid indices
257 0 : if (i <= size(vega_wave)) wv = vega_wave(i)
258 0 : if (i <= size(vega_flux_arr)) fl = vega_flux_arr(i)
259 0 : if (i <= size(conv_flux)) cf = conv_flux(i)
260 0 : if (i <= size(filt_wave)) fwv = filt_wave(i)
261 0 : if (i <= size(filt_trans)) ftr = filt_trans(i)
262 :
263 : ! Write the formatted output
264 : write (line, '(ES14.6, ",", ES14.6, ",", ES14.6, ",", ES14.6, ",", ES14.6)') &
265 0 : wv, fl, cf, fwv, ftr
266 0 : write (10, '(A)') trim(line)
267 : end do
268 :
269 : ! Close the file
270 0 : close (10)
271 : end if
272 :
273 : ! Clean up
274 0 : deallocate (conv_flux, vega_wave, vega_flux_arr)
275 0 : end function calculate_vega_flux
276 :
277 : end module synthetic
|
