"""Check the algebra of section 4 of Janosi & Herein (2026) with computer algebra.

Fact-Checked Science, October 2026, for the page "Rikitake's dynamo gets a physical derivation,
and its reversals come from inside the core" (/rikitake-mean-field/).

NEEDS SymPy (pip install sympy); everything else is the standard library. Takes a few seconds.

The paper (Geophysical Journal International 247(3), ggag390, https://doi.org/10.1093/gji/ggag390,
CC BY 4.0) realizes its canonical Rikitake model with simple sine and cosine patterns, eq. (81),
    b1 = xi1 cos(kz z + phz) e_y
    b2 = xi2 cos(ky y + phy) (e_x - e_z)
    u0 = xi0 [sin(ky y) e_x + sin(kx x) e_y] + eps e_z
with phases phy = -pi/4, phz = -pi/2 (eq. 82), averaged over the box from 0 to a quarter period
in x, y and z (eq. 83), and normalization constants xi0 = 1, xi1 = 2, xi2 = pi/(pi + 2) (eq. 84).
This script computes the coefficients the derivation needs, C_ij (eq. 16), M_ij (eq. 18) and
Lambda_ij (eq. 25), and compares them with the paper's eq. (87). It also checks the
normalization condition (9) and the Laplacian eigenmode property (10).

What it finds: with the paper's basis and constants, every coefficient of eq. (87) is
reproduced exactly, except that the small eps-terms C11 = Lambda11 come out with the opposite
sign (they are dropped in the paper, and its eq. 91 has our sign). The signs and zeros that
produce the Rikitake form hold. The constants of eq. (84), though, give <b1.b1> = 2 and
<b2.b2> = pi/(pi + 2), not 1 as eq. (9) asks; unit norm needs xi1 = sqrt(2) and
xi2 = sqrt(pi/(pi + 2)). That scales every retained coefficient by the same factor,
sqrt(2 (pi + 2) / pi) / 2 = 0.905, about 10%, not the form of the equations.
"""
import sympy as sp

x, y, z = sp.symbols("x y z", real=True)
kx, ky, kz = sp.symbols("k_x k_y k_z", positive=True)
eps = sp.symbols("epsilon", real=True)
phy, phz = -sp.pi / 4, -sp.pi / 2
xi0, xi1, xi2 = sp.Integer(1), sp.Integer(2), sp.pi / (sp.pi + 2)      # the paper's eq. (84)


def vec(*c):
    return sp.Matrix(c)


def curl(F):
    return vec(sp.diff(F[2], y) - sp.diff(F[1], z),
               sp.diff(F[0], z) - sp.diff(F[2], x),
               sp.diff(F[1], x) - sp.diff(F[0], y))


def lap(F):
    return vec(*[sp.diff(f, x, 2) + sp.diff(f, y, 2) + sp.diff(f, z, 2) for f in F])


V = (sp.pi / (2 * kx)) * (sp.pi / (2 * ky)) * (sp.pi / (2 * kz))


def avg(expr):
    """Average over the quarter-period box of eq. (83)."""
    e = sp.integrate(expr, (x, 0, sp.pi / (2 * kx)))
    e = sp.integrate(e, (y, 0, sp.pi / (2 * ky)))
    e = sp.integrate(e, (z, 0, sp.pi / (2 * kz)))
    return sp.simplify(e / V)


b1 = vec(0, xi1 * sp.cos(kz * z + phz), 0)
b2 = vec(xi2 * sp.cos(ky * y + phy), 0, -xi2 * sp.cos(ky * y + phy))
u0 = vec(xi0 * sp.sin(ky * y), xi0 * sp.sin(kx * x), eps)
bs = {1: b1, 2: b2}

print("Normalization, eq. (9), with the paper's constants (each should be 1, or 0 for b1.b2):")
print("  <b1.b1> =", avg(b1.dot(b1)), "  <b2.b2> =", avg(b2.dot(b2)), "  <b1.b2> =", avg(b1.dot(b2)),
      "  <u0.u0> (eps = 0) =", avg(u0.subs(eps, 0).dot(u0.subs(eps, 0))))
print("Laplacian eigenmodes, eq. (10): lap(b1) + kz^2 b1 =", list(sp.simplify(lap(b1) + kz ** 2 * b1)),
      "; lap(b2) + ky^2 b2 =", list(sp.simplify(lap(b2) + ky ** 2 * b2)),
      "; lap(u0) + ky^2 u0 =", list(sp.simplify(lap(u0) + ky ** 2 * u0.subs(eps, 0))), "(zero if kx = ky)")

C = {(i, j): avg(bs[i].dot(curl(u0.cross(bs[j])))) for i in (1, 2) for j in (1, 2)}
M = {(i, j): avg(bs[i].dot(curl(bs[j]))) for i in (1, 2) for j in (1, 2)}
Lam = {(i, j): avg(u0.dot(curl(bs[i]).cross(bs[j]))) for i in (1, 2) for j in (1, 2)}
xi = sp.sqrt(2) / sp.pi * xi2                                           # eq. (88)
paper = {                                                               # eq. (87)
    "C11": eps * 4 / sp.pi * kz, "C12": 16 / sp.pi ** 2 * xi * kx, "C21": sp.sqrt(2) / sp.pi * ky, "C22": 0,
    "M12": 0, "M21": -8 / sp.pi * xi * kz, "Lambda11": eps * 4 / sp.pi * kz, "Lambda22": 0,
    "Lambda12 + Lambda21": -16 / sp.pi ** 2 * xi * kz - (sp.pi - 2) / sp.pi * xi * ky,
}
ours = {"C11": C[1, 1], "C12": C[1, 2], "C21": C[2, 1], "C22": C[2, 2], "M12": M[1, 2], "M21": M[2, 1],
        "Lambda11": Lam[1, 1], "Lambda22": Lam[2, 2], "Lambda12 + Lambda21": Lam[1, 2] + Lam[2, 1]}
print("\nCoefficients, ours against the paper's eq. (87):")
for k in paper:
    o = sp.simplify(ours[k])
    same = sp.simplify(o - sp.sympify(paper[k])) == 0
    print(f"  {k:<20s} ours {str(o):<44s} paper {str(sp.simplify(paper[k])):<44s} {'same' if same else 'DIFFERENT'}")
print("\nThe Rikitake form needs (eqs 30 and 34): C11 = C22 = 0 and Lambda11 = Lambda22 = 0 (up to eps),"
      " C12 = C21 > 0 (with kx chosen by eq. 89), M12 = 0, M21 < 0, Lambda12 + Lambda21 < 0.")
vals = {k: float(sp.N(sp.simplify(v).subs({kx: 1, ky: 1, kz: 1, eps: 0}))) for k, v in ours.items()}
print("  at kx = ky = kz = 1, eps = 0:", ", ".join(f"{k} = {v:.4f}" for k, v in vals.items()))