#!/usr/bin/env python3
"""The curve at 118,000 feet: how much the horizon should bow in a picture, and what the
book's three videos show.

Geometry. From height h above a sphere of radius R the horizon lies a constant angle
    dip = arccos(R / (R + h))
below eye level in every direction: a circle of constant depression. Through a
rectilinear lens of focal length f (pixels) pointed at the horizon, a point of that circle
at azimuth a from the optical axis lands at x = f tan a, y = -f tan(dip)/cos a, so across
a frame of half-width w the edges sit lower than the centre by
    bow = f tan(dip) (1/cos(phi) - 1),   phi = atan(w/f)  = tan(dip) tan(phi/2) w
which barely depends on camera pitch (checked numerically below). On a flat plane dip = 0
and the bow is exactly zero through any lens, at any tilt.

Through-centre rule. Any lens whose distortion is radially symmetric about the optical
axis maps a straight line through the image centre to a straight line. A flat plane's
horizon is a great circle, so wherever it passes through the centre it must image straight;
a globe's horizon is a small circle and still bows there, by the formula above. A spinning
camera sweeps the horizon through the centre in every orientation, which is how the
GoFast frames are tested (measure_through_centre).

    python3 horizon_bow.py             -> the tables
    python3 horizon_bow.py --figure    -> docs/balloon-curve/img/fig-bow-vs-fov.svg
    python3 horizon_bow.py --measure <frames dir>   -> re-run the GoFast through-centre
                                         measurement on 10 fps frames (a0001.jpg ...)
"""
import math
import os
import sys

R = 6371.0                       # km
HERE = os.path.dirname(os.path.abspath(__file__))
OUT = os.path.join(HERE, "..", "docs", "balloon-curve", "img", "fig-bow-vs-fov.svg")


def dip_deg(h_km):
    return math.degrees(math.acos(R / (R + h_km)))


def bow_px(dip, hfov_deg, halfwidth_px):
    """Edge-below-centre bow of the horizon, pixels, level rectilinear lens."""
    return math.tan(math.radians(dip)) * math.tan(math.radians(hfov_deg / 4)) * halfwidth_px


def bow_px_fisheye_level(dip, half_chord_px, f_px):
    """Edge-below-centre bow for a LEVEL camera through an equidistant (fisheye) lens of focal
    length f (r = f * theta): the horizon's image row at x = half_chord minus its row at the
    centre column. This is the level-camera edge figure only; unlike the rectilinear case it
    depends on pitch, and it is NOT the through-centre sagitta (see through_centre_sagitta)."""
    d = math.radians(dip)
    def y(a):
        th = math.acos(math.cos(d) * math.cos(a))
        return f_px * th * math.sin(d) / math.sin(th) if th > 0 else f_px * d
    return y(half_chord_px / f_px) - y(0)


def through_centre_sagitta(dip, chord_px, f_px, exact_equidistant=False):
    """Sagitta of the horizon over a chord of length L through the centre of the frame, for a
    camera pointed at the horizon. Near the centre every radially symmetric projection has the
    same scale f (px per radian), and the horizon is a small circle of geodesic curvature
    tan(dip), so its image curvature is tan(dip)/f and the sagitta is tan(dip) * L^2 / (8 f),
    whatever the projection. exact_equidistant=True instead uses the closed form for an
    equidistant lens: the horizon's angular gap from the tangent great circle at angular
    distance a along it is dip - asin(sin(dip) cos(a)), with a = (L/2)/f."""
    d = math.radians(dip)
    if exact_equidistant:
        a = (chord_px / 2) / f_px
        return f_px * (d - math.asin(math.sin(d) * math.cos(a)))
    return math.tan(d) * chord_px ** 2 / (8 * f_px)


def bow_px_pitched(dip, hfov_deg, pitch_deg, halfwidth_px=640):
    """The same by direct projection with the camera pitched up by pitch_deg: the horizon's
    image row at the frame edge minus its row at the centre. For the numerical check that
    pitch does not matter."""
    import numpy as np
    d = math.radians(dip); phi = math.radians(hfov_deg / 2); t = math.radians(pitch_deg)
    f = halfwidth_px / math.tan(phi)
    xs, ys = [], []
    for a in np.linspace(-1.5, 1.5, 6001):
        v = np.array([math.cos(d) * math.cos(a), math.cos(d) * math.sin(a), -math.sin(d)])
        xc = v[0] * math.cos(t) + v[2] * math.sin(t); zc = -v[0] * math.sin(t) + v[2] * math.cos(t); yc = v[1]
        if xc <= 0:
            continue
        xs.append(f * yc / xc); ys.append(-f * zc / xc)
    xs, ys = np.array(xs), np.array(ys); o = np.argsort(xs)
    y_at = lambda x: float(np.interp(x, xs[o], ys[o]))
    return (y_at(halfwidth_px) + y_at(-halfwidth_px)) / 2 - y_at(0)


