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 bolometric
      21              : 
      22              :    use const_def, only: dp
      23              :    use colors_utils, only: romberg_integration, load_lookup_table
      24              :    use hermite_interp, only: construct_sed_hermite
      25              : 
      26              :    implicit none
      27              : 
      28              :    private
      29              :    public :: calculate_bolometric
      30              : 
      31              : contains
      32              : 
      33              :    !****************************
      34              :    ! Calculate Bolometric Photometry Using Multiple SEDs
      35              :    !****************************
      36            0 :    subroutine calculate_bolometric(teff, log_g, metallicity, R, d, bolometric_magnitude, &
      37              :                                    bolometric_flux, wavelengths, fluxes, sed_filepath)
      38              :       real(dp), intent(in) :: teff, log_g, metallicity, R, d
      39              :       character(len=*), intent(in) :: sed_filepath
      40              :       real(dp), intent(out) :: bolometric_magnitude, bolometric_flux
      41              :       real(dp), dimension(:), allocatable, intent(out) :: wavelengths, fluxes
      42              : 
      43            0 :       real(dp), allocatable :: lu_logg(:), lu_meta(:), lu_teff(:)
      44            0 :       character(len=100), allocatable :: file_names(:)
      45            0 :       REAL, dimension(:, :), allocatable :: lookup_table
      46              :       character(len=256) :: lookup_file
      47              : 
      48            0 :       lookup_file = trim(sed_filepath)//'/lookup_table.csv'
      49              : 
      50              :       call load_lookup_table(lookup_file, lookup_table, file_names, lu_logg, lu_meta, lu_teff)
      51              : 
      52              :       call construct_sed_hermite(teff, log_g, metallicity, R, d, file_names, &
      53              :                                  lu_teff, lu_logg, lu_meta, sed_filepath, &
      54            0 :                                  wavelengths, fluxes)
      55              : 
      56              :       ! Calculate bolometric flux and magnitude
      57            0 :       call calculate_bolometric_phot(wavelengths, fluxes, bolometric_magnitude, bolometric_flux)
      58            0 :    end subroutine calculate_bolometric
      59              : 
      60              :    !****************************
      61              :    ! Calculate Bolometric Magnitude and Flux
      62              :    !****************************
      63            0 :    subroutine calculate_bolometric_phot(wavelengths, fluxes, bolometric_magnitude, bolometric_flux)
      64              :       real(dp), dimension(:), intent(inout) :: wavelengths, fluxes
      65              :       real(dp), intent(out) :: bolometric_magnitude, bolometric_flux
      66              :       integer :: i
      67              : 
      68              :       ! Validate inputs and replace invalid wavelengths with 0
      69            0 :       do i = 1, size(wavelengths) - 1
      70            0 :          if (wavelengths(i) <= 0.0d0 .or. fluxes(i) < 0.0d0) then
      71            0 :             fluxes(i) = 0.0d0  ! Replace invalid wavelength with 0
      72              :          end if
      73              :       end do
      74              : 
      75              :       ! Call Romberg integration
      76            0 :       call romberg_integration(wavelengths, fluxes, bolometric_flux)
      77              : 
      78              :       ! Validate integration result
      79            0 :       if (bolometric_flux <= 0.0d0) then
      80            0 :          print *, "Error: Flux integration resulted in non-positive value."
      81            0 :          bolometric_magnitude = 99.0d0
      82            0 :          return
      83              :       end if
      84              : 
      85              :       ! Calculate bolometric magnitude
      86              :       if (bolometric_flux <= 0.0d0) then
      87              :          print *, "Error: Flux integration resulted in non-positive value."
      88              :          bolometric_magnitude = 99.0d0
      89              :          return
      90            0 :       else if (bolometric_flux < 1.0d-10) then
      91            0 :          print *, "Warning: Flux value is very small, precision might be affected."
      92              :       end if
      93              : 
      94            0 :       bolometric_magnitude = flux_to_magnitude(bolometric_flux)
      95              :    end subroutine calculate_bolometric_phot
      96              : 
      97              :    !****************************
      98              :    ! Convert Flux to Magnitude
      99              :    !****************************
     100            0 :    real(dp) function flux_to_magnitude(flux)
     101              :       real(dp), intent(in) :: flux
     102            0 :       if (flux <= 0.0d0) then
     103            0 :          print *, "Error: Flux must be positive to calculate magnitude."
     104            0 :          flux_to_magnitude = 99.0d0  ! Return an error value
     105              :       else
     106            0 :          flux_to_magnitude = -2.5d0*log10(flux)
     107              :       end if
     108            0 :    end function flux_to_magnitude
     109              : 
     110              : end module bolometric