fix fallback paths, add bolometric flux BB comparison test, add changelog entries and known bugs

Author: Debraheem

colors/defaults/colors.defaults

@@ -91,9 +91,10 @@
    ! 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), so it uses the lowest metallicity in the
-   ! table instead. This is a fallback safeguard and not a dedicated treatment
-   ! for H-free atmospheres.
+   ! 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

colors/private/colors_utils.f90

@@ -487,6 +487,9 @@ function basename(path) result(name)
       name = path(i + 1:)
    end function basename
 
+   ! resolve a path relative to mesa_dir unless it is already explicit.
+   ! legacy colors defaults used '/data/...'; if that absolute-looking path
+   ! does not exist on disk, warn and retry relative to mesa_dir.
    function resolve_path(path) result(full_path)
       use const_def, only: mesa_dir
       character(len=*), intent(in) :: path
@@ -511,10 +514,8 @@ function resolve_path(path) result(full_path)
          if (exists) then
             full_path = p
          else
-            !it might be nice to warn the user that their filepath implies absolute
-            !silently correcting this WILL lead to confusion
-            !warning every step seems a smidge over the top though
-            !write (*, *) trim(p), " not found. Trying ", trim(mesa_dir)//trim(p)
+            write (*, '(5a)') 'colors: path ', trim(p), &
+               ' not found; trying ', trim(mesa_dir)//trim(p), ' instead'
             full_path = trim(mesa_dir)//trim(p)
          end if
 
@@ -559,17 +560,19 @@ subroutine read_strings_from_file(colors_settings, strings, n, ierr)
    end subroutine read_strings_from_file
 
    ! load flux cube from binary file into handle at initialization.
-   ! if the file cannot be opened or the large flux_cube array cannot be
-   ! allocated, sets cube_loaded = .false. so the runtime will fall back to
-   ! loading individual SED files via the lookup table.
-   ! grids and wavelengths are always loaded (small); only the 4-D cube
-   ! allocation is treated as the fallback trigger.
+   ! if the cube is missing, malformed, or too large to allocate, keep
+   ! cube_loaded = .false. and fall back to loading individual SED files.
+   ! validate the header and on-disk size before reading the large 4-D payload
+   ! so a stale or truncated cube reports a clear message instead of failing
+   ! later inside the Fortran runtime.
    subroutine load_flux_cube(rq, stellar_model_dir)
+      use iso_fortran_env, only: int64
       type(Colors_General_Info), intent(inout) :: rq
       character(len=*), intent(in) :: stellar_model_dir
 
       character(len=512) :: bin_filename
       integer :: unit, status, n_teff, n_logg, n_meta, n_lambda
+      integer(int64) :: file_bytes, expected_bytes, header_bytes, real_bytes
       real(dp) :: cube_mb
 
       rq%cube_loaded = .false.
@@ -580,16 +583,40 @@ subroutine load_flux_cube(rq, stellar_model_dir)
             access='STREAM', form='UNFORMATTED', iostat=status)
       if (status /= 0) then
          ! no binary cube available -- will use individual SED files
-         write (*, '(a)') 'colors: no flux_cube.bin found; using per-file SED loading'
+         write (*, '(a)') 'colors: no flux_cube.bin found; using slower per-file SED loading'
          return
       end if
 
       read (unit, iostat=status) n_teff, n_logg, n_meta, n_lambda
       if (status /= 0) then
+         write (*, '(a)') 'colors: could not read flux_cube.bin header; falling back to slower per-file SED loading'
          close (unit)
          return
       end if
 
+      if (n_teff <= 0 .or. n_logg <= 0 .or. n_meta <= 0 .or. n_lambda <= 0) then
+         write (*, '(a)') 'colors: invalid flux_cube.bin header; falling back to slower per-file SED loading'
+         close (unit)
+         return
+      end if
+
+      inquire (file=trim(bin_filename), size=file_bytes, iostat=status)
+      if (status == 0) then
+         header_bytes = 4_int64*int(storage_size(n_teff)/8, int64)
+         real_bytes = int(storage_size(1.0_dp)/8, int64)
+         expected_bytes = header_bytes + real_bytes*( &
+            int(n_teff, int64) + int(n_logg, int64) + int(n_meta, int64) + int(n_lambda, int64) + &
+            int(n_teff, int64)*int(n_logg, int64)*int(n_meta, int64)*int(n_lambda, int64))
+         if (file_bytes /= expected_bytes) then
+            write (*, '(a,i0,a,i0,a)') &
+               'colors: flux_cube.bin size mismatch (', file_bytes, &
+               ' bytes, expected ', expected_bytes, &
+               '); falling back to slower per-file SED loading'
+            close (unit)
+            return
+         end if
+      end if
+
       ! attempt the large allocation first -- this is the one that may fail.
       ! doing it before the small grid allocations avoids partial cleanup.
       allocate (rq%cube_flux(n_teff, n_logg, n_meta, n_lambda), stat=status)
