#!/usr/bin/env python3
"""Independent remeasurement of frozen synthetic radial states, not a replacement Lab engine."""
import argparse,json
from pathlib import Path
import numpy as np
HERE=Path(__file__).resolve().parent
ROOT=HERE.parents[2] if len(HERE.parents)>2 else HERE
REPO_INPUT=ROOT/'docs/astronomy/ac14/radial-results.json'

def main():
    parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--input',type=Path,default=REPO_INPUT if REPO_INPUT.exists() else HERE/'radial-results.json');parser.add_argument('--out',type=Path,default=ROOT/'docs/astronomy/ac14/radial-replay.json' if REPO_INPUT.exists() else HERE/'radial-replay-rerun.json');args=parser.parse_args();source=json.loads(args.input.read_text());results=[]
    for case in source['presets']:
        gas=np.array(case['gas']);dust=np.array(case['dust']);spacing=source['length']/case['cells'];settings=case['settings']
        observation={'gas_mass':float(gas[0].sum()*spacing),'dust_mass':float(dust[0].sum()*spacing),'dust_peak_to_mean':float(dust[0].max()/dust[0].mean()),'gas_peak_to_mean':float(gas[0].max()/gas[0].mean()),'peak_loading':float(np.max(dust[0]/gas[0])),'max_gas_radial_speed':float(np.max(np.abs(gas[1]/gas[0]))),'gas_positive':bool(np.all(gas[0]>0)),'dust_positive':bool(np.all(dust[0]>0))}
        observation['relative_gas_mass_error']=observation['gas_mass']/source['length']-1;observation['relative_dust_mass_error']=observation['dust_mass']/(settings['loading']*source['length'])-1
        results.append({'case':case['preset'],'observation':observation})
    result={'kind':'independent NumPy remeasurement of actual saved synthetic radial finite-volume states, not observed fields or streaming instability','results':results,'spatial_convergence':source['spatialConvergence'],'time_convergence':source['timeConvergence']};args.out.write_text(json.dumps(result,indent=2,allow_nan=False)+'\n');print(json.dumps(results,indent=2))

if __name__=='__main__':main()