#!/usr/bin/env python3
"""How long the Sun and the eclipsed Moon shared a sky, place by place.

The word "selenelion" suggests a special kind of eclipse. It is not one. It is
the overlap of two ordinary regions -- where the Sun's upper limb is up, and
where the Moon's upper limb is up -- and because the Moon during an eclipse sits
within a degree and a half of the point exactly opposite the Sun, and both discs
have width, those two regions always overlap a little. At any single minute the
overlap is a hairline along the day/night line. But the Earth turns, and over
the three and a half hours of umbral phases the hairline sweeps.

This map counts, for every point on Earth, how many minutes it spent inside that
overlap on 10 December 2011. Both of the video's observers are marked. Santa Fe
had two minutes. Cahokia, 1,443 km east, had none -- it sits within a degree of
the footprint's eastern edge, where the window shrinks to nothing.

Sub-solar point, sub-lunar point and both cap radii come from
docs/selenelion/data/eclipse-2011-12-10-3d.json, interpolated to 10-second steps;
nothing is re-derived. Land outlines from the global-land-mask package.

    pip install global_land_mask
    python3 selenelion_band_map.py
"""

import base64, io, json, os
import numpy as np
from PIL import Image
from global_land_mask import globe

HERE = os.path.dirname(os.path.abspath(__file__))
ROOT = os.path.dirname(HERE)
SRC = os.path.join(ROOT, 'docs', 'selenelion', 'data', 'eclipse-2011-12-10-3d.json')
OUT = os.path.join(ROOT, 'docs', 'selenelion', 'img', 'fig-selenelion-band.svg')

INK, DIM, FAINT, RED = '#2E4057', '#6E675B', '#a09a8c', '#c0392b'
SEA, LAND = '#f2f5f8', '#e0dacb'

LAT1, LAT0 = 84.0, -68.0
GRID = 0.25
SAMPLE = 10.0                       # seconds
W, H = 980, 660
ML, MR, MT = 46.0, 16.0, 84.0
MAPH = 236.0                        # world map height
IN_X, IN_Y, IN_W, IN_H = 46.0, 400.0, 330.0, 186.0   # North America inset
IN_LON0, IN_LON1, IN_LAT0, IN_LAT1 = -118.0, -78.0, 26.0, 50.0

HAIRLINE_UT = 14 * 3600 + 3 * 60    # 07:03 MST — the minute the band held Santa Fe

# minutes-of-overlap bands, painted palest to deepest
LEVELS = [(1 / 6, 1.0, '#f3e2b4', 'seconds'),
          (1.0, 2.0, '#ecd08a', '1–2 min'),
          (2.0, 5.0, '#dcae4e', '2–5 min'),
          (5.0, 15.0, '#c1892a', '5–15 min'),
          (15.0, 60.0, '#96631c', '15–60 min'),
          (60.0, 1e9, '#5e3d11', 'over an hour')]

OBS = [('Santa Fe', 35.6870, -105.9378, 'two minutes', 'below'),
       ('Cahokia', 38.6606, -90.0619, 'under six seconds', 'above')]

D = json.load(open(SRC))
ROWS = D['timeline']


def sec(hms):
    h, m, s = (int(v) for v in hms.split(':'))
    return h * 3600 + m * 60 + s


U1, U4 = sec(D['meta']['contacts_ut']['U1']), sec(D['meta']['contacts_ut']['U4'])

# ---- interpolate the two sub-points and the two cap radii to SAMPLE steps ----
t = np.array([r['ut_sec'] for r in ROWS], float)


def col(f):
    return np.array([f(r) for r in ROWS], float)


tt = np.arange(U1, U4 + 1e-6, SAMPLE)


def I(y):
    return np.interp(tt, t, y)


SLA = I(col(lambda r: r['sub_solar'][0]))
SLO = np.degrees(I(np.unwrap(np.radians(col(lambda r: r['sub_solar'][1])))))
MLA = I(col(lambda r: r['sub_lunar'][0]))
MLO = np.degrees(I(np.unwrap(np.radians(col(lambda r: r['sub_lunar'][1])))))
RS = I(col(lambda r: r['cap_radius_sun_deg']))
RM = I(col(lambda r: r['cap_radius_moon_deg']))

lons = np.arange(-180.0, 180.0, GRID)
lats = np.arange(LAT1, LAT0, -GRID)
LON, LAT = np.meshgrid(lons, lats)
CL, SL = np.cos(np.radians(LAT)), np.sin(np.radians(LAT))


def cosarc(la0, lo0):
    return SL * np.sin(np.radians(la0)) + CL * np.cos(np.radians(la0)) * np.cos(np.radians(LON - lo0))


mins = np.zeros(LON.shape, np.float32)
for k in range(len(tt)):
    mins += ((cosarc(SLA[k], SLO[k]) >= np.cos(np.radians(RS[k])))
             & (cosarc(MLA[k], MLO[k]) >= np.cos(np.radians(RM[k]))))
mins *= SAMPLE / 60.0

