app.py

# app.py
# MH370 Flight Simulator — Leaflet map (Folium) + Ping Rings (KML)
# Guided scenarios (ACTIVE):
#   A) Straight heading
#   B) Variable heading & speed  (beam search, BTO-first)
#   C) Variable heading & speed (SE bias after 22:40)  (beam search, BTO-first + tiny east bias only after 22:40)
#   D) Variable heading, CONSTANT southern speed (443 kt)  (beam search, BTO-first)
# Key behavior:
# - Last radar spot FIXED: 18:22:00Z @ 6.578N, 96.341E
# - Turn completion FIXED: 18:40:00Z
# - Guided mode computes a realistic turn between 18:22→18:40 and then southern path.
# - Guided validation timestamps: 18:22, 18:40, 19:40, 20:40, 21:40, 22:40, 00:10, 00:20

import os
import re
import html
import math
import zipfile
import xml.etree.ElementTree as ET
from datetime import datetime, timezone, timedelta

import gradio as gr
import numpy as np
import folium
from folium.plugins import MousePosition
import plotly.graph_objects as go

import flight_plotter_v32 as fp


# ================= CONFIG =================
KML_FILE = "MH370 - 13 nov (3).kml"
SEARCH_AREAS_KMZ = "Search areas.kmz"
MAX_RANGE_KML = "MH370_max_range_6166km_KUL_red.kml"

DEFAULT_LEGS = [
    [300, 500, 10],
    [270, 500, 10],
    [240, 500, 10],
]
NEW_LEG_DEFAULT = [0, 500, 10]

WAYPOINTS = [
    ("IGARI", 6.93611111, 103.585),
    ("VAMPI", 6.18222222, 97.58555556),
    ("NILAM", 6.75638889, 95.97638889),
    ("SANOB", 6.58611111, 95.66916667),
    ("VPG",   5.27963889, 100.26038889),
    ("MEKAR", 6.5, 96.48333333),
]

CAPT_SIMON_HARDY_NAME = "Captain Simon Hardy"
CAPT_SIMON_HARDY_LAT = -(39 + 10 / 60)   # 39°10' S
CAPT_SIMON_HARDY_LON = (88 + 18 / 60)    # 88°18' E

MEASURED_BFO_0020_HZ = 182.0

RING_TIMES = {
    "Ping Ring 1": "18:25:00Z",
    "Ping Ring 2": "18:40:00Z",
    "Ping Ring 3": "19:40:00Z",
    "Ping Ring 4": "20:40:00Z",
    "Ping Ring 5": "21:40:00Z",
    "Ping Ring 6": "22:40:00Z",
    "Ping Ring 7": "24:10:00Z",
}

SATLOC_IDX_MEAS_BFO = 4
SATLOC_IDX_MEAS_BTO = 5

RADAR_FIX_LAT = 6.578
RADAR_FIX_LON = 96.341
RADAR_FIX_TIME_UTC = datetime(2014, 3, 7, 18, 22, 0, tzinfo=timezone.utc)
TURN_COMPLETE_TIME_UTC = datetime(2014, 3, 7, 18, 40, 0, tzinfo=timezone.utc)
RADAR_INITIAL_TRACK_DEG = 296.0

GUIDED_VALIDATION_LABELS = [
    "18:22:00",
    "18:40:00",
    "19:40:00",
    "20:40:00",
    "21:40:00",
    "22:40:00",
    "24:10:00",
    "24:20:00",
]
GUIDED_FINAL_LABEL = "24:20:00"

TURN_BANK_MIN_DEG = 15.0
TURN_BANK_MAX_DEG = 30.0
TURN_BANK_DEFAULT_DEG = 25.0
TURN_STEP_SEC = 30.0

TURN_SPEED_MIN_KT = 470.0
TURN_SPEED_MAX_KT = 520.0
SOUTH_SPEED_MIN_KT = 430.0
SOUTH_SPEED_MAX_KT = 510.0

CONST_SOUTH_SPEED_KT = 443.0

GUIDED_HEADING_MIN = 175
GUIDED_HEADING_MAX = 190
COARSE_HEADING_STEP_DEG = 2.0
REFINE_HEADING_HALF_WINDOW_DEG = 2.0
REFINE_HEADING_STEP_DEG = 0.5

GUIDED_SOUTH_SEG_MIN = [60.0, 60.0, 60.0, 60.0, 90.0, 10.0]
GUIDED_SOUTH_SEG_LABELS_END = ["19:40:00", "20:40:00", "21:40:00", "22:40:00", "24:10:00", "24:20:00"]

HARD_BTO_LABELS = ["18:40:00", "19:40:00", "20:40:00", "21:40:00", "22:40:00", "24:10:00", "24:20:00"]

TOPK_POLISH = 5
ENABLE_POLISH_STRAIGHT = True
POLISH_PASSES = 3
POLISH_WINDOW_KT = 6.0
POLISH_STEP_KT = 2.0
POLISH_OBJ_EPS = 1e-3

VAR_HEADING_MIN = 170.0
VAR_HEADING_MAX = 200.0
VAR_TURN_SPEED_STEP = 10.0
VAR_FIRST_HEADING_STEP = 4.0
VAR_HEADING_DELTA_SET = [-6, -3, 0, 3, 6]
VAR_BEAM_WIDTH = 12
VAR_SMOOTH_PENALTY_US_PER_DEG = 1.5

SE_EAST_DIVERGE_START_SEG_INDEX = 4
SE_VAR_EAST_DELTA_SET = [-8, -4, -2, 0, 2]
SE_EAST_BIAS_PENALTY_US_PER_DEG = 0.8

SEARCH_AREA_COLORS = {
    "Ocean Infinity": "#ff3b3b",
    "Australia": "#1f77b4",
    "Malays": "#2ca02c",
    "Malaysia": "#2ca02c",
    "China": "#ff7f0e",
    "Seabed Constructor": "#9467bd",
    "Fugro": "#8c564b",
}
DEFAULT_SEARCH_COLOR = "#00ffff"


# ================= PATH RESOLUTION =================
def resolve_local_or_mnt(path_like: str) -> str:
    p = (path_like or "").strip()
    if not p:
        return p
    if os.path.exists(p):
        return p
    here = os.path.dirname(os.path.abspath(__file__)) if "__file__" in globals() else os.getcwd()
    p2 = os.path.join(here, p)
    if os.path.exists(p2):
        return p2
    base = os.path.basename(p)
    p3 = os.path.join("/mnt/data", base)
    if os.path.exists(p3):
        return p3
    p4 = os.path.join("/mnt/data", p)
    if os.path.exists(p4):
        return p4
    return p


def color_for_search_area(name: str) -> str:
    lname = (name or "").lower()
    for key, color in SEARCH_AREA_COLORS.items():
        if key.lower() in lname:
            return color
    return DEFAULT_SEARCH_COLOR


# ================= KML (Ping Rings) =================
def load_ping_rings_from_kml(path: str):
    path = resolve_local_or_mnt(path)
    if not os.path.exists(path):
        raise FileNotFoundError(f"KML file not found: '{path}'.")

    tree = ET.parse(path)
    root = tree.getroot()

    m = re.match(r"\{(.*)\}", root.tag)
    ns_uri = m.group(1) if m else "http://www.opengis.net/kml/2.2"
    ns = {"kml": ns_uri}

    rings = {}
    for pm in root.findall(".//kml:Placemark", ns):
        name_el = pm.find("kml:name", ns)
        if name_el is None or name_el.text is None:
            continue
        name = name_el.text.strip()
        if not name or not name.lower().startswith("ping ring"):
            continue

        coords_el = pm.find(".//kml:LineString/kml:coordinates", ns)
        if coords_el is None or coords_el.text is None:
            continue

        pts = []
        for token in re.split(r"\s+", coords_el.text.strip()):
            if not token:
                continue
            parts = token.split(",")
            if len(parts) < 2:
                continue
            lon = float(parts[0])
            lat = float(parts[1])
            pts.append((lat, lon))

        if pts:
            rings[name] = pts

    return rings


try:
    PING_RINGS = load_ping_rings_from_kml(KML_FILE)
except Exception as e:
    PING_RINGS = {}
    KML_LOAD_ERROR = str(e)
else:
    KML_LOAD_ERROR = ""


# ================= KML (Generic features / max range ring) =================
def load_generic_kml_features(path: str):
    path = resolve_local_or_mnt(path)
    if not os.path.exists(path):
        raise FileNotFoundError(f"KML file not found: '{path}'.")

    tree = ET.parse(path)
    root = tree.getroot()

    m = re.match(r"\{(.*)\}", root.tag)
    ns_uri = m.group(1) if m else "http://www.opengis.net/kml/2.2"
    ns = {"kml": ns_uri}

    features = []
    for pm in root.findall(".//kml:Placemark", ns):
        name_el = pm.find("kml:name", ns)
        pm_name = name_el.text.strip() if (name_el is not None and name_el.text) else "KML feature"

        for pt in pm.findall(".//kml:Point", ns):
            coords_el = pt.find(".//kml:coordinates", ns)
            if coords_el is None or coords_el.text is None:
                continue
            coords = coords_el.text.strip().split(",")
            if len(coords) >= 2:
                lon = float(coords[0])
                lat = float(coords[1])
                features.append({
                    "type": "point",
                    "name": pm_name,
                    "lat": lat,
                    "lon": lon,
                })

        for ls in pm.findall(".//kml:LineString", ns):
            coords_el = ls.find(".//kml:coordinates", ns)
            if coords_el is None or coords_el.text is None:
                continue
            pts = []
            for token in re.split(r"\s+", coords_el.text.strip()):
                if not token:
                    continue
                parts = token.split(",")
                if len(parts) < 2:
                    continue
                lon = float(parts[0])
                lat = float(parts[1])
                pts.append((lat, lon))
            if pts:
                features.append({
                    "type": "line",
                    "name": pm_name,
                    "points": pts,
                })

        for poly in pm.findall(".//kml:Polygon", ns):
            outer = poly.find(".//kml:outerBoundaryIs/kml:LinearRing/kml:coordinates", ns)
            if outer is None or outer.text is None:
                continue
            pts = []
            for token in re.split(r"\s+", outer.text.strip()):
                if not token:
                    continue
                parts = token.split(",")
                if len(parts) < 2:
                    continue
                lon = float(parts[0])
                lat = float(parts[1])
                pts.append((lat, lon))
            if pts:
                features.append({
                    "type": "polygon",
                    "name": pm_name,
                    "points": pts,
                })

    return features