@@ -598,7 +625,7 @@ subroutine load_flux_cube(rq, stellar_model_dir)
          cube_mb = real(n_teff, dp)*n_logg*n_meta*n_lambda*8.0_dp/(1024.0_dp**2)
          write (*, '(a,f0.1,a)') &
             'colors: flux cube allocation failed (', cube_mb, &
-            ' MB); falling back to per-file SED loading'
+            ' MB); falling back to slower per-file SED loading'
          close (unit)
          return
       end if
@@ -643,7 +670,7 @@ subroutine load_flux_cube(rq, stellar_model_dir)
 
       ! error cleanup -- deallocate everything that may have been allocated
 900   continue
-      write (*, '(a)') 'colors: error reading flux_cube.bin; falling back to per-file SED loading'
+      write (*, '(a)') 'colors: error reading flux_cube.bin; falling back to slower per-file SED loading'
       if (allocated(rq%cube_flux)) deallocate (rq%cube_flux)
       if (allocated(rq%cube_teff_grid)) deallocate (rq%cube_teff_grid)
       if (allocated(rq%cube_logg_grid)) deallocate (rq%cube_logg_grid)
@@ -1081,4 +1108,4 @@ subroutine interp_linear_internal(x_out, x_in, y_in, y_out)
       end do
 
    end subroutine interp_linear_internal
-end module colors_utils
+end module colors_utils

colors/test/src/test_colors.f90

@@ -44,8 +44,9 @@ program test_colors
       how_many_colors_history_columns, data_for_colors_history_columns, &
       calculate_bolometric
    use colors_def, only: Colors_General_Info
-   use const_def,  only: dp, rsun, mesa_dir
+   use const_def,  only: dp, rsun, boltz_sigma
    use utils_lib,  only: mesa_error
+   use colors_utils, only: resolve_path
 
    implicit none
 
@@ -97,6 +98,11 @@ program test_colors
    ! number of sampled SED points printed for the SED comparison
    integer, parameter :: n_sed_samples = 20
 
+   ! require the integrated SED flux to stay close to sigma*T^4 scaled to
+   ! the test distance. large deficits usually mean the wavelength coverage
+   ! is missing too much UV or IR flux.
+   real(dp), parameter :: bol_flux_rel_tol = 5d-3
+
    character(len=32) :: my_mesa_dir
    integer :: handle, ierr, n_cols, i, j, k
    integer :: model_num
@@ -160,7 +166,7 @@ program test_colors
    ! group 1: representative stellar types
    ! -----------------------------------------------------------------------
 
