"""Fetch daily heliocentric ecliptic-J2000 positions of Mars (barycenter, id 4) and Earth (id 399) 1580-01-01..1606-01-01 from JPL Horizons (DE441).
Writes data/kepler_ephem.csv: jd, mx,my,mz, ex,ey,ez (AU).  Run: python3 kepler_fetch.py"""
import subprocess, urllib.parse, csv, re
def fetch(cmd, a, b):
    q = urllib.parse.urlencode({'format': 'text', 'COMMAND': "'%s'" % cmd, 'OBJ_DATA': "'NO'", 'MAKE_EPHEM': "'YES'", 'EPHEM_TYPE': "'VECTORS'", 'CENTER': "'500@10'",
        'REF_PLANE': "'ECLIPTIC'", 'START_TIME': "'%s'" % a, 'STOP_TIME': "'%s'" % b, 'STEP_SIZE': "'1d'", 'VEC_TABLE': "'1'", 'OUT_UNITS': "'AU-D'", 'CSV_FORMAT': "'YES'"})
    out = subprocess.run(['curl', '-sS', '-m', '120', 'https://ssd.jpl.nasa.gov/api/horizons.api?' + q], capture_output=True, text=True).stdout
    body = out.split('$$SOE')[1].split('$$EOE')[0].strip().splitlines()
    rows = []
    for ln in body:
        f = [x.strip() for x in ln.split(',')]
        rows.append((float(f[0]), float(f[2]), float(f[3]), float(f[4])))
    return rows
ranges = [('1580-01-01', '1587-01-01'), ('1587-01-02', '1594-01-01'), ('1594-01-02', '1601-01-01'), ('1601-01-02', '1606-01-01')]
mars, earth = [], []
for a, b in ranges:
    mars += fetch('4', a, b); earth += fetch('399', a, b); print(a, b, len(mars), len(earth))
assert len(mars) == len(earth)
with open('data/kepler_ephem.csv', 'w') as f:
    w = csv.writer(f); w.writerow(['jd', 'mx', 'my', 'mz', 'ex', 'ey', 'ez'])
    for m, e in zip(mars, earth):
        assert abs(m[0] - e[0]) < 1e-6; w.writerow([m[0], m[1], m[2], m[3], e[1], e[2], e[3]])
print('rows', len(mars))