Updates and fixes before release

.github/workflows/build.yml

@@ -14,7 +14,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        sdk: ["24.7.1", "25.12.1", "26.3.2"]  # pick 2 or 3 most recent.
+        sdk: ["25.12.1", "26.3.2"]  # pick 2 or 3 most recent.
         os: [ubuntu-latest, macos-latest]
     runs-on: ${{ matrix.os }}
     steps:

colors/README.rst

@@ -135,6 +135,27 @@ Plain two-column text files, one row per wavelength point:
 
 No header is required. All SED files within a given atmosphere grid must share the same wavelength grid; if they do not (e.g. BT-Settl), the stencil loader will interpolate non-conforming files onto the canonical wavelength grid of the first file loaded.
 
+The colors data ships with the Kurucz2003 models for solar alpha and alpha = 0.4.
+
+Spectral Grid Variants
+~~~~~~~~~~~~~~~~~~~~~~
+
++-------------------------+--------+-----------+----------------------+-------------------+----------------+------------------------+
+| Name                    | [α/Fe] | N Spectra | Teff (K)             | logg              | [M/H]          | Wavelength (Å)         |
++=========================+========+===========+======================+===================+================+========================+
+| Kurucz2003all           | 0.0    | 3808      | 3500-50000 (76 pt)   | 0.00-5.00 (11 pt) | -2.50-0.50 (8 pt) | 147-1600000 (1199 pt) |
++-------------------------+--------+-----------+----------------------+-------------------+----------------+------------------------+
+| Kurucz2003all__alpha_04 | 0.4    | 4284      | 3500-50000 (76 pt)   | 0.00-5.00 (11 pt) | -4.00-0.50 (9 pt) | 147-1600000 (1199 pt) |
++-------------------------+--------+-----------+----------------------+-------------------+----------------+------------------------+
+
+**[α/Fe]** is the alpha-element enhancement - the abundance of O, Ne, Mg, Si, S, Ar, Ca, and Ti relative to iron compared to solar.
+A value of 0.4 means those elements are enhanced by 0.4 dex above solar relative to iron.
+The alpha-enhanced grid goes down to lower metallicity (-4.00 vs -2.50) because alpha-enriched stars are usually old, metal-poor Pop II stars.
+
+`Castelli & Kurucz 2003 <https://arxiv.org/abs/astro-ph/0405087>`__
+
+`SVO model <https://svo2.cab.inta-csic.es/theory/newov2/index.php?models=Kurucz2003all>`__
+
 flux_cube.bin
 -------------
 

colors/defaults/colors.defaults

@@ -1,43 +1,109 @@
 ! ``colors`` module controls
 ! ==========================
 
-   ! The MESA/colors parameters are given default values here.
-
 
 ! Colors User Parameters
 ! ----------------------
 
    ! ``use_colors``
    ! ~~~~~~~~~~~~~~
+   ! Set to .true. to enable bolometric and synthetic photometry output.
+
    ! ``instrument``
    ! ~~~~~~~~~~~~~~
-   ! ``vega_sed``
-   ! ~~~~~~~~~~~~
+   ! Path to the filter instrument directory (structured as facility/instrument).
+   ! The directory must contain an index file named after the instrument and
+   ! one .dat transmission curve per filter. Each filter becomes a history column.
+
    ! ``stellar_atm``
    ! ~~~~~~~~~~~~~~~
+   ! Path to the directory containing the stellar atmosphere lookup table and
+   ! spectra used by the colors module. The default
+   ! uses the ``Kurucz2003all`` ATLAS9 atmosphere grid shipped with MESA.
+   ! Must contain: lookup_table.csv, SED files, and optionally flux_cube.bin.
+   !Castelli & Kurucz 2003 -- https://arxiv.org/abs/astro-ph/0405087
+   !SVO model -- https://svo2.cab.inta-csic.es/theory/newov2/index.php?models=Kurucz2003all
+
+
+
+   ! ``vega_sed``
+   ! ~~~~~~~~~~~~
+   ! Path to the Vega reference SED file. Required when mag_system = 'Vega'.
+   ! Used to compute photometric zero-points for all filters.
+
    ! ``distance``
    ! ~~~~~~~~~~~~
+   ! Distance to the star in cm. Determines whether output magnitudes are
+   ! absolute (default: 10 pc = 3.0857e19 cm) or apparent (set to source distance).
+
    ! ``make_csv``
    ! ~~~~~~~~~~~~
+   ! If .true., exports the full SED as a CSV file at every profile interval.
+   ! Files are written to colors_results_directory.
+
    ! ``colors_results_directory``
    ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+   ! Directory where CSV outputs (make_csv, sed_per_model) are saved.
