eqsaha.f

EQSAHA is the Yale equation of state for partially ionized gas: it solves the Saha equations for hydrogen, helium (two stages) and ten metals, and returns the density, the ionization fractions and the thermodynamic derivatives. The name "Prather EOS" in EQSTAT2's comments refers to it.

• Who calls it: only EQSTAT2, below the temperature TSCUT. It is called even when the SCV tables are used, since "this EOS replaces the call to EQSAHA, except for derivative purposes"; EQSCVE then overwrites what it can.
• Inputs: FX, the number fractions of the twelve species (set by EQSTAT2 from FXENV in /COMP/); TL, PL; BETA and its relatives; AMU; LATMO (skip the derivatives, for atmospheres); LDERIV (also compute second derivatives).
• Outputs: D, DL, EMU, RMU = $R(\mu_\mathrm{ion}^{-1}+\mu_e^{-1}) = \mathcal{R}/\mu$, the gas constant over the mean molecular weight, so that $P_\mathrm{gas} = \rho\,(\mathcal{R}/\mu)\,T$ (AMU $=\mu_\mathrm{ion}^{-1}$ and EMU $=\mu_e^{-1}$ are its two terms), FXION (H$^+$, He$^+$, He$^{++}$ fractions), QDT, QDP, QCP, DELA and, with LDERIV, their derivatives. KSAHA is the caller's iteration counter, to which this call adds its iterations.
• Modules: only CONSTANTS_MOD (with CMKH the constant term of the Saha equation) and IOUNITS_MOD (the unit ISHORT).

1. Saha constants. For each species $i$, with ionization temperature $\chi_i/k$ (SAHAT, in K: 59 630 for Na up to 250 200 for Ne; SAHET0 = 285 270 and SAHET1 = 631 370 for He and He$^+$) and a constant $u_i$ (SAHAU), $$K_i = \frac{1}{\beta}\exp\!\left[\ln 10\,(u_i + 2.5\log T - \log P + C) - \frac{\chi_i}{kT}\right].$$ The species order is that of ATOMW2 in EQSTAT2: Na, Al, Mg, Fe, Si, C, H, O, N, Ar, Ne, He. The list is in order of increasing ionization temperature. The loop stops at the first species whose exponent is below $-40$: that species and all later ones are treated as neutral (so hydrogen, seventh, counts as ionized only if NZ0 ≥ 7, and FXION(1) is 0 otherwise). A species whose exponent is above $+40$ gets $K_i = 10^{16}$ (fully ionized) and is left out of the derivative sums.

What it leaves out. The solver has no molecules (H$_2$), no H$^-$, no pressure ionization and no Coulomb terms, and treats the electrons as a classical gas. With the code's defaults (LOPALE06 and LSCV false) it nevertheless gives the equation of state of every layer below TSCUT, including the cool, dense envelopes of M dwarfs and brown dwarfs (whether that is adequate there, where molecules and degeneracy would be expected to matter, is the question below). The documentation calls OPAL 2006 "strongly recommended," and 45 of the 51 test-suite runs set LOPALE06 (47 set LSCV), so in practice this solver is used mainly outside those tables and for the second derivatives; whether the default should change is eqstat2.f, question 4.