Skip to content

Szanyi et al. (2025) nucleosynthesis models

Model reader for Szanyiet al. (2025) models.

The models are available to download from here: https://zenodo.org/records/14981333 Set 2 contains the latest nuclear reaction rates.

Szanyi25Reader

Reader for Szanyi et al. (2025) nucleosynthesis models.

Load a given mass, metallicity, and rate set model and provide methods to access the data. This routine will load the closest available model. Please see the example on how to check which model was loaded.

Example usage: 

from rttools.nucmod.szanyi25 import Szanyi25Reader

model = Szanyi25Reader(mass=3.0, metallicity=0.014, rate_set=2)

# print mass, metallicity and rate set of the loaded model
print(f"Loaded model: mass={model.mass}, z={model.z}, rate_set={model.rate_set}")

Source code in src/rttools/nucmod/szanyi25.py 
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
class Szanyi25Reader:
    """Reader for Szanyi et al. (2025) nucleosynthesis models.

    Load a given mass, metallicity, and rate set model and provide methods to access the data.
    This routine will load the closest available model. Please see the example on how to check which model was loaded. 

    Example usage:
    ```python
    from rttools.nucmod.szanyi25 import Szanyi25Reader

    model = Szanyi25Reader(mass=3.0, metallicity=0.014, rate_set=2)

    # print mass, metallicity and rate set of the loaded model
    print(f"Loaded model: mass={model.mass}, z={model.z}, rate_set={model.rate_set}")
    ```
    """

    def __init__(self, mass: float, metallicity: float, rate_set: int = 2):
        """Initialize the reader with the closest model file.

        :param mass: Stellar mass in solar masses.
        :param metallicity: Stellar metallicity.
        :param rate_set: Model set number (0, 1, or 2), defaults to 2 (if not given or invalid).
        """
        self._model_file = get_closest_model(mass, metallicity, rate_set)

        self._mass = None
        self._y = None
        self._z = None
        self._rate_set = None
        self._data = None

        self._tp_cols = None

        self._read_model()

    def _read_model(self):
        """Read in the model data from the file and fill out
        - _mass: Stellar mass (float)
        - _y: Helium mass fraction (float)
        - _z: Metallicity (float)
        - _rate_set: Rate set number (int)
        - _data: Nucleosynthesis data (pandas DataFrame)
        """
        # Read in the first line and parse mass, z, y, and set
        with open(self._model_file, "r") as f:
            hdr = f.readline().strip()

        hdr = hdr.split(",")
        self._mass = float(hdr[0].split("=")[1].strip())
        self._z = float(hdr[1].split("=")[1].strip())
        self._y = float(hdr[2].split("=")[1].strip())
        self._rate_set = float(hdr[3].strip().split(" ")[1].strip())

        # read the data into a pandas DataFrame
        self._data = pd.read_csv(
            self._model_file, skiprows=2, delimiter="\t", index_col=0
        )

        self._tp_cols = self._data.columns[4:]

    def get_co_ratio(self) -> pd.Series:
        """Get the C/O ratio for all thermal pulses (as number ratio)."""
        c12 = self._data.loc["c12"][self._tp_cols]
        c13 = self._data.loc["c13"][self._tp_cols]
        o16 = self._data.loc["o16"][self._tp_cols]
        o17 = self._data.loc["o17"][self._tp_cols]
        o18 = self._data.loc["o18"][self._tp_cols]

        c_total = c12 + c13
        o_total = o16 + o17 + o18

        co_ratio = (c_total / ini.ele["C"].mass) / (o_total / ini.ele["O"].mass)
        return co_ratio

    def get_delta(self, nominator: str, denominator: str) -> pd.Series:
        """Get the delta value for the two isotopes across all thermal pulses."""
        nomin = ini.iso[nominator]
        denom = ini.iso[denominator]

        ind_nomin = "".join(nomin.name.split("-")).lower()
        ind_denom = "".join(denom.name.split("-")).lower()

        ratio = (self._data.loc[ind_nomin][self._tp_cols] / self._data.loc[ind_denom][self._tp_cols])

        delta = ini.iso_delta(nominator, denominator, ratio, mass_fraction=True)
        return delta
    def plot_mod(
        self,
        ax: plt.Axes,
        xisos: [str, str],
        yisos: [str, str],
        marker: str = "o",
        color: str = "tab:blue",
        label=None,
        linestyle: str = "-",
        **kwargs,
    ) -> None:
        """Plot Lugaro Model curve onto a plot.

        Note that **kwargs are passed along to both plotting routines. For C/O < 1, the
        marker is set to 'None'.

        :param ax: Matplotlib axis to plot on.
        :param xisos: Two-element list of isotopes for x-axis (nominator, denominator).
        :param yisos: Two-element list of isotopes for y-axis (nominator, denominator).
        :param marker: Marker of the plot, defaults to 'o'
        :param color: Color of the plot, defaults to 'tab:blue'
        :param linestyle: Linestyle of the plot, defaults to '-'

        :return: Nothing
        """
        xdat_all = self.get_delta(*xisos)
        ydat_all = self.get_delta(*yisos)
        co_ratio = self.get_co_ratio()

        crich_mask = co_ratio >= 1

        # plot all
        ax.plot(
            xdat_all, ydat_all, marker="None", linestyle=linestyle, color=color, **kwargs
        )
        # plot c-rich
        ax.plot(
            xdat_all[crich_mask],
            ydat_all[crich_mask],
            marker=marker,
            color=color,
            linestyle="None",
            label=label,
            **kwargs,
        )

    @property
    def mass(self) -> float:
        """Stellar mass in solar masses."""
        return self._mass

    @property
    def y(self) -> float:
        """Helium mass fraction."""
        return self._y

    @property
    def z(self) -> float:
        """Metallicity."""
        return self._z

    @property
    def rate_set(self) -> int:
        """Nuclear reaction rate set number."""
        return self._rate_set

