The StatStar Python program models the internal structure of a star using the four fundamental differential equations of stellar structure:

1. hydrostatic equilibrium
2. mass continuity
3. energy generation
4. energy transport

It integrates these equations inward from the stellar surface to the center, using physical laws of gas pressure, radiation, and nuclear energy generation.

When you run the program (statstar_run()), it first asks for:

1. the stellar mass \( M \) in solar units
2. the luminosity \( L \) in solar units
3. the effective temperature \( T_{\rm eff} \) in kelvins
4. the chemical composition (fractions of hydrogen \( X \) and metals \( Z \))

From these, the helium fraction is computed as \( Y = 1 - X - Z \). Then the program converts all quantities into cgs units and sets the fundamental physical constants (e.g., \( G, k_B, m_H, c, \sigma \), etc.). It estimates the stellar radius from the Stefan–Boltzmann law:

$$ L = 4 \pi R^2 \sigma T_{\rm eff}^4 $$

which gives

$$ R = \sqrt{\frac{L}{4\pi\sigma}} \, T_{\rm eff}^{-2}. $$

The mean molecular weight is calculated assuming complete ionization:

$$ \mu = \frac{1}{2X + 0.75Y + 0.5Z}. $$

The star is divided into thin spherical shells indexed by radius \( r \). The outermost shell (zone 1) is placed at \( r = R_s \), the surface radius. At the surface, the total mass \( M_r = M_s \) and luminosity \( L_r = L_s \) are known.

The program uses \( T_0 = 0 \) and \( P_0 = 0 \) as mathematical boundary conditions at the outermost layer to simplify the integration start. However, this does not mean that the physical surface temperature is zero—the physical \( T_{\rm eff} \) provided by the user is already used earlier to determine the stellar radius through the Stefan–Boltzmann law. Thus, \( T_0 = 0 \) and \( P_0 = 0 \) simply mark the starting point of numerical integration, not real physical values.

Because the structure equations become numerically unstable near the surface, the code uses approximate analytic expansions for the first few layers. This is done by the function STARTMDL, which estimates new values of pressure \( P \) and temperature \( T \) for small inward steps in radius \( \Delta r \). It assumes the mass and luminosity are constant across these outer zones and uses either:

1. the radiative temperature gradient (if energy is carried mainly by radiation), or
2. the convective gradient (if the region is unstable to convection).

The choice is controlled by a flag irc (0 = radiative, 1 = convective). For the first ten zones, the model is built with these approximate surface solutions.

In each zone, the subroutine EOS computes physical properties of the gas:

1. Density \( \rho \) from the ideal gas law (subtracting radiation pressure)
2. Opacity \( \kappa \) as the sum of bound–free, free–free, and electron–scattering components
3. Energy generation rate \( \epsilon \) from the proton–proton chain and the CNO cycle, using analytic fits to nuclear reaction rates

It also returns an auxiliary factor *tog_bf* (the guillotine-to-Gaunt ratio) used in opacity calculations.

If any nonphysical condition appears (e.g., negative pressure, temperature, or density), the program aborts with an error message.

Once the outer 10 zones are built, the program switches to full numerical integration using the Runge–Kutta method implemented in RUNGE and FUNDEQ.

For each radial step \( \Delta r \) inward, the four differential equations are evaluated:

$$ \frac{dP}{dr} = -\frac{G\rho M_r}{r^2}, \qquad \frac{dM_r}{dr} = 4\pi r^2\rho, \qquad \frac{dL_r}{dr} = 4\pi r^2\rho\epsilon $$

and for the temperature gradient,

$$ \frac{dT}{dr} = \begin{cases} -\frac{3\kappa\rho L_r}{16\pi a c T^3 r^2}, & \text{(radiative)} \\ -\frac{1}{\gamma/(\gamma-1)}\frac{G M_r \mu m_H}{k_B r^2}, & \text{(convective)} \end{cases} $$

These derivatives are combined by a fourth-order Runge–Kutta integrator to estimate new values of \( P, M_r, L_r, T \) at the next inner shell. From the new \( P \) and \( T \), the code recomputes \( \rho, \kappa, \epsilon \) using EOS. The logarithmic pressure–temperature gradient \( d\ln P / d\ln T \) is compared with the adiabatic value \( \gamma/(\gamma-1) \) to decide whether the next zone is radiative or convective. This process repeats zone by zone, moving inward toward the center of the star. The step size \( \Delta r \) is automatically reduced when entering the denser core regions to maintain numerical stability.

The integration stops when any of the following occurs:

1. the radius becomes very small (center reached)
2. the enclosed mass or luminosity becomes non-positive (unphysical)
3. density or temperature behave abnormally
4. the maximum number of shells is reached

At the end, the code estimates the central conditions by extrapolation:

1. central density \( \rho_c = M_r / (\tfrac{4}{3}\pi r^3) \)
2. central pressure from a Taylor expansion
3. central temperature from the ideal gas law
4. central energy generation \( \epsilon_c = L_r / M_r \)

If everything is consistent, the code reports *“CONGRATULATIONS, I THINK YOU FOUND IT!”* Otherwise, it prints diagnostic messages indicating what went wrong (e.g., too high density, negative luminosity, etc.).

Finally, the program writes all computed data to the text file starmodl_py.dat. The header lists:

1. the surface and central conditions
2. total mass, radius, luminosity, and temperature
3. chemical composition
4. central density, temperature, pressure, and energy generation rate

Below that, it prints a table for every computed zone (from center outward) containing:

| r | Qm (1 − M_r/M_total) | L_r | T | P | ρ | κ | ε | dlnP/dlnT |

Each row is tagged with c for convective zones and r for radiative ones. If the pressure–temperature gradient exceeds the allowed limit, a * warning is added.

The file thus records a full radial profile of the modeled star — how its internal temperature, pressure, density, luminosity, and energy generation vary from center to surface — based on the user-supplied mass, luminosity, temperature, and composition.

1. The stellar model is determined entirely by the four structure equations and the input boundary conditions.
2. The Runge–Kutta integration allows the model to evolve smoothly from the surface to the core.
3. Radiative and convective energy transport are treated self-consistently through local temperature gradients.
4. The computed structure directly yields central quantities such as core temperature and density.

1. How does the choice of hydrogen fraction \( X \) affect the core temperature of the resulting model?
2. Why is the Runge–Kutta method particularly suitable for stellar structure equations?
3. Under what conditions would the program detect instability due to excessive radiation pressure?
4. How do radiative and convective gradients differ physically, and how does the code decide between them?