#!/usr/bin/env python3
"""The selenelion band, drawn against the twilight bands it lives inside.

The companion to fig-selenelion-band. Same map, same morning, one change of
frame: longitude is measured from the sub-solar point rather than Greenwich, so
the Sun stays put and the Earth turns underneath. In that frame the terminator
and the three twilight rings are fixed curves, each observer is a horizontal
line travelling out of the night and into the day, and the selenelion band is a
thin strip hugging the sunrise edge.

Which is the point. A twilight ring is the set of places where the Sun is a
given distance below the horizon; the selenelion band is the set of places where
the Sun is up and the eclipsed Moon has not yet set. They are the same kind of
object -- a strip along the day/night line, swept west at fifteen degrees an
hour -- and the selenelion band is the narrower of them. Both men in the video
filmed the wide grey one arriving and called it the Sun.

Reads docs/selenelion/data/eclipse-2011-12-10-3d.json. Land outlines are not
drawn: in a Sun-fixed frame the continents slide, and drawing them at one
instant would imply they stand still.

    python3 twilight_band_map.py
"""

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

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-twilight-bands.svg')

INK, DIM, FAINT, RED = '#2E4057', '#6E675B', '#a09a8c', '#c0392b'
GOLD = '#c1892a'

LAT1, LAT0 = 84.0, -66.0
GRID = 0.2
W, H = 980, 660
ML, MR, MT = 52.0, 16.0, 92.0
MAPH = 226.0
IN_X, IN_Y, IN_W, IN_H = 52.0, 386.0, 372.0, 196.0     # dawn-side zoom
IN_L0, IN_L1, IN_B0, IN_B1 = -100.0, -60.0, 22.0, 52.0

# solar depression bands, outward from the sunrise line
TWILIGHT = [(0.0, '#f7f4ec', 'day'),
            (6.0, '#d9dde4', 'civil twilight'),
            (12.0, '#9ba6b5', 'nautical'),
            (18.0, '#5d6b7d', 'astronomical'),
            (90.0, '#33404f', 'night')]

OBS = [('Santa Fe', 35.6870, -105.9378, -7, [('06:30', 'nautical'), ('06:47', 'the p.169 frame'),
                                             ('07:02', 'sunrise')]),
       ('Cahokia', 38.6606, -90.0619, -6, [('06:30', 'he starts filming'), ('06:55', 'a bite on the Moon'),
                                           ('07:04', '“hardly see anything”')])]

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'])
t = np.array([r['ut_sec'] for r in ROWS], float)


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


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

sub_lat = float(np.mean(SLA))                       # -22.9, near enough fixed all morning

lons = np.arange(-180.0, 180.0, GRID)               # longitude FROM the sub-solar point
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))


arc_sun = np.degrees(np.arccos(np.clip(cosarc(sub_lat, 0.0), -1, 1)))

# the selenelion band, in the same Sun-fixed frame, unioned over the umbral phase
band = np.zeros(LON.shape, bool)
for i, r in enumerate(ROWS):
    if not (U1 <= r['ut_sec'] <= U4):
        continue
    rel_moon_lon = math.degrees(MLO_u[i] - SLO_u[i])
    rel_moon_lon = (rel_moon_lon + 180.0) % 360.0 - 180.0
    band |= ((arc_sun <= RS[i])
             & (np.degrees(np.arccos(np.clip(cosarc(MLA[i], rel_moon_lon), -1, 1))) <= RM[i]))


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


img = np.zeros(LON.shape + (3,))
sunrise_arc = float(np.mean(RS))
for dep, colour, _ in TWILIGHT:
    img[arc_sun >= sunrise_arc + dep] = rgb(colour)
img[arc_sun < sunrise_arc] = rgb(TWILIGHT[0][1])
img[band] = rgb('#7a2d12')

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

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>')


