#!/usr/bin/env python3
"""Reproduce AC11 real-image geometry with Python3 standard library only.

Download comet_orbits.py, svg_figure.py, this file, ison-image.json and
ison-ephemeris.json into one directory. Optional --force-check also needs
relative_release_reference.py and ison-perturbers.json.
"""
from pathlib import Path
import argparse
import hashlib
import json
import math
import statistics
from comet_orbits import AU_METRES, C_AU_DAY, MU, magnitude, subtract, add, scale, kepler, sky_to_native_pixel, inverse_cd, position_angle
from svg_figure import write_figure

HASHES = {'ison-image.json':'db449caff98fcd0644c2b576e7a4891524754e50b200d0c6af37222be158144b',
          'ison-ephemeris.json':'af4ac4b2f6b84de8c7c71f637d4a9c54c33f9daa996ddd3b74c54dbae037e4a4',
          'ison-perturbers.json':'8b2ae1126378c387f7147546ee4a55c2875b63e054b93ffea8a76f39893c0f0c'}

def load_checked(directory, name):
    raw = (directory/name).read_bytes()
    if hashlib.sha256(raw).hexdigest() != HASHES[name]:
        raise ValueError('Frozen input checksum differs: '+name)
    return json.loads(raw)

def morphology_profiles(image, wrong_parity=False):
    values = image['values']
    valid = [(value,x,y) for y,row in enumerate(values) for x,value in enumerate(row)
             if value is not None and image['coverage'][y][x]>=0.95]
    _,x0,y0 = max(valid)
    wcs=image['wcs']
    pixels=[]
    for value,x,y in valid:
        dx,dy=(x-x0)*8,(y-y0)*8
        if wrong_parity:
            dx=-dx
        east=(wcs['CD1_1']*dx+wcs['CD1_2']*dy)*3600
        north=(wcs['CD2_1']*dx+wcs['CD2_2']*dy)*3600
        pixels.append((math.hypot(east,north),math.degrees(math.atan2(east,north))%360,value))
    profiles=[]
    for lower,upper in [(2,6),(4,10),(6,12),(8,15)]:
        annulus=[row for row in pixels if lower<row[0]<upper]
        measurements=[]
        for angle in range(0,360,5):
            sector=[value for _,pa,value in annulus if abs((pa-angle+180)%360-180)<10]
            measurements.append(statistics.median(sector) if sector else None)
        peak=max(range(72),key=lambda index:measurements[index] if measurements[index] is not None else -math.inf)
        profiles.append({'innerArcseconds':lower,'outerArcseconds':upper,'anglesDegrees':list(range(0,360,5)),
                         'medianElectronsPerSecond':measurements,'peakPositionAngleDegrees':peak*5,
                         'interpretation':'Coarse overlapping20degree sectors on5degree grid; morphology comparison, no independent-pixel uncertainty or calibrated flux fit'})
    return [x0,y0],profiles