try:
    MAX_RANGE_FEATURES = load_generic_kml_features(MAX_RANGE_KML)
except Exception as e:
    MAX_RANGE_FEATURES = []
    MAX_RANGE_KML_ERROR = str(e)
else:
    MAX_RANGE_KML_ERROR = ""


# ================= KMZ (Search Areas) =================
def _kml_ns_from_root(root):
    m = re.match(r"\{(.*)\}", root.tag)
    ns_uri = m.group(1) if m else "http://www.opengis.net/kml/2.2"
    return {"kml": ns_uri}


def _parse_kml_coordinates_text(txt: str):
    pts = []
    if not txt:
        return pts
    for token in re.split(r"\s+", txt.strip()):
        if not token:
            continue
        parts = token.split(",")
        if len(parts) < 2:
            continue
        lon = float(parts[0])
        lat = float(parts[1])
        pts.append((lat, lon))
    return pts


def load_search_areas_from_kmz(kmz_path: str):
    kmz_path = resolve_local_or_mnt(kmz_path)
    if not os.path.exists(kmz_path):
        raise FileNotFoundError(f"KMZ file not found: '{kmz_path}'.")

    with zipfile.ZipFile(kmz_path, "r") as zf:
        kml_name = None
        for n in zf.namelist():
            if n.lower().endswith(".kml"):
                kml_name = n
                break
        if not kml_name:
            raise ValueError("KMZ contains no .kml file.")
        root = ET.fromstring(zf.read(kml_name))

    ns = _kml_ns_from_root(root)
    out = []

    for pm in root.findall(".//kml:Placemark", ns):
        name_el = pm.find("kml:name", ns)
        pm_name = (name_el.text.strip() if (name_el is not None and name_el.text) else "Search area")

        points, lines, polys = [], [], []

        for pt in pm.findall(".//kml:Point", ns):
            coords = pt.find(".//kml:coordinates", ns)
            if coords is None or coords.text is None:
                continue
            pts = _parse_kml_coordinates_text(coords.text)
            if pts:
                points.extend(pts[:1])

        for ls in pm.findall(".//kml:LineString", ns):
            coords = ls.find(".//kml:coordinates", ns)
            if coords is None or coords.text is None:
                continue
            pts = _parse_kml_coordinates_text(coords.text)
            if len(pts) >= 2:
                lines.append(pts)

        for poly in pm.findall(".//kml:Polygon", ns):
            outer = poly.find(".//kml:outerBoundaryIs/kml:LinearRing/kml:coordinates", ns)
            if outer is None or outer.text is None:
                continue
            pts = _parse_kml_coordinates_text(outer.text)
            if len(pts) >= 3:
                polys.append(pts)

        if points or lines or polys:
            out.append({"name": pm_name, "points": points, "lines": lines, "polygons": polys})

    if not out:
        raise ValueError("No drawable features found in KMZ/KML.")
    return out


try:
    SEARCH_AREAS = load_search_areas_from_kmz(SEARCH_AREAS_KMZ)
except Exception as e:
    SEARCH_AREAS = []
    KMZ_LOAD_ERROR = str(e)
else:
    KMZ_LOAD_ERROR = ""


# ================= LEGS HELPERS =================
def df_to_rows(df):
    if df is None:
        return []
    if hasattr(df, "to_numpy"):
        return df.to_numpy().tolist()
    return df


def parse_legs_df(df):
    rows = df_to_rows(df)
    legs = []
    for row in rows:
        if row is None or len(row) < 3:
            continue
        h, s, m = row[0], row[1], row[2]
        if h is None or s is None or m is None:
            continue
        if isinstance(h, float) and np.isnan(h):
            continue
        if isinstance(s, float) and np.isnan(s):
            continue
        if isinstance(m, float) and np.isnan(m):
            continue
        legs.append((float(h), float(s), float(m)))
    return legs or [(300, 500, 10), (270, 500, 10), (240, 500, 10)]


def legs_to_rows(legs):
    return [[float(h), float(s), float(m)] for (h, s, m) in legs]


def add_leg(df):
    rows = df_to_rows(df)
    cleaned = []
    for r in rows:
        if r is None or len(r) < 3:
            continue
        a, b, c = r[0], r[1], r[2]
        if (a is None or (isinstance(a, float) and np.isnan(a))) and \
           (b is None or (isinstance(b, float) and np.isnan(b))) and \
           (c is None or (isinstance(c, float) and np.isnan(c))):
            continue
        cleaned.append([a, b, c])
    cleaned.append(list(NEW_LEG_DEFAULT))
    return gr.update(value=cleaned)


def reset_legs():
    return gr.update(value=[row[:] for row in DEFAULT_LEGS])


# ================= GEO HELPERS =================
def _wrap_heading_deg(h: float) -> float:
    return float(h) % 360.0


def _shortest_signed_delta_deg(h_from: float, h_to: float) -> float:
    return ((float(h_to) - float(h_from) + 540.0) % 360.0) - 180.0


def _turn_rate_deg_per_sec(speed_knots: float, bank_deg: float) -> float:
    g = 9.80665
    V = max(10.0, float(speed_knots) * 0.514444)
    bank = math.radians(max(1e-3, float(bank_deg)))
    omega = g * math.tan(bank) / V
    return float(omega * 180.0 / math.pi)


def build_turn_legs(speed_knots: float, method: str, heading_final_deg: float,
                    bank_pref_deg: float = TURN_BANK_DEFAULT_DEG, step_sec: float = TURN_STEP_SEC):
    total_sec = (TURN_COMPLETE_TIME_UTC - RADAR_FIX_TIME_UTC).total_seconds()
    if total_sec <= 0:
        return [], {"error": "Invalid turn window."}

    h0 = float(RADAR_INITIAL_TRACK_DEG)
    h1 = _wrap_heading_deg(float(heading_final_deg))
    delta = _shortest_signed_delta_deg(h0, h1)
    sign = 1.0 if delta >= 0 else -1.0
    delta_abs = abs(delta)

    bank_candidates = [float(bank_pref_deg)] + list(np.linspace(TURN_BANK_MIN_DEG, TURN_BANK_MAX_DEG, 9))
    best_bank = None
    best_turn_sec = None

    for b in bank_candidates:
        rate = _turn_rate_deg_per_sec(speed_knots, b)
        if rate <= 1e-6:
            continue
        t_turn = delta_abs / rate
        if t_turn <= total_sec + 1e-6:
            if best_bank is None:
                best_bank, best_turn_sec = b, t_turn
            else:
                if abs(b - bank_pref_deg) < abs(best_bank - bank_pref_deg):
                    best_bank, best_turn_sec = b, t_turn

    if best_bank is None:
        best_bank = TURN_BANK_MAX_DEG
        best_turn_sec = delta_abs / max(1e-6, _turn_rate_deg_per_sec(speed_knots, best_bank))
        best_turn_sec = min(best_turn_sec, total_sec)

    turn_sec = float(best_turn_sec)
    straight_sec = max(0.0, float(total_sec) - turn_sec)

    legs = []

    if straight_sec > 1.0:
        legs.append((h0, float(speed_knots), straight_sec / 60.0))

    if turn_sec > 1.0 and delta_abs > 1e-6:
        step_sec = max(5.0, float(step_sec))
        n = int(math.ceil(turn_sec / step_sec))
        if n < 1:
            n = 1
        remaining = turn_sec
        for i in range(n):
            this = min(step_sec, remaining)
            remaining -= this
            frac_mid = ((i + 0.5) / n)
            hdg_mid = _wrap_heading_deg(h0 + sign * delta_abs * frac_mid)
            legs.append((hdg_mid, float(speed_knots), this / 60.0))

    info = {
        "h0": h0,
        "h1": h1,
        "delta_deg": delta,
        "bank_deg": float(best_bank),
        "straight_min": straight_sec / 60.0,
        "turn_min": turn_sec / 60.0,
        "total_min": total_sec / 60.0,
        "speed_turn_kt": float(speed_knots),
    }
    return legs, info


# ================= TIME HELPERS =================
def _parse_start_time(s: str) -> datetime:
    s = (s or "").strip()
    if not s:
        return datetime(2014, 3, 7, 18, 22, 0, tzinfo=timezone.utc)

    if re.fullmatch(r"\d{2}:\d{2}:\d{2}Z?", s):
        s2 = s.replace("Z", "")
        hh, mm, ss = map(int, s2.split(":"))
        return datetime(2014, 3, 7, hh, mm, ss, tzinfo=timezone.utc)

    s = s.replace("T", " ").replace("Z", "")
    dt = datetime.fromisoformat(s)
    if dt.tzinfo is None:
        dt = dt.replace(tzinfo=timezone.utc)
    return dt.astimezone(timezone.utc)


def _to_label_utc(dt: datetime) -> str:
    return dt.strftime("%H:%M:%SZ")


def _dt_from_label(label_hms: str) -> datetime:
    s = (label_hms or "").strip().replace("Z", "")
    return fp._time_label_to_dt(s, fp.SAT_EPOCH_UTC)


def dt_0020_utc() -> datetime:
    return fp._time_label_to_dt("24:20:00", fp.SAT_EPOCH_UTC)


# ================= PATCH satLoc =================
def _patch_measured_bfo_0020():
    target_labels = {"24:20:00", "00:20:00"}
    patched = False
    try:
        for row in fp.satLoc:
            try:
                lbl = fp._extract_time_string(row[0])
            except Exception:
                continue
            if lbl in target_labels:
                try:
                    row[SATLOC_IDX_MEAS_BFO] = float(MEASURED_BFO_0020_HZ)
                    patched = True
                except Exception:
                    pass
    except Exception:
        return False
    return patched


_ = _patch_measured_bfo_0020()


# ================= ALTITUDE HELPERS =================
def altitude_ft_at_dt(dt_ping: datetime, start_dt: datetime, alt0_ft: float,
                      vs_fpm: float, descent_start_dt):
    alt0_ft = float(alt0_ft)
    vs_fpm = float(vs_fpm)

    if descent_start_dt is None or vs_fpm == 0.0:
        return max(0.0, alt0_ft)
    if dt_ping < descent_start_dt:
        return max(0.0, alt0_ft)

    mins = (dt_ping - descent_start_dt).total_seconds() / 60.0
    alt = alt0_ft + vs_fpm * mins
    return max(0.0, alt)