band_w = [r['selenelion_band_width_deg'] for r in ROWS if U1 <= r['ut_sec'] <= U4]

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, 28, 'The selenelion band is a twilight band. It is just a narrower one.',
         14, INK, weight='700'),
     txt(W / 2, 48, 'The same morning, with longitude measured from the Sun instead of from '
                    'Greenwich — so the Sun holds still and the Earth turns underneath.', 10.5, FAINT),
     txt(W / 2, 63, 'Every band here is the same kind of thing: a strip along the day/night line, '
                    'sweeping west at fifteen degrees an hour.', 10.5, FAINT),
     f'<image x="{ML}" y="{MT}" width="{MW}" height="{MH}" preserveAspectRatio="none" '
     f'href="data:image/png;base64,{b64}"/>',
     f'<rect x="{ML}" y="{MT}" width="{MW}" height="{MH}" fill="none" stroke="{FAINT}"/>']

for lat in (-60, -30, 0, 30, 60):
    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'))
o.append(txt(X(0), MT + MH + 15, 'noon', 9.5, FAINT))
o.append(txt(X(-90), MT + MH + 15, 'dawn side', 9.5, FAINT))
o.append(txt(X(90), MT + MH + 15, 'dusk side', 9.5, FAINT))
o.append(txt(X(-175), MT + MH + 15, 'midnight', 9.5, FAINT, 'start'))
o.append(txt(X(175), MT + MH + 15, 'midnight', 9.5, FAINT, 'end'))

# each observer is a horizontal line: fixed latitude, longitude-from-the-Sun rising 15deg/hour
def rel_lon(lo, ut):
    i = int(np.argmin(abs(t - ut)))
    return (math.degrees(math.radians(lo) - SLO_u[i]) + 180.0) % 360.0 - 180.0


for name, la, lo, off, ticks in OBS:
    y = Y(la)
    o.append(f'<line x1="{ML}" y1="{y:.1f}" x2="{ML+MW}" y2="{y:.1f}" stroke="{RED}" '
             f'stroke-width="1.1" stroke-dasharray="3 3" stroke-opacity="0.75"/>')
o.append(f'<rect x="{X(IN_L0):.1f}" y="{Y(IN_B1):.1f}" width="{X(IN_L1)-X(IN_L0):.1f}" '
         f'height="{Y(IN_B0)-Y(IN_B1):.1f}" fill="none" stroke="{RED}" stroke-width="1.4"/>')
o.append(txt(X(IN_L0) - 8, Y(IN_B0) + 14, 'both observers', 9.5, RED, 'end', '700', halo=True))

# ---- the dawn-side zoom, where the whole story happens in ten degrees ---------
def IX(l):
    return IN_X + (l - IN_L0) / (IN_L1 - IN_L0) * IN_W


def IY(b):
    return IN_Y + (IN_B1 - b) / (IN_B1 - IN_B0) * IN_H


i0 = int(np.argmin(abs(lats - IN_B1)))
i1 = int(np.argmin(abs(lats - IN_B0)))
j0 = int(np.argmin(abs(lons - IN_L0)))
j1 = int(np.argmin(abs(lons - IN_L1)))
b2 = io.BytesIO()
Image.fromarray((img[i0:i1, j0:j1] * 255).astype(np.uint8)).resize(
    (int(IN_W * 3), int(IN_H * 3)), Image.NEAREST).save(b2, 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(b2.getvalue()).decode()}"/>')
o.append(f'<rect x="{IN_X}" y="{IN_Y}" width="{IN_W}" height="{IN_H}" fill="none" '
         f'stroke="{RED}" stroke-width="1.4"/>')
o.append(txt(IN_X, IN_Y - 8, 'the last forty minutes, magnified', 10.5, INK, 'start', '700'))

