#!/usr/bin/env python3
"""The refraction figure: what the air does, measured against everyday sights.

Two panels, drawn to the same angular scale so they can be compared by eye.

Left  -- every sunset anyone has watched. Horizon refraction is 34.5 arcmin
         and the Sun is 32 arcmin wide, so when the disc appears to sit on
         the horizon it is already wholly below it. The lower limb is lifted
         more than the upper, which squashes the disc by about 6 arcmin --
         the ovalling in everybody's sunset photographs.

Right -- the same air on the morning of 10 December 2011. Sun and Moon are
         both within a degree of the horizon, on opposite sides, and both are
         lifted by very nearly the same amount. This is the point the page
         needs: the selenelion is not bought by bending sunlight especially.
         It is bought twice over, once at each end, by an amount smaller than
         the horizon figure that has been in the tables for a century.

Numbers come from scripts/selenelion_2011.py (JPL DE421) and from Bennett
(1982) for the refraction itself; nothing here is drawn freehand.

    python3 refraction_figure.py
"""

import os, sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))

OUT = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))),
                   'docs', 'selenelion', 'img', 'fig-refraction-lift.svg')

W, H = 880, 520
INK, DIM, FAINT = '#2E4057', '#6E675B', '#a09a8c'
SUN, SUNDK, RED = '#d9a441', '#a8791a', '#c0392b'
SKY, GROUND, RULE = '#eaf1f8', '#dfe6d8', '#cbd6e2'

PXDEG = 112.0          # pixels per degree of altitude
TOP, HY, GND = 62.0, 196.0, 344.0     # panel top, horizon line, ground bottom


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


def txt(x, y, s, size=11, fill=INK, anchor='middle', weight=None):
    w = f' font-weight="{weight}"' if weight else ''
    return (f'<text x="{x:.1f}" y="{y:.1f}" font-size="{size}" text-anchor="{anchor}" '
            f'fill="{fill}"{w}>{esc(s)}</text>')


def panel(x0, w, title, sub):
    """Sky above the horizon, ground below, with room for the true positions."""
    return [f'<rect x="{x0}" y="{TOP}" width="{w}" height="{HY-TOP:.0f}" fill="{SKY}" stroke="{RULE}" rx="4"/>',
            f'<rect x="{x0}" y="{HY:.0f}" width="{w}" height="{GND-HY:.0f}" fill="{GROUND}" stroke="{RULE}" rx="4"/>',
            f'<line x1="{x0}" y1="{HY:.1f}" x2="{x0+w}" y2="{HY:.1f}" stroke="{DIM}" stroke-width="1.6"/>',
            txt(x0 + w / 2, 40, title, 12.5, INK, weight='700'),
            txt(x0 + w / 2, 54, sub, 10.5, FAINT),
            txt(x0 + 8, HY - 7, 'true horizon', 9.5, DIM, 'start'),
            txt(x0 + 8, HY + 15, 'below the horizon', 9.5, FAINT, 'start')]


def disc(cx, cy, r, fill, stroke, squash=1.0):
    return (f'<ellipse cx="{cx:.1f}" cy="{cy:.1f}" rx="{r:.1f}" ry="{r*squash:.1f}" '
            f'fill="{fill}" stroke="{stroke}" stroke-width="1.2"/>')


def ghost(cx, cy, r, stroke, squash=1.0):
    return (f'<ellipse cx="{cx:.1f}" cy="{cy:.1f}" rx="{r:.1f}" ry="{r*squash:.1f}" '
            f'fill="none" stroke="{stroke}" stroke-width="1.1" stroke-dasharray="3 3"/>')


def lift_arrow(x, y_from, y_to, colour, label):
    out = [f'<line x1="{x:.1f}" y1="{y_from:.1f}" x2="{x:.1f}" y2="{y_to+6:.1f}" '
           f'stroke="{colour}" stroke-width="1.5"/>',
           f'<path d="M{x:.1f},{y_to:.1f} L{x-4:.1f},{y_to+8:.1f} L{x+4:.1f},{y_to+8:.1f} Z" fill="{colour}"/>']
    if label:
        out.append(txt(x + 7, (y_from + y_to) / 2 + 4, label, 10.5, colour, 'start', '600'))
    return out


