#!/usr/bin/env python3
"""The Moon's phase and the tilt of its lit side, from angles you measured.

The point of this script is what it does NOT import. There is no ephemeris
here, no orbital elements, no distances, no agency data file and no network
call -- just two directions in the sky and some spherical trigonometry. Feed it
the altitude and azimuth of the Moon and of the Sun, measured however you like,
and it returns the fraction of the disc that should be lit and the direction
its bright limb should face.

That matters for a specific argument. p.168 of Globe Deconstruction? asks
whether heliocentrism "requires us to believe in the dramatic bending of
sunlight". The answer is that the phase needs no bending at all, and the way to
show it is to write the calculation out and see that no bending term appears
anywhere in it -- and, while we are at it, that nothing institutional appears
either.

    The physics, in full:

      A ball lit by a distant source is lit on the hemisphere facing it.
      You see whatever part of that hemisphere faces you.

    which gives, with E the angle between Sun and Moon in your sky:

      phase angle   phi = 180 - E          (the Sun is far, so the third
                                            angle of the triangle is ~0.15 deg)
      lit fraction  f = (1 + cos phi)/2 = (1 - cos E)/2
      bright limb   points along the great circle from Moon toward Sun

    Straight lines throughout. There is nothing for refraction to do, and the
    only place the atmosphere can enter at all is in correcting the two
    altitudes you measured -- about 1-2 arcminutes at any working altitude,
    which moves the answers by hundredths of a degree. That correction is
    offered below as an option so its size can be seen rather than argued
    about.

Usage:

    python3 moon_phase_from_angles.py MOON_ALT MOON_AZ SUN_ALT SUN_AZ
    python3 moon_phase_from_angles.py --check     # the 16 May 2016 evening

How to get the two directions without buying anything: altitude from a
protractor with a thread and a weight on it, or a phone inclinometer; azimuth
from a compass, or from the shadow of a vertical stick at a known time. Fists
at arm's length are about 10 degrees each and will get the elongation to within
about five, which is enough to fix the lit fraction to a few per cent -- the
script prints that sensitivity so you can see what your own measurement is
worth.
"""

import math
import sys

DEG = math.pi / 180.0


def bennett(h_app):
    """Atmospheric refraction, degrees, at apparent altitude h_app (Bennett 1982)."""
    return (1.0 / math.tan((h_app + 7.31 / (h_app + 4.4)) * DEG)) / 60.0


def unit(alt, az):
    """A direction in the sky as a vector. x north, y east, z up."""
    a, z = alt * DEG, az * DEG
    return (math.cos(a) * math.cos(z), math.cos(a) * math.sin(z), math.sin(a))


def elongation(m, s):
    """Angle between the two directions -- the only measurement that sets the phase."""
    d = sum(a * b for a, b in zip(m, s))
    return math.acos(max(-1.0, min(1.0, d))) / DEG


def lit_fraction(E):
    """(1 - cos E)/2. The Sun is far enough that the phase angle is 180 - E."""
    return (1.0 - math.cos(E * DEG)) / 2.0


def limb_angle(m_alt, m_az, s_alt, s_az):
    """Position angle of the bright limb, degrees from straight up at the Moon.

    Positive toward increasing azimuth. This is the direction of the great
    circle from the Moon toward the Sun, which is where the lit side points and
    what the terminator is perpendicular to.
    """
    m, s = unit(m_alt, m_az), unit(s_alt, s_az)
    up = [0.0, 0.0, 1.0]
    dot = sum(a * b for a, b in zip(up, m))
    u = [up[i] - dot * m[i] for i in range(3)]          # "up" in the Moon's own frame
    n = math.sqrt(sum(c * c for c in u))
    u = [c / n for c in u]
    e = [u[1] * m[2] - u[2] * m[1],                     # "toward increasing azimuth"
         u[2] * m[0] - u[0] * m[2],
         u[0] * m[1] - u[1] * m[0]]
    dot = sum(a * b for a, b in zip(s, m))
    d = [s[i] - dot * m[i] for i in range(3)]
    n = math.sqrt(sum(c * c for c in d))
    d = [c / n for c in d]
    return math.atan2(sum(a * b for a, b in zip(d, e)),
                      sum(a * b for a, b in zip(d, u))) / DEG


