Recompute an EMRI frequency evolution against a frozen Teukolsky flux table
===========================================================================

Files
  recompute_frequency_evolution.py  Python 3 standard library only
  inputs.json                       cases, units and the frozen flux table
  expected-output.json              output produced by the command below
  SHA256SUMS                        file identities
  LICENSE.txt                       permission and attribution

Run
  python3 recompute_frequency_evolution.py inputs.json my-output.json
  Compare my-output.json with expected-output.json. Produced with Python
  3.9.6 on macOS; agreement to about 1e-12 relative is expected elsewhere.

Units
  G = c = 1, all lengths and times in units of the large black hole mass M.
  Energy flux is dE/dt per (mu/M)^2. Seconds = M * 4.925490947e-6 s per solar
  mass. Year = 31557600 s. Frequencies in output are hertz for M = 1e6 M_sun.

Model
  Circular, equatorial orbits; exact Kerr circular-orbit energy E(r) and
  frequency Omega(r) = 1/(r^1.5 + a). Adiabatic balance dE/dt = -F, so
  dt/dr = (dE/dr)/F and N_GW = integral of Omega/pi dt. F is either the
  quadrupole flux (32/5) Omega^(10/3) or the Fujita-Tagoshi Table VII flux,
  interpolated through the factor F_table/F_quadrupole with cubic Lagrange
  interpolation on the four nearest tabulated radii. Simpson rule in ln r.

What the table is
  Table VII lists "the gravitational wave luminosity, up through l = 6". We
  read it as the energy flux to infinity per (mu/M)^2. That reading is an
  inference: it reproduces the quadrupole limit and, at r >= 50 M, the
  Tagoshi-Sasaki post-Newtonian flux to infinity within 2e-4. Horizon flux
  is not in the table. Only rows r0/M = 6 to 20 were transcribed.

Checks written into the output
  linearInterpolationCycleChange   linear instead of cubic interpolation
  fourIntervalSimpsonCycleChange   four Simpson intervals instead of 2000
  newtonianEnergyCycleChange       E = -1/(2r) instead of the Kerr energy

Failure cases
  1. Extrapolation: a request below r = 6 M or above 20 M is refused
     (see refusedCases). The a = 0.9 case therefore stops at r = 6 M,
     far outside its ISCO at 2.32 M.
  2. Newtonian energy: using E = -1/(2r) changes the Schwarzschild 10->6
     cycle count by +99 percent. The binding energy matters as much as the flux.
  3. Truncation at l = 6: Table III shows the l = 6 modes carry 0.29 percent
     of the r = 6 M Schwarzschild flux, so omitted l >= 7 modes are likely of
     order 0.1 percent there. This is an extrapolation of the trend.
  4. Horizon flux, eccentricity, inclination, resonances and every
     post-adiabatic (self-force) term are absent. A post-adiabatic phase
     correction is of order one radian and is needed for matched filtering.