+
    ! ``mag_system``
    ! ~~~~~~~~~~~~~~
+   ! Photometric zero-point system. Options: 'Vega', 'AB', 'ST'.
+   !   Vega : zero-point set so Vega has magnitude 0 in all bands.
+   !   AB   : zero-point based on flat f_nu = 3631 Jy.
+   !   ST   : zero-point based on flat f_lambda per unit wavelength.
+
+   ! ``sed_per_model``
+   ! ~~~~~~~~~~~~~~~~~
+   ! Requires make_csv = .true. If .true., each SED output file is suffixed with
+   ! the model number, preserving the full SED history rather than overwriting.
+   ! WARNING: can produce very large numbers of files. Ensure adequate storage.
 
-   ! If ``use_colors`` is true, the colors module is turned on, which will calculate
-   ! bolometric and synthetic magnitudes by interpolating stellar atmosphere model grids and convolving with photometric filter transmission curves.
-   ! Vega SED for Vega photometric system is used for photometric zero points.
-   ! Stellar distance is given in cm.
+
+   ! If ``use_colors`` is true, the colors module computes bolometric and synthetic
+   ! magnitudes by interpolating stellar atmosphere model grids and convolving with
+   ! photometric filter transmission curves. Output is appended to history.data.
    ! ::
 
       use_colors = .false.
-      instrument = '/data/colors_data/filters/Generic/Johnson'
-      stellar_atm = '/data/colors_data/stellar_models/Kurucz2003all/'
-      distance = 3.0857d19
+      instrument = 'data/colors_data/filters/Generic/Johnson'
+      stellar_atm = 'data/colors_data/stellar_models/Kurucz2003all/'
+      vega_sed = 'data/colors_data/stellar_models/vega_flam.csv'
+      distance = 3.0857d19     ! 10 parsecs in cm -> Absolute Magnitudes
       make_csv = .false.
       colors_results_directory = 'SED'
       mag_system = 'Vega'
-      vega_sed = '/data/colors_data/stellar_models/vega_flam.csv'
+      sed_per_model = .false.
+
+
+   ! ``z_over_x_ref``
+   ! ~~~~~~~~~~~~~~~~
+
+   ! ``z_over_x_ref`` is the reference metal-to-hydrogen ratio used to map
+   ! the photospheric composition onto the atmosphere-table metallicity axis:
+   ! [Fe/H] = log10((Z/X)/z_over_x_ref).
+
+   ! The default matches the GS98 solar mixture used by the default
+   ! Kurucz2003all ATLAS9 atmosphere grid and equals
+   ! 0.0169/(1 - 0.0169 - 0.2485) = 2.30057173972d-2.
+
+   ! If the photospheric metallicity, expressed as Z/X, falls outside the
+   ! metallicity range available in the tabulated atmosphere lookup table,
+   ! then MESA uses the nearest metallicity in the table for the atmosphere
+   ! interpolation. If the photospheric hydrogen mass fraction or metal mass
+   ! fraction is not positive, then MESA cannot form
+   ! log10((Z/X)/z_over_x_ref). In that case the module still returns colors
+   ! and SEDs, but it does so by using the lowest metallicity in the table.
+   ! This is a fallback safeguard and not a dedicated treatment for H-free
+   ! atmospheres.
+
+   ! For the shipped Kurucz2003all table in this checkout, the
+   ! metallicity range is [Fe/H] = -2.5 to [Fe/H] = 4.0. With the
+   ! default ``z_over_x_ref``, the lower edge corresponds to
+   ! Z/X ~ 7.275d-5.
+
+   ! ::
+
+      z_over_x_ref = 2.30057173972d-2
 
 
 ! Extra inlist controls
@@ -46,14 +112,16 @@
    ! One can split a colors inlist into pieces using the following parameters.
    ! It works recursively, so the extras can read extras too.
 
+   ! ``read_extra_colors_inlist(1..5)``
+   ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+   ! ``extra_colors_inlist_name(1..5)``
+   ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-         ! ``read_extra_colors_inlist(1..5)``
-         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         ! ``extra_colors_inlist_name(1..5)``
-         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-         ! If ``read_extra_colors_inlist(i)`` is true, then read ``&colors`` from the file ``extra_colors_inlist_name(i)``.
-         ! ::
+   ! If ``read_extra_colors_inlist(i)`` is true, then read ``&colors``
+   ! from the file ``extra_colors_inlist_name(i)``.
+   !
+
+   ! ::
 
       read_extra_colors_inlist(:) = .false.
       extra_colors_inlist_name(:) = 'undefined'