-   write(*,'(a)') '# group1  system=Vega  grid=Kurucz2003  filters=Johnson'
+   call write_section_header('# Group1  system=Vega  grid=Kurucz2003  filters=Johnson')
    do j = 1, n_cases
       model_num = model_num + 1
       call data_for_colors_history_columns( &
@@ -180,7 +186,7 @@ program test_colors
    ! group 2a: vary [M/H]
    ! -----------------------------------------------------------------------
 
-   write(*,'(a)') '# group2a  vary_MH  Teff=5778  logg=4.44'
+   call write_section_header('# Group2a  vary_MH  Teff=5778  logg=4.44')
    do j = 1, n_meta
       model_num = model_num + 1
       call data_for_colors_history_columns( &
@@ -200,7 +206,7 @@ program test_colors
    ! group 2b: vary log g
    ! -----------------------------------------------------------------------
 
-   write(*,'(a)') '# group2b  vary_logg  Teff=5778  MH=0.0'
+   call write_section_header('# Group2b  vary_logg  Teff=5778  MH=0.0')
    do j = 1, n_logg
       model_num = model_num + 1
       call data_for_colors_history_columns( &
@@ -220,7 +226,7 @@ program test_colors
    ! group 2c: vary Teff
    ! -----------------------------------------------------------------------
 
-   write(*,'(a)') '# group2c  vary_Teff  logg=4.0  MH=0.0'
+   call write_section_header('# Group2c  vary_Teff  logg=4.0  MH=0.0')
    do j = 1, n_teff
       model_num = model_num + 1
       call data_for_colors_history_columns( &
@@ -236,23 +242,32 @@ program test_colors
       end do
    end do
 
+   sed_filepath = trim(resolve_path(cs%stellar_atm))
+
    ! -----------------------------------------------------------------------
-   ! SED comparison: solar case, n_sed_samples evenly spaced wavelength/flux
+   ! group 3: solar SED sample plus wavelength-coverage sanity checks
    ! -----------------------------------------------------------------------
 
-   sed_filepath = trim(mesa_dir)//trim(cs%stellar_atm)
+   call write_section_header('# Group3  SED sample + wavelength_coverage_sanity')
+   write(*,'(a)') '# SED sample  case=solar  Teff=5778  logg=4.44  FeH=0.0  columns=wavelength_AA  flux_erg_s_cm2_AA'
    call calculate_bolometric( &
       cs, test_teff(1), test_logg(1), test_meta(1), test_R(1), d_10pc, &
       bol_mag, bol_flux, wavelengths, fluxes, sed_filepath, interp_rad)
 
    n_wav = size(wavelengths)
    stride = max(1, n_wav / n_sed_samples)
 
-   write(*,'(a)') '# SED sample  case=solar  columns=wavelength_AA  flux_erg_s_cm2_AA'
    do i = 1, n_wav, stride
       write(*,'(1pe23.13, 1x, 1pe23.13)') wavelengths(i), fluxes(i)
    end do
 
+   write(*,'(a)') ''
+   write(*,'(a, 1pe10.2)') '# wavelength_coverage_sanity  logg=4.0  FeH=0.0  rel_tol=', bol_flux_rel_tol
+   do j = 1, n_teff
+      call check_bolometric_coverage( &
+         cs, sed_filepath, 'Teff', sweep_teff(j), 4.0d0, 0.0d0, rsun, bol_flux_rel_tol)
+   end do
+
    ! -----------------------------------------------------------------------
    ! cleanup
    ! -----------------------------------------------------------------------
@@ -266,4 +281,46 @@ program test_colors
 
    write(*,*) 'test_colors: passed'
 
+contains
+
+   subroutine write_section_header(title)
+      character(len=*), intent(in) :: title
+
+      write(*,'(a)') ''
+      write(*,'(a)') trim(title)
+      write(*,'(a)') ''
+   end subroutine write_section_header
+
+   subroutine check_bolometric_coverage(rq, sed_path, label, teff, log_g, metallicity, radius, rel_tol)
+      type(Colors_General_Info), intent(inout) :: rq
+      character(len=*), intent(in) :: sed_path, label
+      real(dp), intent(in) :: teff, log_g, metallicity, radius, rel_tol
+
+      real(dp), allocatable :: local_wavelengths(:), local_fluxes(:)
+      real(dp) :: local_mag, local_flux, local_interp_rad
+      real(dp) :: expected_flux, abs_err, rel_err
+
+      call calculate_bolometric( &
+         rq, teff, log_g, metallicity, radius, d_10pc, &
+         local_mag, local_flux, local_wavelengths, local_fluxes, sed_path, local_interp_rad)
+
+      expected_flux = boltz_sigma*teff**4*(radius/d_10pc)**2
+      abs_err = abs(local_flux - expected_flux)
+      rel_err = abs(local_flux - expected_flux)/expected_flux
+
+      write(*,'(a, a, a, f10.1)') '# case: ', trim(label), '=', teff
+      write(*,'(a40, 1pe23.13)') 'Wav_min_AA', minval(local_wavelengths)
+      write(*,'(a40, 1pe23.13)') 'Wav_max_AA', maxval(local_wavelengths)
+      write(*,'(a40, 1pe23.13)') 'Flux_actual', local_flux
+      write(*,'(a40, 1pe23.13)') 'Flux_expected', expected_flux
+      write(*,'(a40, 1pe23.13)') 'Flux_abserr', abs_err
+      write(*,'(a40, 1pe23.13)') 'Flux_relerr', rel_err
+
+      if (rel_err > rel_tol) then
+         write(*,'(a, a, a, 1pe11.3, a, 1pe11.3)') &
+            'wavelength coverage sanity check failed for ', trim(label), &
+            ': rel_err=', rel_err, ' > rel_tol=', rel_tol
+         stop 1
+      end if
+   end subroutine check_bolometric_coverage
 end program test_colors

colors/test/test_output

@@ -1,5 +1,7 @@
 colors: flux cube loaded (76 x 11 x 8 x 1199, 61.2 MB)
-# group1  system=Vega  grid=Kurucz2003  filters=Johnson
+
+# Group1  system=Vega  grid=Kurucz2003  filters=Johnson
+
 # case: solar
                                  Mag_bol    1.6233783742415E+01
                                 Flux_bol    3.2098630468165E-07
@@ -33,7 +35,9 @@ colors: flux cube loaded (76 x 11 x 8 x 1199, 61.2 MB)
                                        R   -2.0538595855424E-01
                                        U    3.9598722643568E+00
                                        V    7.5454159459472E-01
-# group2a  vary_MH  Teff=5778  logg=4.44
+
+# Group2a  vary_MH  Teff=5778  logg=4.44
+
 # MH=  -2.00
                                  Mag_bol    1.6233121880603E+01
                                 Flux_bol    3.2118203670543E-07
@@ -78,7 +82,9 @@ colors: flux cube loaded (76 x 11 x 8 x 1199, 61.2 MB)
                                        R    4.2891305098437E+00
                                        U    5.7511450514938E+00
                                        V    4.7518221763185E+00
-# group2b  vary_logg  Teff=5778  MH=0.0
+
+# Group2b  vary_logg  Teff=5778  MH=0.0
+
 # logg=   1.00
                                  Mag_bol    1.6233777481775E+01
                                 Flux_bol    3.2098815557824E-07
@@ -123,7 +129,9 @@ colors: flux cube loaded (76 x 11 x 8 x 1199, 61.2 MB)
                                        R    4.3286664694765E+00
                                        U    5.5395103794726E+00
                                        V    4.7780378683469E+00
-# group2c  vary_Teff  logg=4.0  MH=0.0
+
+# Group2c  vary_Teff  logg=4.0  MH=0.0
+
 # Teff=     3500.0
                                  Mag_bol    1.8411335119795E+01
                                 Flux_bol    4.3198229964375E-08
@@ -168,7 +176,10 @@ colors: flux cube loaded (76 x 11 x 8 x 1199, 61.2 MB)
                                        R    1.3843911085916E+00
                                        U    1.0937345701750E-01
                                        V    1.2754965059690E+00
-# SED sample  case=solar  columns=wavelength_AA  flux_erg_s_cm2_AA
+
+# Group3  SED sample + wavelength_coverage_sanity
+
+# SED sample  case=solar  Teff=5778  logg=4.44  FeH=0.0  columns=wavelength_AA  flux_erg_s_cm2_AA
     1.4720000000000E+02     2.1385116423739E-41
     5.1200000000000E+02     1.2669686045439E-39
     1.0950000000000E+03     5.5890114496366E-21
@@ -190,4 +201,34 @@ colors: flux cube loaded (76 x 11 x 8 x 1199, 61.2 MB)
     5.6500000000000E+04     5.2911476520239E-14
     7.2600000000000E+04     2.0114488666486E-14
     9.6200000000000E+04     6.6589153670475E-15
+
+# wavelength_coverage_sanity  logg=4.0  FeH=0.0  rel_tol=  5.00E-03
+# case: Teff=    3500.0
+                              Wav_min_AA    1.4720000000000E+02
+                              Wav_max_AA    1.6000000000000E+06
+                             Flux_actual    4.3198229964375E-08
+                           Flux_expected    4.3253426712080E-08
+                             Flux_abserr    5.5196747704715E-11
+                             Flux_relerr    1.2761242726995E-03
+# case: Teff=    6000.0
+                              Wav_min_AA    1.4720000000000E+02
+                              Wav_max_AA    1.6000000000000E+06
+                             Flux_actual    3.7325163050257E-07
+                           Flux_expected    3.7355395930932E-07
+                             Flux_abserr    3.0232880675309E-10
+                             Flux_relerr    8.0933101957231E-04
+# case: Teff=   10000.0
+                              Wav_min_AA    1.4720000000000E+02
+                              Wav_max_AA    1.6000000000000E+06
+                             Flux_actual    2.8809769255562E-06
+                           Flux_expected    2.8823607971398E-06
+                             Flux_abserr    1.3838715836804E-09
+                             Flux_relerr    4.8011740412706E-04
+# case: Teff=   20000.0
+                              Wav_min_AA    1.4720000000000E+02
+                              Wav_max_AA    1.6000000000000E+06
+                             Flux_actual    4.6160109853289E-05
+                           Flux_expected    4.6117772754238E-05
+                             Flux_abserr    4.2337099051924E-08
+                             Flux_relerr    9.1802132938073E-04
  test_colors: passed