fpft.f

FPFT computes, for every shell of a rotating model, the form factors $f_P$ and $f_T$ that correct the stellar structure equations for rotational distortion, together with the mean effective gravity on each equipotential. The equipotentials' shapes come from SHAPE; the angular averages from a Gauss–Legendre quadrature.

• Who calls it: the main program (main.f), CRRECT (the structure correction, iterated), HPOINT and STARIN.
• Inputs: HD, HR, HS ($\log\rho$, $\log r$, $\log m$ in cgs), OMEGA ($\Omega$ per shell, rad s$^{-1}$), M.
• What it calls: SHAPE (the run of $r_0$ and $\eta_2$), TRAPZD (trapezoid rule for the distortion integral), POLINT (polynomial extrapolation) and QGAUSS (the angular quadrature, which calls FUNC).
• Outputs: FP, FT (the form factors), HG (mean gravity $\langle g\rangle$), ETA2 ($\eta_2$, from Radau's equation as computed in SHAPE), R0. It also fills QUADD_MOD with the spherical, rotational and distortion potentials PHISP, PHIROT, PHIDIS, read by CHECKJ, VCIRC and WRTMOD.
• Modules: QUADD_MOD (above); constants from CONSTANTS_MOD.

The method is, in the YREC paper's words (§2.6.2), "a modified version of the R. Kippenhahn & H. C. Thomas (1970) technique," which "assumes uniform rotation on equipotential surfaces" (the paper's §2.6 traces the approach to Endal & Sofia 1976, whom SHAPE's comment cites too): each shell has one $\Omega$ (one OMEGA per shell), and the equipotentials are $r = r_0(1 - aP_2)$ (SHAPE's comment), keeping only the $P_2$ term (how accurate that stays for fast rotators is question 3). For each shell, outward:

1. The distortion amplitude $a = \dfrac{5}{3}\dfrac{\Omega^2 r_0^3}{Gm\,(2+\eta_2)}$.

2. The integral $I(r_0) = \int_0^{r_0}(\ldots)\,dr_0'$ (its integrand, per TRAPZD's comment, is $(\rho/m)\,r_0'^7\,\Omega^2(5+\eta_2)/(2+\eta_2)$) is added shell by shell: the trapezoid rule on the step from the previous shell, refined once and extrapolated to zero step with POLINT (a Romberg step; JMAX = 2, tolerance $10^{-6}$). PHIDIS $= \tfrac{4\pi}{3} I/r^3$ (C4PI*CC13*AINT/RPHI3, line 98↗), the distorted potential of the comment ("needed for the evaluation of the distorted potential," lines 45–47↗); what $I$ is physically is question 2. The interpolation inside TRAPZD that this integral relies on is a suspected bug, trapzd.f, suspected bug 1.

3. QGAUSS integrates over $0 \le \theta \le \pi/2$ the surface element $S$, $gS$ and $S/g$, giving $S_\Phi$, $S_\Phi\langle g\rangle$ and $S_\Phi\langle g^{-1}\rangle$.

4. With $r_\Phi$ the shell's radius (HR), $$f_P = \frac{4\pi r_\Phi^4}{Gm\,S_\Phi\langle g^{-1}\rangle},\qquad f_T = \frac{(4\pi r_\Phi^2)^2}{S_\Phi\langle g\rangle\;S_\Phi\langle g^{-1}\rangle},$$ which is the usual Kippenhahn–Thomas form, $f_T = (4\pi r_\Phi^2/S_\Phi)^2/(\langle g\rangle\langle g^{-1}\rangle)$, provided QGAUSS's outputs are the surface integrals named above (read from its code). HG = $\langle g\rangle$. They enter the structure as $d\ln P/d\ln m = -\dfrac{Gm^2}{4\pi r^4P}f_P$ (COEFFT) and $\nabla_\mathrm{rad} = \dfrac{3\kappa LP}{64\pi\sigma Gm T^4}\dfrac{f_T}{f_P}$ (TPGRAD); both are 1 without rotation.