From de560005e52bc6be3861dc5b403d8ca29e4c6fec Mon Sep 17 00:00:00 2001 From: Ebraheem Farag <63124736+Debraheem@users.noreply.github.com> Date: Sat, 28 Mar 2026 03:32:20 +0100 Subject: [PATCH] bugfix use Fe/H instead of zbase, and allow negative metallicity --- colors/README.rst | 51 ++++++++++++++++++-- colors/defaults/colors.defaults | 75 ++++++++++++++++++++++++------ colors/private/colors_ctrls_io.f90 | 4 ++ colors/private/colors_history.f90 | 5 +- colors/public/colors_def.f90 | 3 +- star/private/history.f90 | 25 +++++++++- 6 files changed, 140 insertions(+), 23 deletions(-) diff --git a/colors/README.rst b/colors/README.rst index d80216147..6222cc4e0 100644 --- a/colors/README.rst +++ b/colors/README.rst @@ -22,7 +22,7 @@ How does the MESA colors module work? The module operates by coupling the stellar structure model with pre-computed grids of stellar atmospheres. -1. **Interpolation**: At each timestep, the module takes the star's current surface parameters—Effective Temperature (:math:`T_{\rm eff}`), Surface Gravity (:math:`\log g`), and Metallicity ([M/H])—and queries a user-specified library of stellar atmospheres (defined in ``stellar_atm``). It interpolates within this grid to construct a specific Spectral Energy Distribution (SED) for the stars current features. +1. **Interpolation**: At each timestep, the module takes the star's current surface parameters, Effective Temperature (:math:`T_{\rm eff}`), Surface Gravity (:math:`\log g`), and Metallicity ([Fe/H]), and queries a user-specified library of stellar atmospheres (defined in ``stellar_atm``). It interpolates within this grid to construct a specific Spectral Energy Distribution (SED) for the stars current features. 2. **Convolution**: This specific SED is then convolved with filter transmission curves (defined in ``instrument``) to calculate the flux passing through each filter. @@ -57,7 +57,7 @@ stellar_atm Specifies the path to the directory containing the grid of stellar atmosphere models. This directory must contain: -1. **lookup_table.csv**: A map linking filenames to physical parameters (:math:`T_{\rm eff}`, :math:`\log g`, [M/H]). +1. **lookup_table.csv**: A map linking filenames to physical parameters (:math:`T_{\rm eff}`, :math:`\log g`, [Fe/H]). 2. **SED files**: The actual spectra (text or binary format). 3. **flux_cube.bin**: (Optional but recommended) A binary cube for rapid interpolation. @@ -87,6 +87,51 @@ The distance to the star in centimeters. distance = 5.1839d20 + +z_over_x_ref +------------ + +**Default:** ``0.0230057173972`` (the GS98 solar :math:`Z/X` ratio) + +Reference metal-to-hydrogen ratio used to convert the model photospheric +composition into the metallicity coordinate passed to the atmosphere grid: + +.. math:: + + [Fe/H] = \log_{10}\left(\frac{(Z/X)_\mathrm{phot}}{(Z/X)_\mathrm{ref}}\right) + +The default is chosen to match the default ``stellar_atm = +'/data/colors_data/stellar_models/Kurucz2003all/'`` grid, which is based on +ATLAS9 models using the Grevesse & Sauval (1998) solar abundance scale. The +numeric value ``0.0230057173972`` is computed from the MESA GS98 constants +``Z_solar = 0.0169`` and ``Y_solar = 0.2485``, so +``(Z/X)_ref = 0.0169 / (1 - 0.0169 - 0.2485)``. + +If you switch to a different atmosphere library, you must set ``z_over_x_ref`` +so that its reference composition matches the metallicity labels in that +library's ``lookup_table.csv``. For custom abundance patterns, making this +choice is the responsibility of the user. + +If the photospheric metallicity, expressed in terms of ``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 +``[Fe/H] = log10((Z/X)/(Z/X)_ref)``, and so it uses the lowest metallicity in +the table instead. This is a fallback safeguard and not a dedicated treatment +for H-free atmospheres. + +For the shipped ``Kurucz2003all`` table in this checkout, the tabulated +metallicity range is ``[Fe/H] = -2.5`` to ``[Fe/H] = 4.0``. With the default +``z_over_x_ref = 2.30057173972d-2``, the lower edge corresponds to an adopted +``Z/X`` of about ``7.275d-5``. + +**Example:** + +.. code-block:: fortran + + z_over_x_ref = 2.30057173972d-2 + make_csv -------- @@ -216,10 +261,10 @@ Below are the default values for the colors module parameters as defined in ``co vega_sed = '/data/colors_data/stellar_models/vega_flam.csv' stellar_atm = '/data/colors_data/stellar_models/Kurucz2003all/' distance = 3.0857d19 ! 10 parsecs in cm (Absolute Magnitude) + z_over_x_ref = 2.30057173972d-2 ! GS98 solar Z/X for Kurucz2003all make_csv = .false. colors_results_directory = 'SED' mag_system = 'Vega' - vega_sed = '/data/colors_data/stellar_models/vega_flam.csv' Visual Summary of Data Flow =========================== diff --git a/colors/defaults/colors.defaults b/colors/defaults/colors.defaults index 0bdb1719a..eb2155252 100644 --- a/colors/defaults/colors.defaults +++ b/colors/defaults/colors.defaults @@ -13,8 +13,6 @@ ! ~~~~~~~~~~~~~~ ! ``vega_sed`` ! ~~~~~~~~~~~~ - ! ``stellar_atm`` - ! ~~~~~~~~~~~~~~~ ! ``distance`` ! ~~~~~~~~~~~~ ! ``make_csv`` @@ -24,15 +22,16 @@ ! ``mag_system`` ! ~~~~~~~~~~~~~~ - ! 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, then MESA will calculate bolometric and + ! synthetic magnitudes by interpolating stellar atmosphere model grids and + ! convolving the resulting spectra with the selected photometric filters. + ! The Vega SED is used for photometric zero points when + ! ``mag_system = 'Vega'``. The stellar distance is given in cm. + ! :: use_colors = .false. instrument = '/data/colors_data/filters/Generic/Johnson' - stellar_atm = '/data/colors_data/stellar_models/Kurucz2003all/' distance = 3.0857d19 make_csv = .false. colors_results_directory = 'SED' @@ -40,6 +39,49 @@ vega_sed = '/data/colors_data/stellar_models/vega_flam.csv' + ! ``stellar_atm`` + ! ~~~~~~~~~~~~~~~ + + ! Path to the directory containing the stellar atmosphere lookup table and + ! spectra used by the colors module. This directory should contain a + ! ``lookup_table.csv`` file and the corresponding SED files. The default + ! uses the ``Kurucz2003all`` ATLAS9 atmosphere grid shipped with MESA. + + ! :: + + stellar_atm = '/data/colors_data/stellar_models/Kurucz2003all/' + + + ! ``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)``, so it uses the lowest metallicity in the + ! table instead. 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 ! --------------------- @@ -47,13 +89,16 @@ ! 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' + read_extra_colors_inlist(:) = .false. + extra_colors_inlist_name(:) = 'undefined' diff --git a/colors/private/colors_ctrls_io.f90 b/colors/private/colors_ctrls_io.f90 index dc560be8e..934905043 100644 --- a/colors/private/colors_ctrls_io.f90 +++ b/colors/private/colors_ctrls_io.f90 @@ -38,6 +38,7 @@ module colors_ctrls_io character(len=32) :: mag_system real(dp) :: distance + real(dp) :: z_over_x_ref logical :: make_csv logical :: use_colors @@ -46,6 +47,7 @@ module colors_ctrls_io vega_sed, & stellar_atm, & distance, & + z_over_x_ref, & make_csv, & mag_system, & colors_results_directory, & @@ -155,6 +157,7 @@ subroutine store_controls(rq, ierr) rq%vega_sed = vega_sed rq%stellar_atm = stellar_atm rq%distance = distance + rq%z_over_x_ref = z_over_x_ref rq%make_csv = make_csv rq%colors_results_directory = colors_results_directory rq%use_colors = use_colors @@ -191,6 +194,7 @@ subroutine set_controls_for_writing(rq) vega_sed = rq%vega_sed stellar_atm = rq%stellar_atm distance = rq%distance + z_over_x_ref = rq%z_over_x_ref make_csv = rq%make_csv colors_results_directory = rq%colors_results_directory use_colors = rq%use_colors diff --git a/colors/private/colors_history.f90 b/colors/private/colors_history.f90 index e65c699e9..dc0421d8a 100644 --- a/colors/private/colors_history.f90 +++ b/colors/private/colors_history.f90 @@ -112,7 +112,8 @@ subroutine data_for_colors_history_columns( & filter_name = trim(remove_dat(color_filter_names(i))) names(i + filter_offset) = filter_name - if (t_eff >= 0 .and. metallicity >= 0) then + ! Negative [M/H] values are valid for metal-poor atmosphere grids. + if (t_eff >= 0) then ! Select precomputed zero-point based on magnitude system select case (trim(cs%mag_system)) case ('VEGA', 'Vega', 'vega') @@ -151,4 +152,4 @@ subroutine data_for_colors_history_columns( & end subroutine data_for_colors_history_columns -end module colors_history \ No newline at end of file +end module colors_history diff --git a/colors/public/colors_def.f90 b/colors/public/colors_def.f90 index ec37f2326..7c59b7238 100644 --- a/colors/public/colors_def.f90 +++ b/colors/public/colors_def.f90 @@ -46,6 +46,7 @@ module colors_def character(len=256) :: mag_system real(dp) :: metallicity real(dp) :: distance + real(dp) :: z_over_x_ref logical :: make_csv logical :: use_colors integer :: handle @@ -206,4 +207,4 @@ subroutine do_free_colors_tables end subroutine do_free_colors_tables -end module colors_def \ No newline at end of file +end module colors_def diff --git a/star/private/history.f90 b/star/private/history.f90 index d6faffd84..e7396e751 100644 --- a/star/private/history.f90 +++ b/star/private/history.f90 @@ -70,7 +70,8 @@ subroutine do_history_info(s, write_flag, ierr) num_extra_cols, num_binary_cols, num_extra_binary_cols, num_colors_cols,num_extra_header_items, n integer, parameter :: num_epsnuc_out = 12 real(dp) :: & - epsnuc_out(num_epsnuc_out), csound_surf, v_surf, envelope_fraction_left + epsnuc_out(num_epsnuc_out), csound_surf, v_surf, envelope_fraction_left, m_div_h, & + min_m_div_h, max_m_div_h integer :: mixing_regions, mix_relr_regions, burning_regions, burn_relr_regions integer, pointer :: mixing_type(:), mix_relr_type(:), burning_type(:), burn_relr_type(:) character (len = maxlen_history_column_name), pointer, dimension(:) :: & @@ -267,7 +268,27 @@ subroutine do_history_info(s, write_flag, ierr) end if colors_col_names(1:num_colors_cols) = 'unknown' colors_col_vals(1:num_colors_cols) = -1d99 - call data_for_colors_history_columns(s%T(1), log10(s%grav(1)), s%R(1), s%kap_rq%Zbase, & + + if (colors_settings% z_over_x_ref <= 0d0) then + write(*, *) 'colors error: z_over_x_ref must be positive' + ierr = -1 + call dealloc + return + end if + + min_m_div_h = minval(colors_settings% lu_meta) + max_m_div_h = maxval(colors_settings% lu_meta) + + ! Map the current photospheric Z/X onto the atmosphere-table metallicity axis. + if (s% X(s% photosphere_cell_k) > 0d0 .and. s% Z(s% photosphere_cell_k) > 0d0) then + m_div_h = log10((s% Z(s% photosphere_cell_k)/s% X(s% photosphere_cell_k)) / & + colors_settings% z_over_x_ref) + m_div_h = max(min_m_div_h, min(max_m_div_h, m_div_h)) + else + m_div_h = min_m_div_h + end if + + call data_for_colors_history_columns(s%T(1), log10(s%grav(1)), s%R(1), m_div_h, & s% colors_handle, num_colors_cols, colors_col_names, colors_col_vals, ierr) if (ierr /= 0) then call dealloc