#!/usr/bin/env python3
"""Freeze the JPL queries used by the AC11 ISON projection tutorial.

Uses curl's normal certificate validation; no authentication is required.
"""
from pathlib import Path
import csv
import io
import json
import subprocess
import argparse

parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('--directory', type=Path, default=Path('/tmp/sandboxphysics-ac11-source'), help='Cache directory for raw queries and parsed ephemerides')
CACHE = parser.parse_args().directory
BASE = 'https://ssd.jpl.nasa.gov/api/horizons.api'
UTC_MJD = (56392.21727116 + 56392.22041947) / 2
COMET = "'90004196'"

def fetch(name, parameters):
    target = CACHE / (name + '.json')
    if target.exists():
        return json.loads(target.read_text())
    arguments = ['curl', '--fail', '--silent', '--show-error', '--get', '--max-time', '45', BASE]
    for key, value in parameters.items():
        arguments += ['--data-urlencode', key + '=' + str(value)]
    raw = subprocess.run(arguments, check=True, capture_output=True).stdout
    result = json.loads(raw)
    if 'error' in result or '$$SOE' not in result.get('result', ''):
        raise ValueError(result)
    target.write_bytes(raw)
    (CACHE / (name + '-query.json')).write_text(json.dumps(parameters, indent=2) + '\n')
    return result

def table(result):
    text = result['result'].split('$$SOE')[1].split('$$EOE')[0].strip()
    return list(csv.reader(io.StringIO(text)))

def main():
    CACHE.mkdir(exist_ok=True)
    observer = fetch('horizons-observer', {
        'format': 'json', 'COMMAND': COMET, 'CENTER': "'500@399'", 'MAKE_EPHEM': 'YES',
        'EPHEM_TYPE': 'OBSERVER', 'TLIST': "'" + repr(UTC_MJD) + "'", 'TLIST_TYPE': 'MJD',
        'TIME_TYPE': 'UT', 'TIME_DIGITS': 'FRACSEC', 'REF_SYSTEM': 'ICRF',
        'QUANTITIES': "'1,19,20,21,24,27,30'", 'CSV_FORMAT': 'YES', 'ANG_FORMAT': 'DEG', 'EXTRA_PREC': 'YES'
    })
    row = table(observer)[0]
    delta_t_seconds = float(row[13])
    observation_tdb = UTC_MJD + 2400000.5 + delta_t_seconds / 86400
    emission_tdb = observation_tdb - float(row[9]) / 1440
    ages = [0, 5, 10, 15, 20, 30, 45, 60, 90, 120, 180, 240]
    common = {'format':'json','MAKE_EPHEM':'YES','EPHEM_TYPE':'VECTORS','TIME_TYPE':'TDB',
              'TLIST_TYPE':'JD','REF_SYSTEM':'ICRF','REF_PLANE':'FRAME','VEC_CORR':'NONE',
              'OUT_UNITS':'AU-D','CSV_FORMAT':'YES','VEC_TABLE':'2'}
    comet = fetch('horizons-comet-release', {**common,'COMMAND':COMET,'CENTER':"'500@10'",
                  'TLIST':' '.join("'"+repr(emission_tdb-age)+"'" for age in ages)})
    earth = fetch('horizons-earth-barycentric', {**common,'COMMAND':"'399'",'CENTER':"'500@0'",'TLIST':"'"+repr(observation_tdb)+"'"})
    sun = fetch('horizons-sun-barycentric', {**common,'COMMAND':"'10'",'CENTER':"'500@0'",'TLIST':"'"+repr(emission_tdb)+"' '"+repr(observation_tdb)+"'"})
    def state(record):
        return {'jdTdb':float(record[0]),'positionAu':[float(x) for x in record[2:5]],'velocityAuPerDay':[float(x) for x in record[5:8]]}
    comet_states = [state(record) for record in table(comet)]
    result = {'observationMjdUtc':UTC_MJD,'observationJdTdb':observation_tdb,'tdbMinusUtcSeconds':delta_t_seconds,
              'emissionJdTdb':emission_tdb,'oneWayLightTimeDays':float(row[9])/1440,
              'frame':'ICRF, heliocentric release states; barycentric observer translated to Sun at emission',
              'cometSolution':'JPL#57, record90004196, solution2025-Nov-21_13:42:48; includes nuclear non-gravitational A1/A2',
              'observerAstrometry':{'raDegrees':float(row[3]),'decDegrees':float(row[4]),'sunDistanceAu':float(row[5]),'observerDistanceAu':float(row[7]),'phaseAngleDegrees':float(row[10]),'antisolarPositionAngleDegrees':float(row[11]),'negativeVelocityPositionAngleDegrees':float(row[12])},
              'releaseStates':[{**record,'ageDays':emission_tdb-record['jdTdb']} for record in comet_states],
              'earthBarycentricAtObservation':state(table(earth)[0]),'sunBarycentric':list(map(state,table(sun)))}
    (CACHE/'ephemeris.json').write_text(json.dumps(result,indent=2)+'\n')
    daily_times = [emission_tdb-day for day in reversed(range(121))]
    daily_options = {'START_TIME':"'JD"+repr(daily_times[0])+"'",'STOP_TIME':"'JD"+repr(daily_times[-1]+1e-6)+"'",'STEP_SIZE':"'1 d'"}
    comet_daily = fetch('horizons-comet-daily', {**common,'COMMAND':COMET,'CENTER':"'500@10'",**daily_options})
    planet_parameters = [
      ('mercury',1,22031.868551),('venus',2,324858.592000),('earth-moon',3,403503.235625),
      ('mars',4,42828.375816),('jupiter',5,126712764.100000),('saturn',6,37940584.841800),
      ('uranus',7,5794556.400000),('neptune',8,6836527.100580)
    ]
    perturbers = []
    for name,identifier,gm in planet_parameters:
        answer = fetch('horizons-'+name+'-daily',{**common,'COMMAND':"'"+str(identifier)+"'",'CENTER':"'500@10'",**daily_options})
        perturbers.append({'name':name,'command':identifier,'gravitationalParameterKm3PerS2':gm,'states':list(map(state,table(answer)))})
    comparison = {'startJdTdb':daily_times[0],'stopJdTdb':daily_times[-1],
      'comet':list(map(state,table(comet_daily))),'planets':perturbers,
      'planetMassSource':'https://ssd.jpl.nasa.gov/astro_par.html (DE440 listed system masses; point-system approximation)',
      'nucleusOutgassing':{'A1AuPerDay2':8.537815093994e-8,'A2AuPerDay2':5.695189237595e-9,'A3AuPerDay2':0,
                         'ALN':0.1112620426,'NK':4.6142,'NM':2.15,'NN':5.093,'R0Au':2.808},
      'purpose':'Independent sensitivity of solar relative-release model to planetary tidal forces and unshared nucleus outgassing, not a full JPL dust integration'}
    (CACHE/'perturber-sample.json').write_text(json.dumps(comparison,separators=(',',':'))+'\n')
    print(json.dumps({'midpointMjdUtc':UTC_MJD,'midpointJdTdb':observation_tdb,'emissionJdTdb':emission_tdb,'tdbMinusUtcSeconds':delta_t_seconds,'releaseStates':len(comet_states),'referencePA':result['observerAstrometry']},indent=2))

if __name__ == '__main__':
    main()