"""Distance models for the latency test.

Every model returns the shortest possible path length (km) between two
lat/lon points *within that model*. On a flat map that is the straight line
across the plane: a cable can bend, detour or loop, but it can never be
shorter than the straight line, so these are hard lower bounds.

Scale conventions (the ones flat-earth maps themselves use):
  globe     sphere, mean radius 6371.0 km, great-circle distance
  ae        azimuthal equidistant centred on the North Pole (the "Gleason"
            map). Pole-to-anywhere distances are true, so 1 degree of
            latitude = 111.2 km everywhere, equator 10,008 km from the pole.
  dualpole  azimuthal equidistant centred on 0N 0E (both poles are points,
            10,008 km above/below the centre). Distances from 0,0 are true.
  mercator  flat Mercator rectangle, true scale on the equator, cut along
            the antimeridian (180E/W) like every wall map.
"""
import math

R = 6371.0
C_VACUUM = 299_792.458          # km/s
N_FIBER = 1.4682                # group index, standard single-mode fibre @1550nm
C_FIBER = C_VACUUM / N_FIBER    # ~204,190 km/s

MODELS = ["globe", "ae", "dualpole", "mercator"]


def globe(lat1, lon1, lat2, lon2):
    p1, p2 = math.radians(lat1), math.radians(lat2)
    dp, dl = p2 - p1, math.radians(lon2 - lon1)
    a = math.sin(dp / 2) ** 2 + math.cos(p1) * math.cos(p2) * math.sin(dl / 2) ** 2
    return 2 * R * math.asin(min(1.0, math.sqrt(a)))


def ae_xy(lat, lon):
    r = R * math.radians(90 - lat)
    t = math.radians(lon)
    return r * math.sin(t), -r * math.cos(t)


def dualpole_xy(lat, lon, lon0=0.0):
    p, l = math.radians(lat), math.radians(lon - lon0)
    cosc = max(-1.0, min(1.0, math.cos(p) * math.cos(l)))
    c = math.acos(cosc)
    k = 1.0 if c < 1e-9 else c / math.sin(c)
    return R * k * math.cos(p) * math.sin(l), R * k * math.sin(p)


def mercator_xy(lat, lon, cut=180.0):
    lat = max(-85.0, min(85.0, lat))
    x = R * math.radians((lon - cut) % 360.0)
    y = R * math.log(math.tan(math.pi / 4 + math.radians(lat) / 2))
    return x, y


def _plane(f, a, b, **kw):
    x1, y1 = f(*a, **kw)
    x2, y2 = f(*b, **kw)
    return math.hypot(x2 - x1, y2 - y1)


def distance(model, a, b):
    """a, b = (lat, lon). Returns km."""
    if model == "globe":
        return globe(*a, *b)
    if model == "ae":
        return _plane(ae_xy, a, b)
    if model == "dualpole":
        return _plane(dualpole_xy, a, b)
    if model == "mercator":
        return _plane(mercator_xy, a, b)
    raise ValueError(model)


def floor_ms(km, speed=C_FIBER):
    """Minimum possible round-trip time in ms for a path of `km` each way."""
    return 2 * km / speed * 1000


if __name__ == "__main__":
    syd, akl, la, sel = (-33.87, 151.21), (-36.85, 174.76), (34.05, -118.24), (37.57, 126.98)
    for name, a, b in [("SYD-AKL", syd, akl), ("SYD-LAX", syd, la), ("SEL-SYD", sel, syd)]:
        row = [f"{m}={distance(m, a, b):8.0f}km/{floor_ms(distance(m, a, b)):6.1f}ms" for m in MODELS]
        print(name, *row)