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 colors_utils
21 : use const_def, only: dp, strlen, mesa_dir
22 : use colors_def, only: Colors_General_Info
23 :
24 : implicit none
25 :
26 : public :: dilute_flux, trapezoidal_integration, romberg_integration, &
27 : simpson_integration, load_sed, load_filter, load_vega_sed, &
28 : load_lookup_table, remove_dat
29 : contains
30 :
31 : !---------------------------------------------------------------------------
32 : ! Apply dilution factor to convert surface flux to observed flux
33 : !---------------------------------------------------------------------------
34 0 : subroutine dilute_flux(surface_flux, R, d, calibrated_flux)
35 : real(dp), intent(in) :: surface_flux(:)
36 : real(dp), intent(in) :: R, d ! R = stellar radius, d = distance (both in the same units, e.g., cm)
37 : real(dp), intent(out) :: calibrated_flux(:)
38 :
39 : ! Check that the output array has the same size as the input
40 0 : if (size(calibrated_flux) /= size(surface_flux)) then
41 0 : print *, "Error in dilute_flux: Output array must have the same size as input array."
42 0 : stop 1
43 : end if
44 :
45 : ! Apply the dilution factor (R/d)^2 to each element
46 0 : calibrated_flux = surface_flux*((R/d)**2)
47 0 : end subroutine dilute_flux
48 :
49 : !###########################################################
50 : !## MATHS
51 : !###########################################################
52 :
53 : !****************************
54 : !Trapezoidal and Simpson Integration For Flux Calculation
55 : !****************************
56 :
57 0 : subroutine trapezoidal_integration(x, y, result)
58 : real(dp), dimension(:), intent(in) :: x, y
59 : real(dp), intent(out) :: result
60 :
61 : integer :: i, n
62 0 : REAL :: sum
63 :
64 0 : n = size(x)
65 0 : sum = 0.0
66 :
67 : ! Validate input sizes
68 0 : if (size(x) /= size(y)) then
69 0 : print *, "Error: x and y arrays must have the same size."
70 0 : stop
71 : end if
72 :
73 0 : if (size(x) < 2) then
74 0 : print *, "Error: x and y arrays must have at least 2 elements."
75 0 : stop
76 : end if
77 :
78 : ! Perform trapezoidal integration
79 0 : do i = 1, n - 1
80 0 : sum = sum + 0.5*(x(i + 1) - x(i))*(y(i + 1) + y(i))
81 : end do
82 :
83 0 : result = sum
84 0 : end subroutine trapezoidal_integration
85 :
86 0 : subroutine simpson_integration(x, y, result)
87 : real(dp), dimension(:), intent(in) :: x, y
88 : real(dp), intent(out) :: result
89 :
90 : integer :: i, n
91 0 : real(dp) :: sum, h1, h2, f1, f2, f0
92 :
93 0 : n = size(x)
94 0 : sum = 0.0_dp
95 :
96 : ! Validate input sizes
97 0 : if (size(x) /= size(y)) then
98 0 : print *, "Error: x and y arrays must have the same size."
99 0 : stop
100 : end if
101 :
102 0 : if (size(x) < 2) then
103 0 : print *, "Error: x and y arrays must have at least 2 elements."
104 0 : stop
105 : end if
106 :
107 : ! Perform adaptive Simpson’s rule
108 0 : do i = 1, n - 2, 2
109 0 : h1 = x(i + 1) - x(i) ! Step size for first interval
110 0 : h2 = x(i + 2) - x(i + 1) ! Step size for second interval
111 :
112 0 : f0 = y(i)
113 0 : f1 = y(i + 1)
114 0 : f2 = y(i + 2)
115 :
116 : ! Simpson's rule: (h/3) * (f0 + 4f1 + f2)
117 0 : sum = sum + (h1 + h2)/6.0_dp*(f0 + 4.0_dp*f1 + f2)
118 : end do
119 :
120 : ! Handle the case where n is odd (last interval)
121 0 : if (MOD(n, 2) == 0) then
122 0 : sum = sum + 0.5_dp*(x(n) - x(n - 1))*(y(n) + y(n - 1))
123 : end if
124 :
125 0 : result = sum
126 0 : end subroutine simpson_integration
127 :
128 0 : subroutine romberg_integration(x, y, result)
129 : real(dp), dimension(:), intent(in) :: x, y
130 : real(dp), intent(out) :: result
131 :
132 : integer :: i, j, k, n, m
133 0 : real(dp), dimension(:), allocatable :: R
134 0 : real(dp) :: h, sum, factor
135 :
136 0 : n = size(x)
137 0 : m = int(log(real(n, DP))/log(2.0_dp)) + 1 ! Number of refinement levels
138 :
139 : ! Validate input sizes
140 0 : if (size(x) /= size(y)) then
141 0 : print *, "Error: x and y arrays must have the same size."
142 0 : stop
143 : end if
144 :
145 0 : if (n < 2) then
146 0 : print *, "Error: x and y arrays must have at least 2 elements."
147 0 : stop
148 : end if
149 :
150 0 : allocate (R(m))
151 :
152 : ! Compute initial trapezoidal rule estimate
153 0 : h = x(n) - x(1)
154 0 : R(1) = 0.5_dp*h*(y(1) + y(n))
155 :
156 : ! Refinement using Romberg's method
157 0 : do j = 2, m
158 0 : sum = 0.0_dp
159 0 : do i = 1, 2**(j - 2)
160 0 : sum = sum + y(1 + (2*i - 1)*(n - 1)/(2**(j - 1)))
161 : end do
162 :
163 0 : h = h/2.0_dp
164 0 : R(j) = 0.5_dp*R(j - 1) + h*sum
165 :
166 : ! Richardson extrapolation
167 0 : factor = 4.0_dp
168 0 : do k = j, 2, -1
169 0 : R(k - 1) = (factor*R(k) - R(k - 1))/(factor - 1.0_dp)
170 0 : factor = factor*4.0_dp
171 : end do
172 : end do
173 :
174 0 : result = R(1)
175 0 : end subroutine romberg_integration
176 :
177 : !-----------------------------------------------------------------------
178 : ! File I/O functions
179 : !-----------------------------------------------------------------------
180 :
181 : !****************************
182 : ! Load Vega SED for Zero Point Calculation
183 : !****************************
184 0 : subroutine load_vega_sed(filepath, wavelengths, flux)
185 : character(len=*), intent(in) :: filepath
186 : real(dp), dimension(:), allocatable, intent(out) :: wavelengths, flux
187 : character(len=512) :: line
188 : integer :: unit, n_rows, status, i
189 0 : real(dp) :: temp_wave, temp_flux
190 :
191 0 : unit = 20
192 0 : open (unit, file=trim(filepath), status='OLD', action='READ', iostat=status)
193 0 : if (status /= 0) then
194 0 : print *, "Error: Could not open Vega SED file ", trim(filepath)
195 0 : stop
196 : end if
197 :
198 : ! Skip header line
199 0 : read (unit, '(A)', iostat=status) line
200 0 : if (status /= 0) then
201 0 : print *, "Error: Could not read header from Vega SED file ", trim(filepath)
202 0 : stop
203 : end if
204 :
205 : ! Count the number of data lines
206 : n_rows = 0
207 0 : do
208 0 : read (unit, '(A)', iostat=status) line
209 0 : if (status /= 0) exit
210 0 : n_rows = n_rows + 1
211 : end do
212 :
213 0 : rewind (unit)
214 0 : read (unit, '(A)', iostat=status) line ! Skip header again
215 :
216 0 : allocate (wavelengths(n_rows))
217 0 : allocate (flux(n_rows))
218 :
219 0 : i = 0
220 0 : do
221 0 : read (unit, *, iostat=status) temp_wave, temp_flux ! Ignore any extra columns
222 0 : if (status /= 0) exit
223 0 : i = i + 1
224 0 : wavelengths(i) = temp_wave
225 0 : flux(i) = temp_flux
226 : end do
227 :
228 0 : close (unit)
229 0 : end subroutine load_vega_sed
230 :
231 : !****************************
232 : ! Load Filter File
233 : !****************************
234 0 : subroutine load_filter(directory, filter_wavelengths, filter_trans)
235 : character(len=*), intent(in) :: directory
236 : real(dp), dimension(:), allocatable, intent(out) :: filter_wavelengths, filter_trans
237 :
238 : character(len=512) :: line
239 : integer :: unit, n_rows, status, i
240 0 : REAL :: temp_wavelength, temp_trans
241 :
242 : ! Open the file
243 0 : unit = 20
244 0 : open (unit, file=trim(directory), status='OLD', action='READ', iostat=status)
245 0 : if (status /= 0) then
246 0 : print *, "Error: Could not open file ", trim(directory)
247 0 : stop
248 : end if
249 :
250 : ! Skip header line
251 0 : read (unit, '(A)', iostat=status) line
252 0 : if (status /= 0) then
253 0 : print *, "Error: Could not read the file", trim(directory)
254 0 : stop
255 : end if
256 :
257 : ! Count rows in the file
258 : n_rows = 0
259 0 : do
260 0 : read (unit, '(A)', iostat=status) line
261 0 : if (status /= 0) exit
262 0 : n_rows = n_rows + 1
263 : end do
264 :
265 : ! Allocate arrays
266 0 : allocate (filter_wavelengths(n_rows))
267 0 : allocate (filter_trans(n_rows))
268 :
269 : ! Rewind to the first non-comment line
270 0 : rewind (unit)
271 : do
272 0 : read (unit, '(A)', iostat=status) line
273 0 : if (status /= 0) then
274 0 : print *, "Error: Could not rewind file", trim(directory)
275 0 : stop
276 : end if
277 0 : if (line(1:1) /= "#") exit
278 : end do
279 :
280 : ! Read and parse data
281 : i = 0
282 0 : do
283 0 : read (unit, *, iostat=status) temp_wavelength, temp_trans
284 0 : if (status /= 0) exit
285 0 : i = i + 1
286 :
287 0 : filter_wavelengths(i) = temp_wavelength
288 0 : filter_trans(i) = temp_trans
289 : end do
290 :
291 0 : close (unit)
292 0 : end subroutine load_filter
293 :
294 : !****************************
295 : ! Load Lookup Table For Identifying Stellar Atmosphere Models
296 : !****************************
297 0 : subroutine load_lookup_table(lookup_file, lookup_table, out_file_names, &
298 : out_logg, out_meta, out_teff)
299 :
300 : character(len=*), intent(in) :: lookup_file
301 : REAL, dimension(:, :), allocatable, intent(out) :: lookup_table
302 : character(len=100), allocatable, intent(inout) :: out_file_names(:)
303 : real(dp), allocatable, intent(inout) :: out_logg(:), out_meta(:), out_teff(:)
304 :
305 : integer :: i, n_rows, status, unit
306 : character(len=512) :: line
307 : character(len=*), parameter :: delimiter = ","
308 0 : character(len=100), allocatable :: columns(:), headers(:)
309 : integer :: logg_col, meta_col, teff_col
310 :
311 : ! Open the file
312 0 : unit = 10
313 0 : open (unit, file=lookup_file, status='old', action='read', iostat=status)
314 0 : if (status /= 0) then
315 0 : print *, "Error: Could not open file", lookup_file
316 0 : stop
317 : end if
318 :
319 : ! Read header line
320 0 : read (unit, '(A)', iostat=status) line
321 0 : if (status /= 0) then
322 0 : print *, "Error: Could not read header line"
323 0 : stop
324 : end if
325 :
326 : call split_line(line, delimiter, headers)
327 :
328 : ! Determine column indices for logg, meta, and teff
329 0 : logg_col = get_column_index(headers, "logg")
330 0 : teff_col = get_column_index(headers, "teff")
331 :
332 0 : meta_col = get_column_index(headers, "meta")
333 0 : if (meta_col < 0) then
334 0 : meta_col = get_column_index(headers, "feh")
335 : end if
336 :
337 0 : n_rows = 0
338 0 : do
339 0 : read (unit, '(A)', iostat=status) line
340 0 : if (status /= 0) exit
341 0 : n_rows = n_rows + 1
342 : end do
343 0 : rewind (unit)
344 :
345 : ! Skip header
346 0 : read (unit, '(A)', iostat=status) line
347 :
348 : ! Allocate output arrays
349 0 : allocate (out_file_names(n_rows))
350 0 : allocate (out_logg(n_rows), out_meta(n_rows), out_teff(n_rows))
351 :
352 : ! Read and parse the file
353 0 : i = 0
354 0 : do
355 0 : read (unit, '(A)', iostat=status) line
356 0 : if (status /= 0) exit
357 0 : i = i + 1
358 :
359 0 : call split_line(line, delimiter, columns)
360 :
361 : ! Populate arrays
362 0 : out_file_names(i) = columns(1)
363 :
364 0 : if (logg_col > 0) then
365 0 : if (columns(logg_col) /= "") then
366 0 : read (columns(logg_col), *) out_logg(i)
367 : else
368 0 : out_logg(i) = -999.0
369 : end if
370 : else
371 0 : out_logg(i) = -999.0
372 : end if
373 :
374 0 : if (meta_col > 0) then
375 0 : if (columns(meta_col) /= "") then
376 0 : read (columns(meta_col), *) out_meta(i)
377 : else
378 0 : out_meta(i) = 0.0
379 : end if
380 : else
381 0 : out_meta(i) = 0.0
382 : end if
383 :
384 0 : if (teff_col > 0) then
385 0 : if (columns(teff_col) /= "") then
386 0 : read (columns(teff_col), *) out_teff(i)
387 : else
388 0 : out_teff(i) = 0.0
389 : end if
390 : else
391 0 : out_teff(i) = 0.0
392 : end if
393 : end do
394 :
395 0 : close (unit)
396 :
397 : contains
398 :
399 0 : function get_column_index(headers, target) result(index)
400 : character(len=100), intent(in) :: headers(:)
401 : character(len=*), intent(in) :: target
402 : integer :: index, i
403 : character(len=100) :: clean_header, clean_target
404 :
405 0 : index = -1
406 0 : clean_target = trim(adjustl(target)) ! Clean the target string
407 :
408 0 : do i = 1, size(headers)
409 0 : clean_header = trim(adjustl(headers(i))) ! Clean each header
410 0 : if (clean_header == clean_target) then
411 : index = i
412 : exit
413 : end if
414 : end do
415 0 : end function get_column_index
416 :
417 0 : subroutine split_line(line, delimiter, tokens)
418 : character(len=*), intent(in) :: line, delimiter
419 : character(len=100), allocatable, intent(out) :: tokens(:)
420 : integer :: num_tokens, pos, start, len_delim
421 :
422 0 : len_delim = len_trim(delimiter)
423 0 : start = 1
424 0 : num_tokens = 0
425 0 : if (allocated(tokens)) deallocate (tokens)
426 :
427 0 : do
428 0 : pos = index(line(start:), delimiter)
429 :
430 0 : if (pos == 0) exit
431 0 : num_tokens = num_tokens + 1
432 0 : call append_token(tokens, line(start:start + pos - 2))
433 0 : start = start + pos + len_delim - 1
434 : end do
435 :
436 0 : num_tokens = num_tokens + 1
437 0 : call append_token(tokens, line(start:))
438 0 : end subroutine split_line
439 :
440 0 : subroutine append_token(tokens, token)
441 : character(len=*), intent(in) :: token
442 : character(len=100), allocatable, intent(inout) :: tokens(:)
443 0 : character(len=100), allocatable :: temp(:)
444 : integer :: n
445 :
446 0 : if (.not. allocated(tokens)) then
447 0 : allocate (tokens(1))
448 0 : tokens(1) = token
449 : else
450 0 : n = size(tokens)
451 0 : allocate (temp(n))
452 0 : temp = tokens ! Backup the current tokens
453 0 : deallocate (tokens) ! Deallocate the old array
454 0 : allocate (tokens(n + 1)) ! Allocate with one extra space
455 0 : tokens(1:n) = temp ! Restore old tokens
456 0 : tokens(n + 1) = token ! Add the new token
457 0 : deallocate (temp) ! unsure if this is till needed.
458 : end if
459 0 : end subroutine append_token
460 :
461 : end subroutine load_lookup_table
462 :
463 0 : subroutine load_sed(directory, index, wavelengths, flux)
464 : character(len=*), intent(in) :: directory
465 : integer, intent(in) :: index
466 : real(dp), dimension(:), allocatable, intent(out) :: wavelengths, flux
467 :
468 : character(len=512) :: line
469 : integer :: unit, n_rows, status, i
470 0 : real(dp) :: temp_wavelength, temp_flux
471 :
472 : ! Open the file
473 0 : unit = 20
474 0 : open (unit, file=trim(directory), status='OLD', action='READ', iostat=status)
475 0 : if (status /= 0) then
476 0 : print *, "Error: Could not open file ", trim(directory)
477 0 : stop
478 : end if
479 :
480 : ! Skip header lines
481 : do
482 0 : read (unit, '(A)', iostat=status) line
483 0 : if (status /= 0) then
484 0 : print *, "Error: Could not read the file", trim(directory)
485 0 : stop
486 : end if
487 0 : if (line(1:1) /= "#") exit
488 : end do
489 :
490 : ! Count rows in the file
491 : n_rows = 0
492 0 : do
493 0 : read (unit, '(A)', iostat=status) line
494 0 : if (status /= 0) exit
495 0 : n_rows = n_rows + 1
496 : end do
497 :
498 : ! Allocate arrays
499 0 : allocate (wavelengths(n_rows))
500 0 : allocate (flux(n_rows))
501 :
502 : ! Rewind to the first non-comment line
503 0 : rewind (unit)
504 : do
505 0 : read (unit, '(A)', iostat=status) line
506 0 : if (status /= 0) then
507 0 : print *, "Error: Could not rewind file", trim(directory)
508 0 : stop
509 : end if
510 0 : if (line(1:1) /= "#") exit
511 : end do
512 :
513 : ! Read and parse data
514 : i = 0
515 0 : do
516 0 : read (unit, *, iostat=status) temp_wavelength, temp_flux
517 0 : if (status /= 0) exit
518 0 : i = i + 1
519 : ! Convert f_lambda to f_nu
520 0 : wavelengths(i) = temp_wavelength
521 0 : flux(i) = temp_flux
522 : end do
523 :
524 0 : close (unit)
525 :
526 0 : end subroutine load_sed
527 :
528 : !-----------------------------------------------------------------------
529 : ! Helper function for file names
530 : !-----------------------------------------------------------------------
531 :
532 0 : function remove_dat(path) result(base)
533 : ! Extracts the portion of the string before the first dot
534 : character(len=*), intent(in) :: path
535 : character(len=strlen) :: base
536 : integer :: first_dot
537 :
538 : ! Find the position of the first dot
539 0 : first_dot = 0
540 0 : do while (first_dot < len_trim(path) .and. path(first_dot + 1:first_dot + 1) /= '.')
541 0 : first_dot = first_dot + 1
542 : end do
543 :
544 : ! Check if a dot was found
545 0 : if (first_dot < len_trim(path)) then
546 : ! Extract the part before the dot
547 0 : base = path(:first_dot)
548 : else
549 : ! No dot found, return the input string
550 0 : base = path
551 : end if
552 0 : end function remove_dat
553 :
554 1 : function basename(path) result(name)
555 : character(len=*), intent(in) :: path
556 : character(len=strlen) :: name
557 : integer :: i
558 1 : if (len_trim(path) == 0) then
559 0 : name = ''
560 0 : return
561 : end if
562 1 : i = index(path, '/', back=.true.)
563 1 : name = path(i + 1:)
564 : end function basename
565 :
566 1 : subroutine read_strings_from_file(colors_settings, strings, n, ierr)
567 : character(len=512) :: filename
568 : character(len=100), allocatable, intent(out) :: strings(:)
569 : integer, intent(out) :: n, ierr
570 : integer :: unit, i, status
571 : character(len=100) :: line
572 : type(Colors_General_Info), pointer :: colors_settings
573 :
574 1 : ierr = 0
575 :
576 : filename = trim(mesa_dir)//trim(colors_settings%instrument)//"/"// &
577 1 : trim(basename(colors_settings%instrument))
578 1 : n = 0
579 1 : unit = 10
580 1 : open (unit, file=filename, status='old', action='read', iostat=status)
581 1 : if (status /= 0) then
582 0 : ierr = -1
583 0 : print *, "Error: Could not open file", filename
584 0 : stop
585 : end if
586 :
587 6 : do
588 7 : read (unit, '(A)', iostat=status) line
589 7 : if (status /= 0) exit
590 6 : n = n + 1
591 : end do
592 1 : rewind (unit)
593 1 : if (allocated(strings)) deallocate (strings)
594 1 : allocate (strings(n))
595 7 : do i = 1, n
596 7 : read (unit, '(A)') strings(i)
597 : end do
598 1 : close (unit)
599 1 : end subroutine read_strings_from_file
600 :
601 : end module colors_utils
|