def build_altitude_plot(start_dt: datetime, total_time_min: float, alt0_ft: float,
                        vs_fpm: float, descent_start_dt):
    t = np.linspace(0.0, max(1.0, float(total_time_min)), 120)
    alts = []
    for mins in t:
        dt = start_dt + timedelta(minutes=float(mins))
        alts.append(altitude_ft_at_dt(dt, start_dt, float(alt0_ft), float(vs_fpm), descent_start_dt))

    fig = go.Figure()
    fig.add_trace(go.Scatter(x=t, y=alts, mode="lines", name="Altitude (ft)"))
    fig.update_layout(
        title="Altitude vs Time",
        xaxis_title="Minutes since start",
        yaxis_title="Altitude (ft)",
        height=320,
        margin=dict(l=60, r=30, t=60, b=50),
    )
    return fig


# ================= SPEED PROFILE PLOT =================
def build_speed_profile_plot(start_dt: datetime, legs):
    xs = [0.0]
    ys = []
    cum_min = 0.0
    for (_hdg, spd, mins) in legs:
        mins = float(mins)
        spd = float(spd)
        if mins <= 0:
            continue
        if not ys:
            ys.append(spd)
        else:
            ys.append(ys[-1])
        xs.append(cum_min)
        ys.append(spd)
        cum_min += mins
        xs.append(cum_min)
        ys.append(spd)

    fig = go.Figure()
    if ys:
        fig.add_trace(go.Scatter(x=xs[1:], y=ys, mode="lines", name="Speed (kt)"))
    fig.update_layout(
        title="Speed Profile",
        xaxis_title="Minutes since start",
        yaxis_title="Speed (kt)",
        height=320,
        margin=dict(l=60, r=30, t=60, b=50),
    )
    return fig


# ================= SIM CORE =================
def simulate_path(start_lat, start_lon, legs, method):
    lat, lon = float(start_lat), float(start_lon)
    pts = [(lat, lon)]
    seg_dists_m = []
    total_time_min = 0.0

    for heading, speed_knots, minutes in legs:
        dist_m = float(speed_knots) * 1852.0 * (float(minutes) / 60.0)
        if method == "VINCENTY":
            (lat, lon), _ = fp.vincenty_direct((lat, lon), float(heading), dist_m)
        else:
            lat, lon = fp.rhumb_direct((lat, lon), float(heading), dist_m)

        pts.append((lat, lon))
        seg_dists_m.append(dist_m)
        total_time_min += float(minutes)

    return pts, seg_dists_m, total_time_min


def compute_stats(path_pts, seg_dists_m, total_time_min, legs):
    total_dist_m = float(np.sum(seg_dists_m)) if seg_dists_m else 0.0
    total_dist_km = total_dist_m / 1000.0
    total_dist_nm = total_dist_m / 1852.0
    total_time_hr = total_time_min / 60.0
    avg_gs_knots = (total_dist_nm / total_time_hr) if total_time_hr > 0 else 0.0
    end_lat, end_lon = path_pts[-1]
    return "\n\n".join([
        f"**Legs:** {len(legs)}",
        f"**Total time:** {total_time_min:.1f} min ({total_time_hr:.2f} hr)",
        f"**Total distance:** {total_dist_nm:.1f} nm ({total_dist_km:.1f} km)",
        f"**Average ground speed:** {avg_gs_knots:.1f} kt",
        f"**Final position (polyline end):** lat {end_lat:.4f}, lon {end_lon:.4f}",
    ])


# ================= EXACT STATE HELPERS =================
def _predict_state_from_leg_start(method, leg_start_lat, leg_start_lon, hdg_deg, spd_knots, t_inside_leg_sec):
    dist_m = float(spd_knots) * 1852.0 * (float(t_inside_leg_sec) / 3600.0)
    if method == "VINCENTY":
        (lat, lon), _ = fp.vincenty_direct((float(leg_start_lat), float(leg_start_lon)), float(hdg_deg), dist_m)
    else:
        lat, lon = fp.rhumb_direct((float(leg_start_lat), float(leg_start_lon)), float(hdg_deg), dist_m)
    return float(lat), float(lon)


def _precompute_leg_starts(start_lat, start_lon, legs, method):
    cur_lat = float(start_lat)
    cur_lon = float(start_lon)
    leg_start_positions = [(cur_lat, cur_lon)]
    for (hdg, spd, minutes) in legs[:-1]:
        dist_m = float(spd) * 1852.0 * (float(minutes) / 60.0)
        if method == "VINCENTY":
            (cur_lat, cur_lon), _ = fp.vincenty_direct((cur_lat, cur_lon), float(hdg), dist_m)
        else:
            cur_lat, cur_lon = fp.rhumb_direct((cur_lat, cur_lon), float(hdg), dist_m)
        leg_start_positions.append((cur_lat, cur_lon))
    return leg_start_positions


def _state_at_dt(start_lat, start_lon, legs, method, start_dt, dt_k):
    if dt_k < start_dt:
        return None

    leg_secs = [float(m) * 60.0 for _, _, m in legs]
    boundaries = np.cumsum([0.0] + leg_secs)
    total_end = float(boundaries[-1])

    t_sec = float((dt_k - start_dt).total_seconds())
    if t_sec > total_end:
        t_sec = total_end

    i = int(np.searchsorted(boundaries, t_sec, side="right") - 1)
    i = max(0, min(i, len(legs) - 1))
    t_inside = t_sec - float(boundaries[i])

    leg_start_positions = _precompute_leg_starts(start_lat, start_lon, legs, method)
    leg_start_lat, leg_start_lon = leg_start_positions[i]
    hdg, spd, _m = legs[i]
    lat, lon = _predict_state_from_leg_start(method, leg_start_lat, leg_start_lon, hdg, spd, t_inside)

    dist_km = 0.0
    for j in range(i):
        _hj, spdj, mj = legs[j]
        dist_km += float(spdj) * 1.852 * (float(mj) / 60.0)
    dist_km += float(spd) * 1.852 * (float(t_inside) / 3600.0)

    return {
        "dt": dt_k,
        "elapsed_min": t_sec / 60.0,
        "lat": float(lat),
        "lon": float(lon),
        "speed_kt": float(spd),
        "heading_deg": float(hdg),
        "dist_km": float(dist_km),
    }


def _map_marker_times(start_dt, legs):
    if start_dt == RADAR_FIX_TIME_UTC:
        times = [_dt_from_label(x) for x in GUIDED_VALIDATION_LABELS]
        return sorted(list({t: t for t in times}.values()))

    times = [start_dt]
    cum = 0.0
    for (_h, _s, m) in legs:
        cum += float(m)
        times.append(start_dt + timedelta(minutes=cum))
    return times


# ================= CLIP TRAJECTORY =================
def clip_legs_to_end_dt(legs, start_dt: datetime, end_dt: datetime):
    total_target_sec = (end_dt - start_dt).total_seconds()
    if total_target_sec <= 0:
        return []

    out = []
    remaining = float(total_target_sec)

    for (hdg, spd, minutes) in legs:
        sec = float(minutes) * 60.0
        if remaining <= 0:
            break
        if sec <= remaining + 1e-9:
            out.append((float(hdg), float(spd), float(minutes)))
            remaining -= sec
        else:
            trunc_min = max(0.0, remaining / 60.0)
            if trunc_min > 1e-6:
                out.append((float(hdg), float(spd), float(trunc_min)))
            remaining = 0.0
            break

    return out


# ================= BFO ANALYSIS =================
def _ping_points_after(start_dt: datetime):
    pts = []
    for row in fp.satLoc:
        label = fp._extract_time_string(row[0])
        meas_bfo = row[SATLOC_IDX_MEAS_BFO]
        meas_bto = row[SATLOC_IDX_MEAS_BTO]
        if meas_bfo is None:
            continue
        try:
            dt_ping = fp._time_label_to_dt(label, fp.SAT_EPOCH_UTC)
        except Exception:
            continue
        if dt_ping < start_dt:
            continue
        pts.append((label, dt_ping, float(meas_bfo), None if meas_bto is None else float(meas_bto)))
    return pts


def _descent_start_dt_from_choice(descent_start_choice: str):
    ch = (descent_start_choice or "").strip()
    if ch.startswith("At 00:20Z only"):
        return dt_0020_utc()
    if ch.startswith("None"):
        return None
    return None


def bfo_analysis_with_legs(start_lat, start_lon, legs, method, start_dt, alt_ft, vs_fpm, descent_start_choice):
    pts = _ping_points_after(start_dt)
    descent_start_dt = _descent_start_dt_from_choice(descent_start_choice)

    leg_secs = [float(m) * 60.0 for _, _, m in legs]
    boundaries = np.cumsum([0.0] + leg_secs)
    total_end = float(boundaries[-1])
    leg_start_positions = _precompute_leg_starts(start_lat, start_lon, legs, method)

    x_labels, measured, calculated = [], [], []

    dt_0020 = dt_0020_utc()
    have_0020 = any((lbl.replace("Z", "") in {"24:20:00", "00:20:00"}) for (lbl, *_rest) in pts)
    if (start_dt <= dt_0020 <= (start_dt + timedelta(seconds=total_end))) and not have_0020:
        pts.append(("24:20:00", dt_0020, float(MEASURED_BFO_0020_HZ), None))

    pts = sorted(pts, key=lambda x: x[1])

    for label, dt_ping, meas_bfo, _meas_bto in pts:
        t_sec = (dt_ping - start_dt).total_seconds()
        if t_sec < 0 or t_sec > total_end:
            continue

        i = int(np.searchsorted(boundaries, t_sec, side="right") - 1)
        i = max(0, min(i, len(legs) - 1))
        t_inside = float(t_sec) - float(boundaries[i])

        leg_start_lat, leg_start_lon = leg_start_positions[i]
        hdg, spd, _m = legs[i]
        lat, lon = _predict_state_from_leg_start(method, leg_start_lat, leg_start_lon, hdg, spd, t_inside)

        alt_now_ft = altitude_ft_at_dt(dt_ping, start_dt, float(alt_ft), float(vs_fpm), descent_start_dt)
        alt_m = alt_now_ft * 0.3048

        vs_mps = float(vs_fpm) * 0.00508
        if descent_start_dt is not None and dt_ping < descent_start_dt:
            vs_mps = 0.0

        spd_mps = float(spd) * 1852.0 / 3600.0
        pred_bfo = float(fp.calculateBFO_at_time(spd_mps, float(hdg), float(lat), float(lon), alt_m, vs_mps, dt_ping))

        x_labels.append(label if label.endswith("Z") else (label + "Z"))
        measured.append(float(meas_bfo))
        calculated.append(float(pred_bfo))

    fig = go.Figure()
    if x_labels:
        fig.add_trace(go.Scatter(
            x=x_labels, y=measured,
            mode="markers",
            name="Measured BFO",
            marker=dict(size=11, symbol="circle")
        ))
        fig.add_trace(go.Scatter(
            x=x_labels, y=calculated,
            mode="lines+markers",
            name="Calculated BFO",
            line=dict(width=3),
            marker=dict(size=10, symbol="diamond")
        ))
    fig.update_layout(
        title="BFO Analysis",
        xaxis_title="Time",
        yaxis_title="Hz",
        legend=dict(x=1.02, y=1.0),
        margin=dict(l=60, r=160, t=80, b=60),
        height=520,
    )
    return fig, ""


