Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2026 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 : ! ***********************************************************************
11 :
12 : module hermite_interp_bounded
13 : use const_def, only: dp
14 : use colors_def, only: Colors_General_Info
15 : use colors_utils, only: simpson_integration
16 : use hermite_interp, only: construct_sed_hermite
17 : use linear_interp, only: construct_sed_linear
18 : use utils_lib, only: is_inf, is_nan, mesa_error
19 :
20 : implicit none
21 :
22 : private
23 : public :: construct_sed_hermite_bounded
24 :
25 : ! Accept negligible Hermite undershoot when both its integrated magnitude
26 : ! and its peak amplitude are below 0.1% of the positive SED.
27 : real(dp), parameter :: negative_flux_fraction_tol = 1.0d-3
28 :
29 : contains
30 :
31 7 : subroutine construct_sed_hermite_bounded(rq, teff, log_g, metallicity, R, d, &
32 : stellar_model_dir, wavelengths, fluxes)
33 : type(Colors_General_Info), intent(inout) :: rq
34 : real(dp), intent(in) :: teff, log_g, metallicity, R, d
35 : character(len=*), intent(in) :: stellar_model_dir
36 : real(dp), allocatable, intent(out) :: wavelengths(:), fluxes(:)
37 :
38 : real(dp) :: integrated_flux
39 : logical :: interp_ok
40 :
41 : call construct_sed_hermite(rq, teff, log_g, metallicity, R, d, &
42 7 : stellar_model_dir, wavelengths, fluxes)
43 :
44 : ! Accept Hermite whenever it is finite and physically usable. Negligible
45 : ! negative undershoots are projected onto the physical flux bound.
46 7 : interp_ok = valid_sed(wavelengths, fluxes, .false.)
47 7 : if (interp_ok) call repair_negative_flux(wavelengths, fluxes, interp_ok)
48 :
49 7 : if (interp_ok) then
50 6 : call simpson_integration(wavelengths, fluxes, integrated_flux)
51 : interp_ok = .not. is_nan(integrated_flux) .and. &
52 6 : .not. is_inf(integrated_flux) .and. integrated_flux > 0.0_dp
53 : end if
54 :
55 7 : if (interp_ok) return
56 :
57 : ! Hermite was unusable. Fall back to the positivity-preserving linear
58 : ! interpolant so a Colors failure does not terminate the MESA run.
59 : call construct_sed_linear(rq, teff, log_g, metallicity, R, d, &
60 1 : stellar_model_dir, wavelengths, fluxes)
61 :
62 1 : interp_ok = valid_sed(wavelengths, fluxes, .true.)
63 1 : if (interp_ok) then
64 1 : call simpson_integration(wavelengths, fluxes, integrated_flux)
65 : interp_ok = .not. is_nan(integrated_flux) .and. &
66 1 : .not. is_inf(integrated_flux) .and. integrated_flux > 0.0_dp
67 : end if
68 :
69 1 : if (.not. interp_ok) then
70 0 : write (*, *) 'colors: Hermite and linear interpolation returned invalid SEDs'
71 0 : call mesa_error(__FILE__, __LINE__)
72 : end if
73 : end subroutine construct_sed_hermite_bounded
74 :
75 7 : subroutine repair_negative_flux(wavelengths, fluxes, repair_ok)
76 : real(dp), intent(in) :: wavelengths(:)
77 : real(dp), intent(inout) :: fluxes(:)
78 : logical, intent(out) :: repair_ok
79 :
80 7 : real(dp), allocatable :: flux_part(:)
81 : real(dp) :: negative_integral, positive_integral
82 : real(dp) :: negative_peak, positive_peak
83 : integer :: n
84 :
85 7 : repair_ok = .false.
86 7 : n = size(fluxes)
87 :
88 7 : if (n < 2 .or. size(wavelengths) /= n) return
89 :
90 : ! Nothing needs repairing. Zero flux is already physically safe.
91 2411 : if (.not. any(fluxes < 0.0_dp)) then
92 12 : repair_ok = any(fluxes > 0.0_dp)
93 3 : return
94 : end if
95 :
96 : ! A useful SED must retain a positive interior.
97 62 : if (.not. any(fluxes > 0.0_dp)) return
98 :
99 9 : allocate (flux_part(n))
100 :
101 2407 : flux_part = max(0.0_dp, -fluxes)
102 3 : call simpson_integration(wavelengths, flux_part, negative_integral)
103 2404 : negative_peak = maxval(flux_part)
104 :
105 2407 : flux_part = max(0.0_dp, fluxes)
106 3 : call simpson_integration(wavelengths, flux_part, positive_integral)
107 2404 : positive_peak = maxval(flux_part)
108 :
109 3 : if (is_nan(negative_integral) .or. is_inf(negative_integral)) return
110 3 : if (is_nan(positive_integral) .or. is_inf(positive_integral)) return
111 3 : if (is_nan(negative_peak) .or. is_inf(negative_peak)) return
112 3 : if (is_nan(positive_peak) .or. is_inf(positive_peak)) return
113 3 : if (positive_integral <= 0.0_dp .or. positive_peak <= 0.0_dp) return
114 :
115 3 : if (negative_integral/positive_integral > negative_flux_fraction_tol) return
116 3 : if (negative_peak/positive_peak > negative_flux_fraction_tol) return
117 :
118 : ! The undershoot is negligible; project only the unphysical samples onto
119 : ! the physical constraint f_lambda >= 0.
120 2404 : where (fluxes < 0.0_dp) fluxes = 0.0_dp
121 :
122 2462 : repair_ok = all(fluxes >= 0.0_dp) .and. any(fluxes > 0.0_dp)
123 7 : end subroutine repair_negative_flux
124 :
125 8 : logical function valid_sed(wavelengths, fluxes, require_positive_flux)
126 : real(dp), intent(in) :: wavelengths(:), fluxes(:)
127 : logical, intent(in) :: require_positive_flux
128 :
129 8 : valid_sed = size(wavelengths) >= 2 .and. size(wavelengths) == size(fluxes)
130 : if (.not. valid_sed) return
131 :
132 : valid_sed = .not. any(is_nan(wavelengths)) .and. .not. any(is_inf(wavelengths)) .and. &
133 19264 : .not. any(is_nan(fluxes)) .and. .not. any(is_inf(fluxes))
134 : if (.not. valid_sed) return
135 :
136 : valid_sed = all(wavelengths > 0.0_dp) .and. &
137 9616 : all(wavelengths(2:) > wavelengths(:size(wavelengths) - 1))
138 8 : if (.not. valid_sed) return
139 :
140 11 : if (require_positive_flux) valid_sed = all(fluxes > 0.0_dp)
141 : end function valid_sed
142 :
143 : end module hermite_interp_bounded
|
