#!/usr/bin/env python3
"""What a wings-level, hands-off leg does, at the equator and away from it.

Section 4 claims that an aircraft holding zero bank with no lateral reference
follows a great circle, and that a great circle coincides with a parallel only
at the equator. This computes the consequence: how far south of its starting
parallel such an aircraft ends up, and how much its true heading has swung.

The result is the thing a home flight simulator can be made to reproduce --
X-Plane will write latitude, longitude and true heading to a file -- which is
a check on this arithmetic rather than on the Earth, since the simulator is a
globe by construction.

    python3 wings_level_drift.py
"""

import numpy as np
D,G=np.radians,np.degrees
R=6371.0088; KT=560*1.609344   # 560 mph in km/h

def unit(lat,lon):
    f,l=D(lat),D(lon)
    return np.array([np.cos(f)*np.cos(l),np.cos(f)*np.sin(l),np.sin(f)])

def gc_point(lat0,lon0,az,d_km):
    """Great-circle propagation: wings level, no lateral input."""
    f1,l1,a=D(lat0),D(lon0),D(az); s=d_km/R
    f2=np.arcsin(np.sin(f1)*np.cos(s)+np.cos(f1)*np.sin(s)*np.cos(a))
    l2=l1+np.arctan2(np.sin(a)*np.sin(s)*np.cos(f1), np.cos(s)-np.sin(f1)*np.sin(f2))
    return G(f2),(G(l2)+540)%360-180

print("Wings level, due east, no lateral reference, 560 mph, zero wind.")
print("Departure from the parallel you started on:\n")
print(f"{'start lat':>9} " + "".join(f"{h:>10}" for h in ['1 h','2 h','4 h','8 h']))
for lat0 in [0.0,10.0,30.0,45.0,60.0]:
    row=[]
    for hours in [1,2,4,8]:
        f2,l2=gc_point(lat0,0.0,90.0,KT*hours)
        # distance south of the starting parallel, along a meridian
        row.append((lat0-f2)*np.pi/180*R*0.621371)   # statute miles
    print(f"{lat0:9.1f} " + "".join(f"{v:9.1f} " for v in row) + " mi south of the parallel")

print("\nSame legs, as a latitude readout:")
for lat0 in [0.0,45.0]:
    print(f"  start {lat0:.1f} deg:", ", ".join(
        f"{hours}h -> {gc_point(lat0,0.0,90.0,KT*hours)[0]:8.4f}" for hours in [1,2,4,8]))

print("\nAnd the heading a wings-level aircraft ends up on (true):")
for lat0 in [0.0,45.0]:
    for hours in [4]:
        f2,l2=gc_point(lat0,0.0,90.0,KT*hours)
        # final azimuth of the great circle at that point
        f1,l1=D(lat0),0.0; F2,L2=D(f2),D(l2); dl=L2-l1
        y=np.sin(-dl)*np.cos(f1); x=np.cos(F2)*np.sin(f1)-np.sin(F2)*np.cos(f1)*np.cos(dl)
        back=(G(np.arctan2(y,x))+180)%360
        print(f"  start {lat0:.1f} deg, after {hours} h: heading {back:.2f} deg true")


# --- and what a flat map predicts for the same hands-off legs ---------------
#
# On a north-polar azimuthal-equidistant map a parallel is a circle of radius
# r0 = R*(90-lat) in degrees-as-distance, and flying straight is a straight
# line in the plane. A straight line is never a circle, so the departure
#     sqrt(r0^2 + d^2) - r0
# is positive at EVERY latitude. There is no radius at which it vanishes.
# The globe has exactly one such line, and it is the equator.

def flat_departure(lat0, d_km):
    r0 = R * np.pi/180 * (90.0 - lat0)
    return (np.hypot(r0, d_km) - r0) * 0.621371     # statute miles

if __name__ == '__main__':
    print("\n\nTHE SAME LEGS ON A POLE-CENTRED FLAT MAP\n")
    print(f"{'start lat':>9} " + "".join(f"{h:>10}" for h in ['1 h','2 h','4 h','8 h'])
          + "     globe at 4 h")
    for lat0 in [0.0, 30.0, 45.0, 60.0]:
        row = [flat_departure(lat0, KT*h) for h in (1,2,4,8)]
        g4 = (lat0 - gc_point(lat0,0.0,90.0,KT*4)[0]) * np.pi/180 * R * 0.621371
        print(f"{lat0:9.1f} " + "".join(f"{v:9.1f} " for v in row) + f"{g4:14.1f}")
    print("\n  The flat column has no zero at any latitude; the globe column has")
    print("  exactly one, and it is the line the chapter chose to fly.")