# ================= VALIDATION TABLE =================
def _lookup_measured_at_dt(dt_ping: datetime):
    if hasattr(fp, "nearest_measured_for_time"):
        idx = fp.nearest_measured_for_time(dt_ping)[0]
        row = fp.satLoc[idx]
        return row[SATLOC_IDX_MEAS_BFO], row[SATLOC_IDX_MEAS_BTO]
    return None, None


def build_validation_table(start_lat, start_lon, method, start_dt, legs, alt_ft, vs_fpm, descent_start_choice):
    descent_start_dt = _descent_start_dt_from_choice(descent_start_choice)
    is_guided = (start_dt == RADAR_FIX_TIME_UTC)

    if is_guided:
        times = [_dt_from_label(x) for x in GUIDED_VALIDATION_LABELS]
        times = sorted(list({t: t for t in times}.values()))
    else:
        times = [start_dt]
        cum = 0.0
        for (_h, _s, m) in legs:
            cum += float(m)
            times.append(start_dt + timedelta(minutes=cum))

    rows = []
    for dt_k in times:
        st = _state_at_dt(start_lat, start_lon, legs, method, start_dt, dt_k)
        if st is None:
            continue

        lat = st["lat"]
        lon = st["lon"]
        spd = st["speed_kt"]
        hdg = st["heading_deg"]

        alt_now_ft = altitude_ft_at_dt(dt_k, start_dt, float(alt_ft), float(vs_fpm), descent_start_dt)
        alt_m = alt_now_ft * 0.3048

        vs_mps = float(vs_fpm) * 0.00508
        if descent_start_dt is not None and dt_k < descent_start_dt:
            vs_mps = 0.0

        spd_mps = float(spd) * 1852.0 / 3600.0

        bfo_calc = float(fp.calculateBFO_at_time(spd_mps, float(hdg), float(lat), float(lon), alt_m, vs_mps, dt_k))
        bto_calc = float(fp.calculateBTO_at_time(float(lat), float(lon), alt_m, dt_k))

        bfo_meas, bto_meas = _lookup_measured_at_dt(dt_k)

        try:
            if dt_k == dt_0020_utc():
                bfo_meas = float(MEASURED_BFO_0020_HZ)
        except Exception:
            pass

        bfo_err = None if bfo_meas is None else (float(bfo_calc) - float(bfo_meas))
        bto_err = None if bto_meas is None else (float(bto_calc) - float(bto_meas))

        rows.append([
            _to_label_utc(dt_k),
            f"{lat:.5f}",
            f"{lon:.5f}",
            f"{float(spd):.1f}",
            f"{float(hdg):.1f}",
            f"{bfo_calc:.2f}",
            "" if bfo_meas is None else f"{float(bfo_meas):.2f}",
            "" if bfo_err is None else f"{float(bfo_err):.2f}",
            f"{bto_calc:.1f}",
            "" if bto_meas is None else f"{float(bto_meas):.1f}",
            "" if bto_err is None else f"{float(bto_err):.1f}",
        ])

    headers = [
        "time_utc", "lat", "lon", "speed_kt", "heading_deg",
        "bfo_calc_hz", "bfo_meas_hz", "bfo_error_hz",
        "bto_calc_us", "bto_meas_us", "bto_error_us",
    ]
    return headers, rows


# ================= MAP HELPERS =================
def _all_bounds(points_list):
    lats, lons = [], []
    for pts in points_list:
        for (la, lo) in pts:
            if la is None or lo is None:
                continue
            lats.append(float(la))
            lons.append(float(lo))
    if not lats:
        return None
    return [[min(lats), min(lons)], [max(lats), max(lons)]]


def add_permanent_label(m, name, lat, lon):
    tip = f"{html.escape(name)}<br>Lat: {float(lat):.5f}<br>Lon: {float(lon):.5f}"
    pop = folium.Popup(f"<b>{html.escape(name)}</b><br>Lat: {float(lat):.5f}<br>Lon: {float(lon):.5f}", max_width=260)

    folium.CircleMarker(
        location=(lat, lon),
        radius=4,
        color="black",
        weight=2,
        fill=True,
        fill_color="white",
        fill_opacity=1.0,
        tooltip=folium.Tooltip(tip, sticky=True),
        popup=pop,
    ).add_to(m)

    label_html = f"""
    <div style="
        font-size: 13px;
        font-weight: 700;
        color: black;
        text-shadow: -1px -1px 0 #fff, 1px -1px 0 #fff, -1px 1px 0 #fff, 1px 1px 0 #fff;
        white-space: nowrap;
        transform: translate(8px, -12px);
        pointer-events: none;
    ">{html.escape(name)}</div>
    """
    folium.Marker(location=(lat, lon), icon=folium.DivIcon(html=label_html)).add_to(m)


def add_named_marker(m, name: str, lat: float, lon: float):
    tip = f"{html.escape(name)}<br>Lat: {float(lat):.5f}<br>Lon: {float(lon):.5f}"
    pop = folium.Popup(f"<b>{html.escape(name)}</b><br>Lat: {float(lat):.5f}<br>Lon: {float(lon):.5f}", max_width=260)
    folium.CircleMarker(
        location=(float(lat), float(lon)),
        radius=6,
        color="white",
        weight=2,
        fill=True,
        fill_color="cyan",
        fill_opacity=0.9,
        tooltip=folium.Tooltip(tip, sticky=True),
        popup=pop,
    ).add_to(m)


def add_max_range_features(m, max_range_features):
    bounds_pts = []

    for feat in max_range_features:
        ftype = feat.get("type")
        name = feat.get("name", "Maximum range")

        if ftype in ("line", "polygon"):
            pts = feat.get("points", []) or []
            if not pts:
                continue
            bounds_pts.extend(pts)
            folium.PolyLine(
                locations=pts,
                color="red",
                weight=3,
                opacity=0.95,
                tooltip=name,
                popup=folium.Popup(f"<b>{html.escape(name)}</b>", max_width=260),
            ).add_to(m)

        elif ftype == "point":
            lat = float(feat["lat"])
            lon = float(feat["lon"])
            bounds_pts.append((lat, lon))
            add_named_marker(m, name, lat, lon)

    return bounds_pts


