func.f

FUNC evaluates, at one colatitude $\theta$ on an equipotential surface of a rotating star, the effective gravity $g = |\nabla\Phi|$ and the surface element. It is the integrand for the angular averages $\langle g\rangle$ and $\langle g^{-1}\rangle$ that the rotation form factors need.

• Who calls it: only QGAUSS, a ten-point Gauss–Legendre quadrature over $0 \le \theta \le \pi/2$, called by FPFT.
• Inputs: X ($\theta$), R0 (mean radii of the equipotentials), HS ($\log m$), AINT and Q (the distortion integral and its integrand, from TRAPZD), W2 ($\Omega^2$), A (the distortion amplitude), I (the shell).
• Outputs: G ($g$) and S (the surface element per unit $\theta$).
• Modules: constants only (CONSTANTS_MOD: $4\pi$, $1/3$, $\ln 10$, $\log G$).

With $P_2 = 1 - \tfrac32\sin^2\theta$, the equipotential is $r = r_0(1 - aP_2)$. Then $$\frac{\partial\Phi}{\partial r} = \frac{Gm}{r^2} - \frac{4\pi}{3r^4}P_2\left(3I - rQ\,\frac{\partial r_0}{\partial r}\right) - \Omega^2 r\sin^2\theta,$$ $$\frac{\partial\Phi}{\partial\theta} = \frac{4\pi}{3r^4}\left(P_2\,Q\,\frac{\partial r_0}{\partial\theta} - 3I\cos\theta\sin\theta\right) - \Omega^2 r\cos\theta\sin\theta,$$ with $I$ = AINT, $\partial r_0/\partial r = 1/(1-4aP_2)$ and $\partial r_0/\partial\theta = 3a\cos\theta\sin\theta\,(\partial r_0/\partial r)/(1-aP_2)$, and $$g = \sqrt{(\partial_r\Phi)^2 + (\partial_\theta\Phi)^2},\qquad S = 4\pi r\sin\theta\, r_0\sqrt{(1-aP_2)^2 + (3a\cos\theta\sin\theta)^2}.$$

S carries $4\pi$ rather than the $2\pi$ of a ring, which counts both hemispheres: FPFT calls QGAUSS with the upper limit $\pi/2$ (B = 0.125D0*C4PI, fpft.f, line 39↗), so that the integral of S is the whole area of the equipotential (for a sphere, $\int_0^{\pi/2} 4\pi r^2\sin\theta\,d\theta = 4\pi r^2$) (settled from the code, Nov 9, 2026: lines 46–47↗).