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