mass property

Stellar mass in solar masses.

rate_set property

Nuclear reaction rate set number.

y property

Helium mass fraction.

z property

Metallicity.

__init__(mass, metallicity, rate_set=2)

Initialize the reader with the closest model file.

Parameters:

Name Type Description Default
mass float

Stellar mass in solar masses.

 required
metallicity float

Stellar metallicity.

 required
rate_set int

Model set number (0, 1, or 2), defaults to 2 (if not given or invalid).

 2
Source code in src/rttools/nucmod/szanyi25.py 
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
def __init__(self, mass: float, metallicity: float, rate_set: int = 2):
    """Initialize the reader with the closest model file.

    :param mass: Stellar mass in solar masses.
    :param metallicity: Stellar metallicity.
    :param rate_set: Model set number (0, 1, or 2), defaults to 2 (if not given or invalid).
    """
    self._model_file = get_closest_model(mass, metallicity, rate_set)

    self._mass = None
    self._y = None
    self._z = None
    self._rate_set = None
    self._data = None

    self._tp_cols = None

    self._read_model()

get_co_ratio()

Get the C/O ratio for all thermal pulses (as number ratio).

Source code in src/rttools/nucmod/szanyi25.py 
113
114
115
116
117
118
119
120
121
122
123
124
125
def get_co_ratio(self) -> pd.Series:
    """Get the C/O ratio for all thermal pulses (as number ratio)."""
    c12 = self._data.loc["c12"][self._tp_cols]
    c13 = self._data.loc["c13"][self._tp_cols]
    o16 = self._data.loc["o16"][self._tp_cols]
    o17 = self._data.loc["o17"][self._tp_cols]
    o18 = self._data.loc["o18"][self._tp_cols]

    c_total = c12 + c13
    o_total = o16 + o17 + o18

    co_ratio = (c_total / ini.ele["C"].mass) / (o_total / ini.ele["O"].mass)
    return co_ratio

get_delta(nominator, denominator)

Get the delta value for the two isotopes across all thermal pulses.

Source code in src/rttools/nucmod/szanyi25.py 
127
128
129
130
131
132
133
134
135
136
137
138
def get_delta(self, nominator: str, denominator: str) -> pd.Series:
    """Get the delta value for the two isotopes across all thermal pulses."""
    nomin = ini.iso[nominator]
    denom = ini.iso[denominator]

    ind_nomin = "".join(nomin.name.split("-")).lower()
    ind_denom = "".join(denom.name.split("-")).lower()

    ratio = (self._data.loc[ind_nomin][self._tp_cols] / self._data.loc[ind_denom][self._tp_cols])

    delta = ini.iso_delta(nominator, denominator, ratio, mass_fraction=True)
    return delta