def build_leaflet_map(path_pts, legs, method, start_dt, show_max_range_ring=False):
    start_lat, start_lon = path_pts[0]
    is_guided = (start_dt == RADAR_FIX_TIME_UTC)

    m = folium.Map(
        location=[start_lat, start_lon],
        zoom_start=4,
        tiles=None,
        control_scale=True,
    )

    folium.TileLayer("OpenStreetMap", name="OpenStreetMap", overlay=False, control=True).add_to(m)
    folium.TileLayer(
        tiles="https://server.arcgisonline.com/ArcGIS/rest/services/Ocean/World_Ocean_Base/MapServer/tile/{z}/{y}/{x}",
        attr="Esri",
        name="Esri Ocean",
        overlay=False,
        control=True,
    ).add_to(m)
    folium.TileLayer(
        tiles="https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}",
        attr="Esri",
        name="Esri Satellite",
        overlay=False,
        control=True,
    ).add_to(m)

    MousePosition(
        position="topright",
        separator=" | ",
        prefix="Lat/Lon:",
        lat_formatter="function(num) {return L.Util.formatNum(num, 5);}",
        lng_formatter="function(num) {return L.Util.formatNum(num, 5);}",
    ).add_to(m)

    search_area_pts_for_bounds = []
    for area in SEARCH_AREAS:
        nm = area.get("name", "Search area")
        col = color_for_search_area(nm)

        for (la, lo) in area.get("points", []) or []:
            search_area_pts_for_bounds.append((la, lo))
            folium.CircleMarker(
                location=(la, lo),
                radius=6,
                color=col,
                weight=2,
                fill=True,
                fill_color=col,
                fill_opacity=0.85,
                tooltip=nm,
            ).add_to(m)

        for line in area.get("lines", []) or []:
            if not line:
                continue
            search_area_pts_for_bounds.extend(line)
            folium.PolyLine(locations=line, color=col, weight=3, opacity=0.9, tooltip=nm).add_to(m)

        for poly in area.get("polygons", []) or []:
            if not poly:
                continue
            search_area_pts_for_bounds.extend(poly)
            folium.Polygon(
                locations=poly,
                color=col,
                weight=2,
                opacity=0.9,
                fill=True,
                fill_color=col,
                fill_opacity=0.18,
                tooltip=nm,
            ).add_to(m)

    for ring_name in sorted(PING_RINGS.keys()):
        ring_pts = PING_RINGS.get(ring_name)
        if not ring_pts:
            continue
        ring_time = RING_TIMES.get(ring_name)
        tooltip = f"{ring_name} ({ring_time})" if ring_time else ring_name
        folium.PolyLine(ring_pts, weight=2, opacity=0.6, color="yellow", tooltip=tooltip).add_to(m)

    max_range_pts_for_bounds = []
    if show_max_range_ring:
        max_range_pts_for_bounds = add_max_range_features(m, MAX_RANGE_FEATURES)

    folium.PolyLine(path_pts, color="white", weight=8, opacity=0.3).add_to(m)
    folium.PolyLine(path_pts, color="red", weight=4, opacity=1.0, tooltip="Flight path").add_to(m)

    folium.CircleMarker(
        location=path_pts[0],
        radius=7,
        fill=True,
        color="green",
        fill_opacity=1.0,
        popup=folium.Popup(
            f"<b>Start</b><br>Lat: {path_pts[0][0]:.5f}<br>Lon: {path_pts[0][1]:.5f}",
            max_width=320,
        ),
    ).add_to(m)

    if not is_guided:
        folium.CircleMarker(
            location=path_pts[-1],
            radius=7,
            fill=True,
            color="purple",
            fill_opacity=1.0,
            popup=folium.Popup(
                f"<b>End</b><br>Lat: {path_pts[-1][0]:.5f}<br>Lon: {path_pts[-1][1]:.5f}",
                max_width=320,
            ),
        ).add_to(m)

    marker_times = _map_marker_times(start_dt, legs)
    marker_pts_for_bounds = []

    for dt_k in marker_times:
        st = _state_at_dt(path_pts[0][0], path_pts[0][1], legs, method, start_dt, dt_k)
        if st is None:
            continue

        lat = st["lat"]
        lon = st["lon"]
        marker_pts_for_bounds.append((lat, lon))

        tlabel = _to_label_utc(dt_k)
        elapsed_min = st["elapsed_min"]
        hdg = st["heading_deg"]
        spd = st["speed_kt"]
        dist_km = st["dist_km"]

        is_final = (start_dt == RADAR_FIX_TIME_UTC) and (dt_k == _dt_from_label(GUIDED_FINAL_LABEL))

        popup_lines = [
            f"<b>{'FINAL RESTING POINT' if is_final else 'Validation point'}</b>",
            f"UTC time: {tlabel}",
            f"Lat: {lat:.5f}",
            f"Lon: {lon:.5f}",
            f"Time elapsed: {elapsed_min:.1f} min",
            f"Heading: {hdg:.1f}°",
            f"Speed: {spd:.1f} kt",
            f"Distance flown: {dist_km:.2f} km",
        ]

        if is_final:
            folium.CircleMarker(
                location=(lat, lon),
                radius=9,
                fill=True,
                color="magenta",
                fill_color="magenta",
                fill_opacity=0.95,
                tooltip=f"{tlabel} (FINAL)",
                popup=folium.Popup("<br>".join(popup_lines), max_width=360),
            ).add_to(m)
        else:
            folium.CircleMarker(
                location=(lat, lon),
                radius=6,
                fill=True,
                color="white",
                fill_color="white",
                fill_opacity=1.0,
                tooltip=tlabel,
                popup=folium.Popup("<br>".join(popup_lines), max_width=360),
            ).add_to(m)

    for name, lat, lon in WAYPOINTS:
        add_permanent_label(m, float(name) if False else name, float(lat), float(lon))

    add_named_marker(m, CAPT_SIMON_HARDY_NAME, CAPT_SIMON_HARDY_LAT, CAPT_SIMON_HARDY_LON)

    waypoint_pts = [(float(lat), float(lon)) for _, lat, lon in WAYPOINTS] + [(CAPT_SIMON_HARDY_LAT, CAPT_SIMON_HARDY_LON)]
    all_ring_pts = []
    for k in PING_RINGS:
        all_ring_pts.extend(PING_RINGS.get(k, []))

    bounds = _all_bounds(
        [path_pts]
        + [waypoint_pts]
        + [all_ring_pts]
        + ([search_area_pts_for_bounds] if search_area_pts_for_bounds else [])
        + ([marker_pts_for_bounds] if marker_pts_for_bounds else [])
        + ([max_range_pts_for_bounds] if max_range_pts_for_bounds else [])
    )
    if bounds:
        m.fit_bounds(bounds, padding=(20, 20))

    map_html = m.get_root().render()
    escaped = html.escape(map_html, quote=True)
    return f"""
    <iframe
        srcdoc="{escaped}"
        style="width: 100%; height: 650px; border: 0; border-radius: 12px;"
        loading="lazy"
        referrerpolicy="no-referrer"
    ></iframe>
    """


# ================= Guided helpers =================
def _clamp(x, lo, hi):
    return max(lo, min(hi, x))


def _frange(lo: float, hi: float, step: float):
    step = float(step)
    if step <= 0:
        return []
    n = int(math.floor((hi - lo) / step + 1e-9)) + 1
    out = [lo + i * step for i in range(max(0, n))]
    return [x for x in out if x <= hi + 1e-9]


def _refine_range(center: float, half_window: float, step: float, lo: float, hi: float):
    lo2 = _clamp(center - half_window, lo, hi)
    hi2 = _clamp(center + half_window, lo, hi)
    return _frange(lo2, hi2, step)


def _guided_start_dt():
    return RADAR_FIX_TIME_UTC


def _lookup_measured_bto(dt_ping: datetime):
    _bfo, bto = _lookup_measured_at_dt(dt_ping)
    return None if bto is None else float(bto)


def _available_hard_bto_labels():
    out = []
    for lbl in HARD_BTO_LABELS:
        dt = _dt_from_label(lbl)
        if _lookup_measured_bto(dt) is not None:
            out.append(lbl)
    return out


def _simulate_one_leg_end(method: str, lat0: float, lon0: float, hdg: float, spd_kt: float, minutes: float):
    dist_m = float(spd_kt) * 1852.0 * (float(minutes) / 60.0)
    if method == "VINCENTY":
        (lat1, lon1), _ = fp.vincenty_direct((float(lat0), float(lon0)), float(hdg), dist_m)
    else:
        lat1, lon1 = fp.rhumb_direct((float(lat0), float(lon0)), float(hdg), dist_m)
    return float(lat1), float(lon1)


def _bto_calc(method: str, lat: float, lon: float, dt_ping: datetime, alt_ft: float,
              vs_fpm: float, start_dt: datetime, descent_start_dt):
    alt_now_ft = altitude_ft_at_dt(dt_ping, start_dt, float(alt_ft), float(vs_fpm), descent_start_dt)
    alt_m = alt_now_ft * 0.3048
    return float(fp.calculateBTO_at_time(float(lat), float(lon), alt_m, dt_ping))


def _golden_min_abs(fn, lo: float, hi: float, iters: int = 24):
    grc = (math.sqrt(5.0) - 1.0) / 2.0
    a, b = float(lo), float(hi)
    c = b - grc * (b - a)
    d = a + grc * (b - a)
    fc = abs(float(fn(c)))
    fd = abs(float(fn(d)))
    for _ in range(int(iters)):
        if fc <= fd:
            b, d, fd = d, c, fc
            c = b - grc * (b - a)
            fc = abs(float(fn(c)))
        else:
            a, c, fc = c, d, fd
            d = a + grc * (b - a)
            fd = abs(float(fn(d)))
    x_best = (a + b) / 2.0
    return float(x_best), abs(float(fn(x_best)))


def _fit_segment_speed(method: str, lat0: float, lon0: float, hdg: float, minutes: float,
                       dt_end: datetime, alt_ft: float, vs_fpm: float, start_dt: datetime,
                       descent_start_dt, vmin: float, vmax: float):
    meas = _lookup_measured_bto(dt_end)
    if meas is None:
        v = 0.5 * (vmin + vmax)
        lat1, lon1 = _simulate_one_leg_end(method, lat0, lon0, hdg, v, minutes)
        return float(v), float("nan"), (lat1, lon1)

    def err(v):
        lat1, lon1 = _simulate_one_leg_end(method, lat0, lon0, hdg, v, minutes)
        bto = _bto_calc(method, lat1, lon1, dt_end, alt_ft, vs_fpm, start_dt, descent_start_dt)
        return float(bto - meas)

    e_lo = float(err(vmin))
    e_hi = float(err(vmax))

    if (e_lo < 0.0 and e_hi > 0.0) or (e_lo > 0.0 and e_hi < 0.0):
        a, b = float(vmin), float(vmax)
        fa, fb = e_lo, e_hi
        for _ in range(26):
            m = 0.5 * (a + b)
            fm = float(err(m))
            if abs(fm) < 0.2:
                a = b = m
                break
            if (fa < 0 and fm > 0) or (fa > 0 and fm < 0):
                b, fb = m, fm
            else:
                a, fa = m, fm
        v_best = 0.5 * (a + b)
        lat1, lon1 = _simulate_one_leg_end(method, lat0, lon0, hdg, v_best, minutes)
        return float(v_best), float(err(v_best)), (lat1, lon1)

    v_best, _ = _golden_min_abs(err, vmin, vmax, iters=18)
    lat1, lon1 = _simulate_one_leg_end(method, lat0, lon0, hdg, v_best, minutes)
    return float(v_best), float(err(v_best)), (lat1, lon1)


def _compute_bto_errors_at_enforced(method: str, start_dt: datetime, alt_ft: float, vs_fpm: float,
                                    descent_start_dt, full_legs, enforced_dts):
    errs = []
    for dt_k in enforced_dts:
        meas = _lookup_measured_bto(dt_k)
        if meas is None:
            continue
        st = _state_at_dt(RADAR_FIX_LAT, RADAR_FIX_LON, full_legs, method, start_dt, dt_k)
        if st is None:
            continue
        bto_calc = _bto_calc(method, st["lat"], st["lon"], dt_k, alt_ft, vs_fpm, start_dt, descent_start_dt)
        errs.append(float(bto_calc - meas))
    return errs


def _build_south_legs_piecewise_constant_heading(h_south: float, speeds, minutes):
    legs = []
    for v, m in zip(speeds, minutes):
        if float(m) <= 0:
            continue
        legs.append((float(h_south), float(v), float(m)))
    return legs


def _build_south_legs_variable(headings, speeds, minutes):
    legs = []
    for h, v, m in zip(headings, speeds, minutes):
        if float(m) <= 0:
            continue
        legs.append((float(h), float(v), float(m)))
    return legs


