"""Independent Green extraction and inner harmonic tests of the new solar fields."""
import json
import numpy as np
from scipy.special import eval_legendre
from dark_cmc_model import BASE, OUT, CMCResponse, digest
from dark_cassini_followup import load_field, composite_gauss, field_derivatives
from dark_solar_green_check import normalization_check
import dark_mercury_check as mercury


def green(field,law,order):
    spec=field["spec"]
    radius,angular,spline=load_field(BASE/field["field"])
    y,wy=composite_gauss(order,np.log(radius))
    u,wu=composite_gauss(order,np.r_[-1.,angular,1.])
    f=field_derivatives(spline,y,u,spec["external"]/spec["scale"])
    offset=law.offsets(f["g"])[0]
    integrand=offset*((3*u*u-1)/2*f["gr"]-u*(1-u*u)*f["gu"]/f["r"])
    length=np.sqrt(6.67430e-11*1.98847e30/spec["scale"])
    measure=-4.5*spec["scale"]/length*(wy/np.exp(y))[:,None]*integrand*wu
    q=float(np.sum(measure))
    assert np.isfinite(q)
    multipoles={}
    for degree in [3,4]:
        p=eval_legendre(degree,u)
        derivative=degree*(eval_legendre(degree-1,u)-u*p)/(1-u*u)
        source=offset*((degree+1)*p*f["gr"]-(1-u*u)*derivative*f["gu"]/f["r"])
        coefficient=.5*np.einsum("i,ij,j->",wy*np.exp((1-degree)*y),source,wu)
        multipoles[str(degree)]=float(coefficient)
    return {"Q2_s_minus2":q,"positive_integral_s_minus2":float(measure[measure>0].sum()),
            "negative_integral_s_minus2":float(measure[measure<0].sum()),
            "absolute_integral_over_abs_Q2":float(np.sum(abs(measure))/abs(q)),
            "dimensionless_potential_multipole_coefficients":multipoles}


def inner_profile(field):
    spec=field["spec"]
    data=np.genfromtxt(BASE/field["field"],delimiter=",",names=True)
    radius=np.unique(data["r"])
    angular=np.unique(data["u"])
    potential=data["correction"].reshape(len(radius),len(angular))
    p2=(3*angular**2-1)/2
    p2-=p2.mean()
    centered=potential-potential.mean(axis=1)[:,None]
    coefficient=np.einsum("ij,j->i",centered,p2)/np.sum(p2*p2)
    length=np.sqrt(6.67430e-11*1.98847e30/spec["scale"])
    local=-3*coefficient/radius**2*spec["scale"]/length
    # A conservative floating-point projection scale, not a rigorous PDE
    # error estimate. The quadrupole is much smaller than the monopole.
    roundoff=30*np.finfo(float).eps*np.max(abs(potential),axis=1)
    q_roundoff=3*roundoff/radius**2*spec["scale"]/length
    rows=[]
    for lower,upper in [(0.3,0.5),(1.,3.),(10.,30.),(60.,120.)]:
        mask=(radius*length/mercury.AU>=lower)&(radius*length/mercury.AU<=upper)
        if np.any(mask):
            rows.append({"radius_interval_au":[lower,upper],"radial_cells":int(mask.sum()),
                "minimum_Q2_from_shell_s_minus2":float(local[mask].min()),
                "maximum_Q2_from_shell_s_minus2":float(local[mask].max()),
                "median_Q2_from_shell_s_minus2":float(np.median(local[mask])),
                "maximum_roundoff_projection_scale_s_minus2":float(q_roundoff[mask].max())})
    return {"minimum_grid_radius_au":float(radius[0]*length/mercury.AU),"shells":rows,
            "scope":"Diagnostic shell P2 projections; inner harmonic continuation remains the controlled Q2 estimator"}


def main():
    fields=json.loads((OUT/"fields_summary.json").read_text())
    controls=json.loads((OUT/"controls_summary.json").read_text())
    for report in [fields,controls]:
        assert report["source_sha256"]==digest(BASE/"dark_cmc_model.py")
    law=CMCResponse()
    central=min(fields["cases"],key=lambda c:abs(c["solar"]["spec"]["external"]-2.32e-10))["solar"]
    largest=next(c["solar"] for c in controls["cases"] if c["solar"]["spec"]["nr"]==512 and c["solar"]["spec"]["rmin"]==.0003)
    extended=next(c["solar"] for c in controls["cases"] if c["solar"]["spec"]["rmin"]==.00001)
    green_rows=[]
    for field in [central,largest]:
        integrals=[green(field,law,order) for order in [2,4]]
        values=[item["Q2_s_minus2"] for item in integrals]
        direct=field["result"]["Q2_s_minus2"][1]
        difference=abs(values[-1]-direct)
        change=abs(values[-1]-values[0])
        assert difference<.04e-27 and change<.01e-27,(values,direct)
        green_rows.append({"nr":field["spec"]["nr"],"Green_Q2_two_and_four_nodes_s_minus2":values,
            "inner_projection_Q2_s_minus2":direct,"difference_s_minus2":difference,
            "quadrature_change_s_minus2":change,"integrals":integrals})
    external_integrals=[]
    for item in fields["cases"]:
        field=item["solar"]
        integral=(green_rows[0]["integrals"][0] if field is central else green(field,law,2))
        difference=abs(integral["Q2_s_minus2"]-field["result"]["Q2_s_minus2"][1])
        assert difference<.06e-27,difference
        row={"external_acceleration_m_s2":field["spec"]["external"],"integral":integral,
             "difference_from_inner_projection_s_minus2":difference}
        external_integrals.append(row)
        print(json.dumps(row),flush=True)
    result={"source_sha256":digest(__file__),"normalization_check":normalization_check(),
        "green_checks":green_rows,"external_integrals":external_integrals,
        "inner_extended_grid":inner_profile(extended),
        "scope":"Green extraction of the same AQUAL fields, not a second nonlinear solver; shell profiles are numerical diagnostics",
        "input_sha256":{str(p.relative_to(BASE)):digest(p) for p in [
            OUT/"fields_summary.json",OUT/"controls_summary.json",BASE/"dark_cmc_model.py",
            BASE/"dark_cassini_followup.py",BASE/"dark_solar_green_check.py",BASE/"dark_mercury_check.py"]}}
    (OUT/"solar_checks_summary.json").write_text(json.dumps(result,indent=2,allow_nan=False)+"\n")
    print(json.dumps(result,indent=2),flush=True)


if __name__=="__main__":
    main()