kh = int(np.argmin(abs(tt - HAIRLINE_UT)))
hair = ((cosarc(SLA[kh], SLO[kh]) >= np.cos(np.radians(RS[kh])))
        & (cosarc(MLA[kh], MLO[kh]) >= np.cos(np.radians(RM[kh]))))
night = cosarc(SLA[kh], SLO[kh]) < np.cos(np.radians(RS[kh]))

land = globe.is_land(np.clip(LAT, -89.99, 89.99), np.clip(LON, -179.99, 179.99))


# ---- paint -------------------------------------------------------------------
def rgb(h):
    return np.array([int(h[i:i + 2], 16) / 255 for i in (1, 3, 5)])


img = np.ones(LON.shape + (3,)) * rgb(SEA)
img[land] = rgb(LAND)


def tint(mask, colour, a=1.0):
    img[mask] = img[mask] * (1 - a) + rgb(colour) * a


tint(night, '#2E4057', 0.16)
for lo, hi, colour, _ in LEVELS:
    tint((mins >= lo) & (mins < hi), colour, 0.92)
tint(hair, '#7a2d12', 0.95)

buf = io.BytesIO()
Image.fromarray((img * 255).astype(np.uint8)).save(buf, format='PNG', optimize=True)
b64 = base64.b64encode(buf.getvalue()).decode()

# ---- svg ---------------------------------------------------------------------
MW, MH = W - ML - MR, MAPH


def X(lon):
    return ML + (lon + 180.0) / 360.0 * MW


def Y(lat):
    return MT + (LAT1 - lat) / (LAT1 - LAT0) * MH


def esc(s):
    return s.replace('&', '&amp;').replace('<', '&lt;').replace('>', '&gt;')


def txt(x, y, s, size=11, fill=INK, anchor='middle', weight=None, halo=False):
    w = f' font-weight="{weight}"' if weight else ''
    h = (' stroke="#fff" stroke-width="3" stroke-linejoin="round" paint-order="stroke"'
         if halo else '')
    return (f'<text x="{x:.1f}" y="{y:.1f}" font-size="{size}" text-anchor="{anchor}" '
            f'fill="{fill}"{w}{h}>{esc(s)}</text>')


w_area = np.cos(np.radians(LAT))
w_area = w_area / w_area.sum()
frac_any = float(((mins > SAMPLE / 60.0) * w_area).sum())
median_min = float(np.median(mins[mins > 0]))

o = [f'<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 {W} {H}" width="{W}" height="{H}">',
     f'<rect width="{W}" height="{H}" fill="#fff"/>',
     txt(W / 2, 26, 'A quarter of the Earth was crossed by the selenelion band that morning. '
                    'Almost all of it for about five minutes — 2% of the eclipse.',
         13.5, INK, weight='700'),
     txt(W / 2, 45, 'Minutes with the risen Sun and the eclipsed Moon above opposite horizons, '
                    '10 December 2011, between first and last umbral contact.', 10.5, FAINT),
     txt(W / 2, 60, 'The dark line is the same band at a single minute, 14:03 UT — a hairline '
                    'along the day/night line. The shading is where that hairline swept.', 10.5, FAINT),
     f'<image x="{ML}" y="{MT}" width="{MW}" height="{MH}" '
     f'preserveAspectRatio="none" href="data:image/png;base64,{b64}"/>',
     f'<rect x="{ML}" y="{MT}" width="{MW}" height="{MH}" fill="none" stroke="{FAINT}"/>']

for lon in range(-180, 181, 60):
    o.append(f'<line x1="{X(lon):.1f}" y1="{MT}" x2="{X(lon):.1f}" y2="{MT+MH}" '
             f'stroke="#fff" stroke-opacity=".3" stroke-dasharray="2 5"/>')
    o.append(txt(X(lon), MT + MH + 15,
                 f'{abs(lon)}°{"" if lon == 0 else ("E" if lon > 0 else "W")}', 9, FAINT))
for lat in (-60, -30, 0, 30, 60):
    o.append(f'<line x1="{ML}" y1="{Y(lat):.1f}" x2="{ML+MW}" y2="{Y(lat):.1f}" '
             f'stroke="#fff" stroke-opacity=".3" stroke-dasharray="2 5"/>')
    o.append(txt(ML - 7, Y(lat) + 3.5,
                 f'{abs(lat)}°{"" if lat == 0 else ("N" if lat > 0 else "S")}', 9, FAINT, 'end'))

for name, la, lo, note, side in OBS:
    x, y = X(lo), Y(la)
    o.append(f'<circle cx="{x:.1f}" cy="{y:.1f}" r="3.6" fill="#fff" stroke="{RED}" stroke-width="1.8"/>')
o.append(f'<rect x="{X(IN_LON0):.1f}" y="{Y(IN_LAT1):.1f}" width="{X(IN_LON1)-X(IN_LON0):.1f}" '
         f'height="{Y(IN_LAT0)-Y(IN_LAT1):.1f}" fill="none" stroke="{RED}" stroke-width="1.4"/>')