# ================= Top-K polish =================
def _polish_south_speeds_topk(method: str, start_dt: datetime, end_dt: datetime, alt_ft: float,
                              vs_fpm: float, descent_start_dt, h_south: float, turn_legs,
                              enforced_dts, init_speeds, passes: int = POLISH_PASSES,
                              window_kt: float = POLISH_WINDOW_KT, step_kt: float = POLISH_STEP_KT):
    speeds = [float(v) for v in init_speeds]

    pts_turn, _, _ = simulate_path(RADAR_FIX_LAT, RADAR_FIX_LON, turn_legs, method)
    turn_end_lat, turn_end_lon = float(pts_turn[-1][0]), float(pts_turn[-1][1])

    def build_full(speeds_now):
        lat_cur, lon_cur = turn_end_lat, turn_end_lon
        seg_errs = []
        for mins, end_lbl, vseg in zip(GUIDED_SOUTH_SEG_MIN, GUIDED_SOUTH_SEG_LABELS_END, speeds_now):
            dt_end = _dt_from_label(end_lbl)
            lat_cur, lon_cur = _simulate_one_leg_end(method, lat_cur, lon_cur, h_south, vseg, mins)
            meas = _lookup_measured_bto(dt_end)
            if meas is None:
                seg_errs.append(float("nan"))
            else:
                bto_here = _bto_calc(method, lat_cur, lon_cur, dt_end, alt_ft, vs_fpm, start_dt, descent_start_dt)
                seg_errs.append(float(bto_here - meas))

        south_legs = _build_south_legs_piecewise_constant_heading(h_south, speeds_now, GUIDED_SOUTH_SEG_MIN)
        full_legs = clip_legs_to_end_dt(turn_legs + south_legs, start_dt, end_dt)
        return full_legs, seg_errs

    def objective(full_legs):
        bto_errs = _compute_bto_errors_at_enforced(method, start_dt, alt_ft, vs_fpm, descent_start_dt, full_legs, enforced_dts)
        if not bto_errs:
            return 1e9
        abs_arr = np.abs(np.array(bto_errs, dtype=float))
        max_abs = float(np.max(abs_arr))
        sum_abs = float(np.sum(abs_arr))
        return float(max_abs + POLISH_OBJ_EPS * sum_abs)

    full_legs, seg_errs = build_full(speeds)
    best_obj = objective(full_legs)
    best_speeds = speeds[:]
    best_full = full_legs
    best_seg_errs = seg_errs

    for _p in range(int(passes)):
        improved_any = False
        for i in range(len(best_speeds)):
            v0 = float(best_speeds[i])
            v_lo = _clamp(v0 - float(window_kt), SOUTH_SPEED_MIN_KT, SOUTH_SPEED_MAX_KT)
            v_hi = _clamp(v0 + float(window_kt), SOUTH_SPEED_MIN_KT, SOUTH_SPEED_MAX_KT)

            cand = list(np.arange(v_lo, v_hi + 1e-9, float(step_kt)))
            if not cand:
                continue

            local_best_obj = best_obj
            local_best_speeds = best_speeds[:]
            local_best_full = best_full
            local_best_seg_errs = best_seg_errs

            for v_try in cand:
                sp = best_speeds[:]
                sp[i] = float(v_try)
                full_try, seg_try = build_full(sp)
                obj_try = objective(full_try)
                if obj_try < local_best_obj - 1e-9:
                    local_best_obj = obj_try
                    local_best_speeds = sp
                    local_best_full = full_try
                    local_best_seg_errs = seg_try

            if local_best_obj < best_obj - 1e-9:
                best_obj = local_best_obj
                best_speeds = local_best_speeds
                best_full = local_best_full
                best_seg_errs = local_best_seg_errs
                improved_any = True

        if not improved_any:
            break

    return best_speeds, best_full, best_seg_errs


# ================= Guided Solver: Straight heading =================
def solve_guided_straight(method, alt_ft, vs_fpm, descent_start_choice):
    start_dt = _guided_start_dt()
    end_dt = _dt_from_label(GUIDED_FINAL_LABEL)
    descent_start_dt = _descent_start_dt_from_choice(descent_start_choice)

    pts_all = _ping_points_after(start_dt)
    pts_window = [(lbl, dt, bfo, bto) for (lbl, dt, bfo, bto) in pts_all if start_dt <= dt <= end_dt]
    if len(pts_window) < 2:
        return None, "Guided solve: not enough ping points to score.", None

    enforced_labels = _available_hard_bto_labels()
    if not enforced_labels:
        return None, "Guided solve: no measured BTO available at any enforced ring times.", None

    enforced_dts = [_dt_from_label(lbl) for lbl in enforced_labels]

    def _evaluate(h_south: float, v_turn: float, do_polish: bool):
        h_south = float(h_south)
        v_turn = float(v_turn)

        turn_legs, _turn_info = build_turn_legs(speed_knots=float(v_turn), method=method, heading_final_deg=float(h_south))
        pts_turn, _, _ = simulate_path(RADAR_FIX_LAT, RADAR_FIX_LON, turn_legs, method)
        lat_cur, lon_cur = float(pts_turn[-1][0]), float(pts_turn[-1][1])

        south_speeds = []
        south_errs = []
        for mins, end_lbl in zip(GUIDED_SOUTH_SEG_MIN, GUIDED_SOUTH_SEG_LABELS_END):
            dt_end = _dt_from_label(end_lbl)
            v_seg, err_seg, (lat1, lon1) = _fit_segment_speed(
                method=method,
                lat0=lat_cur,
                lon0=lon_cur,
                hdg=float(h_south),
                minutes=float(mins),
                dt_end=dt_end,
                alt_ft=alt_ft,
                vs_fpm=vs_fpm,
                start_dt=start_dt,
                descent_start_dt=descent_start_dt,
                vmin=SOUTH_SPEED_MIN_KT,
                vmax=SOUTH_SPEED_MAX_KT,
            )
            south_speeds.append(float(v_seg))
            south_errs.append(float(err_seg))
            lat_cur, lon_cur = float(lat1), float(lon1)

        south_legs = _build_south_legs_piecewise_constant_heading(float(h_south), south_speeds, GUIDED_SOUTH_SEG_MIN)
        full_legs = clip_legs_to_end_dt(turn_legs + south_legs, start_dt, end_dt)

        if do_polish and ENABLE_POLISH_STRAIGHT:
            south_speeds, full_legs, south_errs = _polish_south_speeds_topk(
                method=method,
                start_dt=start_dt,
                end_dt=end_dt,
                alt_ft=alt_ft,
                vs_fpm=vs_fpm,
                descent_start_dt=descent_start_dt,
                h_south=float(h_south),
                turn_legs=turn_legs,
                enforced_dts=enforced_dts,
                init_speeds=south_speeds,
            )

        bto_errs = _compute_bto_errors_at_enforced(method, start_dt, alt_ft, vs_fpm, descent_start_dt, full_legs, enforced_dts)
        if not bto_errs:
            return None

        abs_arr = np.abs(np.array(bto_errs, dtype=float))
        max_abs = float(np.max(abs_arr))
        sum_abs = float(np.sum(abs_arr))
        rmse = float(np.sqrt(np.mean(np.array(bto_errs, dtype=float) ** 2)))

        details = {
            "h_south": float(h_south),
            "v_turn": float(v_turn),
            "south_speeds": [float(v) for v in south_speeds],
            "south_errs_us": [float(e) for e in south_errs],
            "enforced_labels": enforced_labels,
            "max_abs_err_us": max_abs,
            "sum_abs_err_us": sum_abs,
            "rmse_bto_us": rmse,
        }
        return (max_abs, sum_abs, rmse), full_legs, details

    def best_for_heading(h_south: float):
        coarse_turns = list(np.arange(TURN_SPEED_MIN_KT, TURN_SPEED_MAX_KT + 1e-9, 5.0))
        best_key = None
        best_legs = None
        best_det = None

        for v in coarse_turns:
            out = _evaluate(h_south, v, do_polish=False)
            if out is None:
                continue
            key, legs, det = out
            if best_key is None or key < best_key:
                best_key, best_legs, best_det = key, legs, det

        if best_key is None:
            return None

        v0 = float(best_det["v_turn"])
        v_lo = _clamp(v0 - 7.5, TURN_SPEED_MIN_KT, TURN_SPEED_MAX_KT)
        v_hi = _clamp(v0 + 7.5, TURN_SPEED_MIN_KT, TURN_SPEED_MAX_KT)

        def obj(v):
            out2 = _evaluate(h_south, v, do_polish=False)
            if out2 is None:
                return 1e9
            (max_abs, sum_abs, _rmse), _legs, _det = out2
            return float(max_abs) + POLISH_OBJ_EPS * float(sum_abs)

        v_best, _ = _golden_min_abs(obj, v_lo, v_hi, iters=16)
        out3 = _evaluate(h_south, v_best, do_polish=False)
        if out3 is None:
            return best_key, best_legs, best_det
        return out3

    coarse_headings = _frange(GUIDED_HEADING_MIN, GUIDED_HEADING_MAX, COARSE_HEADING_STEP_DEG)
    candidates = []
    best = None
    best_legs = None
    best_det = None

    for h in coarse_headings:
        out = best_for_heading(h)
        if out is None:
            continue
        key, legs, det = out
        candidates.append((key, float(det["h_south"]), float(det["v_turn"])))
        if best is None or key < best:
            best, best_legs, best_det = key, legs, det

    if best_legs is None:
        return None, "Guided solve: no solution found within bounds.", None

    h0 = float(best_det["h_south"])
    refine_headings = _refine_range(h0, REFINE_HEADING_HALF_WINDOW_DEG, REFINE_HEADING_STEP_DEG, GUIDED_HEADING_MIN, GUIDED_HEADING_MAX)
    for h in refine_headings:
        out = best_for_heading(h)
        if out is None:
            continue
        key, legs, det = out
        candidates.append((key, float(det["h_south"]), float(det["v_turn"])))
        if best is None or key < best:
            best, best_legs, best_det = key, legs, det

    candidates_sorted = sorted(candidates, key=lambda x: x[0])
    uniq = []
    seen = set()
    for key, hh, vv in candidates_sorted:
        k2 = (round(hh, 4), round(vv, 4))
        if k2 in seen:
            continue
        seen.add(k2)
        uniq.append((key, hh, vv))
        if len(uniq) >= max(1, int(TOPK_POLISH)):
            break

    best_final_key, best_final_legs, best_final_det = best, best_legs, best_det
    if ENABLE_POLISH_STRAIGHT and uniq:
        for _key, hh, vv in uniq:
            outp = _evaluate(hh, vv, do_polish=True)
            if outp is None:
                continue
            keyp, legsp, detp = outp
            if best_final_key is None or keyp < best_final_key:
                best_final_key, best_final_legs, best_final_det = keyp, legsp, detp

    return best_final_legs, "", start_dt


