#!/usr/bin/env python3
"""Which route carries the most information, for the chapter's own experiment.

The chapter picks the equator, and picks it for the turn rate: one circuit is
360 degrees of turning on a flat surface and 0 on a globe. That is the right
instinct applied to the wrong quantity, because a turn rate is a rate, and a
flat map with a free scale parameter can always be rescaled to produce a
different one. What cannot be rescaled is a quantity that changes sign.

This script computes, for a set of real scheduled routes, the two things a
gyro reports at each end:

    the vertical component of earth rate,   15.041 * sin(phi)  deg/hr
    the tilt of the spin axis to local vertical,  the co-latitude, 90 - phi

against the three models actually on the table:

    globe          15.041 * sin(phi),  tilt = co-latitude
    spinning disc  15.041 everywhere,  tilt = 0 everywhere
                   (the disc's axis is normal to the plane at every point,
                   so the whole rate is vertical wherever you stand)
    the book       0 everywhere; it calls itself geocentric and does not spin

and alongside them the great-circle heading change, which is the curvature
discriminator from the closed-circuit argument -- so the trade-off between the
two is visible in one table rather than asserted.

No external data. Spherical earth, R = 6371 km, sidereal rate 15.041 deg/hr.

    python3 equator_route_choice.py
"""

import numpy as np

D, G = np.radians, np.degrees
R = 6371.0
OMEGA = 15.041  # deg/hr, sidereal

# ICAO code -> (latitude, longitude), airport reference points
AP = {
    'LHR': (51.4775, -0.4614),   'JFK': (40.6413, -73.7781),
    'ATL': (33.6407, -84.4277),  'SCL': (-33.3930, -70.7858),
    'EZE': (-34.8222, -58.5358), 'JNB': (-26.1367, 28.2411),
    'CPT': (-33.9715, 18.6021),  'SIN': (1.3644, 103.9915),
    'HEL': (60.3172, 24.9633),   'SYD': (-33.9399, 151.1753),
    'UIO': (-0.1292, -78.3576),  'PNK': (-0.1507, 109.4039),
}


def great_circle(p, q):
    """Initial azimuth, forward azimuth on arrival, and distance in km."""
    f1, l1, f2, l2 = D(p[0]), D(p[1]), D(q[0]), D(q[1])
    dl = l2 - l1
    a1 = np.arctan2(np.sin(dl) * np.cos(f2),
                    np.cos(f1) * np.sin(f2) - np.sin(f1) * np.cos(f2) * np.cos(dl))
    a2 = np.arctan2(np.sin(-dl) * np.cos(f1),
                    np.cos(f2) * np.sin(f1) - np.sin(f2) * np.cos(f1) * np.cos(dl)) + np.pi
    d = 2 * np.arcsin(np.sqrt(np.sin((f2 - f1) / 2) ** 2
                              + np.cos(f1) * np.cos(f2) * np.sin(dl / 2) ** 2)) * R
    return G(a1) % 360, G(a2) % 360, d


def earth_rate(lat):
    """Vertical component of earth rate at a latitude, deg/hr."""
    return OMEGA * np.sin(D(lat))


def report(routes):
    print(f"{'route':>9} {'lat A':>7} {'lat B':>7} {'rate A':>7} {'rate B':>7} "
          f"{'change':>7} {'tilt A':>7} {'tilt B':>7} {'km':>7} {'d hdg':>7}")
    for a, b in routes:
        p, q = AP[a], AP[b]
        r1, r2 = earth_rate(p[0]), earth_rate(q[0])
        h1, h2, d = great_circle(p, q)
        dh = (h2 - h1 + 540) % 360 - 180
        flip = ' <- sign reversal' if r1 * r2 < 0 else ''
        print(f"{a + '-' + b:>9} {p[0]:7.1f} {q[0]:7.1f} {r1:7.2f} {r2:7.2f} "
              f"{r2 - r1:7.2f} {90 - p[0]:7.1f} {90 - q[0]:7.1f} {d:7.0f} {dh:7.1f}{flip}")


def models():
    print(f"\n{'':>6} {'lat':>7} {'globe':>8} {'disc':>8} {'book':>7}")
    for code in ['HEL', 'LHR', 'JFK', 'SIN', 'JNB', 'SCL', 'CPT']:
        lat = AP[code][0]
        print(f"{code:>6} {lat:7.1f} {earth_rate(lat):8.2f} {OMEGA:8.2f} {0.0:7.2f}")


def margin():
    """The floor is bias stability, not quantisation."""
    lsb = 0.0039 * 3600          # deg/hr, one count of a 0.0039 deg/s output
    print(f"\n  one output count      {lsb:8.3f} deg/hr")
    for hours in (1, 4, 11):
        n = 50 * 3600 * hours    # 50 Hz
        print(f"  averaged over {hours:2d} h   {lsb / np.sqrt(n):8.4f} deg/hr"
              f"   ({n:,} samples at 50 Hz)")
    print(f"  navigation-grade bias stability   0.003 to 0.010 deg/hr"
          f"  <- the real floor; it does not average down")
    print(f"  signal on LHR-CPT                {abs(earth_rate(AP['CPT'][0]) - earth_rate(AP['LHR'][0])):7.2f} deg/hr")


if __name__ == '__main__':
    print(__doc__.split('\n\n')[0])
    print("\nROUTES\n")
    report([('JFK', 'LHR'), ('UIO', 'PNK'), ('ATL', 'SCL'),
            ('JFK', 'EZE'), ('LHR', 'JNB'), ('LHR', 'CPT'), ('HEL', 'SIN')])
    print("\nWHAT EACH MODEL PREDICTS FOR THE VERTICAL-AXIS READING (deg/hr)")
    models()
    print("\nMARGIN")
    margin()