def report(m_alt, m_az, s_alt, s_az, label=None):
    if label:
        print(f'\n{label}')
    print(f'  measured   Moon {m_alt:6.2f} deg up, azimuth {m_az:6.2f}')
    print(f'             Sun  {s_alt:6.2f} deg up, azimuth {s_az:6.2f}')

    E = elongation(unit(m_alt, m_az), unit(s_alt, s_az))
    f = lit_fraction(E)
    pa = limb_angle(m_alt, m_az, s_alt, s_az)
    side = 'right' if pa > 0 else 'left'
    print(f'\n  elongation           {E:7.2f} deg   <- the one angle the phase needs')
    print(f'  lit fraction         {f * 100:7.1f} %     = (1 - cos E)/2')
    print(f'  bright limb points   {abs(pa):7.2f} deg {side} of straight up')
    print(f'  so the terminator runs at right angles to that.')

    # what the atmosphere is worth here, as a number rather than an argument
    dm, ds = bennett(max(m_alt, -0.5)), bennett(max(s_alt, -0.5))
    E2 = elongation(unit(m_alt - dm, m_az), unit(s_alt - ds, s_az))
    pa2 = limb_angle(m_alt - dm, m_az, s_alt - ds, s_az)
    print(f'\n  refraction on the Moon {dm * 60:5.2f} arcmin, on the Sun {ds * 60:5.2f} arcmin')
    print(f'  undo both and the answers move by {abs(f - lit_fraction(E2)) * 100:.3f} '
          f'points and {abs(pa - pa2):.3f} deg')
    print(f'  -- which is the whole of the atmosphere\'s involvement in a phase.')

    # what a rough measurement is worth
    d = 5.0
    lo, hi = lit_fraction(E - d), lit_fraction(E + d)
    print(f'\n  measure E by fists (about +/-{d:.0f} deg) and the lit fraction is '
          f'{min(lo, hi) * 100:.0f}-{max(lo, hi) * 100:.0f}%')
    return E, f, pa


def check():
    """The evening of Example 1, 16 May 2016, central Scotland.

    The altitudes and azimuths below are the ones the page's table carries.
    They came from an ephemeris -- but only as a convenience, and the point of
    running them through this script is that it makes no use of one. Anyone who
    stood in that field with a protractor would have measured the same four
    numbers and got the same answers out.
    """
    print('Check: the evening of Example 1 (16 May 2016, central Scotland)')
    print('The four inputs are directions in the sky. Nothing else is used.')
    rows = [('16:00 BST', 4.1, 95.3, 41.3, 238.1),
            ('18:00 BST', 19.1, 121.0, 25.6, 266.6),
            ('19:00 BST', 25.3, 135.3, 17.4, 278.9),
            ('21:00 BST', 32.7, 168.0, 2.0, 303.1)]
    print(f"\n{'BST':>10} {'elongation':>11} {'lit':>7} {'bright limb':>14}")
    for label, ma, maz, sa, saz in rows:
        E = elongation(unit(ma, maz), unit(sa, saz))
        f = lit_fraction(E)
        pa = limb_angle(ma, maz, sa, saz)
        print(f'{label:>10} {E:10.1f}  {f * 100:6.1f}% {abs(pa):9.1f} '
              f'{"right" if pa > 0 else "left"}')
    print('\nThe page says the Moon was 78% lit and about 123 deg from the Sun')
    print('all evening. Both fall out of the four measured directions alone.')
    report(*rows[-1][1:], label='In full, for the last of those:')


if __name__ == '__main__':
    if len(sys.argv) == 2 and sys.argv[1] == '--check':
        check()
    elif len(sys.argv) == 5:
        report(*[float(x) for x in sys.argv[1:5]], label='From your measurements:')
    else:
        print(__doc__)
        sys.exit(1)