colors/private/bolometric.f90

@@ -1,6 +1,6 @@
 ! ***********************************************************************
 !
-!   Copyright (C) 2025  Niall Miller & The MESA Team
+!   Copyright (C) 2026  Niall Miller & The MESA Team
 !
 !   This program is free software: you can redistribute it and/or modify
 !   it under the terms of the GNU Lesser General Public License
@@ -19,7 +19,9 @@
 
 module bolometric
 
-   use const_def, only: dp, boltz_sigma
+   use const_def, only: dp
+   use colors_def, only: Colors_General_Info
+   use colors_utils, only: simpson_integration
    use hermite_interp, only: construct_sed_hermite
    use linear_interp, only: construct_sed_linear
    use knn_interp, only: construct_sed_knn
@@ -31,64 +33,81 @@ module bolometric
 
 contains
 
-subroutine calculate_bolometric(teff, log_g, metallicity, R, d, bolometric_magnitude, &
-                                   bolometric_flux, wavelengths, fluxes, sed_filepath, interpolation_radius, &
-                                   lu_file_names, lu_teff, lu_logg, lu_meta)
+   ! rq carries cached lookup table data and the preloaded flux cube (if available)
+   subroutine calculate_bolometric(rq, teff, log_g, metallicity, R, d, bolometric_magnitude, &
+                                   bolometric_flux, wavelengths, fluxes, sed_filepath, interpolation_radius)
+      type(Colors_General_Info), intent(inout) :: rq
       real(dp), intent(in) :: teff, log_g, metallicity, R, d
       character(len=*), intent(in) :: sed_filepath
       real(dp), intent(out) :: bolometric_magnitude, bolometric_flux, interpolation_radius
       real(dp), dimension(:), allocatable, intent(out) :: wavelengths, fluxes
-      character(len=100), intent(in) :: lu_file_names(:)
-      real(dp), intent(in) :: lu_teff(:), lu_logg(:), lu_meta(:)
+
       character(len=32) :: interpolation_method
 
-      interpolation_method = 'Hermite'
+      interpolation_method = 'Hermite'   ! or 'Linear' / 'KNN' later
 
+      ! how far (teff, log_g, metallicity) is from the nearest grid point
       interpolation_radius = compute_interp_radius(teff, log_g, metallicity, &
-                                                   lu_teff, lu_logg, lu_meta)
+                                                   rq%lu_teff, rq%lu_logg, rq%lu_meta)
 
       select case (interpolation_method)
       case ('Hermite', 'hermite', 'HERMITE')
-         call construct_sed_hermite(teff, log_g, metallicity, R, d, lu_file_names, &
-                                    lu_teff, lu_logg, lu_meta, sed_filepath, &
-                                    wavelengths, fluxes)
+         call construct_sed_hermite(rq, teff, log_g, metallicity, R, d, &
+                                    sed_filepath, wavelengths, fluxes)
+
       case ('Linear', 'linear', 'LINEAR')
-         call construct_sed_linear(teff, log_g, metallicity, R, d, lu_file_names, &
-                                   lu_teff, lu_logg, lu_meta, sed_filepath, &
-                                   wavelengths, fluxes)
+         call construct_sed_linear(rq, teff, log_g, metallicity, R, d, &
+                                   sed_filepath, wavelengths, fluxes)
+
       case ('KNN', 'knn', 'Knn')
-         call construct_sed_knn(teff, log_g, metallicity, R, d, lu_file_names, &
-                                lu_teff, lu_logg, lu_meta, sed_filepath, &
-                                wavelengths, fluxes)
+         call construct_sed_knn(rq, teff, log_g, metallicity, R, d, &
+                                sed_filepath, wavelengths, fluxes)
+
       case default
-         call construct_sed_hermite(teff, log_g, metallicity, R, d, lu_file_names, &
-                                    lu_teff, lu_logg, lu_meta, sed_filepath, &
-                                    wavelengths, fluxes)
+         ! fallback: hermite
+         call construct_sed_hermite(rq, teff, log_g, metallicity, R, d, &
+                                    sed_filepath, wavelengths, fluxes)
       end select
 
-      ! bolometric flux is sigma*Teff^4 by definition, diluted by (R/d)^2
-      bolometric_flux = boltz_sigma * teff**4 * (R/d)**2
-      bolometric_magnitude = flux_to_magnitude(bolometric_flux)
-
+      call calculate_bolometric_phot(wavelengths, fluxes, bolometric_magnitude, bolometric_flux)
    end subroutine calculate_bolometric
 
