"""Schwarzschild vacuum regression: excised interior mass, zero-weight test particles."""
import numpy as np
from pbhgr.spherical_ev import Grid, Particles, RunConfig, Run, solve_metric
from pbhgr.analytic import circular_orbit_L, schwarzschild_radial_infall, eta_of_t

M = 1.0


def _grid(K=600):
    return Grid(r_out=40.0, K=K, r_in=2.5 * M, m_in=M)


def test_exterior_metric_is_schwarzschild():
    g = _grid()
    P = Particles(np.array([10.0]), np.array([0.0]), np.array([0.0]), np.array([0.0]), np.array([1.0]))
    met = solve_metric(g, P)
    assert np.allclose(met.m_edge, M)
    exact = 0.5 * np.log(1 - 2 * M / g.r_edge)
    e1 = np.max(np.abs(met.mu_edge - exact))
    g2 = _grid(K=1200)
    e2 = np.max(np.abs(solve_metric(g2, P).mu_edge - 0.5 * np.log(1 - 2 * M / g2.r_edge)))
    assert e1 < 1e-3                       # midpoint-rule error, steepest near the excision
    assert 3.0 < e1 / e2 < 5.0             # second-order convergence


def test_circular_orbit_stays_circular():
    g = _grid()
    r0 = 8.0 * M
    P = Particles(np.array([r0]), np.array([0.0]), np.array([circular_orbit_L(M, r0)]),
                  np.array([0.0]), np.array([1.0]))
    errs = []
    for K in (600, 1200):
        P = Particles(np.array([r0]), np.array([0.0]), np.array([circular_orbit_L(M, r0)]),
                      np.array([0.0]), np.array([1.0]))
        run = Run(_grid(K), P, RunConfig(t_end=2 * np.pi * np.sqrt(r0**3 / M) * 3, cfl=0.4, diag_every=100000))
        run.run()
        errs.append(abs(run.P.r[0] - r0))
        assert abs(run.P.w[0]) < 1e-4
    assert errs[0] < 1e-3 * r0
    assert 3.0 < errs[0] / errs[1] < 5.0   # second-order convergence of the orbit radius


def test_radial_infall_coordinate_time():
    g = _grid(K=1200)
    R0 = 15.0 * M
    P = Particles(np.array([R0]), np.array([0.0]), np.array([0.0]), np.array([0.0]), np.array([1.0]))
    t_end = 60.0
    run = Run(g, P, RunConfig(t_end=t_end, cfl=0.4, diag_every=100000))
    run.run()
    eta = eta_of_t(M, R0, t_end)
    r_exact, tau_exact, _ = schwarzschild_radial_infall(M, R0, eta)
    assert abs(run.P.r[0] - r_exact) < 2e-3 * R0
    assert abs(run.P.tau[0] - tau_exact) < 2e-3 * tau_exact