Chris773 commited on
Commit
cbd02a3
·
verified ·
1 Parent(s): 749227c

Update app.py

Browse files
Files changed (1) hide show
  1. app.py +468 -436
app.py CHANGED
@@ -1,8 +1,9 @@
1
  # app.py
2
  # MH370 Flight Simulator — Leaflet map (Folium) + Ping Rings (KML)
3
  #
4
- # Guided scenario (ACTIVE):
5
  # A) Straight (single heading)
 
6
  #
7
  # Key behavior:
8
  # - Last radar spot FIXED: 18:22:00Z @ 6.578N, 96.341E
@@ -10,32 +11,24 @@
10
  # - Guided mode computes a realistic turn between 18:22→18:40 and then southern path.
11
  # - Guided validation timestamps: 18:22, 18:40, 19:40, 20:40, 21:40, 22:40, 00:10, 00:20
12
  #
13
- # FIXES IMPLEMENTED:
14
  # 1) Remove trajectory (red line) AFTER 00:20
15
  # - In guided mode, we CLIP the simulated legs so the polyline ends exactly at 00:20Z.
16
- # 2) Restore BFO plot
17
- # - Full BFO analysis computation + chart + MAE/RMSE markdown.
18
  #
19
- # TURN MODEL:
20
  # 3) Plausible final major turn (18:22→18:40) can adjust BOTH heading and speed
 
 
 
 
21
  #
22
- # SPEED MODEL:
23
- # 4) Southern trajectory speed is NOT constant:
24
- # - Solver chooses 2 speeds:
25
- # v1 from 18:40 → 19:40
26
- # v2 from 19:40 → 00:10
27
- # - Then 00:10 → 00:20 uses v2 (for display), but 00:20 is ignored in scoring.
28
- #
29
- # OPTIMIZER MIRROR + PERFORMANCE (THIS UPDATE):
30
- # - Uses bounded optimizer (L-BFGS-B) on [heading, turn speed, v1, v2]
31
- # - FAST seed pruning: pre-score coarse seeds using only 3 pings, keep top K
32
- # - Strict scoring pings: only 18:40,19:40,20:40,21:40,22:40,00:10 (00:20 excluded)
33
- # - Aggressive memoization cache for BFO/BTO calculations (big speed-up)
34
- #
35
- # SCORING NOTE:
36
- # 5) 00:20 (24:20) is IGNORED for scoring:
37
- # - Excluded from guided scoring set.
38
- # - Validation table shows calc at 00:20, but measured BFO/BTO are blanked there.
39
 
40
  import os
41
  import re
@@ -52,9 +45,6 @@ import plotly.graph_objects as go
52
 
53
  import flight_plotter_v32 as fp
54
 
55
- # Requires SciPy installed in your HF Space
56
- from scipy.optimize import minimize
57
-
58
 
59
  # ================= CONFIG =================
60
  KML_FILE = "MH370 - 13 nov (3).kml"
@@ -86,9 +76,9 @@ RING_TIMES = {
86
  "Ping Ring 7": "24:10:00Z",
87
  }
88
 
89
- # Normalization scales (used to normalize residuals)
90
  BFO_SCALE_HZ = 5.0
91
- BTO_SCALE_US = 25.0 # tighter BTO enforcement
92
 
93
  # EXPLICIT WEIGHTS
94
  DEFAULT_WEIGHT_BFO = 1.0
@@ -102,12 +92,15 @@ SATLOC_IDX_MEAS_BTO = 5
102
  RADAR_FIX_LAT = 6.578
103
  RADAR_FIX_LON = 96.341
104
  RADAR_FIX_TIME_UTC = datetime(2014, 3, 7, 18, 22, 0, tzinfo=timezone.utc)
 
105
  TURN_COMPLETE_TIME_UTC = datetime(2014, 3, 7, 18, 40, 0, tzinfo=timezone.utc)
 
106
 
107
  # Initial track shown in your radar annotation
108
  RADAR_INITIAL_TRACK_DEG = 296.0
109
 
110
  # ===== Guided validation timestamps of interest =====
 