+   subroutine calculate_bolometric_phot(wavelengths, fluxes, bolometric_magnitude, bolometric_flux)
+      real(dp), dimension(:), intent(inout) :: wavelengths, fluxes
+      real(dp), intent(out) :: bolometric_magnitude, bolometric_flux
+      integer :: i
+
+      ! zero out any invalid flux/wavelength values
+      do i = 1, size(wavelengths) - 1
+         if (wavelengths(i) <= 0.0d0 .or. fluxes(i) < 0.0d0) then
+            fluxes(i) = 0.0d0
+         end if
+      end do
+
+      call simpson_integration(wavelengths, fluxes, bolometric_flux)
+
+      if (bolometric_flux <= 0.0d0) then
+         print *, "Error: Flux integration resulted in non-positive value."
+         bolometric_magnitude = 99.0d0
+         return
+      else if (bolometric_flux < 1.0d-10) then
+         print *, "Warning: Flux value is very small, precision might be affected."
+      end if
+
+      bolometric_magnitude = flux_to_magnitude(bolometric_flux)
+   end subroutine calculate_bolometric_phot
 
-   !****************************
-   ! Convert Flux to Magnitude
-   !****************************
    real(dp) function flux_to_magnitude(flux)
       real(dp), intent(in) :: flux
       if (flux <= 0.0d0) then
          print *, "Error: Flux must be positive to calculate magnitude."
          flux_to_magnitude = 99.0d0
       else
-         flux_to_magnitude = -2.5d0 * log10(flux)
+         flux_to_magnitude = -2.5d0*log10(flux)
       end if
    end function flux_to_magnitude
 
-   !--------------------------------------------------------------------
-   ! Scalar metric: distance to nearest grid point in normalized space
-   !--------------------------------------------------------------------
+   ! scalar metric: distance to nearest grid point in normalized space
    real(dp) function compute_interp_radius(teff, log_g, metallicity, &
                                            lu_teff, lu_logg, lu_meta)
 
@@ -105,7 +124,7 @@ real(dp) function compute_interp_radius(teff, log_g, metallicity, &
       logical :: use_teff, use_logg, use_meta
       real(dp), parameter :: eps = 1.0d-12
 
-      ! Detect dummy columns (entire axis is 0 or ±999)
+      ! detect dummy columns (entire axis is 0 or ±999)
       use_teff = .not. (all(lu_teff == 0.0d0) .or. &
                         all(lu_teff == 999.0d0) .or. &
                         all(lu_teff == -999.0d0))
@@ -118,47 +137,46 @@ real(dp) function compute_interp_radius(teff, log_g, metallicity, &
                         all(lu_meta == 999.0d0) .or. &
                         all(lu_meta == -999.0d0))
 
-      ! Compute min/max for valid axes
       if (use_teff) then
          teff_min = minval(lu_teff)
          teff_max = maxval(lu_teff)
          teff_range = max(teff_max - teff_min, eps)
-         norm_teff = (teff - teff_min) / teff_range
+         norm_teff = (teff - teff_min)/teff_range
       end if
 
       if (use_logg) then
          logg_min = minval(lu_logg)
          logg_max = maxval(lu_logg)
          logg_range = max(logg_max - logg_min, eps)
-         norm_logg = (log_g - logg_min) / logg_range
+         norm_logg = (log_g - logg_min)/logg_range
       end if
 
       if (use_meta) then
          meta_min = minval(lu_meta)
          meta_max = maxval(lu_meta)
          meta_range = max(meta_max - meta_min, eps)
-         norm_meta = (metallicity - meta_min) / meta_range
+         norm_meta = (metallicity - meta_min)/meta_range
       end if
 
-      ! Find minimum distance to any grid point
+      ! find minimum distance to any grid point
       d_min = huge(1.0d0)
       n = size(lu_teff)
 
       do i = 1, n
          d = 0.0d0
 
          if (use_teff) then
-            grid_teff = (lu_teff(i) - teff_min) / teff_range
+            grid_teff = (lu_teff(i) - teff_min)/teff_range
             d = d + (norm_teff - grid_teff)**2
          end if
 
          if (use_logg) then
-            grid_logg = (lu_logg(i) - logg_min) / logg_range
+            grid_logg = (lu_logg(i) - logg_min)/logg_range
             d = d + (norm_logg - grid_logg)**2
          end if
 
          if (use_meta) then
-            grid_meta = (lu_meta(i) - meta_min) / meta_range
+            grid_meta = (lu_meta(i) - meta_min)/meta_range
             d = d + (norm_meta - grid_meta)**2
          end if