"""Find oppositions of Mars (heliocentric longitude of Earth == that of Mars) in data/kepler_ephem.csv; print date, JD, Mars heliocentric ecliptic longitude (J2000)."""
import csv, math, numpy as np
rows = np.array([[float(x) for x in r] for r in list(csv.reader(open('data/kepler_ephem.csv')))[1:]])
jd = rows[:, 0]; lm = np.arctan2(rows[:, 2], rows[:, 1]); le = np.arctan2(rows[:, 5], rows[:, 4])
def wrap(x): return (x + math.pi) % (2 * math.pi) - math.pi
d = wrap(lm - le)
def cal(j):  # Julian-calendar date for JD (Tycho used the Julian calendar; Denmark adopted Gregorian in 1700)
    j = j + 0.5; Z = int(j); F = j - Z; A = Z; B = A + 1524; C = int((B - 122.1) / 365.25); D = int(365.25 * C); E = int((B - D) / 30.6001)
    day = B - D - int(30.6001 * E) + F; month = E - 1 if E < 14 else E - 13; year = C - 4716 if month > 2 else C - 4715
    return year, month, day
opp = []
for i in range(len(jd) - 1):
    if d[i] > 0 and d[i + 1] <= 0 and abs(d[i]) < 0.5:  # Earth overtakes Mars: sign change + to -
        t = d[i] / (d[i] - d[i + 1]); j = jd[i] + t * (jd[i + 1] - jd[i])
        # linear interpolation of positions
        pm = rows[i, 1:3] + t * (rows[i + 1, 1:3] - rows[i, 1:3]); lam = math.degrees(math.atan2(pm[1], pm[0])) % 360
        rm = float(np.hypot(*pm)); opp.append((j, lam, rm))
if __name__ == '__main__':
    print('n oppositions', len(opp))
    for j, lam, rm in opp:
        y, m, dd = cal(j); print('%4d-%02d-%05.2f (Julian)  JD %.2f  Mars helio lon %.4f deg  r=%.4f AU' % (y, m, dd, j, lam, rm))