# ================= Guided Solver: Variable heading & speed =================
def _solve_guided_variable_core(method, alt_ft, vs_fpm, descent_start_choice, se_bias_after_2240: bool):
    start_dt = _guided_start_dt()
    end_dt = _dt_from_label(GUIDED_FINAL_LABEL)
    descent_start_dt = _descent_start_dt_from_choice(descent_start_choice)

    enforced_labels = _available_hard_bto_labels()
    if not enforced_labels:
        return None, "Guided solve: no measured BTO available at any enforced ring times.", None
    enforced_dts = [_dt_from_label(lbl) for lbl in enforced_labels]

    turn_speeds = list(np.arange(TURN_SPEED_MIN_KT, TURN_SPEED_MAX_KT + 1e-9, VAR_TURN_SPEED_STEP))
    first_headings = _frange(VAR_HEADING_MIN, VAR_HEADING_MAX, VAR_FIRST_HEADING_STEP)

    best_global = None
    best_global_legs = None

    def score_from_errors(errs):
        if not errs:
            return None
        abs_arr = np.abs(np.array(errs, dtype=float))
        return float(np.max(abs_arr)), float(np.sum(abs_arr)), float(np.sqrt(np.mean(abs_arr * abs_arr)))

    for v_turn in turn_speeds:
        for h1 in first_headings:
            turn_legs, _turn_info = build_turn_legs(speed_knots=float(v_turn), method=method, heading_final_deg=float(h1))
            pts_turn, _, _ = simulate_path(RADAR_FIX_LAT, RADAR_FIX_LON, turn_legs, method)
            lat0, lon0 = float(pts_turn[-1][0]), float(pts_turn[-1][1])

            beam = [{
                "lat": lat0,
                "lon": lon0,
                "headings": [],
                "speeds": [],
                "h_prev": float(h1),
                "v_turn": float(v_turn),
                "turn_legs": turn_legs,
                "max_abs": 0.0,
                "sum_abs": 0.0,
                "smooth_sum": 0.0,
                "east_bias_sum": 0.0,
            }]

            seg_i = -1
            for mins, end_lbl in zip(GUIDED_SOUTH_SEG_MIN, GUIDED_SOUTH_SEG_LABELS_END):
                seg_i += 1
                dt_end = _dt_from_label(end_lbl)
                meas = _lookup_measured_bto(dt_end)
                next_beam = []

                use_se_bias = bool(se_bias_after_2240 and seg_i >= SE_EAST_DIVERGE_START_SEG_INDEX)
                delta_set = SE_VAR_EAST_DELTA_SET if use_se_bias else VAR_HEADING_DELTA_SET

                for st in beam:
                    h_prev = float(st["h_prev"])

                    cand_h = []
                    for d in delta_set:
                        h_try = h_prev + float(d)
                        h_try = _clamp(h_try, VAR_HEADING_MIN, VAR_HEADING_MAX)
                        cand_h.append(h_try)
                    cand_h = sorted(list({round(x, 6): x for x in cand_h}.values()))

                    for h_try in cand_h:
                        v_seg, err_seg, (lat1, lon1) = _fit_segment_speed(
                            method=method,
                            lat0=float(st["lat"]),
                            lon0=float(st["lon"]),
                            hdg=float(h_try),
                            minutes=float(mins),
                            dt_end=dt_end,
                            alt_ft=alt_ft,
                            vs_fpm=vs_fpm,
                            start_dt=start_dt,
                            descent_start_dt=descent_start_dt,
                            vmin=SOUTH_SPEED_MIN_KT,
                            vmax=SOUTH_SPEED_MAX_KT,
                        )

                        abs_err = 0.0 if meas is None or (err_seg != err_seg) else abs(float(err_seg))
                        smooth_pen = abs(float(h_try) - float(h_prev)) * VAR_SMOOTH_PENALTY_US_PER_DEG

                        east_bias_pen = 0.0
                        if use_se_bias:
                            dh = float(h_try) - float(h_prev)
                            if dh > 0.0:
                                east_bias_pen = dh * float(SE_EAST_BIAS_PENALTY_US_PER_DEG)

                        max_abs = max(float(st["max_abs"]), float(abs_err))
                        sum_abs = float(st["sum_abs"]) + float(abs_err)
                        smooth_sum = float(st["smooth_sum"]) + float(smooth_pen)
                        east_bias_sum = float(st.get("east_bias_sum", 0.0)) + float(east_bias_pen)

                        next_beam.append({
                            "lat": float(lat1),
                            "lon": float(lon1),
                            "headings": st["headings"] + [float(h_try)],
                            "speeds": st["speeds"] + [float(v_seg)],
                            "h_prev": float(h_try),
                            "v_turn": float(st["v_turn"]),
                            "turn_legs": st["turn_legs"],
                            "max_abs": max_abs,
                            "sum_abs": sum_abs,
                            "smooth_sum": smooth_sum,
                            "east_bias_sum": east_bias_sum,
                        })

                if not next_beam:
                    beam = []
                    break

                next_beam.sort(key=lambda x: (x["max_abs"], x["sum_abs"], x.get("east_bias_sum", 0.0), x["smooth_sum"]))
                beam = next_beam[:max(1, int(VAR_BEAM_WIDTH))]

            if not beam:
                continue

            for st in beam:
                south_legs = _build_south_legs_variable(st["headings"], st["speeds"], GUIDED_SOUTH_SEG_MIN)
                full_legs = clip_legs_to_end_dt(st["turn_legs"] + south_legs, start_dt, end_dt)

                bto_errs = _compute_bto_errors_at_enforced(method, start_dt, alt_ft, vs_fpm, descent_start_dt, full_legs, enforced_dts)
                if not bto_errs:
                    continue

                max_abs, sum_abs, rmse = score_from_errors(bto_errs)
                key = (max_abs, sum_abs, float(st.get("east_bias_sum", 0.0)), float(st["smooth_sum"]), rmse)

                if best_global is None or key < best_global:
                    best_global = key
                    best_global_legs = full_legs

    if best_global_legs is None:
        return None, "Guided solve: Variable mode found no solution within bounds.", None

    return best_global_legs, "", start_dt


def solve_guided_variable(method, alt_ft, vs_fpm, descent_start_choice):
    return _solve_guided_variable_core(method, alt_ft, vs_fpm, descent_start_choice, se_bias_after_2240=False)


def solve_guided_variable_se(method, alt_ft, vs_fpm, descent_start_choice):
    return _solve_guided_variable_core(method, alt_ft, vs_fpm, descent_start_choice, se_bias_after_2240=True)


# ================= Guided Solver: Variable heading, CONSTANT southern speed =================
def solve_guided_constant_speed(method, alt_ft, vs_fpm, descent_start_choice):
    start_dt = _guided_start_dt()
    end_dt = _dt_from_label(GUIDED_FINAL_LABEL)
    descent_start_dt = _descent_start_dt_from_choice(descent_start_choice)

    enforced_labels = _available_hard_bto_labels()
    if not enforced_labels:
        return None, "Guided solve: no measured BTO available at any enforced ring times.", None
    enforced_dts = [_dt_from_label(lbl) for lbl in enforced_labels]

    turn_speeds = list(np.arange(TURN_SPEED_MIN_KT, TURN_SPEED_MAX_KT + 1e-9, VAR_TURN_SPEED_STEP))
    first_headings = _frange(VAR_HEADING_MIN, VAR_HEADING_MAX, VAR_FIRST_HEADING_STEP)

    best_global = None
    best_global_legs = None

    def score_from_errors(errs):
        if not errs:
            return None
        abs_arr = np.abs(np.array(errs, dtype=float))
        return float(np.max(abs_arr)), float(np.sum(abs_arr)), float(np.sqrt(np.mean(abs_arr * abs_arr)))

    v_const = float(CONST_SOUTH_SPEED_KT)

    for v_turn in turn_speeds:
        for h1 in first_headings:
            turn_legs, _turn_info = build_turn_legs(speed_knots=float(v_turn), method=method, heading_final_deg=float(h1))
            pts_turn, _, _ = simulate_path(RADAR_FIX_LAT, RADAR_FIX_LON, turn_legs, method)
            lat0, lon0 = float(pts_turn[-1][0]), float(pts_turn[-1][1])

            beam = [{
                "lat": lat0,
                "lon": lon0,
                "headings": [],
                "h_prev": float(h1),
                "v_turn": float(v_turn),
                "turn_legs": turn_legs,
                "max_abs": 0.0,
                "sum_abs": 0.0,
                "smooth_sum": 0.0,
            }]

            for mins, end_lbl in zip(GUIDED_SOUTH_SEG_MIN, GUIDED_SOUTH_SEG_LABELS_END):
                dt_end = _dt_from_label(end_lbl)
                meas = _lookup_measured_bto(dt_end)
                next_beam = []

                for st in beam:
                    h_prev = float(st["h_prev"])
                    cand_h = []
                    for d in VAR_HEADING_DELTA_SET:
                        h_try = _clamp(h_prev + float(d), VAR_HEADING_MIN, VAR_HEADING_MAX)
                        cand_h.append(h_try)
                    cand_h = sorted(list({round(x, 6): x for x in cand_h}.values()))

                    for h_try in cand_h:
                        lat1, lon1 = _simulate_one_leg_end(method, float(st["lat"]), float(st["lon"]), float(h_try), v_const, float(mins))

                        if meas is None:
                            abs_err = 0.0
                        else:
                            bto_here = _bto_calc(method, lat1, lon1, dt_end, alt_ft, vs_fpm, start_dt, descent_start_dt)
                            abs_err = abs(float(bto_here - float(meas)))

                        smooth_pen = abs(float(h_try) - float(h_prev)) * VAR_SMOOTH_PENALTY_US_PER_DEG

                        max_abs = max(float(st["max_abs"]), float(abs_err))
                        sum_abs = float(st["sum_abs"]) + float(abs_err)
                        smooth_sum = float(st["smooth_sum"]) + float(smooth_pen)

                        next_beam.append({
                            "lat": float(lat1),
                            "lon": float(lon1),
                            "headings": st["headings"] + [float(h_try)],
                            "h_prev": float(h_try),
                            "v_turn": float(st["v_turn"]),
                            "turn_legs": st["turn_legs"],
                            "max_abs": max_abs,
                            "sum_abs": sum_abs,
                            "smooth_sum": smooth_sum,
                        })

                if not next_beam:
                    beam = []
                    break

                next_beam.sort(key=lambda x: (x["max_abs"], x["sum_abs"], x["smooth_sum"]))
                beam = next_beam[:max(1, int(VAR_BEAM_WIDTH))]

            if not beam:
                continue

            for st in beam:
                speeds = [v_const] * len(GUIDED_SOUTH_SEG_MIN)
                south_legs = _build_south_legs_variable(st["headings"], speeds, GUIDED_SOUTH_SEG_MIN)
                full_legs = clip_legs_to_end_dt(st["turn_legs"] + south_legs, start_dt, end_dt)

                bto_errs = _compute_bto_errors_at_enforced(method, start_dt, alt_ft, vs_fpm, descent_start_dt, full_legs, enforced_dts)
                if not bto_errs:
                    continue

                max_abs, sum_abs, rmse = score_from_errors(bto_errs)
                key = (max_abs, sum_abs, float(st["smooth_sum"]), rmse)

                if best_global is None or key < best_global:
                    best_global = key
                    best_global_legs = full_legs

    if best_global_legs is None:
        return None, "Guided solve: Constant-speed mode found no solution within bounds.", None

    return best_global_legs, "", start_dt