for name, la, lo, off, ticks in OBS:
    y = IY(la)
    up = name == 'Cahokia'
    o.append(f'<line x1="{IN_X}" y1="{y:.1f}" x2="{IN_X+IN_W}" y2="{y:.1f}" stroke="{RED}" '
             f'stroke-width="1.3" stroke-dasharray="4 3"/>')
    o.append(txt(IN_X + 5, y + (-7 if up else 13), name, 10.5, RED, 'start', '700', halo=True))
    for k, (hhmm, note) in enumerate(ticks):
        hh, mm = (int(v) for v in hhmm.split(':'))
        x = IX(rel_lon(lo, (hh - off) * 3600 + mm * 60))
        o.append(f'<circle cx="{x:.1f}" cy="{y:.1f}" r="3.6" fill="#fff" stroke="{RED}" '
                 f'stroke-width="1.7"/>')
        ty = y - 17 - (k % 2) * 22 if up else y + 21 + (k % 2) * 22
        o.append(f'<line x1="{x:.1f}" y1="{y:.1f}" x2="{x:.1f}" y2="{ty + (4 if up else -5):.1f}" '
                 f'stroke="{RED}" stroke-width="0.9" stroke-opacity="0.55"/>')
        o.append(txt(x, ty, hhmm, 9.5, RED, 'middle', '700', halo=True))
        if note:
            o.append(txt(x, ty + (-10 if up else 10), note, 8.5, DIM, 'middle', halo=True))

y = MT + MH + 34
lx = ML + 2
for dep, colour, lab in TWILIGHT:
    o.append(f'<rect x="{lx}" y="{y-9}" width="15" height="11" fill="{colour}" '
             f'stroke="{FAINT}" stroke-width="0.5" rx="2"/>')
    o.append(txt(lx + 20, y, lab, 9.5, DIM, 'start'))
    lx += 118
o.append(f'<rect x="{lx}" y="{y-9}" width="15" height="11" fill="#7a2d12" rx="2"/>')
o.append(txt(lx + 20, y, 'the selenelion band', 9.5, DIM, 'start'))

cx, cy = IN_X + IN_W + 34, IN_Y + 4
for line, colour, wt in (
    (f'Civil twilight is a ring 6° of arc wide.', INK, '700'),
    (f'The selenelion band, the same morning, runs', INK, '700'),
    (f'{min(band_w):.1f}° to {max(band_w):.1f}° — a quarter to a half of it,', INK, '700'),
    ('tucked against the sunrise edge.', INK, '700'),
    ('', DIM, None),
    ('Both men spent their whole footage in the wide', DIM, None),
    ('grey band, walking toward the narrow dark one.', DIM, None),
    ('That is what “the sky is getting brighter” is:', DIM, None),
    ('civil twilight arriving. The Sun stays below the', DIM, None),
    ('horizon the entire time, which is why neither of', DIM, None),
    ('them ever films it.', DIM, None),
    ('', DIM, None),
    ('And the narrow band is not a different kind of', DIM, None),
    ('thing. It is the same geometry with the Moon', DIM, None),
    ('added — the strip where the Sun has cleared one', DIM, None),
    ('horizon before the eclipsed Moon has left the', DIM, None),
    ('other. The seam on the dusk side is the same', DIM, None),
    ('band with the roles swapped.', DIM, None),
):
    if line:
        o.append(txt(cx, cy, line, 10.5, colour, 'start', wt))
    cy += 14

o.append(txt(W / 2, H - 10, 'Twilight rings are solar depression of 0°, 6°, 12° and 18° from the '
                            'sub-solar point; the band is upper limbs with 34.5′ refraction and '
                            'lunar parallax removed. JPL DE421. scripts/twilight_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'  civil ring 6.0° wide; selenelion band {min(band_w):.2f}°–{max(band_w):.2f}°')
for name, la, lo, off, ticks in OBS:
    for hhmm, _ in ticks:
        hh, mm = (int(v) for v in hhmm.split(':'))
        ut = (hh - off) * 3600 + mm * 60
        i = int(np.argmin(abs(t - ut)))
        rel = (math.degrees(math.radians(lo) - SLO_u[i]) + 180) % 360 - 180
        a = math.degrees(math.acos(max(-1, min(1,
            math.sin(math.radians(la)) * math.sin(math.radians(sub_lat))
            + math.cos(math.radians(la)) * math.cos(math.radians(sub_lat))
            * math.cos(math.radians(rel - 0.0))))))
        print(f'  {name:9s} {hhmm}  rel lon {rel:+7.2f}°  arc from sub-solar {a:6.2f}°  '
              f'(sun {a - sunrise_arc:+.2f}° below the sunrise ring)')