#!/usr/bin/env python3
"""Replay the real DSHARP Elias24 published brightness conversion, separate from synthetic models."""
import argparse,hashlib,json,math
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_DATA=ROOT/'src/features/labs/protoplanetary-dust/data'
DEFAULT_DATA=REPO_DATA if REPO_DATA.exists() else HERE

def header_value(raw,key):
    for index in range(0,len(raw),80):
        card=raw[index:index+80].decode('ascii')
        if card[:8].strip()==key:return float(card[10:].split('/')[0].replace('D','E'))
    raise ValueError('Missing FITS header '+key)

def main():
    parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--data',type=Path,default=DEFAULT_DATA);parser.add_argument('--out',type=Path,default=ROOT/'docs/astronomy/ac14/brightness-replay.json' if REPO_DATA.exists() else HERE/'brightness-replay-rerun.json');args=parser.parse_args()
    raw=(args.data/'Elias24.profile.txt').read_bytes();header=(args.data/'Elias24_continuum.header').read_bytes()
    data=np.loadtxt(args.data/'Elias24.profile.txt');frequency=header_value(header,'CRVAL3');major=math.radians(header_value(header,'BMAJ'));minor=math.radians(header_value(header,'BMIN'))
    area=math.pi*major*minor/(4*math.log(2));c=299792458.;k=1.380649e-23;h=6.62607015e-34
    intensity=data[:,2]*1e-26/area
    rj=intensity*c*c/(2*k*frequency*frequency)
    planck=np.full(len(data),np.nan);valid=intensity>0
    planck[valid]=h*frequency/k/np.log1p(2*h*frequency**3/(c*c*intensity[valid]))
    rings=data[(data[:,0]>=65)&(data[:,0]<=95)]
    gap=data[(data[:,0]>=45)&(data[:,0]<=65)]
    peak=rings[np.argmax(rings[:,2])];minimum=gap[np.argmin(gap[:,2])]
    points=[{'radius_au':float(row[0]),'radius_arcsec':float(row[1]),'brightness_jy_beam':float(row[2]),'pointwise_error_jy_beam':float(row[3]),'rayleigh_jeans_k':float(rj[index]),'planck_k':float(planck[index])if np.isfinite(planck[index])else None}for index,row in enumerate(data)]
    result={'kind':'observed public DSHARP Elias24 profile replay; no pressure, mass or planet inference','source_sha256':hashlib.sha256(raw).hexdigest(),'header_sha256':hashlib.sha256(header).hexdigest(),'rows':len(data),'frequency_hz':frequency,'beam_major_arcsec':math.degrees(major)*3600,'beam_minor_arcsec':math.degrees(minor)*3600,'beam_solid_angle_sr':area,'max_rj_difference_from_rounded_archive_k':float(np.max(np.abs(rj-data[:,4]))),'max_planck_difference_from_rounded_archive_k':float(np.max(np.abs(planck[valid]-data[valid,6]))),'nonpositive_intensity_rows':int((~valid).sum()),'declared_radial_search_windows_au':{'ring':[65,95],'gap':[45,65]},'window_peak_radius_au':float(peak[0]),'window_gap_radius_au':float(minimum[0]),'window_brightness_contrast':float(peak[2]/minimum[2]),'uncertainty':'published pointwise errors retained; beam-correlated radial bins and calibration uncertainty preclude independent-bin confidence intervals','points':points}
    args.out.parent.mkdir(parents=True,exist_ok=True);args.out.write_text(json.dumps(result,indent=2,allow_nan=False)+'\n');print(json.dumps({key:value for key,value in result.items()if key!='points'},indent=2))

if __name__=='__main__':main()