# ================= RUNNERS =================
def _run_common(start_lat, start_lon, legs, method, start_dt, alt_ft, vs_fpm, descent_start_choice, show_max_range_ring):
    if start_dt == RADAR_FIX_TIME_UTC:
        end_dt = _dt_from_label(GUIDED_FINAL_LABEL)
        legs = clip_legs_to_end_dt(legs, start_dt, end_dt)

    path_pts, seg_dists_m, total_time_min = simulate_path(start_lat, start_lon, legs, method)

    html_map = build_leaflet_map(path_pts, legs, method, start_dt, show_max_range_ring=show_max_range_ring)
    stats_md = compute_stats(path_pts, seg_dists_m, total_time_min, legs)

    bfo_fig, _bfo_md = bfo_analysis_with_legs(start_lat, start_lon, legs, method, start_dt, alt_ft, vs_fpm, descent_start_choice)

    descent_start_dt = _descent_start_dt_from_choice(descent_start_choice)
    alt_fig = build_altitude_plot(start_dt, total_time_min, alt_ft, vs_fpm, descent_start_dt)
    speed_fig = build_speed_profile_plot(start_dt, legs)

    headers, rows = build_validation_table(start_lat, start_lon, method, start_dt, legs, alt_ft, vs_fpm, descent_start_choice)
    return html_map, stats_md, bfo_fig, alt_fig, speed_fig, headers, rows


def run_sim_manual(start_lat, start_lon, legs_df, method, start_time_str, alt_ft, vs_fpm, descent_start_choice, show_max_range_ring):
    legs = parse_legs_df(legs_df)
    start_dt = _parse_start_time(start_time_str)

    html_map, stats_md, bfo_fig, alt_fig, speed_fig, headers, rows = _run_common(
        start_lat, start_lon, legs, method, start_dt, alt_ft, vs_fpm, descent_start_choice, show_max_range_ring
    )

    return gr.update(), "", html_map, stats_md, bfo_fig, alt_fig, speed_fig, gr.update(headers=headers, value=rows)


def run_guided(guided_mode, start_lat, start_lon, method, alt_ft, vs_fpm, descent_start_choice, show_max_range_ring):
    if guided_mode == "Straight heading":
        legs, _report_md, start_dt = solve_guided_straight(method, alt_ft, vs_fpm, descent_start_choice)
    elif guided_mode == "Variable heading & speed":
        legs, _report_md, start_dt = solve_guided_variable(method, alt_ft, vs_fpm, descent_start_choice)
    elif guided_mode == "Variable heading & speed (SE after 22:40)":
        legs, _report_md, start_dt = solve_guided_variable_se(method, alt_ft, vs_fpm, descent_start_choice)
    elif guided_mode.startswith("Variable heading (constant"):
        legs, _report_md, start_dt = solve_guided_constant_speed(method, alt_ft, vs_fpm, descent_start_choice)
    else:
        legs, _report_md, start_dt = solve_guided_straight(method, alt_ft, vs_fpm, descent_start_choice)

    if legs is None:
        fig = go.Figure()
        fig.update_layout(title="BFO Analysis", xaxis_title="Time", yaxis_title="Hz", height=520)
        alt_fig = go.Figure()
        alt_fig.update_layout(title="Altitude vs Time", xaxis_title="Minutes since start", yaxis_title="Altitude (ft)", height=320)
        speed_fig = go.Figure()
        speed_fig.update_layout(title="Speed Profile", xaxis_title="Minutes since start", yaxis_title="Speed (kt)", height=320)
        empty_table = gr.update(headers=[
            "time_utc","lat","lon","speed_kt","heading_deg",
            "bfo_calc_hz","bfo_meas_hz","bfo_error_hz",
            "bto_calc_us","bto_meas_us","bto_error_us"
        ], value=[])
        return gr.update(), "", "", "", fig, alt_fig, speed_fig, empty_table

    html_map, stats_md, bfo_fig, alt_fig, speed_fig, headers, rows = _run_common(
        RADAR_FIX_LAT, RADAR_FIX_LON, legs, method, start_dt, alt_ft, vs_fpm, descent_start_choice, show_max_range_ring
    )

    return gr.update(value=legs_to_rows(legs)), "", html_map, stats_md, bfo_fig, alt_fig, speed_fig, gr.update(headers=headers, value=rows)


# ================= UI =================
with gr.Blocks() as demo:
    gr.Markdown("## ✈️ MH370 Flight Simulator — Leaflet + Esri Ocean/Satellite")

    if KML_LOAD_ERROR:
        gr.Markdown(
            "### ⚠️ KML not loaded\n"
            f"**Error:** {KML_LOAD_ERROR}\n\n"
            "Make sure your KML file is uploaded and `KML_FILE` matches exactly."
        )

    if KMZ_LOAD_ERROR:
        gr.Markdown(
            "### ⚠️ Search areas KMZ not loaded\n"
            f"**Error:** {KMZ_LOAD_ERROR}\n\n"
            "Make sure `Search areas.kmz` is uploaded and `SEARCH_AREAS_KMZ` matches exactly."
        )

    if MAX_RANGE_KML_ERROR:
        gr.Markdown(
            "### ⚠️ Maximum range KML not loaded\n"
            f"**Error:** {MAX_RANGE_KML_ERROR}\n\n"
            "Make sure `MH370_max_range_6166km_KUL_red.kml` is uploaded and `MAX_RANGE_KML` matches exactly."
        )

    guided_mode = gr.Dropdown(
        [
            "Straight heading",
            "Variable heading & speed",
            "Variable heading & speed (SE after 22:40)",
            f"Variable heading (constant {CONST_SOUTH_SPEED_KT:.0f} kt)",
        ],
        value="Straight heading",
        label="Guided scenario",
    )

    with gr.Row():
        start_lat = gr.Slider(-10, 15, value=6.76, label="Start Latitude (Manual only)")
        start_lon = gr.Slider(70, 115, value=95.69, label="Start Longitude (Manual only)")
        method = gr.Dropdown(["VINCENTY", "RHUMB"], value="VINCENTY", label="Method")

    gr.Markdown("### Map layers")
    show_max_range_ring = gr.Checkbox(
        value=False,
        label="Show maximum range ring (6166 km from KUL)"
    )

    gr.Markdown("### BFO/BTO Settings")
    with gr.Row():
        start_time_str = gr.Textbox(value="18:22:00Z", label="Simulation start time (UTC) — Manual mode")
        alt_ft = gr.Number(value=35000, label="Altitude at start (ft)")
        vs_fpm = gr.Number(value=0, label="Vertical speed (ft/min) (negative = descent)")
        descent_start_choice = gr.Dropdown(
            ["None (constant VS from start)", "At 00:20Z only"],
            value="None (constant VS from start)",
            label="When does descent start? (VS applies after this time)",
        )

    with gr.Row():
        guided_btn = gr.Button("Generate: Guided scenario")
        manual_btn = gr.Button("Simulate: Manual legs")

    status_md = gr.Markdown(visible=False)

    out_map = gr.HTML()
    stats = gr.Markdown()

    gr.Markdown("### BFO Analysis")
    bfo_plot = gr.Plot()

    gr.Markdown("### Altitude vs Time")
    alt_plot = gr.Plot()

    gr.Markdown("### Speed Profile")
    speed_plot = gr.Plot()

    gr.Markdown("### Validation table")
    validation_df = gr.Dataframe(
        headers=[
            "time_utc","lat","lon","speed_kt","heading_deg",
            "bfo_calc_hz","bfo_meas_hz","bfo_error_hz",
            "bto_calc_us","bto_meas_us","bto_error_us"
        ],
        value=[],
        interactive=False,
        row_count=(40, "dynamic"),
        column_count=(11, "fixed"),
        wrap=True,
        label="Validation"
    )

    gr.Markdown("### Multi-leg plan (Manual)")
    legs_df = gr.Dataframe(
        headers=["heading_deg", "speed_knots", "minutes"],
        datatype=["number", "number", "number"],
        row_count=(200, "dynamic"),
        column_count=(3, "fixed"),
        value=[row[:] for row in DEFAULT_LEGS],
        label="Legs (manual inputs)",
    )

    with gr.Row():
        add_btn = gr.Button("➕ Add leg")
        reset_btn = gr.Button("♻️ Reset legs")

    guided_btn.click(
        run_guided,
        inputs=[guided_mode, start_lat, start_lon, method, alt_ft, vs_fpm, descent_start_choice, show_max_range_ring],
        outputs=[legs_df, status_md, out_map, stats, bfo_plot, alt_plot, speed_plot, validation_df],
    )

    manual_btn.click(
        run_sim_manual,
        inputs=[start_lat, start_lon, legs_df, method, start_time_str, alt_ft, vs_fpm, descent_start_choice, show_max_range_ring],
        outputs=[legs_df, status_md, out_map, stats, bfo_plot, alt_plot, speed_plot, validation_df],
    )

    add_btn.click(add_leg, inputs=[legs_df], outputs=[legs_df])
    reset_btn.click(reset_legs, outputs=[legs_df])

demo.launch()