eqstat.f

EQSTAT is the equation of state that the rest of YREC calls. Its header says it was added ("LLP 10-22-06") as a "new front end" so that the derivatives of the equation of state can be computed numerically, by calling the older routine, now EQSTAT2, at neighbouring points. All the physics is in EQSTAT2.

• Who calls it: ten routines across the code: COEFFT, ENVINT, HSUBP, MASSLOSS, PHYSIC, QATM, QENV, SCONVEC, STARIN and WRTOUT.
• Inputs: TL, PL ($\log_{10}$ of $T$ in K and $P$ in dyn cm$^{-2}$), X, Z, and the flags LDERIV (second derivatives wanted) and LATMO (atmosphere). DL serves as a starting guess.
• What it calls: only EQSTAT2.
• Outputs: those of EQSTAT2: T, P, D, DL, BETA, FXION, RMU, AMU, EMU, ETA, QDT $=\partial\ln\rho/\partial\ln T|_P$, QDP $=\partial\ln\rho/\partial\ln P|_T$, QCP ($c_P$), DELA ($\nabla_\mathrm{ad}$), and the derivatives QDTT, QDTP, QAT, QAP, QCPT, QCPP.
• Modules: it reads only CLN ($\ln 10$, from CONSTANTS_MOD) and LNumDeriv (from NMLPAR_MOD, set in the PHYSICS namelist read by PARMIN).

If LDERIV and LNumDeriv are both set, EQSTAT2 is called four times without derivatives, at $\log T \pm h_T$ and $\log P \pm h_P$, with $h_T = 3\times10^{-4}$ and $h_P = 1.5\times10^{-3}$ (the "approximate table intervals," 0.03 and 0.15, scaled down by 100: "empirically, scaling the intervals down by 100 seemed to work best"). Central differences give the second derivatives, for example $$\frac{\partial Q_{\rho T}}{\partial\ln T} = \frac{Q_{\rho T}^{+} - Q_{\rho T}^{-}}{2h_T\ln 10},\qquad \frac{\partial\log c_P}{\partial\log T} = \frac{\log c_P^{+} - \log c_P^{-}}{2h_T},\qquad \frac{\partial\log\nabla_\mathrm{ad}}{\partial\log T} = \frac{\log\nabla_\mathrm{ad}^{+} - \log\nabla_\mathrm{ad}^{-}}{2h_T},$$ and the same in $\log P$. A last call at the point itself, without derivatives, gives the values; its derivative outputs go to dummy variables so they cannot overwrite the numerical ones. Otherwise EQSTAT2 is called once, with analytic derivatives if LDERIV is set. The second derivatives are needed because the Henyey method linearizes the structure equations in $\ln P$ and $\ln T$: TPGRAD differentiates the temperature gradient (QAT, QAP $=\partial\ln\nabla_\mathrm{ad}/\partial\ln T, \partial\ln P$, and QCPT, QDTT, QDTP for the convective gradient), and COEFFT differentiates $\epsilon_\mathrm{grav}$, which contains $Q_{\rho T}$ and $\nabla_\mathrm{ad}$ (QDTT, QDTP, QAT, QAP). The symbols are in the overview's Notation: $Q_{\rho T}$ = QDT $= -\delta$ and $Q_{\rho P}$ = QDP $= \alpha$ in the notation of Kippenhahn, Weigert & Weiss.

Microphysics in YREC: LNUMDERIV (default false) replaces YREC's semi-analytic EOS derivatives by numerical ones from the OPAL tables, at about three times the cost and with little effect on the solution.