111
  GUIDED_VALIDATION_LABELS = [
112
  "18:22:00",
113
  "18:40:00",
@@ -116,68 +109,49 @@ GUIDED_VALIDATION_LABELS = [
116
  "21:40:00",
117
  "22:40:00",
118
  "24:10:00", # 00:10Z
119
- "24:20:00", # 00:20Z (final)
120
  ]
121
  GUIDED_FINAL_LABEL = "24:20:00" # 00:20Z FINAL RESTING POINT
122
 
123
- # Score ONLY these pings (and NEVER 00:20)
124
- GUIDED_SCORE_LABELS = [
125
- "18:40:00",
126
- "19:40:00",
127
- "20:40:00",
128
- "21:40:00",
129
- "22:40:00",
130
- "24:10:00", # 00:10Z
131
- ]
132
- GUIDED_FAST_SEED_LABELS = [
133
- "18:40:00",
134
- "19:40:00",
135
- "22:40:00",
136
- ]
137
-
138
  # Turn modeling parameters
139
  TURN_BANK_MIN_DEG = 15.0
140
  TURN_BANK_MAX_DEG = 30.0
141
  TURN_BANK_DEFAULT_DEG = 25.0
142
  TURN_STEP_SEC = 30.0
143
 
144
- # Turn speed bounds
145
- TURN_SPEED_MIN_ABS_KT = 470.0
146
- TURN_SPEED_MAX_ABS_KT = 500.0
147
-
148
- # Southern speed bounds (2 segments)
149
- SOUTH_SPEED_CENTER_GUESS_KT = 445.0
150
- SOUTH_SPEED_MIN_ABS_KT = 400.0
151
- SOUTH_SPEED_MAX_ABS_KT = 500.0
152
-
153
-
154
- # Southern schedule:
155
- # - scoring window ends at 00:10 (24:10)
156
- # - solve v1: 18:40->19:40 (60 min), v2: 19:40->00:10 (270 min)
157
- # - final 00:10->00:20 (10 min) uses v2 for display; 00:20 ignored in scoring
158
- SOUTH_SCHEDULE_SCORE_MIN = [60.0, 270.0]
159
- SOUTH_SCHEDULE_FULL_MIN = [60.0, 270.0, 10.0]
160
-
161
- # Regularization to avoid overfitting crazy speed behavior
162
- SPEED_SMOOTH_WEIGHT = 0.05
163
- SPEED_CENTER_WEIGHT = 0.01
164
 
165
  # Guided heading search (south-ish)
166
  GUIDED_HEADING_MIN = 150.0
167
  GUIDED_HEADING_MAX = 210.0
 
168
 
169
- # ===== Seed grid + pruning =====
170
- COARSE_HEADING_STEP_DEG = 2.0
171
- COARSE_TURN_SPEED_STEP_KT = 5.0
172
- SEED_PRUNE_TOPK = 8 # BIG speed-up: only optimize from top K seeds
 
 
 
 
 
173
 
174
- # Optimizer iterations
175
- OPT_MAXITER = 160 # lower than before (faster), usually enough with good seeds
 
 
 
176
 
177
- # Cache controls
178
- CACHE_MAX_ITEMS = 50000 # safety cap
 
 
179
 
180
- # Search area colors
181
  SEARCH_AREA_COLORS = {
182
  "Ocean Infinity": "#ff3b3b",
183
  "Australia": "#1f77b4",
@@ -475,9 +449,11 @@ def build_turn_legs(
475
 
476
  legs = []
477
 
 
478
  if straight_sec > 1.0:
479
  legs.append((h0, float(speed_knots), straight_sec / 60.0))
480
 
 
481
  if turn_sec > 1.0 and delta_abs > 1e-6:
482
  step_sec = max(5.0, float(step_sec))
483
  n = int(math.ceil(turn_sec / step_sec))
@@ -692,43 +668,12 @@ def _map_marker_times(start_dt, legs):
692
  return times
693
 
694
 
695
- def _make_state_evaluator(start_lat: float, start_lon: float, legs, method: str, start_dt: datetime):
696
- """
697
- Precompute boundaries + leg start positions once, then evaluate lat/lon/spd/hdg at dt quickly.
698
- """
699
- legs = [(float(h), float(s), float(m)) for (h, s, m) in legs]
700
- leg_secs = np.array([m * 60.0 for (_h, _s, m) in legs], dtype=float)
701
- boundaries = np.cumsum(np.concatenate([[0.0], leg_secs]))
702
- total_end = float(boundaries[-1])
703
-
704
- leg_starts = _precompute_leg_starts(start_lat, start_lon, legs, method)
705
-
706
- def state_at(dt_k: datetime):
707
- if dt_k < start_dt:
708
- return None
709
- t_sec = float((dt_k - start_dt).total_seconds())
710
- if t_sec > total_end:
711
- t_sec = total_end
712
-
713
- i = int(np.searchsorted(boundaries, t_sec, side="right") - 1)
714
- i = max(0, min(i, len(legs) - 1))
715
- t_inside = t_sec - float(boundaries[i])
716
-
717
- lat0, lon0 = leg_starts[i]
718
- hdg, spd, _m = legs[i]
719
- lat, lon = _predict_state_from_leg_start(method, lat0, lon0, hdg, spd, t_inside)
720
- return {
721
- "lat": float(lat),
722
- "lon": float(lon),
723
- "speed_kt": float(spd),
724
- "heading_deg": float(hdg),
725
- }
726
-
727
- return state_at
728
-
729
-
730
  # ================= CLIP TRAJECTORY TO 00:20 (GUIDED) =================
731
  def clip_legs_to_end_dt(legs, start_dt: datetime, end_dt: datetime):
 
 
 
 
732
  total_target_sec = (end_dt - start_dt).total_seconds()
733
  if total_target_sec <= 0:
734
  return []
@@ -744,6 +689,7 @@ def clip_legs_to_end_dt(legs, start_dt: datetime, end_dt: datetime):
744
  out.append((float(hdg), float(spd), float(minutes)))
745
  remaining -= sec
746
  else:
 
747
  trunc_min = max(0.0, remaining / 60.0)
748
  if trunc_min > 1e-6:
749
  out.append((float(hdg), float(spd), float(trunc_min)))
@@ -860,50 +806,12 @@ def bfo_analysis_with_legs(start_lat, start_lon, legs, method, start_dt, alt_ft,
860
 
861
 
862
  # ================= VALIDATION TABLE =================
863
- def _lookup_measured_at_dt(dt_ping: datetime, prefer_bto: bool = True, bto_search_radius: int = 4):
864
- """
865
- Returns (meas_bfo, meas_bto).
866
-
867
- IMPORTANT:
868
- Sometimes the nearest satLoc sample has BFO but BTO is None (your 18:40 case).
869
- If prefer_bto=True, we search nearby indices for the closest non-None BTO.
870
- """
871
- if not hasattr(fp, "nearest_measured_for_time"):
872
- return None, None
873
-
874
- idx0 = int(fp.nearest_measured_for_time(dt_ping)[0])
875
- row0 = fp.satLoc[idx0]
876
- bfo0 = row0[SATLOC_IDX_MEAS_BFO]
877
- bto0 = row0[SATLOC_IDX_MEAS_BTO]
878
-
879
- if (not prefer_bto) or (bto0 is not None):
880
- return bfo0, bto0
881
-
882
- # Search neighbors for a valid BTO
883
- best = (None, None, None) # (abs_dt_sec, bto, idx)
884
- for d in range(1, int(bto_search_radius) + 1):
885
- for idx in (idx0 - d, idx0 + d):
886
- if idx < 0 or idx >= len(fp.satLoc):
887
- continue
888
- row = fp.satLoc[idx]
889
- bto = row[SATLOC_IDX_MEAS_BTO]
890
- if bto is None:
891
- continue
892
- # Approx time distance using satLoc label time
893
- try:
894
- lab = fp._extract_time_string(row[0])
895
- dt_row = fp._time_label_to_dt(lab, fp.SAT_EPOCH_UTC)
896
- abs_sec = abs((dt_row - dt_ping).total_seconds())
897
- except Exception:
898
- abs_sec = float(d) * 60.0
899
-
900
- if best[0] is None or abs_sec < best[0]:
901
- best = (abs_sec, bto, idx)
902
-
903
- if best[1] is None:
904
- return bfo0, bto0
905
-
906
- return bfo0, best[1]
907
 
908
 
909
  def build_validation_table(
@@ -951,11 +859,7 @@ def build_validation_table(
951
  bfo_calc = float(fp.calculateBFO_at_time(spd_mps, float(hdg), float(lat), float(lon), alt_m, vs_mps, dt_k))
952
  bto_calc = float(fp.calculateBTO_at_time(float(lat), float(lon), alt_m, dt_k))
953
 
954
- bfo_meas, bto_meas = _lookup_measured_at_dt(dt_k, prefer_bto=True)
955
-
956
- # remove measured values at 00:20 (final)
957
- if is_guided and dt_k == _dt_from_label(GUIDED_FINAL_LABEL):
958
- bfo_meas, bto_meas = None, None
959
 
960
  bfo_err = None if bfo_meas is None else (float(bfo_calc) - float(bfo_meas))
961
  bto_err = None if bto_meas is None else (float(bto_calc) - float(bto_meas))
@@ -1022,6 +926,13 @@ def add_permanent_label(m, name, lat, lon):
1022
 
1023
 
1024
  def build_leaflet_map(path_pts, legs, method, start_dt):
 
 
 
 
 
 
 
1025
  start_lat, start_lon = path_pts[0]
1026
  is_guided = (start_dt == RADAR_FIX_TIME_UTC)
1027
 
@@ -1166,6 +1077,7 @@ def build_leaflet_map(path_pts, legs, method, start_dt):
1166
  f"Distance flown: {dist_km:.2f} km",
1167
  ]
1168
 
 
1169
  if is_final:
1170
  folium.CircleMarker(
1171
  location=(lat, lon),
@@ -1220,319 +1132,404 @@ def build_leaflet_map(path_pts, legs, method, start_dt):
1220
  """
1221
 
1222
 
1223
- # ================= Guided scoring point selection =================
1224
- def _score_points_for_labels(labels: list[str]):
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1225
  """
1226
- Return [(label, dt, meas_bfo, meas_bto)] using nearest measured sample for each label time.
1227
- Ensures unique rows (if two labels snap to same satLoc index).
1228
  """
1229
- out = []
1230
- used_idx = set()
1231
- for lab in labels:
1232
- dt = _dt_from_label(lab)
1233
- if hasattr(fp, "nearest_measured_for_time"):
1234
- idx = fp.nearest_measured_for_time(dt)[0]
1235
- if idx in used_idx:
1236
- continue
1237
- used_idx.add(idx)
1238
- row = fp.satLoc[idx]
1239
- bfo = row[SATLOC_IDX_MEAS_BFO]
1240
- bto = row[SATLOC_IDX_MEAS_BTO]
1241
- if bfo is None:
1242
- continue
1243
- out.append((lab, dt, float(bfo), None if bto is None else float(bto)))
1244
- else:
1245
- # fallback: scan satLoc text label match
1246
- best = None
1247
- for r in fp.satLoc:
1248
- label_r = fp._extract_time_string(r[0])
1249
- if label_r.replace("Z", "") == lab.replace("Z", ""):
1250
- best = r
1251
- break
1252
- if best is None or best[SATLOC_IDX_MEAS_BFO] is None:
1253
- continue
1254
- out.append((lab, dt, float(best[SATLOC_IDX_MEAS_BFO]),
1255
- None if best[SATLOC_IDX_MEAS_BTO] is None else float(best[SATLOC_IDX_MEAS_BTO])))
1256
- out.sort(key=lambda x: x[1])
1257
  return out
1258
 
1259
 
1260
- # ================= Guided Solver (OPTIMIZER MIRROR + FAST) =================
1261
- def solve_guided_straight(
1262
- method, alt_ft, vs_fpm,
1263
- descent_start_choice, weight_bfo, weight_bto
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1264
  ):
1265
- start_dt = RADAR_FIX_TIME_UTC
1266
- end_dt = _dt_from_label(GUIDED_FINAL_LABEL)
 
 
 
 
1267
 
1268
- descent_start_dt = None
1269
- if (descent_start_choice or "").startswith("After Ping Ring 6"):
1270
- descent_start_dt = dt_ping6_utc()
1271
- elif (descent_start_choice or "").startswith("After Ping Ring 7"):
1272
- descent_start_dt = dt_ping7_utc()
1273
 
1274
- # Strict scoring pings (6 points)
1275
- pts_score = _score_points_for_labels(GUIDED_SCORE_LABELS)
1276
- if len(pts_score) < 2:
1277
- return None, "Guided solve: not enough scoring pings (check satLoc).", None
1278
-
1279
- pts_fast = _score_points_for_labels(GUIDED_FAST_SEED_LABELS)
1280
- if len(pts_fast) < 2:
1281
- pts_fast = pts_score[:3]
1282
-
1283
- # --- Cache for expensive BFO/BTO computations ---
1284
- cache = {}
1285
-
1286
- def _cache_get(key):
1287
- return cache.get(key, None)
1288
-
1289
- def _cache_put(key, val):
1290
- if len(cache) >= int(CACHE_MAX_ITEMS):
1291
- # simple eviction: clear (safe + fast)
1292
- cache.clear()
1293
- cache[key] = val
1294
-
1295
- def _compute_residuals_at_point_cached(lat, lon, hdg_deg, spd_knots, dt_ping):
1296
- """
1297
- Returns (r_bfo, r_bto) where residual = calc - meas.
1298
- Uses memoization with rounding to speed repeated optimizer calls.
1299
- """
1300
- # altitude & VS depend on time
1301
- alt_now_ft = altitude_ft_at_dt(dt_ping, RADAR_FIX_TIME_UTC, float(alt_ft), float(vs_fpm), descent_start_dt)
1302
  alt_m = alt_now_ft * 0.3048
1303
 
1304
  vs_mps = float(vs_fpm) * 0.00508
1305
  if descent_start_dt is not None and dt_ping < descent_start_dt:
1306
  vs_mps = 0.0
1307
 
1308
- spd_mps = float(spd_knots) * 1852.0 / 3600.0
 
 
 
1309
 
1310
- key = (
1311
- round(float(lat), 4),
1312
- round(float(lon), 4),
1313
- round(float(hdg_deg), 1),
1314
- round(float(spd_knots), 1),
1315
- int(dt_ping.timestamp()),
1316
- round(float(alt_m), 1),
1317
- round(float(vs_mps), 3),
1318
- )
1319
- got = _cache_get(key)
1320
- if got is not None:
1321
- return got
1322
 
1323
- bfo_calc = float(fp.calculateBFO_at_time(spd_mps, float(hdg_deg), float(lat), float(lon), alt_m, vs_mps, dt_ping))
1324
- bto_calc = float(fp.calculateBTO_at_time(float(lat), float(lon), alt_m, dt_ping))
1325
- bfo_meas, bto_meas = _lookup_measured_at_dt(dt_ping, prefer_bto=True)
1326
 
1327
- r_bfo = None if bfo_meas is None else (bfo_calc - float(bfo_meas))
1328
- r_bto = None if bto_meas is None else (bto_calc - float(bto_meas))
 
 
1329
 
1330
- _cache_put(key, (r_bfo, r_bto))
1331
- return r_bfo, r_bto
 
1332
 
1333
- def _regularization_speed(v1: float, v2: float) -> float:
1334
- pen = 0.0
1335
- pen += float(SPEED_CENTER_WEIGHT) * ((float(v1) - float(SOUTH_SPEED_CENTER_GUESS_KT)) / 10.0) ** 2
1336
- pen += float(SPEED_CENTER_WEIGHT) * ((float(v2) - float(SOUTH_SPEED_CENTER_GUESS_KT)) / 10.0) ** 2
1337
- pen += float(SPEED_SMOOTH_WEIGHT) * ((float(v2) - float(v1)) / 10.0) ** 2
1338
- return float(pen)
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1339
 
1340
- def _build_legs_for_params(h_south, spd_turn, v1, v2, full=False):
1341
- turn_legs, _ = build_turn_legs(
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1342
  speed_knots=float(spd_turn),
1343
  method=method,
1344
  heading_final_deg=float(h_south),
1345
  )
1346
- if not turn_legs:
1347
- return None
1348
- if full:
1349
- south = [
1350
- (float(h_south), float(v1), float(SOUTH_SCHEDULE_FULL_MIN[0])),
1351
- (float(h_south), float(v2), float(SOUTH_SCHEDULE_FULL_MIN[1])),
1352
- (float(h_south), float(v2), float(SOUTH_SCHEDULE_FULL_MIN[2])),
1353
- ]
1354
- else:
1355
- south = [
1356
- (float(h_south), float(v1), float(SOUTH_SCHEDULE_SCORE_MIN[0])),
1357
- (float(h_south), float(v2), float(SOUTH_SCHEDULE_SCORE_MIN[1])),
1358
- ]
1359
- return turn_legs + south
1360
-
1361
- # Objective on a chosen ping set
1362
- def _objective_on_pts(x, pts):
1363
- h_south = float(x[0])
1364
- spd_turn = float(x[1])
1365
- v1 = float(x[2])
1366
- v2 = float(x[3])
1367
-
1368
- legs_score = _build_legs_for_params(h_south, spd_turn, v1, v2, full=False)
1369
- if legs_score is None:
1370
- return 1e12
1371
-
1372
- state_at = _make_state_evaluator(RADAR_FIX_LAT, RADAR_FIX_LON, legs_score, method, start_dt)
1373
-
1374
- sse = 0.0
1375
- n_used = 0
1376
-
1377
- for (_lbl, dt_ping, _mbfo, _mbto) in pts:
1378
- st = state_at(dt_ping)
1379
- if st is None:
1380
- continue
1381
-
1382
- r_bfo, r_bto = _compute_residuals_at_point_cached(
1383
- st["lat"], st["lon"], st["heading_deg"], st["speed_kt"], dt_ping
1384
- )
1385
 
1386
- sse_bto = 0.0
1387
- sse_bfo = 0.0
1388
- n_bto = 0
1389
- n_bfo = 0
 
1390
 
1391
- for (_lbl, dt_ping, _mbfo, _mbto) in pts:
1392
- st = state_at(dt_ping)
1393
- if st is None:
1394
- continue
1395
 
1396
- r_bfo, r_bto = _compute_residuals_at_point_cached(
1397
- st["lat"], st["lon"], st["heading_deg"], st["speed_kt"], dt_ping
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1398
  )
1399
 
1400
- # PRIMARY: BTO
1401
- if r_bto is not None:
1402
- sse_bto += (float(r_bto) / float(BTO_SCALE_US)) ** 2
1403
- n_bto += 1
 
 
 
 
1404
 
1405
- # SECONDARY: BFO
1406
- if r_bfo is not None:
1407
- sse_bfo += (float(r_bfo) / float(BFO_SCALE_HZ)) ** 2
1408
- n_bfo += 1
1409
-
1410
- # If no BTO points, reject (we're optimizing for BTO)
1411
- if n_bto == 0:
1412
- return 1e12
 
 
1413
 
1414
- rmse_bto_n = math.sqrt(sse_bto / max(1, n_bto))
1415
- rmse_bfo_n = math.sqrt(sse_bfo / max(1, n_bfo)) if n_bfo > 0 else 0.0
1416
 
1417
- # Objective: BTO-first, BFO is a tie-breaker/regularizer
1418
- # weight_bto acts like "how strongly you want BTO", but BTO is always primary here.
1419
- obj = float(weight_bto) * rmse_bto_n + float(weight_bfo) * rmse_bfo_n
1420
 
1421
- # Keep your speed regularization (helps prevent crazy solutions)
1422
- obj += _regularization_speed(v1, v2)
 
 
 
 
 
1423
 
1424
- return float(obj)
 
 
 
 
 
1425
 
1426
- # Full objective (strict scoring pings)
1427
- def objective(x):
1428
- return _objective_on_pts(x, pts_score)
 
 
1429
 
1430
- # Bounds (mirrors bounded optimizer style)
1431
- bounds = [
1432
- (GUIDED_HEADING_MIN, GUIDED_HEADING_MAX), # heading
1433
- (TURN_SPEED_MIN_ABS_KT, TURN_SPEED_MAX_ABS_KT), # turn speed
1434
- (SOUTH_SPEED_MIN_ABS_KT, SOUTH_SPEED_MAX_ABS_KT), # v1
1435
- (SOUTH_SPEED_MIN_ABS_KT, SOUTH_SPEED_MAX_ABS_KT), # v2
1436
- ]
1437
 
1438
- # ===== Seed pruning (FAST) =====
1439
- coarse_headings = np.arange(GUIDED_HEADING_MIN, GUIDED_HEADING_MAX + 1e-9, COARSE_HEADING_STEP_DEG)
1440
- coarse_turn_speeds = np.arange(TURN_SPEED_MIN_ABS_KT, TURN_SPEED_MAX_ABS_KT + 1e-9, COARSE_TURN_SPEED_STEP_KT)
1441
 
1442
- seed_scores = []
1443
- for h in coarse_headings:
1444
- for spd in coarse_turn_speeds:
1445
- x0 = np.array([float(h), float(spd), float(SOUTH_SPEED_CENTER_GUESS_KT), float(SOUTH_SPEED_CENTER_GUESS_KT)], dtype=float)
1446
- f0 = _objective_on_pts(x0, pts_fast)
1447
- seed_scores.append((float(f0), float(h), float(spd)))
1448
-
1449
- seed_scores.sort(key=lambda t: t[0])
1450
- top = seed_scores[:max(1, int(SEED_PRUNE_TOPK))]
1451
-
1452
- # ===== Optimize only top K seeds =====
1453
- best_res = None
1454
- for (_f0, h0, spd0) in top:
1455
- x0 = np.array([float(h0), float(spd0), float(SOUTH_SPEED_CENTER_GUESS_KT), float(SOUTH_SPEED_CENTER_GUESS_KT)], dtype=float)
1456
- res = minimize(
1457
- objective,
1458
- x0=x0,
1459
- method="L-BFGS-B",
1460
- bounds=bounds,
1461
- options={"maxiter": int(OPT_MAXITER), "ftol": 1e-9},
1462
  )
1463
- if not np.isfinite(res.fun):
1464
- continue
1465
- if (best_res is None) or (float(res.fun) < float(best_res.fun)):
1466
- best_res = res
1467
-
1468
- if best_res is None:
1469
- return None, "Guided solve: no valid solution found (optimizer failed).", None
1470
-
1471
- h_best = float(best_res.x[0])
1472
- spd_turn_best = float(best_res.x[1])
1473
- v1_best = float(best_res.x[2])
1474
- v2_best = float(best_res.x[3])
1475
-
1476
- # Build final legs for display through 00:20
1477
- full_legs = _build_legs_for_params(h_best, spd_turn_best, v1_best, v2_best, full=True)
1478
- full_legs = clip_legs_to_end_dt(full_legs, start_dt, end_dt)
1479
-
1480
- # Compute RMSEs for reporting (strict scoring pings)
1481
- def final_metrics():
1482
- legs_score = _build_legs_for_params(h_best, spd_turn_best, v1_best, v2_best, full=False)
1483
- state_at = _make_state_evaluator(RADAR_FIX_LAT, RADAR_FIX_LON, legs_score, method, start_dt)
1484
-
1485
- sum_bfo2 = 0.0
1486
- sum_bto2 = 0.0
1487
- n_bfo = 0
1488
- n_bto = 0
1489
-
1490
- for (_lbl, dt_ping, _mbfo, _mbto) in pts_score:
1491
- st = state_at(dt_ping)
1492
- if st is None:
1493
- continue
1494
- r_bfo, r_bto = _compute_residuals_at_point_cached(
1495
- st["lat"], st["lon"], st["heading_deg"], st["speed_kt"], dt_ping
1496
- )
1497
- if r_bfo is not None:
1498
- sum_bfo2 += float(r_bfo) ** 2
1499
- n_bfo += 1
1500
- if r_bto is not None:
1501
- sum_bto2 += float(r_bto) ** 2
1502
- n_bto += 1
1503
-
1504
- rmse_bfo = math.sqrt(sum_bfo2 / max(1, n_bfo)) if n_bfo > 0 else 0.0
1505
- rmse_bto = math.sqrt(sum_bto2 / max(1, n_bto)) if n_bto > 0 else None
1506
- return n_bfo, n_bto, rmse_bfo, rmse_bto
1507
-
1508
- n_bfo, n_bto, rmse_bfo, rmse_bto = final_metrics()
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1509
 
1510
  report = "\n".join([
1511
- "**Guided mode: Straight + computed turn (18:22→18:40), clipped to 00:20**",
1512
- "- Solver: **L-BFGS-B** on **[heading, turn_speed, v1, v2]**",
1513
- f"- Seed pruning: coarse grid → **top {SEED_PRUNE_TOPK}** seeds (fast pre-score on {len(pts_fast)} pings)",
1514
- f"- Strict scoring pings: **{', '.join(GUIDED_SCORE_LABELS)}**",
1515
- f"- Best southern heading: **{h_best:.2f}°**",
1516
- f"- Best turn speed (18:22→18:40): **{spd_turn_best:.2f} kt**",
1517
- f"- Southern speed segments:",
1518
- f" - 18:40 → 19:40: **{v1_best:.2f} kt**",
1519
- f" - 19:40 → 00:10: **{v2_best:.2f} kt**",
1520
- f" - 00:10 → 00:20 (display only): **{v2_best:.2f} kt**",
1521
- f"- Objective (normalized): **{float(best_res.fun):.6f}**",
1522
- f"- Optimizer status: **{best_res.message}**",
1523
- f"- Function evals: **{getattr(best_res, 'nfev', 'n/a')}** | Iters: **{getattr(best_res, 'nit', 'n/a')}**",
1524
- f"- Cache size (end of solve): **{len(cache)}** items (cap {CACHE_MAX_ITEMS})",
1525
- f"- Regularization: smooth={SPEED_SMOOTH_WEIGHT:.2f}, center={SPEED_CENTER_WEIGHT:.2f}",
1526
- f"- BFO points: **{n_bfo}** | RMSE_BFO: **{rmse_bfo:.2f} Hz**",
1527
- f"- BTO points: **{n_bto}** | RMSE_BTO: **{rmse_bto:.1f} µs**" if rmse_bto is not None else
1528
- f"- BTO points: **{n_bto}** | RMSE_BTO: n/a",
1529
  ])
1530
 
1531
- return full_legs, report, start_dt
1532
 
1533
 
1534
  # ================= RUNNERS =================
1535
  def _run_common(start_lat, start_lon, legs, method, start_dt, alt_ft, vs_fpm, descent_start_choice):
 
1536
  if start_dt == RADAR_FIX_TIME_UTC:
1537
  end_dt = _dt_from_label(GUIDED_FINAL_LABEL)
1538
  legs = clip_legs_to_end_dt(legs, start_dt, end_dt)
@@ -1542,6 +1539,7 @@ def _run_common(start_lat, start_lon, legs, method, start_dt, alt_ft, vs_fpm, de
1542
  html_map = build_leaflet_map(path_pts, legs, method, start_dt)
1543
  stats_md = compute_stats(path_pts, seg_dists_m, total_time_min, legs)
1544
 
 
1545
  bfo_fig, bfo_md = bfo_analysis_with_legs(
1546
  start_lat, start_lon, legs, method, start_dt, alt_ft, vs_fpm, descent_start_choice
1547
  )
@@ -1572,13 +1570,26 @@ def run_sim_manual(start_lat, start_lon, legs_df, method, start_time_str, alt_ft
1572
 
1573
  def run_guided(
1574
  start_lat, start_lon, method, alt_ft, vs_fpm, descent_start_choice,
1575
- weight_bfo, weight_bto
 
 
1576
  ):
1577
- legs, report_md, start_dt = solve_guided_straight(
1578
- method,
1579
- alt_ft, vs_fpm,
1580
- descent_start_choice, weight_bfo, weight_bto
1581
- )
 
 
 
 
 
 
 
 
 
 
 
1582
 
1583
  if legs is None:
1584
  fig = go.Figure()
@@ -1592,6 +1603,7 @@ def run_guided(
1592
  ], value=[])
1593
  return gr.update(), f"### ❌ {report_md}", "", "", fig, "", alt_fig, empty_table
1594
 
 
1595
  html_map, stats_md, bfo_fig, bfo_md, alt_fig, headers, rows = _run_common(
1596
  RADAR_FIX_LAT, RADAR_FIX_LON, legs, method, start_dt, alt_ft, vs_fpm, descent_start_choice
1597
  )
@@ -1623,10 +1635,8 @@ with gr.Blocks() as demo:
1623
  f"- **Turn completion (fixed):** 18:40:00Z\n"
1624
  f"- **Initial track (fixed):** {RADAR_INITIAL_TRACK_DEG:.1f}°\n"
1625
  f"- **Final resting point:** **00:20:00Z** (trajectory clipped exactly at this time)\n\n"
1626
- f"**Turn speed bounds:** **{TURN_SPEED_MIN_ABS_KT:.0f}–{TURN_SPEED_MAX_ABS_KT:.0f} kt**\n"
1627
- f"**Southern speed bounds:** **{SOUTH_SPEED_MIN_ABS_KT:.0f}–{SOUTH_SPEED_MAX_ABS_KT:.0f} kt** (center {SOUTH_SPEED_CENTER_GUESS_KT:.0f})\n"
1628
- f"**Scoring pings:** **{', '.join(GUIDED_SCORE_LABELS)}** (00:20 excluded)\n"
1629
- f"**Performance:** seed_prune_topk={SEED_PRUNE_TOPK}, opt_maxiter={OPT_MAXITER}, cache_cap={CACHE_MAX_ITEMS}"
1630
  )
1631
 
1632
  with gr.Row():
@@ -1645,8 +1655,28 @@ with gr.Blocks() as demo:
1645
  label="When does descent start? (VS applies after this time)",
1646
  )
1647
 
1648
- gr.Markdown("### Guided mode (Straight only)")
1649
- gr.Markdown("- Optimizer-based (L-BFGS-B) with fast seed pruning + caching.")
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1650
 
1651
  gr.Markdown("### Guided scoring weights")
1652
  with gr.Row():
@@ -1654,7 +1684,7 @@ with gr.Blocks() as demo:
1654
  weight_bto = gr.Slider(0.0, 10.0, value=DEFAULT_WEIGHT_BTO, step=0.1, label="Weight: BTO (increase to emphasize BTO)")
1655
 
1656
  with gr.Row():
1657
- guided_btn = gr.Button("Generate: Guided (Straight)")
1658
  manual_btn = gr.Button("Simulate: Manual legs")
1659
 
1660
  status_md = gr.Markdown()
@@ -1701,7 +1731,9 @@ with gr.Blocks() as demo:
1701
  run_guided,
1702
  inputs=[
1703
  start_lat, start_lon, method, alt_ft, vs_fpm, descent_start_choice,
1704
- weight_bfo, weight_bto
 
 
1705
  ],
1706
  outputs=[legs_df, status_md, out_map, stats, bfo_plot, bfo_score, alt_plot, validation_df],
1707
  )
 
1
  # app.py
2
  # MH370 Flight Simulator — Leaflet map (Folium) + Ping Rings (KML)
3
  #
4
+ # Guided scenarios (ACTIVE):
5
  # A) Straight (single heading)
6
+ # A2) Straight + drift (per-leg ±deg, penalized)
7
  #
8
  # Key behavior:
9
  # - Last radar spot FIXED: 18:22:00Z @ 6.578N, 96.341E
 
11
  # - Guided mode computes a realistic turn between 18:22→18:40 and then southern path.
12
  # - Guided validation timestamps: 18:22, 18:40, 19:40, 20:40, 21:40, 22:40, 00:10, 00:20
13
  #
14
+ # FIXES IMPLEMENTED IN THIS VERSION:
15
  # 1) Remove trajectory (red line) AFTER 00:20
16
  # - In guided mode, we CLIP the simulated legs so the polyline ends exactly at 00:20Z.
17
+ # 2) Restore BFO plot (it disappeared previously)
18
+ # - Re-added full BFO analysis computation + chart + MAE/RMSE markdown.
19
  #
20
+ # NEW (PREVIOUS UPDATE):
21
  # 3) Plausible final major turn (18:22→18:40) can adjust BOTH heading and speed
22
+ # - Turn speed is no longer fixed to the southern cruise speed.
23
+ # - The solver searches a plausible TURN speed range and chooses the best-scoring solution.
24
+ # - After 18:40, the southern trajectory continues “as per usual” using the selected speed profile.
25
+ # - Applies to BOTH scenarios: Straight and Straight + drift.
26
  #
27
+ # NEW (THIS UPDATE):
28
+ # 4) Two-stage guided search (COARSE → REFINE) to make generation much faster
29
+ # - Stage 1: coarse grid search to find the region of the best solution
30
+ # - Stage 2: refined grid around the best candidate only
31
+ # - Applies to BOTH scenarios: Straight and Straight + drift
 
 
 
 
 
 
 
 
 
 
 
 
32
 
33
  import os
34
  import re
 
45
 
46
  import flight_plotter_v32 as fp
47
 
 
 
 
48
 
49
  # ================= CONFIG =================
50
  KML_FILE = "MH370 - 13 nov (3).kml"
 
76
  "Ping Ring 7": "24:10:00Z",
77
  }
78
 
79
+ # Normalization (used for combining BFO + BTO into a single score in guided solve)
80
  BFO_SCALE_HZ = 5.0
81
+ BTO_SCALE_US = 50.0
82
 
83
  # EXPLICIT WEIGHTS
84
  DEFAULT_WEIGHT_BFO = 1.0
 
92
  RADAR_FIX_LAT = 6.578
93
  RADAR_FIX_LON = 96.341
94
  RADAR_FIX_TIME_UTC = datetime(2014, 3, 7, 18, 22, 0, tzinfo=timezone.utc)
95
+
96
  TURN_COMPLETE_TIME_UTC = datetime(2014, 3, 7, 18, 40, 0, tzinfo=timezone.utc)
97
+ TURN_WINDOW_MIN = (TURN_COMPLETE_TIME_UTC - RADAR_FIX_TIME_UTC).total_seconds() / 60.0
98
 
99
  # Initial track shown in your radar annotation
100
  RADAR_INITIAL_TRACK_DEG = 296.0
101
 
102
  # ===== Guided validation timestamps of interest =====
103
+ # (00:10 and 00:20 are represented as 24:10 and 24:20 in your SAT epoch helper)
104
  GUIDED_VALIDATION_LABELS = [
105
  "18:22:00",
106
  "18:40:00",
 
109
  "21:40:00",
110
  "22:40:00",
111
  "24:10:00", # 00:10Z
112
+ "24:20:00", # 00:20Z
113
  ]
114
  GUIDED_FINAL_LABEL = "24:20:00" # 00:20Z FINAL RESTING POINT
115
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
116
  # Turn modeling parameters
117
  TURN_BANK_MIN_DEG = 15.0
118
  TURN_BANK_MAX_DEG = 30.0
119
  TURN_BANK_DEFAULT_DEG = 25.0
120
  TURN_STEP_SEC = 30.0
121
 
122
+ # Turn-speed search parameters (turn only; southern uses speed profile)
123
+ TURN_SPEED_SEARCH_DELTA_KT = 40.0 # search ± around the southern cruise speed
124
+ TURN_SPEED_STEP_KT = 5.0 # refined candidate spacing
125
+ TURN_SPEED_MIN_ABS_KT = 350.0 # hard min plausible
126
+ TURN_SPEED_MAX_ABS_KT = 550.0 # hard max plausible
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
127
 
128
  # Guided heading search (south-ish)
129
  GUIDED_HEADING_MIN = 150.0
130
  GUIDED_HEADING_MAX = 210.0
131
+ STRAIGHT_STEP_DEG = 0.5 # refined step
132
 
133
+ # Southern schedule (kept)
134
+ GUIDED_LEG_SCHEDULE_MIN = [12.0, 60.0, 60.0, 60.0, 60.0, 90.0, 10.0]
135
+
136
+ # Drift
137
+ DEFAULT_DRIFT_MAX_DEG = 2.0
138
+ DRIFT_STEP_DEG = 0.5 # refined step
139
+ DEFAULT_WEIGHT_DRIFT = 0.30
140
+ DEFAULT_DRIFT_ITERS = 2
141
+ DRIFT_SCALE_DEG = 1.0
142
 
143
+ # ===== NEW: Two-stage search knobs =====
144
+ # Stage 1 (coarse)
145
+ COARSE_HEADING_STEP_DEG = 2.0
146
+ COARSE_TURN_SPEED_STEP_KT = 10.0
147
+ COARSE_DRIFT_STEP_DEG = 1.0
148
 
149
+ # Stage 2 (refine around best)
150
+ REFINE_HEADING_HALF_WINDOW_DEG = 2.0 # +/- around best
151
+ REFINE_TURN_SPEED_HALF_WINDOW_KT = 10.0
152
+ REFINE_DRIFT_HALF_WINDOW_DEG = 1.0
153
 
154
+ # Search area colors (optional)
155
  SEARCH_AREA_COLORS = {
156
  "Ocean Infinity": "#ff3b3b",
157
  "Australia": "#1f77b4",
 
449
 
450
  legs = []
451
 
452
+ # initial straight on radar track
453
  if straight_sec > 1.0:
454
  legs.append((h0, float(speed_knots), straight_sec / 60.0))
455
 
456
+ # discretized turn segments
457
  if turn_sec > 1.0 and delta_abs > 1e-6:
458
  step_sec = max(5.0, float(step_sec))
459
  n = int(math.ceil(turn_sec / step_sec))
 
668
  return times
669
 
670
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
671
  # ================= CLIP TRAJECTORY TO 00:20 (GUIDED) =================
672
  def clip_legs_to_end_dt(legs, start_dt: datetime, end_dt: datetime):
673
+ """
674
+ Returns a new legs list truncated so total duration ends exactly at end_dt.
675
+ If legs already shorter than end_dt-start_dt, returns legs unchanged.
676
+ """
677
  total_target_sec = (end_dt - start_dt).total_seconds()
678
  if total_target_sec <= 0:
679
  return []
 
689
  out.append((float(hdg), float(spd), float(minutes)))
690
  remaining -= sec
691
  else:
692
+ # truncate this leg
693
  trunc_min = max(0.0, remaining / 60.0)
694
  if trunc_min > 1e-6:
695
  out.append((float(hdg), float(spd), float(trunc_min)))
 
806
 
807
 
808
  # ================= VALIDATION TABLE =================
809
+ def _lookup_measured_at_dt(dt_ping: datetime):
810
+ if hasattr(fp, "nearest_measured_for_time"):
811
+ idx = fp.nearest_measured_for_time(dt_ping)[0]
812
+ row = fp.satLoc[idx]
813
+ return row[SATLOC_IDX_MEAS_BFO], row[SATLOC_IDX_MEAS_BTO]
814
+ return None, None
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
815
 
816
 
817
  def build_validation_table(
 
859
  bfo_calc = float(fp.calculateBFO_at_time(spd_mps, float(hdg), float(lat), float(lon), alt_m, vs_mps, dt_k))
860
  bto_calc = float(fp.calculateBTO_at_time(float(lat), float(lon), alt_m, dt_k))
861
 
862
+ bfo_meas, bto_meas = _lookup_measured_at_dt(dt_k)
 
 
 
 
863
 
864
  bfo_err = None if bfo_meas is None else (float(bfo_calc) - float(bfo_meas))
865
  bto_err = None if bto_meas is None else (float(bto_calc) - float(bto_meas))
 
926
 
927
 
928
  def build_leaflet_map(path_pts, legs, method, start_dt):
929
+ """
930
+ - Ping rings always drawn
931
+ - Search areas always drawn (color-coded)
932
+ - Validation markers are the authoritative markers
933
+ - Guided mode does NOT draw a purple End marker
934
+ - Trajectory is clipped to 00:20 in guided mode by clipping legs before simulate_path
935
+ """
936
  start_lat, start_lon = path_pts[0]
937
  is_guided = (start_dt == RADAR_FIX_TIME_UTC)
938
 
 
1077
  f"Distance flown: {dist_km:.2f} km",
1078
  ]
1079
 
1080
+ # Make final point visually distinct
1081
  if is_final:
1082
  folium.CircleMarker(
1083
  location=(lat, lon),
 
1132
  """
1133
 
1134
 
1135
+ # ================= Guided Speed Profile =================
1136
+ def isa_speed_of_sound_mps(alt_m: float) -> float:
1137
+ T0 = 288.15
1138
+ L = -0.0065
1139
+ h_trop = 11000.0
1140
+ if alt_m <= h_trop:
1141
+ T = T0 + L * alt_m
1142
+ else:
1143
+ T = 216.65
1144
+ gamma = 1.4
1145
+ R = 287.05
1146
+ return math.sqrt(gamma * R * T)
1147
+
1148
+
1149
+ def mach_to_tas_knots(mach: float, alt_ft: float) -> float:
1150
+ alt_m = float(alt_ft) * 0.3048
1151
+ a = isa_speed_of_sound_mps(alt_m)
1152
+ tas_mps = float(mach) * a
1153
+ return tas_mps / 0.514444
1154
+
1155
+
1156
+ def speed_from_profile(profile: str, alt_ft: float, custom_mach: float, custom_knots: float) -> float:
1157
+ profile = (profile or "").strip()
1158
+ if profile.startswith("MRC"):
1159
+ return mach_to_tas_knots(0.82, alt_ft)
1160
+ if profile.startswith("LRC"):
1161
+ return mach_to_tas_knots(0.84, alt_ft)
1162
+ if profile.startswith("ECON"):
1163
+ return mach_to_tas_knots(0.80, alt_ft)
1164
+ if profile.startswith("Custom Mach"):
1165
+ return mach_to_tas_knots(float(custom_mach), alt_ft)
1166
+ if profile.startswith("Custom knots"):
1167
+ return float(custom_knots)
1168
+ return mach_to_tas_knots(0.82, alt_ft)
1169
+
1170
+
1171
+ def _turn_speed_candidates(south_speed_knots: float, step_kt: float):
1172
  """
1173
+ Plausible turn-speed candidates around the southern cruise speed.
1174
+ Hard-clipped to [TURN_SPEED_MIN_ABS_KT, TURN_SPEED_MAX_ABS_KT].
1175
  """
1176
+ base = float(south_speed_knots)
1177
+ lo = max(TURN_SPEED_MIN_ABS_KT, base - float(TURN_SPEED_SEARCH_DELTA_KT))
1178
+ hi = min(TURN_SPEED_MAX_ABS_KT, base + float(TURN_SPEED_SEARCH_DELTA_KT))
1179
+ step = max(1.0, float(step_kt))
1180
+ if hi < lo:
1181
+ lo, hi = hi, lo
1182
+ cand = list(np.arange(lo, hi + 1e-9, step))
1183
+ # Ensure base is included
1184
+ if all(abs(c - base) > 1e-9 for c in cand):
1185
+ cand.append(base)
1186
+ cand = sorted({round(float(c), 6) for c in cand})
1187
+ return cand
1188
+
1189
+
1190
+ def _clamp(x, lo, hi):
1191
+ return max(lo, min(hi, x))
1192
+
1193
+
1194
+ def _frange(lo: float, hi: float, step: float):
1195
+ # inclusive hi, stable float stepping
1196
+ step = float(step)
1197
+ if step <= 0:
1198
+ return []
1199
+ n = int(math.floor((hi - lo) / step + 1e-9)) + 1
1200
+ out = [lo + i * step for i in range(max(0, n))]
1201
+ # trim due to floating errors
1202
+ out = [x for x in out if x <= hi + 1e-9]
 
1203
  return out
1204
 
1205
 
1206
+ def _refine_range(center: float, half_window: float, step: float, lo: float, hi: float):
1207
+ lo2 = _clamp(center - half_window, lo, hi)
1208
+ hi2 = _clamp(center + half_window, lo, hi)
1209
+ return _frange(lo2, hi2, step)
1210
+
1211
+
1212
+ # ================= Guided Legs Creation =================
1213
+ def make_legs_from_schedule(headings, speed_knots, schedule_minutes):
1214
+ legs = []
1215
+ if isinstance(headings, (int, float)):
1216
+ headings = [float(headings)] * len(schedule_minutes)
1217
+ for hdg, m in zip(headings, schedule_minutes):
1218
+ m = float(m)
1219
+ if m <= 0:
1220
+ continue
1221
+ legs.append((float(hdg), float(speed_knots), float(m)))
1222
+ return legs
1223
+
1224
+
1225
+ # ================= Scoring helper (used in guided solvers) =================
1226
+ def score_legs_bfo_bto(
1227
+ start_lat, start_lon, method, start_dt, end_dt, legs,
1228
+ alt_ft, vs_fpm, descent_start_dt,
1229
+ weight_bfo: float, weight_bto: float,
1230
+ pts_window=None
1231
  ):
1232
+ # pts_window: optional precomputed list [(lbl, dt, meas_bfo, meas_bto), ...] already clipped to [start_dt, end_dt]
1233
+ if pts_window is None:
1234
+ pts = _ping_points_after(start_dt)
1235
+ pts = [(lbl, dt, bfo, bto) for (lbl, dt, bfo, bto) in pts if start_dt <= dt <= end_dt]
1236
+ else:
1237
+ pts = pts_window
1238
 
1239
+ if len(pts) < 2:
1240
+ return None
 
 
 
1241
 
1242
+ leg_secs = [float(m) * 60.0 for _, _, m in legs]
1243
+ boundaries = np.cumsum([0.0] + leg_secs)
1244
+ total_end = float(boundaries[-1])
1245
+
1246
+ leg_start_positions = _precompute_leg_starts(start_lat, start_lon, legs, method)
1247
+
1248
+ bfo_preds, bfo_meas = [], []
1249
+ bto_preds, bto_meas = [], []
1250
+
1251
+ for _lbl, dt_ping, meas_bfo, meas_bto in pts:
1252
+ t_sec = (dt_ping - start_dt).total_seconds()
1253
+ if t_sec < 0 or t_sec > total_end:
1254
+ continue
1255
+
1256
+ i = int(np.searchsorted(boundaries, t_sec, side="right") - 1)
1257
+ i = max(0, min(i, len(legs) - 1))
1258
+ t_inside = float(t_sec) - float(boundaries[i])
1259
+
1260
+ leg_start_lat, leg_start_lon = leg_start_positions[i]
1261
+ hdg, spd, _ = legs[i]
1262
+ lat, lon = _predict_state_from_leg_start(method, leg_start_lat, leg_start_lon, hdg, spd, t_inside)
1263
+
1264
+ alt_now_ft = altitude_ft_at_dt(dt_ping, start_dt, float(alt_ft), float(vs_fpm), descent_start_dt)
 
 
 
 
 
1265
  alt_m = alt_now_ft * 0.3048
1266
 
1267
  vs_mps = float(vs_fpm) * 0.00508
1268
  if descent_start_dt is not None and dt_ping < descent_start_dt:
1269
  vs_mps = 0.0
1270
 
1271
+ spd_mps = float(spd) * 1852.0 / 3600.0
1272
+ pred_bfo = float(fp.calculateBFO_at_time(spd_mps, float(hdg), float(lat), float(lon), alt_m, vs_mps, dt_ping))
1273
+ bfo_preds.append(pred_bfo)
1274
+ bfo_meas.append(float(meas_bfo))
1275
 
1276
+ if meas_bto is not None:
1277
+ pred_bto = float(fp.calculateBTO_at_time(float(lat), float(lon), alt_m, dt_ping))
1278
+ bto_preds.append(pred_bto)
1279
+ bto_meas.append(float(meas_bto))
1280
+
1281
+ if len(bfo_meas) < 2:
1282
+ return None
 
 
 
 
 
1283
 
1284
+ e_bfo = np.array(bfo_preds) - np.array(bfo_meas)
1285
+ rmse_bfo = float(np.sqrt(np.mean(e_bfo * e_bfo)))
 
1286
 
1287
+ rmse_bto = None
1288
+ if len(bto_meas) >= 2:
1289
+ e_bto = np.array(bto_preds) - np.array(bto_meas)
1290
+ rmse_bto = float(np.sqrt(np.mean(e_bto * e_bto)))
1291
 
1292
+ combined = float(weight_bfo) * (rmse_bfo / BFO_SCALE_HZ)
1293
+ if rmse_bto is not None:
1294
+ combined += float(weight_bto) * (rmse_bto / BTO_SCALE_US)
1295
 
1296
+ return {
1297
+ "n_bfo": len(bfo_meas),
1298
+ "n_bto": len(bto_meas),
1299
+ "rmse_bfo": rmse_bfo,
1300
+ "rmse_bto": rmse_bto,
1301
+ "combined": combined,
1302
+ }
1303
+
1304
+
1305
+ # ================= Guided Solvers =================
1306
+ def _guided_start_dt():
1307
+ return RADAR_FIX_TIME_UTC
1308
+
1309
+
1310
+ def solve_guided_straight(
1311
+ method, south_speed_knots, alt_ft, vs_fpm,
1312
+ descent_start_choice, weight_bfo, weight_bto
1313
+ ):
1314
+ # We clip everything to final 00:20Z; end_dt is fixed by label.
1315
+ start_dt = _guided_start_dt()
1316
+ end_dt = _dt_from_label(GUIDED_FINAL_LABEL)
1317
 
1318
+ # descent start
1319
+ descent_start_dt = None
1320
+ if (descent_start_choice or "").startswith("After Ping Ring 6"):
1321
+ descent_start_dt = dt_ping6_utc()
1322
+ elif (descent_start_choice or "").startswith("After Ping Ring 7"):
1323
+ descent_start_dt = dt_ping7_utc()
1324
+
1325
+ # Precompute ping points once (windowed)
1326
+ pts_all = _ping_points_after(start_dt)
1327
+ pts_window = [(lbl, dt, bfo, bto) for (lbl, dt, bfo, bto) in pts_all if start_dt <= dt <= end_dt]
1328
+ if len(pts_window) < 2:
1329
+ return None, "Guided solve: not enough ping points to score.", None
1330
+
1331
+ def try_candidate(h_south: float, spd_turn: float):
1332
+ southern_legs = make_legs_from_schedule(h_south, float(south_speed_knots), GUIDED_LEG_SCHEDULE_MIN)
1333
+ preturn_legs, _turn_info = build_turn_legs(
1334
  speed_knots=float(spd_turn),
1335
  method=method,
1336
  heading_final_deg=float(h_south),
1337
  )
1338
+ full_legs = preturn_legs + southern_legs
1339
+ full_legs = clip_legs_to_end_dt(full_legs, start_dt, end_dt)
1340
+ s = score_legs_bfo_bto(
1341
+ RADAR_FIX_LAT, RADAR_FIX_LON, method, start_dt, end_dt, full_legs,
1342
+ alt_ft, vs_fpm, descent_start_dt,
1343
+ weight_bfo, weight_bto,
1344
+ pts_window=pts_window
1345
+ )
1346
+ return s, full_legs
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1347
 
1348
+ # -------- Stage 1: COARSE search --------
1349
+ best_h = None
1350
+ best_turn_spd = None
1351
+ best_full_legs = None
1352
+ best_score = None
1353
 
1354
+ coarse_headings = _frange(GUIDED_HEADING_MIN, GUIDED_HEADING_MAX, COARSE_HEADING_STEP_DEG)
1355
+ coarse_turn_speeds = _turn_speed_candidates(float(south_speed_knots), COARSE_TURN_SPEED_STEP_KT)
 
 
1356
 
1357
+ for h in coarse_headings:
1358
+ for spd_turn in coarse_turn_speeds:
1359
+ s, legs = try_candidate(h, spd_turn)
1360
+ if s is not None and (best_score is None or s["combined"] < best_score["combined"]):
1361
+ best_score = s
1362
+ best_h = float(h)
1363
+ best_turn_spd = float(spd_turn)
1364
+ best_full_legs = legs
1365
+
1366
+ if best_full_legs is None:
1367
+ return None, "Guided solve: not enough ping points to score.", None
1368
+
1369
+ # -------- Stage 2: REFINE around best --------
1370
+ refine_headings = _refine_range(best_h, REFINE_HEADING_HALF_WINDOW_DEG, STRAIGHT_STEP_DEG, GUIDED_HEADING_MIN, GUIDED_HEADING_MAX)
1371
+ refine_turn_speeds = _refine_range(
1372
+ best_turn_spd, REFINE_TURN_SPEED_HALF_WINDOW_KT, TURN_SPEED_STEP_KT,
1373
+ max(TURN_SPEED_MIN_ABS_KT, float(south_speed_knots) - TURN_SPEED_SEARCH_DELTA_KT),
1374
+ min(TURN_SPEED_MAX_ABS_KT, float(south_speed_knots) + TURN_SPEED_SEARCH_DELTA_KT),
1375
  )
1376
 
1377
+ for h in refine_headings:
1378
+ for spd_turn in refine_turn_speeds:
1379
+ s, legs = try_candidate(h, spd_turn)
1380
+ if s is not None and (best_score is None or s["combined"] < best_score["combined"]):
1381
+ best_score = s
1382
+ best_h = float(h)
1383
+ best_turn_spd = float(spd_turn)
1384
+ best_full_legs = legs
1385
 
1386
+ report = "\n".join([
1387
+ "**Guided mode: Straight + computed turn (18:22→18:40), clipped to 00:20**",
1388
+ f"- Search mode: **two-stage (coarse → refine)**",
1389
+ f"- Southern speed (speed profile): **{float(south_speed_knots):.1f} kt**",
1390
+ f"- Best turn speed (18:22→18:40): **{float(best_turn_spd):.1f} kt**",
1391
+ f"- Best southern heading: **{best_h:.1f}°**",
1392
+ f"- BFO points: **{best_score['n_bfo']}** | RMSE_BFO: **{best_score['rmse_bfo']:.2f} Hz**",
1393
+ f"- BTO points: **{best_score['n_bto']}** | RMSE_BTO: **{best_score['rmse_bto']:.1f} µs**" if best_score.get("rmse_bto") is not None else
1394
+ f"- BTO points: **{best_score['n_bto']}** | RMSE_BTO: n/a",
1395
+ ])
1396
 
1397
+ return best_full_legs, report, start_dt
 
1398
 
 
 
 
1399
 
1400
+ def solve_guided_straight_with_drift(
1401
+ method, south_speed_knots, alt_ft, vs_fpm,
1402
+ descent_start_choice, weight_bfo, weight_bto,
1403
+ drift_max_deg, drift_weight
1404
+ ):
1405
+ start_dt = _guided_start_dt()
1406
+ end_dt = _dt_from_label(GUIDED_FINAL_LABEL)
1407
 
1408
+ # descent start
1409
+ descent_start_dt = None
1410
+ if (descent_start_choice or "").startswith("After Ping Ring 6"):
1411
+ descent_start_dt = dt_ping6_utc()
1412
+ elif (descent_start_choice or "").startswith("After Ping Ring 7"):
1413
+ descent_start_dt = dt_ping7_utc()
1414
 
1415
+ # Precompute ping points once (windowed)
1416
+ pts_all = _ping_points_after(start_dt)
1417
+ pts_window = [(lbl, dt, bfo, bto) for (lbl, dt, bfo, bto) in pts_all if start_dt <= dt <= end_dt]
1418
+ if len(pts_window) < 2:
1419
+ return None, "Guided drift solve: not enough ping points to score.", None
1420
 
1421
+ drift_max_deg = float(drift_max_deg)
 
 
 
 
 
 
1422
 
1423
+ def try_candidate(h0_nominal: float, drift_deg: float, spd_turn: float):
1424
+ h_eff = _wrap_heading_deg(float(h0_nominal) + float(drift_deg))
 
1425
 
1426
+ southern_legs = make_legs_from_schedule(h_eff, float(south_speed_knots), GUIDED_LEG_SCHEDULE_MIN)
1427
+ preturn_legs, _turn_info = build_turn_legs(
1428
+ speed_knots=float(spd_turn),
1429
+ method=method,
1430
+ heading_final_deg=float(h_eff),
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1431
  )
1432
+ full_legs = preturn_legs + southern_legs
1433
+ full_legs = clip_legs_to_end_dt(full_legs, start_dt, end_dt)
1434
+
1435
+ s = score_legs_bfo_bto(
1436
+ RADAR_FIX_LAT, RADAR_FIX_LON, method, start_dt, end_dt, full_legs,
1437
+ alt_ft, vs_fpm, descent_start_dt,
1438
+ weight_bfo, weight_bto,
1439
+ pts_window=pts_window
1440
+ )
1441
+ if s is None:
1442
+ return None, None, None
1443
+
1444
+ combined_plus = float(s["combined"]) + float(drift_weight) * (abs(float(drift_deg)) / DRIFT_SCALE_DEG)
1445
+ return s, combined_plus, full_legs
1446
+
1447
+ # -------- Stage 1: COARSE search --------
1448
+ best_h0 = None
1449
+ best_d = None
1450
+ best_turn_spd = None
1451
+ best_full_legs = None
1452
+ best_score = None
1453
+ best_combined_plus = None
1454
+
1455
+ coarse_headings = _frange(GUIDED_HEADING_MIN, GUIDED_HEADING_MAX, COARSE_HEADING_STEP_DEG)
1456
+
1457
+ # Drift candidates (coarse)
1458
+ if drift_max_deg <= 1e-9:
1459
+ coarse_drifts = [0.0]
1460
+ else:
1461
+ coarse_drifts = list(np.arange(-drift_max_deg, drift_max_deg + 1e-9, COARSE_DRIFT_STEP_DEG))
1462
+ if 0.0 not in coarse_drifts:
1463
+ coarse_drifts.append(0.0)
1464
+ coarse_drifts = sorted({round(float(x), 6) for x in coarse_drifts})
1465
+
1466
+ coarse_turn_speeds = _turn_speed_candidates(float(south_speed_knots), COARSE_TURN_SPEED_STEP_KT)
1467
+
1468
+ for h0 in coarse_headings:
1469
+ for d in coarse_drifts:
1470
+ for spd_turn in coarse_turn_speeds:
1471
+ s, combined_plus, legs = try_candidate(h0, d, spd_turn)
1472
+ if s is None:
1473
+ continue
1474
+ if best_combined_plus is None or float(combined_plus) < float(best_combined_plus):
1475
+ best_score = dict(s)
1476
+ best_combined_plus = float(combined_plus)
1477
+ best_h0 = float(h0)
1478
+ best_d = float(d)
1479
+ best_turn_spd = float(spd_turn)
1480
+ best_full_legs = legs
1481
+
1482
+ if best_full_legs is None:
1483
+ return None, "Guided drift solve: not enough ping points to score.", None
1484
+
1485
+ # -------- Stage 2: REFINE around best --------
1486
+ refine_headings = _refine_range(best_h0, REFINE_HEADING_HALF_WINDOW_DEG, 1.0, GUIDED_HEADING_MIN, GUIDED_HEADING_MAX)
1487
+
1488
+ # Refine drift around best_d with DRIFT_STEP_DEG (0.5)
1489
+ d_lo = max(-drift_max_deg, best_d - REFINE_DRIFT_HALF_WINDOW_DEG)
1490
+ d_hi = min(+drift_max_deg, best_d + REFINE_DRIFT_HALF_WINDOW_DEG)
1491
+ refine_drifts = _frange(d_lo, d_hi, DRIFT_STEP_DEG)
1492
+ if not refine_drifts:
1493
+ refine_drifts = [best_d]
1494
+
1495
+ refine_turn_speeds = _refine_range(
1496
+ best_turn_spd, REFINE_TURN_SPEED_HALF_WINDOW_KT, TURN_SPEED_STEP_KT,
1497
+ max(TURN_SPEED_MIN_ABS_KT, float(south_speed_knots) - TURN_SPEED_SEARCH_DELTA_KT),
1498
+ min(TURN_SPEED_MAX_ABS_KT, float(south_speed_knots) + TURN_SPEED_SEARCH_DELTA_KT),
1499
+ )
1500
+
1501
+ for h0 in refine_headings:
1502
+ for d in refine_drifts:
1503
+ for spd_turn in refine_turn_speeds:
1504
+ s, combined_plus, legs = try_candidate(h0, d, spd_turn)
1505
+ if s is None:
1506
+ continue
1507
+ if best_combined_plus is None or float(combined_plus) < float(best_combined_plus):
1508
+ best_score = dict(s)
1509
+ best_combined_plus = float(combined_plus)
1510
+ best_h0 = float(h0)
1511
+ best_d = float(d)
1512
+ best_turn_spd = float(spd_turn)
1513
+ best_full_legs = legs
1514
 
1515
  report = "\n".join([
1516
+ "**Guided mode: Straight + drift + computed turn (18:22→18:40), clipped to 00:20**",
1517
+ f"- Search mode: **two-stage (coarse → refine)**",
1518
+ f"- Southern speed (speed profile): **{float(south_speed_knots):.1f} kt**",
1519
+ f"- Best turn speed (18:22→18:40): **{float(best_turn_spd):.1f} kt**",
1520
+ f"- Best nominal heading: **{best_h0:.1f}°**",
1521
+ f"- Drift used: **{best_d:+.1f}°** | Drift penalty weight: **{float(drift_weight):.2f}**",
1522
+ f"- BFO points: **{best_score['n_bfo']}** | RMSE_BFO: **{best_score['rmse_bfo']:.2f} Hz**",
1523
+ f"- BTO points: **{best_score['n_bto']}** | RMSE_BTO: **{best_score['rmse_bto']:.1f} µs**" if best_score.get("rmse_bto") is not None else
1524
+ f"- BTO points: **{best_score['n_bto']}** | RMSE_BTO: n/a",
 
 
 
 
 
 
 
 
 
1525
  ])
1526
 
1527
+ return best_full_legs, report, start_dt
1528
 
1529
 
1530
  # ================= RUNNERS =================
1531
  def _run_common(start_lat, start_lon, legs, method, start_dt, alt_ft, vs_fpm, descent_start_choice):
1532
+ # In guided mode, ensure the polyline ends at 00:20 by clipping legs to end_dt first.
1533
  if start_dt == RADAR_FIX_TIME_UTC:
1534
  end_dt = _dt_from_label(GUIDED_FINAL_LABEL)
1535
  legs = clip_legs_to_end_dt(legs, start_dt, end_dt)
 
1539
  html_map = build_leaflet_map(path_pts, legs, method, start_dt)
1540
  stats_md = compute_stats(path_pts, seg_dists_m, total_time_min, legs)
1541
 
1542
+ # BFO plot (restored)
1543
  bfo_fig, bfo_md = bfo_analysis_with_legs(
1544
  start_lat, start_lon, legs, method, start_dt, alt_ft, vs_fpm, descent_start_choice
1545
  )
 
1570
 
1571
  def run_guided(
1572
  start_lat, start_lon, method, alt_ft, vs_fpm, descent_start_choice,
1573
+ speed_profile, custom_mach, custom_knots, guided_scenario,
1574
+ weight_bfo, weight_bto,
1575
+ drift_max_deg, drift_weight
1576
  ):
1577
+ # This is the SOUTHERN cruise speed (post-18:40). Turn speed is searched independently.
1578
+ south_spd_knots = speed_from_profile(speed_profile, alt_ft, custom_mach, custom_knots)
1579
+
1580
+ if (guided_scenario or "").startswith("Straight + drift"):
1581
+ legs, report_md, start_dt = solve_guided_straight_with_drift(
1582
+ method, south_spd_knots,
1583
+ alt_ft, vs_fpm,
1584
+ descent_start_choice, weight_bfo, weight_bto,
1585
+ drift_max_deg, drift_weight
1586
+ )
1587
+ else:
1588
+ legs, report_md, start_dt = solve_guided_straight(
1589
+ method, south_spd_knots,
1590
+ alt_ft, vs_fpm,
1591
+ descent_start_choice, weight_bfo, weight_bto
1592
+ )
1593
 
1594
  if legs is None:
1595
  fig = go.Figure()
 
1603
  ], value=[])
1604
  return gr.update(), f"### ❌ {report_md}", "", "", fig, "", alt_fig, empty_table
1605
 
1606
+ # Guided always starts at radar fix
1607
  html_map, stats_md, bfo_fig, bfo_md, alt_fig, headers, rows = _run_common(
1608
  RADAR_FIX_LAT, RADAR_FIX_LON, legs, method, start_dt, alt_ft, vs_fpm, descent_start_choice
1609
  )
 
1635
  f"- **Turn completion (fixed):** 18:40:00Z\n"
1636
  f"- **Initial track (fixed):** {RADAR_INITIAL_TRACK_DEG:.1f}°\n"
1637
  f"- **Final resting point:** **00:20:00Z** (trajectory clipped exactly at this time)\n\n"
1638
+ f"**Turn behavior:** 18:22→18:40 turn is solved with **adjustable turn speed**.\n"
1639
+ f"**Performance:** guided solver uses **two-stage search (coarse → refine)** to run faster."
 
 
1640
  )
1641
 
1642
  with gr.Row():
 
1655
  label="When does descent start? (VS applies after this time)",
1656
  )
1657
 
1658
+ gr.Markdown("### Guided mode (Recommended)")
1659
+ with gr.Row():
1660
+ guided_scenario = gr.Dropdown(
1661
+ [
1662
+ "Straight (single heading)",
1663
+ "Straight + drift (per-leg ±deg, penalized)",
1664
+ ],
1665
+ value="Straight (single heading)",
1666
+ label="Scenario",
1667
+ )
1668
+ speed_profile = gr.Dropdown(
1669
+ ["MRC Mach 0.82", "LRC Mach 0.84", "ECON Mach 0.80", "Custom Mach", "Custom knots"],
1670
+ value="MRC Mach 0.82",
1671
+ label="Speed profile (sets SOUTHERN speed after 18:40)",
1672
+ )
1673
+ custom_mach = gr.Number(value=0.82, label="Custom Mach (if selected)")
1674
+ custom_knots = gr.Number(value=470, label="Custom knots (if selected)")
1675
+
1676
+ gr.Markdown("### Heading drift settings (for “Straight + drift” scenario)")
1677
+ with gr.Row():
1678
+ drift_max_deg = gr.Slider(0.0, 5.0, value=DEFAULT_DRIFT_MAX_DEG, step=0.1, label="Max drift per leg (±deg)")
1679
+ drift_weight = gr.Slider(0.0, 2.0, value=DEFAULT_WEIGHT_DRIFT, step=0.05, label="Drift penalty weight (higher = less drift)")
1680
 
1681
  gr.Markdown("### Guided scoring weights")
1682
  with gr.Row():
 
1684
  weight_bto = gr.Slider(0.0, 10.0, value=DEFAULT_WEIGHT_BTO, step=0.1, label="Weight: BTO (increase to emphasize BTO)")
1685
 
1686
  with gr.Row():
1687
+ guided_btn = gr.Button("Generate: Guided scenario")
1688
  manual_btn = gr.Button("Simulate: Manual legs")
1689
 
1690
  status_md = gr.Markdown()
 
1731
  run_guided,
1732
  inputs=[
1733
  start_lat, start_lon, method, alt_ft, vs_fpm, descent_start_choice,
1734
+ speed_profile, custom_mach, custom_knots, guided_scenario,
1735
+ weight_bfo, weight_bto,
1736
+ drift_max_deg, drift_weight
1737
  ],
1738
  outputs=[legs_df, status_md, out_map, stats, bfo_plot, bfo_score, alt_plot, validation_df],
1739
  )