#!/usr/bin/env python3
"""Optional independent analytic Parker reference. Python3.10+, SciPy1.13.1."""
import argparse
import json
import math
from pathlib import Path
import numpy as np
from scipy.special import lambertw
from scipy import __version__ as scipy_version
import importlib.util
spec=importlib.util.spec_from_file_location('independent_wind',Path(__file__).with_name('replay-wind.py'))
wind=importlib.util.module_from_spec(spec)
spec.loader.exec_module(wind)


def main():
    parser=argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--out',type=Path,default=Path(__file__).with_name('lambert-reference-rerun.json'))
    args=parser.parse_args()
    comparisons=[]
    for x in [.03,.1,.3,.8,.9999,1,1.0001,1.2,2,5,20]:
        d=math.exp(3-4/x)/x**4
        mach=1.0 if x==1 else math.sqrt(-float(lambertw(-d,0 if x<1 else -1).real))
        ode=math.exp(wind.transonic_log_mach(x))
        comparisons.append({'radius_over_sonic':x,'lambert_mach':mach,'ode_mach':ode,'relative_difference':abs(ode/mach-1)})
    x=2;d=math.exp(3-4/x)/x**4
    wrong=math.sqrt(-float(lambertw(-d,0).real))
    result={'kind':'synthetic analytic branch benchmark; not observed wind data','scipy_version':scipy_version,'numpy_version':np.__version__,'comparisons':comparisons,'wrong_branch_at_twice_sonic':{'principal_branch_mach':wrong,'accelerating_branch_mach':comparisons[8]['lambert_mach']},'convention':'W0 inside sonic point, W-1 outside; exact sonic point set by its regular limit'}
    args.out.write_text(json.dumps(result,indent=2,allow_nan=False)+'\n')
    print(json.dumps(result,indent=2))

if __name__=='__main__':main()