"""Evaluate the RHS of Turyshev & Toth 2020 Eq153, not a mission forecast."""
from pathlib import Path
import csv
import json
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

OUT=Path(__file__).resolve().parent
YEAR=365.25*86400

def exposure(side=100, diameter_m=1, output_snr=7, z_au=650, target_pc=30):
    """Reproduce the displayed fixed scenario, not arbitrary mission designs.

    Supported: integer 25–200 pixels per side; aperture 1 or 2 m;
    reconstructed SNR 7, observer 650 AU and target 30 pc.
    """
    if isinstance(side, bool) or not isinstance(side, int) or not 25 <= side <= 200:
        raise ValueError('Choose an integer from 25 to 200 pixels per side.')
    if isinstance(diameter_m, bool) or diameter_m not in (1, 2):
        raise ValueError('The displayed comparison supports a 1 or 2 metre aperture.')
    if (output_snr, z_au, target_pc) != (7, 650, 30):
        raise ValueError('This reproduction fixes SNR=7, observer=650 AU and target=30 pc.')
    n=side**2
    pixel=3.54e-3*n*output_snr**2*(1+.79*(650/z_au)**5.1+.05*(z_au/650)**2.65)*diameter_m**-3*(target_pc/30)**2*(650/z_au)**3.4
    return {'side_pixels':side,'total_pixels':n,'aperture_m':diameter_m,'output_snr':output_snr,'observer_AU':z_au,'target_pc':target_pc,'pixel_rhs_seconds':pixel,'total_rhs_seconds':pixel*n,'total_rhs_Julian_years':pixel*n/YEAR}

def main():
    rows=[exposure(s,d) for d in [1,2] for s in [25,50,75,100,150,200]]
    with (OUT/'exposure-scaling.csv').open('w',newline='',encoding='utf-8') as f:
        w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)
    methods={'classification':'Arithmetic evaluation of the right-hand side of a published approximate inequality; not executed image recovery or operational upper guarantee',
             'source':'https://arxiv.org/pdf/2002.06492v2','doi':'10.1103/PhysRevD.102.024038','locator':'Eqs152-153, manuscript p28; NIAC PhaseII section2.4.4',
             'formula':'tpix <= 3.54e-3*N*SNR_R^2*[1+0.79*(650/z)^5.1+0.05*(z/650)^2.65]*d^-3*(z0/30)^2*(650/z)^3.4 seconds; displayed value is RHS, total=N*RHS',
             'units':{'z':'AU','z0':'parsec','d':'metre','N':'total pixels (side squared)','year':'365.25 days'},
             'fixed_assumptions':'Full-phase Lambertian Earth/Sun analogue, broadband albedo0.3 and incident irradiance1366.83W/m2; solar-corona-dominated shot noise assuming mean corona removal and an annular coronagraph; empirical filtering factor a=10 retained from the paper; SNR_R=7 after deconvolution, target30pc, observer650AU, one sequential sampler.',
             'exclusions':['No travel time, slews, downlink, downtime, navigation or mission lifetime.', 'Not a prediction for the synthetic 6m optical benchmark, whose noise and source are different.', 'Do not equate N with side length, or combine these fixed-input scalings with different published megapixel scenarios.', 'Aperture/noise/filter/source assumptions can change results; the paper inequality is not an operational guarantee.'],
             'reference':exposure(),'two_metre_comparison':exposure(diameter_m=2)}
    (OUT/'exposure-methods.json').write_text(json.dumps(methods,indent=2)+'\n',encoding='utf-8')
    plt.rcParams.update({'font.family':'sans-serif','font.sans-serif':['Arial','DejaVu Sans','sans-serif'],'font.size':13,'svg.fonttype':'none'})
    for mobile in [False,True]:
        plt.rcParams.update({'font.size':22 if mobile else 13})
        fig,ax=plt.subplots(figsize=(6,6) if mobile else (10,6),layout='constrained')
        sides=list(range(25,201))
        for d,c in [(1,'#9b288f'),(2,'#176a8a')]:
            ax.plot(sides,[exposure(s,d)['total_rhs_Julian_years'] for s in sides],color=c,lw=2,label=f'{d} m aperture')
        ax.scatter([100],[exposure()['total_rhs_Julian_years']],color='#9b288f')
        ax.set(xlabel='Pixels along each side\n(N = side²)' if mobile else 'Pixels along each side (N = side²)',ylabel='Integration RHS (years)' if mobile else 'Calculated integration RHS (Julian years)',yscale='log')
        if mobile:ax.set_xticks([25,100,200])
        ax.grid(alpha=.25,which='both');ax.legend();ax.set_title('Published exposure scaling\nSNR 7 · 30 pc · 650 AU')
        stem='exposure-mobile' if mobile else 'exposure'
        fig.savefig(OUT/(stem+'.svg'));fig.savefig(OUT/(stem+'.png'),dpi=130)
        if not mobile:fig.savefig(OUT/'exposure.pdf')
        plt.close(fig)
    print(json.dumps(methods['reference']))

if __name__=='__main__':main()