plot_mod(ax, xisos, yisos, marker='o', color='tab:blue', label=None, linestyle='-', **kwargs)

Plot Lugaro Model curve onto a plot.

Note that **kwargs are passed along to both plotting routines. For C/O < 1, the marker is set to 'None'.

Parameters:

Name Type Description Default
ax Axes

Matplotlib axis to plot on.

 required
xisos [str, str]

Two-element list of isotopes for x-axis (nominator, denominator).

 required
yisos [str, str]

Two-element list of isotopes for y-axis (nominator, denominator).

 required
marker str

Marker of the plot, defaults to 'o'

 'o'
color str

Color of the plot, defaults to 'tab:blue'

 'tab:blue'
linestyle str

Linestyle of the plot, defaults to '-'

 '-'

Returns:

Type Description
None

Nothing
Source code in src/rttools/nucmod/szanyi25.py 
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
def plot_mod(
    self,
    ax: plt.Axes,
    xisos: [str, str],
    yisos: [str, str],
    marker: str = "o",
    color: str = "tab:blue",
    label=None,
    linestyle: str = "-",
    **kwargs,
) -> None:
    """Plot Lugaro Model curve onto a plot.

    Note that **kwargs are passed along to both plotting routines. For C/O < 1, the
    marker is set to 'None'.

    :param ax: Matplotlib axis to plot on.
    :param xisos: Two-element list of isotopes for x-axis (nominator, denominator).
    :param yisos: Two-element list of isotopes for y-axis (nominator, denominator).
    :param marker: Marker of the plot, defaults to 'o'
    :param color: Color of the plot, defaults to 'tab:blue'
    :param linestyle: Linestyle of the plot, defaults to '-'

    :return: Nothing
    """
    xdat_all = self.get_delta(*xisos)
    ydat_all = self.get_delta(*yisos)
    co_ratio = self.get_co_ratio()

    crich_mask = co_ratio >= 1

    # plot all
    ax.plot(
        xdat_all, ydat_all, marker="None", linestyle=linestyle, color=color, **kwargs
    )
    # plot c-rich
    ax.plot(
        xdat_all[crich_mask],
        ydat_all[crich_mask],
        marker=marker,
        color=color,
        linestyle="None",
        label=label,
        **kwargs,
    )

get_closest_model(mass, metallicity, rate_set=2)

Get the closest model file for a given mass, metallicity, and set.

Parameters:

Name Type Description Default
mass float

Stellar mass in solar masses.

 required
metallicity float

Stellar metallicity.

 required
rate_set int

Model set number (0, 1, or 2), defaults to 2 (if not given or invalid).

 2

Returns:

Type Description
Path

Path to the closest model file.
Source code in src/rttools/nucmod/szanyi25.py 
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
def get_closest_model(mass: float, metallicity: float, rate_set: int = 2) -> Path:
    """Get the closest model file for a given mass, metallicity, and set.

    :param mass: Stellar mass in solar masses.
    :param metallicity: Stellar metallicity.
    :param rate_set: Model set number (0, 1, or 2), defaults to 2 (if not given or invalid).

    :return: Path to the closest model file.
    """
    if rate_set not in [0, 1, 2]:
        rate_set = 2

    if mass < 2.5:
        mass_str = "m2"
    elif mass < 3.5:
        mass_str = "m3"
    else:
        mass_str = "m4"

    if mass_str == "m2":
        metall_str = "z014"  # only one available
    else:  # now we have 0.007, 0.014, 0.03
        if metallicity < 0.0105:
            metall_str = "z007"
        elif metallicity < 0.022:
            metall_str = "z014"
        else:
            metall_str = "z030"

    fname = MODEL_DIR.joinpath(f"srf_{mass_str}{metall_str}_set{rate_set}.dat")
    assert fname.exists(), f"Model file {fname} does not exist."

    return fname