def build():
    import selenelion_2011 as S

    R0 = S.bennett(0.0) * 60.0                  # horizon refraction, arcmin
    sd_sun = 16.0 / 60.0
    R_up = S.bennett(2 * sd_sun) * 60.0
    squash = R0 - R_up                          # arcmin of flattening
    sq_ratio = 1.0 - (squash / 60.0) / (2 * sd_sun)

    r = S.state(S.at_mst(7, 3))
    m_geo, m_app, s_geo, s_app = r['m_alt'], r['m_app'], r['s_alt'], r['s_app']
    sdm, sds = r['sd_moon'], r['sd_sun']

    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, 22, 'What the air does, and where you have already watched it doing it',
             13.5, INK, weight='700')]

    # ---- left: the ordinary sunset ---------------------------------------
    o += panel(28, 384, 'Every sunset you have watched',
               'the disc appears to rest on the horizon')
    cx, rp = 172.0, sd_sun * PXDEG
    ty = HY + (R0 / 60.0) * PXDEG
    o.append(ghost(cx, ty, rp, SUNDK))
    o.append(txt(cx, ty + rp + 14, 'where it actually is', 10, SUNDK))
    o.append(disc(cx, HY - rp * sq_ratio, rp, SUN, SUNDK, squash=sq_ratio))
    o.append(txt(cx, HY - rp * sq_ratio * 2 - 10, 'where you see it', 10.5, SUNDK, weight='600'))
    o += lift_arrow(cx + rp + 30, ty, HY - rp * sq_ratio, RED, f'{R0:.1f}′')

    # ---- right: the morning in question ----------------------------------
    o += panel(468, 384, '10 December 2011, 07:03 MST, Santa Fe',
               'Sun rising in the east, eclipsed Moon setting in the west')
    for cxx, geo, app, sd, fill, edge, name in [
            (596.0, s_geo, s_app, sds, SUN, SUNDK, 'Sun, east'),
            (776.0, m_geo, m_app, sdm, '#6b3a33', '#8d8677', 'Moon, west')]:
        rr = sd * PXDEG
        o.append(ghost(cxx, HY - geo * PXDEG, rr, edge))
        o.append(disc(cxx, HY - app * PXDEG, rr, fill, edge))
        o.append(txt(cxx, HY - app * PXDEG - rr - 9, name, 10.5, edge, weight='600'))
        o += lift_arrow(cxx + rr + 26, HY - geo * PXDEG, HY - app * PXDEG, RED,
                        f'{(app-geo)*60:.0f}′')

    # ---- the reading, below the panels rather than over them -------------
    y = GND + 26
    for line in [f'The Sun is {2*sd_sun*60:.0f}′ across and the lift at the horizon is {R0:.1f}′.',
                 'So the disc you watch touch the horizon is already wholly below it,',
                 f'and its lower limb is lifted {squash:.1f}′ more than its upper — the',
                 'flattening in every sunset photograph anyone has taken.']:
        o.append(txt(220, y, line, 10.8, DIM)); y += 15
    y = GND + 26
    for line in ['Dashed is where each body really is; solid is where it is seen.',
                 f'Both are within a degree of the horizon and both are lifted — {(s_app-s_geo)*60:.0f}′',
                 f'and {(m_app-m_geo)*60:.0f}′. The same air, doing the same thing at each end,',
                 'in the amount the tables have carried for a century.']:
        o.append(txt(660, y, line, 10.8, DIM)); y += 15

    o.append(txt(W / 2, H - 30, 'Both panels are drawn to the same angular scale, so the lifts '
                                'can be compared by eye.', 10, FAINT))
    o.append(txt(W / 2, H - 14, f'Refraction from Bennett (1982); positions from JPL DE421 for '
                                f'Santa Fe at 2,130 m.', 10, FAINT))
    o.append('</svg>')

    with open(OUT, 'w') as f:
        f.write('\n'.join(o))
    print(f'wrote {OUT}')
    print(f'  horizon refraction {R0:.2f}′, upper limb {R_up:.2f}′, squash {squash:.2f}′')
    print(f'  07:03  Sun {s_geo:+.3f} -> {s_app:+.3f} ({(s_app-s_geo)*60:.1f}′);'
          f' Moon {m_geo:+.3f} -> {m_app:+.3f} ({(m_app-m_geo)*60:.1f}′)')


if __name__ == '__main__':
    build()