#!/usr/bin/env python3
"""Draw the ball test: what the ball and the Moon actually look like side by side.

Not an illustration. Every disc here is drawn from the geometry -- the lit
fraction from the elongation, the terminator as the projected half-ellipse, the
tilt as the position angle of the great circle toward the Sun -- for the evening
of Example 1, 16 May 2016 in central Scotland, at 21:00 BST.

Three discs, because the third is the objection:

  1. the Moon, elongation 125.2 deg
  2. a ball held IN LINE with it. Its elongation differs by 0.15 deg, the
     Earth-Moon distance seen from the Sun, so the two are the same picture
  3. a ball held 20 deg off, moved ALONG the Moon-Sun great circle so that its
     tilt is unchanged and only its elongation moves. It shows a visibly
     different phase -- which is what someone reports when they say the ball and
     the Moon do not match, and it is a way of running the test wrong. It also
     separates the two halves of the test: the tilt can match while the phase
     does not.

    python3 ball_test_figure.py
"""

import math
import os

DEG = math.pi / 180.0
HERE = os.path.dirname(os.path.abspath(__file__))
OUT = os.path.join(os.path.dirname(HERE), 'docs', 'moon-tilt', 'img', 'fig-ball-test.svg')

# 21:00 BST, 16 May 2016, central Scotland -- the last row of the page's table.
ELONG = 125.23          # Moon-Sun angle in the sky, degrees
PA = 59.73              # bright limb, degrees from straight up, toward the right
OFFSET = 0.1472         # Earth-Moon distance seen from the Sun
OFF_BALL = 20.0         # how far off the Moon the third ball is held

INK, DIM, RULE = '#181715', '#6E675B', '#E4DDCB'
LIT, DARK = '#f2e4b8', '#2b2b30'


def lit_fraction(elong):
    return (1.0 - math.cos(elong * DEG)) / 2.0


def disc(cx, cy, r, elong, pa, n=96):
    """A Moon-like disc, drawn as two arcs met end to end.

    Work in a basis (u, v) at the disc: u toward the Sun, v across it. In that
    basis the limb is the circle r*(cos t, sin t) and the terminator is the
    ellipse r*(s*cos t, sin t), where s = cos(phase angle) is SIGNED -- positive
    for a gibbous phase, when the terminator bulges away from the Sun, negative
    for a crescent, when it bulges toward it. The dark region is the anti-Sun
    half of the limb closed by the matching half of that ellipse, so it is drawn
    as one polygon and no arc flags have to be reasoned about.
    """
    phi = 180.0 - elong                       # phase angle at the body
    s = math.cos(phi * DEG)                   # signed: + gibbous, - crescent
    a = pa * DEG                              # bright limb, from up, toward the right
    ux, uy = math.sin(a), -math.cos(a)        # toward the Sun, screen coords (y down)
    vx, vy = -uy, ux                          # across it

    def pt(t, ku, kv):
        return (cx + r * (ku * math.cos(t) * ux + kv * math.sin(t) * vx),
                cy + r * (ku * math.cos(t) * uy + kv * math.sin(t) * vy))

    ts = [math.pi / 2 + i * math.pi / n for i in range(n + 1)]      # 90 deg -> 270 deg
    limb = [pt(t, 1.0, 1.0) for t in ts]
    term = [pt(t, s, 1.0) for t in reversed(ts)]
    poly = ' '.join(f'{x:.2f},{y:.2f}' for x, y in limb + term)

    out = [f'<circle cx="{cx:.2f}" cy="{cy:.2f}" r="{r:.2f}" fill="{LIT}"/>',
           f'<polygon points="{poly}" fill="{DARK}"/>',
           f'<circle cx="{cx:.2f}" cy="{cy:.2f}" r="{r:.2f}" fill="none" '
           f'stroke="{DIM}" stroke-width="0.8"/>',
           f'<line x1="{cx + r * ux * 1.10:.2f}" y1="{cy + r * uy * 1.10:.2f}" '
           f'x2="{cx + r * ux * 1.62:.2f}" y2="{cy + r * uy * 1.62:.2f}" '
           f'stroke="{INK}" stroke-width="1.6" marker-end="url(#ah)"/>']
    return '\n'.join(out)


def main():
    W, H, R = 760, 272, 62
    xs = [140, 380, 620]
    labels = [
        ('The Moon', f'{ELONG:.1f}&#176; from the Sun', f'{lit_fraction(ELONG)*100:.1f}% lit'),
        ('A ball held in line with it', f'{ELONG + OFFSET:.1f}&#176; from the Sun',
         f'{lit_fraction(ELONG + OFFSET)*100:.1f}% lit'),
        (f'A ball held {OFF_BALL:.0f}&#176; toward the Sun', f'{ELONG - OFF_BALL:.1f}&#176; from the Sun',
         f'{lit_fraction(ELONG - OFF_BALL)*100:.1f}% lit'),
    ]
    elongs = [ELONG, ELONG + OFFSET, ELONG - OFF_BALL]

    o = [f'<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 {W} {H}" width="{W}" height="{H}" '
         f'font-family="Inter,system-ui,sans-serif">',
         '<defs><marker id="ah" viewBox="0 0 8 8" refX="7" refY="4" markerWidth="6" '
         f'markerHeight="6" orient="auto"><path d="M0 0 L8 4 L0 8 z" fill="{INK}"/></marker></defs>',
         f'<rect width="{W}" height="{H}" fill="#fff"/>']

    for x, e, (t1, t2, t3) in zip(xs, elongs, labels):
        o.append(disc(x, 118, R, e, PA))
        o.append(f'<text x="{x}" y="215" text-anchor="middle" font-size="14" '
                 f'font-weight="700" fill="{INK}">{t1}</text>')
        o.append(f'<text x="{x}" y="234" text-anchor="middle" font-size="12.5" fill="{DIM}">{t2}</text>')
        o.append(f'<text x="{x}" y="251" text-anchor="middle" font-size="12.5" fill="{DIM}">{t3}</text>')

    o.append(f'<line x1="255" y1="40" x2="255" y2="196" stroke="{RULE}" stroke-width="1.5"/>')
    o.append(f'<line x1="500" y1="40" x2="500" y2="196" stroke="{RULE}" stroke-width="1.5"/>')
    o.append(f'<text x="{W/2:.0f}" y="26" text-anchor="middle" font-size="12.5" fill="{DIM}">'
             f'16 May 2016, 21:00 BST &#8212; arrows point along the sky toward the Sun</text>')
    o.append('</svg>')

    os.makedirs(os.path.dirname(OUT), exist_ok=True)
    with open(OUT, 'w') as f:
        f.write('\n'.join(o) + '\n')
    print(f'wrote {OUT}')
    for e, (t1, _, _) in zip(elongs, labels):
        print(f'  {t1:32s} elongation {e:7.2f}   lit {lit_fraction(e)*100:5.1f}%')


if __name__ == '__main__':
    main()