def main():
    parser=argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--directory',type=Path,default=Path(__file__).parent)
    parser.add_argument('--output',type=Path,default=Path('ison-results.json'))
    parser.add_argument('--figure',type=Path,default=Path('ison-projection.svg'))
    parser.add_argument('--force-check',action='store_true')
    args=parser.parse_args()
    image=load_checked(args.directory,'ison-image.json')
    eph=load_checked(args.directory,'ison-ephemeris.json')
    nucleus,profiles=morphology_profiles(image)
    _,mirrored_profiles=morphology_profiles(image,True)
    current=min(eph['releaseStates'],key=lambda row:row['ageDays'])
    sun=min(eph['sunBarycentric'],key=lambda row:abs(row['jdTdb']-eph['emissionJdTdb']))
    observer=subtract(eph['earthBarycentricAtObservation']['positionAu'],sun['positionAu'])
    line=subtract(current['positionAu'],observer)
    ra=math.degrees(math.atan2(line[1],line[0]))%360
    dec=math.degrees(math.asin(line[2]/magnitude(line)))
    nucleus_pixel=sky_to_native_pixel(line,image['wcs'])
    expected=eph['observerAstrometry']
    assert abs(ra-expected['raDegrees'])<1e-8 and abs(dec-expected['decDegrees'])<1e-8
    antisolar=position_angle(current['positionAu'],ra,dec)
    negative_velocity=position_angle(scale(current['velocityAuPerDay'],-1),ra,dec)
    assert abs(antisolar-expected['antisolarPositionAngleDegrees'])<0.001
    assert abs(negative_velocity-expected['negativeVelocityPositionAngleDegrees'])<0.001
    y_angle=math.degrees(math.atan2(image['wcs']['CD1_2'],image['wcs']['CD2_2']))
    assert abs(y_angle-image['wcs']['ORIENTAT'])<1e-10
    points=[]
    bare_parent_residuals=[]
    for release in eph['releaseStates']:
        age=round(release['ageDays'])
        if age==0:
            continue
        parent,parent_velocity=kepler(release['positionAu'],release['velocityAuPerDay'],age)
        bare=subtract(parent,current['positionAu'])
        bare_parent_residuals.append({'ageDays':age,'bareParentMinusJplNucleusKm':magnitude(bare)*AU_METRES/1000,
          'falseZeroBetaPositionAngleDegrees':position_angle(bare,ra,dec),
          'interpretation':'Solar parent is a separate model. Subtracting the JPL nucleus directly would create a spurious zero-beta tail.'})
        if age>120:
            continue
        for beta in [0.005,0.01,0.03,0.1,0.3,0.5]:
            dust,dust_velocity=kepler(release['positionAu'],release['velocityAuPerDay'],age,MU*(1-beta))
            relative=subtract(dust,parent)
            line_dust=add(line,relative)
            dust_pixel=sky_to_native_pixel(line_dust,image['wcs'])
            offset=[(a-b)/8 for a,b in zip(dust_pixel,nucleus_pixel)]
            extra_lighttime=abs(magnitude(line_dust)-magnitude(line))/C_AU_DAY
            velocity_registered=add(current['velocityAuPerDay'],subtract(dust_velocity,parent_velocity))
            lighttime_shift=magnitude(velocity_registered)*extra_lighttime*AU_METRES/1000
            points.append({'ageDays':age,'beta':beta,'relativePositionAu':relative,'blockOffset':offset,
               'positionAngleDegrees':position_angle(relative,ra,dec),
               'separationInTangentArcseconds':math.hypot((image['wcs']['CD1_1']*offset[0]+image['wcs']['CD1_2']*offset[1])*8*3600,(image['wcs']['CD2_1']*offset[0]+image['wcs']['CD2_2']*offset[1])*8*3600),
               'differentialLightTimePositionShiftEstimateKm':lighttime_shift})
    results={'kind':'Real HST median morphology and frozen JPL geometry compared with a declared solar relative-release model',
       'imageNucleusBlock':nucleus,'registeredJplNucleusNativePixel':nucleus_pixel,
       'reference':{'rightAscensionDegrees':ra,'declinationDegrees':dec,'antisolarPaDegrees':antisolar,'negativeVelocityPaDegrees':negative_velocity,'imagePositiveYAxisPaDegrees':y_angle},
       'profiles':profiles,'modelPoints':points,'bareSolarParentResiduals':bare_parent_residuals,
       'failureCases':{'mirroredCompassPeakPositionAnglesDegrees':[row['peakPositionAngleDegrees'] for row in mirrored_profiles],
         'naiveBrightestNativePixel':image['naiveBrightestNativePixel'],'correctedZeroBetaOffsetAu':[0,0,0]},
       'scope':'No ion-tail image synthesis, calibrated brightness or unique grain-size/age fit. Source2026 calibration is distinct from2013 observation. Same-emission-epoch dust approximation; differential light-time estimate saved per point.'}
    if args.force_check:
        from relative_release_reference import check_relative_model
        snapshot=load_checked(args.directory,'ison-perturbers.json')
        results['relativeForceSensitivity']=check_relative_model(snapshot,eph)
    args.output.write_text(json.dumps(results,indent=2,allow_nan=False)+'\n')
    write_figure(args.figure,image,results)
    print(json.dumps({'sectorPeakPAs':[row['peakPositionAngleDegrees'] for row in profiles],'antisolarPA':antisolar,
                      'negativeVelocityPA':negative_velocity,'modelPoints':len(points),'output':str(args.output),'figure':str(args.figure)}))

if __name__=='__main__':
    main()