o.append(txt(X(IN_LON0) - 6, Y(IN_LAT1) - 5, 'both observers', 9, RED, 'end', '700'))

# ---- inset: North America, same data, room for the labels --------------------
def IX(lon):
    return IN_X + (lon - IN_LON0) / (IN_LON1 - IN_LON0) * IN_W


def IY(lat):
    return IN_Y + (IN_LAT1 - lat) / (IN_LAT1 - IN_LAT0) * IN_H


i0 = int(np.argmin(abs(lats - IN_LAT1)))
i1 = int(np.argmin(abs(lats - IN_LAT0)))
j0 = int(np.argmin(abs(lons - IN_LON0)))
j1 = int(np.argmin(abs(lons - IN_LON1)))
buf2 = io.BytesIO()
Image.fromarray((img[i0:i1, j0:j1] * 255).astype(np.uint8)).resize(
    (int(IN_W * 3), int(IN_H * 3)), Image.LANCZOS).save(buf2, format='PNG', optimize=True)
o.append(f'<image x="{IN_X}" y="{IN_Y}" width="{IN_W}" height="{IN_H}" preserveAspectRatio="none" '
         f'href="data:image/png;base64,{base64.b64encode(buf2.getvalue()).decode()}"/>')
o.append(f'<rect x="{IN_X}" y="{IN_Y}" width="{IN_W}" height="{IN_H}" fill="none" stroke="{RED}" stroke-width="1.4"/>')
o.append(txt(IN_X, IN_Y - 8, 'The edge of the footprint runs between them', 10.5, INK, 'start', '700'))
for name, la, lo, note, side in OBS:
    x, y = IX(lo), IY(la)
    o.append(f'<circle cx="{x:.1f}" cy="{y:.1f}" r="5" fill="#fff" stroke="{RED}" stroke-width="2"/>')
    if side == 'below':
        o.append(txt(x, y + 18, name, 11.5, RED, 'middle', '700', halo=True))
        o.append(txt(x, y + 30, note, 10, RED, halo=True))
    else:
        o.append(txt(x, y - 20, name, 11.5, RED, 'middle', '700', halo=True))
        o.append(txt(x, y - 9, note, 10, RED, halo=True))

y = MT + MH + 40
lx = ML + 4
for _, _, colour, lab in LEVELS:
    o.append(f'<rect x="{lx}" y="{y-9}" width="15" height="11" fill="{colour}" rx="2"/>')
    o.append(txt(lx + 20, y, lab, 9.5, DIM, 'start'))
    lx += 128
o.append(f'<rect x="{lx}" y="{y-9}" width="15" height="11" fill="#7a2d12" rx="2"/>')
o.append(txt(lx + 20, y, 'the band at 14:03 UT', 9.5, DIM, 'start'))

cx = IN_X + IN_W + 30
cy = IN_Y + 6
for line, colour, wt in (
    (f'{frac_any*100:.0f} per cent of the planet spent some time in the band.', INK, '700'),
    (f'Where it happened at all the median was {median_min:.0f} minutes; the longest', DIM, None),
    (f'anywhere was {mins.max():.0f} minutes, on the polar rim. Santa Fe had two.', DIM, None),
    ('', DIM, None),
    ('The band is not a feature of this eclipse. Every umbral', DIM, None),
    ('eclipse has one, for the same reason: the Moon sits within', DIM, None),
    ('a degree and a half of the anti-solar point, both discs have', DIM, None),
    ('width, and the air lifts every image by half a degree.', DIM, None),
    ('', DIM, None),
    ('What is uncommon is the band falling somewhere with', INK, '700'),
    ('people, clear air and a camera — rarity of witnesses,', INK, '700'),
    ('not rarity of physics.', INK, '700'),
    ('', DIM, None),
    ('The window lengthens westward. The two men in the video', DIM, None),
    ('are 1,443 km apart and straddle the footprint’s eastern', DIM, None),
    ('edge — which is why one filmed a selenelion, and the', DIM, None),
    ('other filmed a dawn.', DIM, None),
):
    if line:
        o.append(txt(cx, cy, line, 10.5, colour, 'start', wt))
    cy += 14

o.append(txt(W / 2, H - 8, 'Upper limbs, 34.5′ horizon refraction, lunar parallax removed, sea-level '
                           'horizon, no terrain. Positions from JPL DE421 at 10-second steps; '
                           'scripts/selenelion_band_map.py', 9, FAINT))
o.append('</svg>')

os.makedirs(os.path.dirname(OUT), exist_ok=True)
open(OUT, 'w').write('\n'.join(o))
print(f'wrote {OUT} ({os.path.getsize(OUT)/1e3:.0f} kB)')
print(f'  any overlap: {frac_any*100:.2f}% of Earth by area; median {median_min:.1f} min; '
      f'max {mins.max():.0f} min')
for name, la, lo, _, _ in OBS:
    i, j = int(np.argmin(abs(lats - la))), int(np.argmin(abs(lons - lo)))
    print(f'  {name}: {mins[i, j]*60:.0f} s')