# --- the book's three sources, as measured (frames at 1280 px wide unless stated) ---
RENDER = dict(alt_km=36.0, width_px=1210, eye_row=337, horizon_centre_row=390, horizon_edge_row=422,
              dip_px=53, bow_px=32)            # p. 235 turned upright (1210 x 677); the render's own
                                               # "Eye level" line, and the ground/glow edge at the centre and at both edges
RENDER_PX_PER_DEG = RENDER["dip_px"] / dip_deg(RENDER["alt_km"])
RENDER_F = RENDER_PX_PER_DEG * 180 / math.pi
RENDER_HFOV = 2 * math.degrees(math.atan((RENDER["width_px"] / 2) / RENDER_F))

# Focozz compilation: top edge of the limb glow, quadratic fit across 32 column bins,
# edge-below-centre in px on 1280-wide frames (negative = edges higher)
FOCOZZ_BOW = {"clip 1, t=60 s": -3.2, "clip 2, t=700 s": -8.8, "the p. 234 frame, t=921 s": 4.2}
GLOW_TOP_TANGENT_KM = 22.0        # the glow-top proxy sits about this high; its depression is smaller than the surface's

# GoFast apogee sequence (172-293 s of the re-upload), frames where the limb passes within
# 40 px of the frame centre against black space: parabola sagitta toward space, px, over the
# chord, and the chord's orientation. Six shown on the page; the clean set is 26 frames.
GOFAST_CENTRE = [(213.2, 16.5, 73), (216.4, 21.2, 59), (226.4, 20.5, 75), (241.0, 10.0, 59), (262.7, 6.0, -62), (289.7, 15.0, 49)]
# through_centre.csv / through_centre_v2.json: the limb is traced as the boundary of the region
# brighter than a threshold, on frames where it passes within 40 px of the centre, with a
# dark side under luminance 20 and a parabola rms under 6 px. Two thresholds trace two edges:
#   60  = the faint OUTER edge of the atmospheric glow (a tangent height well above the surface,
#         so a smaller depression than the surface horizon): 21 frames, median 8.4 px
#   120 = the bright limb itself, close to the surface horizon: 12 frames, median 26 px
# Five threshold-60 frames (206-225 s) returned sagittas of -245 to -354 px, which no horizon
# gives (a parabola that deep over a 900-px chord is a circle of radius ~300 px): the detector
# had run along the edge of a sun-glare band. They are excluded from both sets.
GOFAST_SAGITTAS = sorted([-8.0, 0.6, 1.5, 3.2, 4.4, 6.0, 6.2, 7.1, 7.2, 7.4, 8.4, 8.8, 10.0, 11.8, 15.0, 16.5, 18.7, 20.2, 20.5, 21.2, 26.2])
GOFAST_GLARE_ARTEFACTS = 5
GOFAST_CLEAN_N = len(GOFAST_SAGITTAS)                                   # 21 (outer glow edge)
GOFAST_CLEAN_MEDIAN = GOFAST_SAGITTAS[len(GOFAST_SAGITTAS) // 2]        # 8.4
GOFAST_CLEAN_POSITIVE = sum(1 for v in GOFAST_SAGITTAS if v > 0)        # 20
GOFAST_CLEAN_ANGLES = (-62, 83)
GOFAST_LIMB = [(202.7, 13.8, 734), (206.0, 17.0, 735), (206.5, 16.8, 730), (207.0, 22.0, 731), (207.5, 28.8, 738), (209.5, 29.9, 765),
               (213.7, 35.4, 736), (216.4, 24.9, 824), (217.1, 35.2, 811), (226.4, 32.2, 724), (226.9, 8.5, 710), (248.2, 27.2, 1401)]   # (t, sagitta, chord)
GOFAST_LIMB_SAGITTAS = sorted(v for _, v, _ in GOFAST_LIMB)
GOFAST_LIMB_N = len(GOFAST_LIMB)                                          # 12 (bright limb)
GOFAST_LIMB_MEDIAN = (GOFAST_LIMB_SAGITTAS[5] + GOFAST_LIMB_SAGITTAS[6]) / 2   # 26.1
GOFAST_LIMB_POSITIVE = sum(1 for v in GOFAST_LIMB_SAGITTAS if v > 0)     # 12
GOFAST_APOGEE_KM = 117.6          # 385,800 ft, CSXT's reported apogee for the 2014 flight (the book says 116 km, the 2004 figure)
GOFAST_CHORD = 780.0              # typical chord through the centre, px (710-830 in the limb set)
GOFAST_SPIN_STOP_S = 227          # in the re-upload: roll of roughly a turn a second to nil within a second
GOPRO_HFOV = 118.0                # GoPro-class wide setting, degrees, treated as equidistant: f = 640 / 1.03 rad = 622 px/rad
GOPRO_F = 640 / math.radians(GOPRO_HFOV / 2)


# --- MAPHEUS-5 (the book's QR 72), VSB-30 nominal flight sequence and what it implies ---
# Nominal main events of a VSB-30 flight (Palmerio et al., "VSB-30 sounding rocket: history of
# flight performance", J. Aerosp. Technol. Manag., Table 1): time s, altitude km. DLR's flown
# figures for MAPHEUS-5 (30 June 2015): apogee 253 km, microgravity from 74 s for over six
# minutes, despin at about 70 km (DLR press release; MORABA MAPHEUS page).
VSB30_EVENTS = [("S31 ignition", 0.0, 0.051), ("S31 burnout / separation", 13.5, 3.617), ("S30 ignition", 15.0, 4.323),
                ("S30 burnout", 44.0, 43.1), ("nosecone ejection", 55.0, 64.3), ("yo-yo despin", 56.0, 66.1),
                ("payload separation", 59.0, 71.7), ("apogee", 259.0, 252.7)]
MAPHEUS5_APOGEE_KM, MAPHEUS5_MICROG_START_S = 253.0, 74.0
# US Standard Atmosphere 1976 density, kg/m^3, at geometric altitude km
US76_RHO = {0: 1.225, 253: 6.0e-11, 4: 0.8194, 10: 0.4135, 20: 0.08891, 30: 0.01841, 40: 0.003996, 43: 0.002620, 50: 0.001027,
            60: 3.097e-4, 66: 1.35e-4, 70: 8.283e-5, 100: 5.6e-7, 117: 3.8e-8}
G0 = 9.81


def density_pct(alt_km):
    """Air density as a percentage of sea level, US76, linear-in-log between tabulated points."""
    ks = sorted(US76_RHO)
    for a, b in zip(ks, ks[1:]):
        if a <= alt_km <= b:
            la, lb = math.log(US76_RHO[a]), math.log(US76_RHO[b])
            return 100 * math.exp(la + (lb - la) * (alt_km - a) / (b - a)) / US76_RHO[0]
    raise ValueError(alt_km)


def coast(h0_km, target_km, inverse_square=True, dt=0.01):
    """Vertical velocity at h0 that just reaches target under g(h) (inverse-square about R, or
    constant G0), and the time it takes: a ballistic coast with nothing else acting."""
    def g(h_m):
        return G0 * (R * 1000 / (R * 1000 + h_m)) ** 2 if inverse_square else G0
    def apogee(v0):
        h, v, t = h0_km * 1000, v0, 0.0
        while v > 0:
            v -= g(h) * dt; h += v * dt; t += dt
        return h / 1000, t
    lo, hi = 500.0, 4000.0
    for _ in range(40):
        mid = (lo + hi) / 2
        if apogee(mid)[0] < target_km:
            lo = mid
        else:
            hi = mid
    return mid, apogee(mid)[1]


def mapheus():
    ev = dict((e[0], (e[1], e[2])) for e in VSB30_EVENTS)
    t_bo, h_bo = ev["S30 burnout"]; t_ap, h_ap = ev["apogee"]
    v, t = coast(h_bo, h_ap); v_c, t_c = coast(h_bo, h_ap, inverse_square=False)
    print(f"MAPHEUS-5 / VSB-30: S30 burns {ev['S30 ignition'][0]:.0f}-{t_bo:.0f} s, {ev['S30 ignition'][1]:.1f}-{h_bo:.1f} km; "
          f"air density {density_pct(ev['S30 ignition'][1]):.0f} % of sea level at ignition, {density_pct(h_bo):.2f} % at burnout, "
          f"{density_pct(ev['yo-yo despin'][1]):.3f} % at despin ({ev['yo-yo despin'][1]:.0f} km), {density_pct(117):.0e} % at 117 km")
    print(f"  coast {h_bo:.1f} -> {h_ap:.1f} km needs {v:.0f} m/s vertical and takes {t:.0f} s with inverse-square g "
          f"(table: {t_ap - t_bo:.0f} s); with constant 9.81: {v_c:.0f} m/s, {t_c:.0f} s")
    v_go = 3580 * 0.44704
    print(f"  GoFast: 3,580 mph = {v_go:.0f} m/s; a ballistic rise of {v_go**2 / (2 * 9.7) / 1000:.0f} km from burnout, apogee reported {GOFAST_APOGEE_KM} km")


def gofast_expected(hfov_deg, alt_km=GOFAST_APOGEE_KM, chord=GOFAST_CHORD, width_px=1280):
    """Expected through-centre sagitta over a chord, rectilinear lens of the given field."""
    f = (width_px / 2) / math.tan(math.radians(hfov_deg / 2))
    # exact for a rectilinear lens (pitch-independent, so the level-camera formula applies): f tan(d) (1/cos(phi) - 1)
    return bow_px(dip_deg(alt_km), 2 * math.degrees(math.atan((chord / 2) / f)), chord / 2)


def gofast_expected_fisheye(alt_km=GOFAST_APOGEE_KM, chord=GOFAST_CHORD, f=GOPRO_F):
    """Expected through-centre sagitta for the GoPro-class fisheye (exact equidistant form)."""
    return through_centre_sagitta(dip_deg(alt_km), chord, f, exact_equidistant=True)


def tables():
    print(f"dip at 36 km {dip_deg(36):.2f} deg; at {GOFAST_APOGEE_KM} km {dip_deg(GOFAST_APOGEE_KM):.2f} deg; at 80 km {dip_deg(80):.2f} deg")
    print("bow on a 1280-px frame, level rectilinear lens, 36 km:", "  ".join(f"{h} deg: {bow_px(dip_deg(36), h, 640):.0f} px" for h in (40, 60, 90, 120)))
    print("  pitch check at 90 deg: " + ", ".join(f"{p:+d} deg -> {bow_px_pitched(dip_deg(36), 90, p):.1f}" for p in (-20, 0, 20)))
    print(f"render: {RENDER['dip_px']} px for {dip_deg(36):.2f} deg -> {RENDER_PX_PER_DEG:.1f} px/deg, f {RENDER_F:.0f} px, HFOV {RENDER_HFOV:.0f} deg; predicted bow {bow_px(dip_deg(36), RENDER_HFOV, RENDER['width_px']/2):.0f} px vs read {RENDER['bow_px']}")
    glow = math.degrees(math.acos((R + GLOW_TOP_TANGENT_KM) / (R + 36)))
    print(f"Focozz glow-top depression ~{glow:.1f} deg; expected glow-top edge bow, level camera, at 60 / {RENDER_HFOV:.0f} deg rectilinear: {bow_px(glow, 60, 640):.0f} / {bow_px(glow, RENDER_HFOV, 640):.0f} px, fisheye {bow_px_fisheye_level(glow, 640, GOPRO_F):.0f}; measured {FOCOZZ_BOW}")
    print(f"GoFast expected through-centre sagitta over a {GOFAST_CHORD:.0f}-px chord at {GOFAST_APOGEE_KM} km, rectilinear:", "  ".join(f"{h} deg: {gofast_expected(h):.0f} px" for h in (60, 90, 120)), "; at 60 km, 90 deg:", f"{gofast_expected(90, 60):.0f} px")
    print(f"  GoPro-class fisheye ({GOPRO_HFOV:.0f} deg, f {GOPRO_F:.0f} px/rad): {gofast_expected_fisheye():.0f} px at {GOFAST_APOGEE_KM} km, {gofast_expected_fisheye(80):.0f} at 80 km, {gofast_expected_fisheye(60):.0f} at 60 km; over 710-830 px chords at apogee {gofast_expected_fisheye(chord=710):.0f}-{gofast_expected_fisheye(chord=830):.0f}; flat plane: 0")
    print(f"  level-camera fisheye edge bow (not through-centre), full 1210 width at {GOFAST_APOGEE_KM} km: {bow_px_fisheye_level(dip_deg(GOFAST_APOGEE_KM), 605, GOPRO_F):.0f} px; through-centre over the same 1210 chord: {gofast_expected_fisheye(chord=1210):.0f} px")
    mapheus()
    print("GoFast measured, outer glow edge (thr 60):", f"n={GOFAST_CLEAN_N} (+{GOFAST_GLARE_ARTEFACTS} glare artefacts), median {GOFAST_CLEAN_MEDIAN} px, positive {GOFAST_CLEAN_POSITIVE}/{GOFAST_CLEAN_N}, orientations {GOFAST_CLEAN_ANGLES}")
    print("GoFast measured, bright limb (thr 120):", f"n={GOFAST_LIMB_N}, median {GOFAST_LIMB_MEDIAN:.0f} px, positive {GOFAST_LIMB_POSITIVE}/{GOFAST_LIMB_N}, range {GOFAST_LIMB_SAGITTAS[0]}-{GOFAST_LIMB_SAGITTAS[-1]}; spin stops at {GOFAST_SPIN_STOP_S} s")


def measure_through_centre(frames_dir, w=1280, h=720, thr=60, near=40):
    """Re-run the GoFast measurement: limb = boundary of the bright (Earth) region; frames
    where it passes within `near` px of the frame centre; sagitta of a parabola fitted along
    the longest chord, signed toward the dark side. Needs OpenCV and numpy."""
    import glob
    import re
    import cv2
    import numpy as np
    cx, cy = w / 2, h / 2; rows = []
    for f in sorted(glob.glob(os.path.join(frames_dir, "a*.jpg"))):
        t = 172 + (int(re.search(r"a(\d+)\.jpg", f).group(1)) - 1) / 10
        im = cv2.imread(f); g = cv2.GaussianBlur(cv2.cvtColor(im, cv2.COLOR_BGR2GRAY), (5, 5), 0)
        m = (g > thr).astype(np.uint8); frac = m.mean()
        if frac < 0.2 or frac > 0.8:
            continue
        cnts, _ = cv2.findContours(m, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
        c = max(cnts, key=cv2.contourArea)[:, 0, :].astype(float)
        c = c[(c[:, 0] > 6) & (c[:, 0] < w - 7) & (c[:, 1] > 6) & (c[:, 1] < h - 7)]
        if len(c) < 300 or np.hypot(c[:, 0] - cx, c[:, 1] - cy).min() > near:
            continue
        P = c[[np.argmin(c[:, 0]), np.argmax(c[:, 0]), np.argmin(c[:, 1]), np.argmax(c[:, 1])]]; best = None
        for a in range(4):
            for b in range(a + 1, 4):
                L = np.hypot(*(P[b] - P[a]))
                if best is None or L > best[0]:
                    best = (L, P[a], P[b])
        L, p0, p1 = best; v = (p1 - p0) / L; n = np.array([-v[1], v[0]])
        q1 = (cx + 80 * n[0], cy + 80 * n[1]); q2 = (cx - 80 * n[0], cy - 80 * n[1])
        b1 = g[int(np.clip(q1[1], 0, h - 1)), int(np.clip(q1[0], 0, w - 1))]; b2 = g[int(np.clip(q2[1], 0, h - 1)), int(np.clip(q2[0], 0, w - 1))]
        if b1 > b2:
            n = -n; b1, b2 = b2, b1
        s = (c - p0) @ v; dev = (c - p0) @ n; p = np.polyfit(s, dev, 2)
        rows.append((round(t, 1), round(L), round(float(-p[0] * (L / 2) ** 2), 1), int(b1), round(math.degrees(math.atan2(v[1], v[0])))))
    return rows


def figure(out=OUT):
    d36 = dip_deg(36); glow = math.degrees(math.acos((R + GLOW_TOP_TANGENT_KM) / (R + 36)))
    W, H = 900, 420; L, Rr, T, B = 70, 860, 40, 350
    X = lambda h: L + (h - 20) / (140 - 20) * (Rr - L)
    Y = lambda b: B - (b + 20) / (90) * (B - T)
    o = [f'<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 {W} {H}" width="{W}" height="{H}" font-family="Georgia, serif">', f'<rect width="{W}" height="{H}" fill="#fff"/>',
         f'<text x="{W/2}" y="24" font-size="14.5" text-anchor="middle" fill="#2E4057" font-weight="700">How much the horizon should bow at 36 km, by lens (1,280-px-wide frame)</text>']
    for b in range(-20, 71, 10):
        o.append(f'<line x1="{L}" x2="{Rr}" y1="{Y(b):.1f}" y2="{Y(b):.1f}" stroke="#eee"/><text x="{L-8}" y="{Y(b)+4:.1f}" font-size="10" text-anchor="end" fill="#6E675B">{b}</text>')
    for hh in range(20, 141, 20):
        o.append(f'<text x="{X(hh):.1f}" y="{B+16}" font-size="10" text-anchor="middle" fill="#6E675B">{hh}&#176;</text>')
    o.append(f'<line x1="{L}" x2="{Rr}" y1="{B}" y2="{B}" stroke="#6E675B"/><line x1="{L}" x2="{L}" y1="{T}" y2="{B}" stroke="#6E675B"/>')
    o.append(f'<text x="{(L+Rr)/2}" y="{B+34}" font-size="11.5" text-anchor="middle" fill="#3A352D">horizontal field of view of the lens</text>')
    o.append(f'<text x="18" y="{(T+B)/2}" font-size="11.5" text-anchor="middle" fill="#3A352D" transform="rotate(-90 18 {(T+B)/2})">edge below centre, px</text>')
    # Focozz band
    o.append(f'<rect x="{L}" y="{Y(4.2):.1f}" width="{Rr-L}" height="{Y(-9)-Y(3):.1f}" fill="#f2c14e" opacity=".35"/>')
    o.append(f'<text x="{Rr-6}" y="{Y(-11):.1f}" font-size="10" text-anchor="end" fill="#8a5710">balloon compilation frames, limb-glow top: &#8722;9 to +4 px (lens, height unknown)</text>')
    for dep, col, dash, lab in ((d36, "#2c5a94", "", f"surface horizon, dip {d36:.2f}&#176;"), (glow, "#2c5a94", ' stroke-dasharray="5 4"', f"top of the limb glow, ~{glow:.1f}&#176; down")):
        pts = " ".join(f"{X(hh):.1f},{Y(bow_px(dep, hh, 640)):.1f}" for hh in range(20, 141, 2))
        o.append(f'<polyline points="{pts}" fill="none" stroke="{col}" stroke-width="{2 if not dash else 1.3}"{dash}/>')
        o.append(f'<text x="{X(140)-4:.1f}" y="{Y(bow_px(dep,140,640))-6:.1f}" font-size="10" text-anchor="end" fill="{col}">{lab}</text>')
    o.append(f'<line x1="{L}" x2="{Rr}" y1="{Y(0):.1f}" y2="{Y(0):.1f}" stroke="#992e26" stroke-width="2"/><text x="{L+6}" y="{Y(0)-6:.1f}" font-size="10.5" fill="#992e26" font-weight="700">flat plane: zero bow, any lens, any tilt</text>')
    rb = RENDER["bow_px"] * 640 / (RENDER["width_px"] / 2)
    o.append(f'<circle cx="{X(RENDER_HFOV):.1f}" cy="{Y(rb):.1f}" r="6" fill="#276b2b"/><text x="{X(RENDER_HFOV)-10:.1f}" y="{Y(rb)-10:.1f}" font-size="10.5" text-anchor="end" fill="#276b2b" font-weight="700">the p. 235 render: ~{RENDER_HFOV:.0f}&#176;, {RENDER["bow_px"]} px on 1,210</text>')
    o.append("</svg>")
    os.makedirs(os.path.dirname(out), exist_ok=True); open(out, "w").write("\n".join(o)); print("wrote", os.path.normpath(out))


if __name__ == "__main__":
    if "--figure" in sys.argv:
        figure()
    elif "--measure" in sys.argv:
        for r in measure_through_centre(sys.argv[sys.argv.index("--measure") + 1]):
            print(r)
    else:
        tables()