open-flight · jewbetcha · Sep 23, 2026 · Sep 14, 2026 · Sep 16, 2026 · Sep 23, 2026
@@ -25,14 +25,17 @@ The coefficients depend on the **spin parameter** $S_p = r\omega / v$ — the
 ratio of surface speed to translational speed:
 
 $$
-C_d = C_{d,\text{base}} + k_d S_p
+C_d = a + b\,S_p + c\,S_p^2
 \qquad
-C_l = \frac{C_{l,\text{sat}} \, S_p}{C_{l,\text{half}} + S_p}
+C_l = \max\!\left(0,\; d + e\,S_p + f\,S_p^2\right)
 $$
 
-Drag rises linearly with spin. Lift follows a Hill-type saturating form: it
-approaches a ceiling as spin increases rather than growing without bound, which
-is what the wind-tunnel data shows.
+Both are second-order polynomials in $S_p$, the form **Ferguson, McNally &
+McPhee (2022)** fitted to 1040 measured shots; the coefficients are their
+published values. Lift rises with spin, peaks near $S_p \approx 0.52$ and then
+falls, which is what the measured data shows and what a saturating form cannot
+represent. Above $S_p = 0.75$, the top of the fitted range, both curves are held
+at their end value rather than extrapolated.
 
 These parametric forms are consistent with **Bearman & Harvey (1976)** and
 **Kensrud & Smith (2018)** for dimpled balls past the drag crisis
@@ -54,10 +57,9 @@ matters across a six-second flight.
 | `BALL_MASS_KG` | 0.04593 | USGA **maximum**-conforming ball, 45.93 g |
 | `BALL_RADIUS_M` | 0.02135 | 42.7 mm diameter |
 | `AIR_DENSITY_STD` | 1.225 kg/m³ | Sea level, 15 °C ISA |
-| `CD_BASE` | 0.205 | Drag at zero spin |
-| `CD_SPIN_COEFF` | 0.18 | Linear drag rise with $S_p$ |
-| `CL_SATURATION` | 0.32 | Lift ceiling |
-| `CL_HALF_SP` | 0.15 | $S_p$ at half the lift ceiling |
+| `CD_POLY` | (0.1304, 0.9287, -0.8259) | $C_d$ polynomial $(a, b, c)$ |
+| `CL_POLY` | (0.0504, 1.2031, -1.1490) | $C_l$ polynomial $(d, e, f)$ |
+| `SP_FIT_MAX` | 0.75 | Top of the fitted $S_p$ range; curves held beyond it |
 | `SPIN_DECAY_RATE` | 0.04 /s | ≈4 %/s |
 | `GRAVITY` | 9.81 m/s² | |
 | `DT_SECONDS` | 0.002 | 500 Hz integration |

@@ -91,10 +91,9 @@ references.
 | `BALL_MASS_KG` | 0.04593 |
 | `BALL_RADIUS_M` | 0.02135 |
 | `AIR_DENSITY_STD` | 1.225 |
-| `CD_BASE` | 0.205 |
-| `CD_SPIN_COEFF` | 0.18 |
-| `CL_SATURATION` | 0.32 |
-| `CL_HALF_SP` | 0.15 |
+| `CD_POLY` | (0.1304, 0.9287, -0.8259) |
+| `CL_POLY` | (0.0504, 1.2031, -1.1490) |
+| `SP_FIT_MAX` | 0.75 |
 | `SPIN_DECAY_RATE` | 0.04 |
 | `GRAVITY` | 9.81 |
 | `DT_SECONDS` | 0.002 |

@@ -1,11 +1,9 @@
-"""Tune the four aerodynamic coefficients in
+"""Tune the six polynomial aerodynamic coefficients in
 :mod:`openflight.ballistics` against TrackMan-measured carry.
 
-Parameters swept:
-    CD_BASE          — drag coefficient at zero spin
-    CD_SPIN_COEFF    — slope of Cd(Sp) (linear in spin parameter)
-    CL_SATURATION    — Cl asymptote at high Sp
-    CL_HALF_SP       — Sp at which Cl reaches CL_SATURATION/2
+Parameters swept (Cd = a + b*Sp + c*Sp^2, Cl = d + e*Sp + f*Sp^2):
+    CD_POLY = (a, b, c)
+    CL_POLY = (d, e, f)
 
 Method: scipy.optimize.differential_evolution (global) followed by
 Nelder-Mead refinement. Loss is overall RMSE on the TrackMan-inputs
@@ -77,23 +75,21 @@
 # dimpled golf balls in the post-drag-crisis regime; narrow enough that
 # differential evolution converges in a few minutes.
 PARAM_BOUNDS: List[Tuple[float, float]] = [
-    (0.16, 0.28),   # CD_BASE
-    (0.00, 0.40),   # CD_SPIN_COEFF
-    (0.18, 0.42),   # CL_SATURATION
-    (0.03, 0.35),   # CL_HALF_SP
+    (0.08, 0.30),  # Cd a: drag at zero spin
+    (0.00, 1.80),  # Cd b
+    (-1.80, 0.00),  # Cd c
+    (0.00, 0.15),  # Cl d
+    (0.40, 2.20),  # Cl e
+    (-2.20, -0.30),  # Cl f: negative so lift peaks and turns over
 ]
-PARAM_NAMES = ["CD_BASE", "CD_SPIN_COEFF", "CL_SATURATION", "CL_HALF_SP"]
-DEFAULT_COEFFS = (
-    bl.CD_BASE,
-    bl.CD_SPIN_COEFF,
-    bl.CL_SATURATION,
-    bl.CL_HALF_SP,
-)
+PARAM_NAMES = ["cd_a", "cd_b", "cd_c", "cl_d", "cl_e", "cl_f"]
+Coeffs = Tuple[float, ...]
+DEFAULT_COEFFS: Coeffs = (*bl.CD_POLY, *bl.CL_POLY)
 
 
 @dataclass
 class FitResult:
-    coeffs: Tuple[float, float, float, float]
+    coeffs: Coeffs
     rmse: float
     preds: List[float]
 
@@ -123,28 +119,28 @@ def _filter_shots(shots: List[TMShot]) -> List[TMShot]:
 
 def simulate_with_coeffs(
     shots: List[TMShot],
-    coeffs: Tuple[float, float, float, float],
+    coeffs: Coeffs,
 ) -> List[float]:
     """Monkey-patch the ballistics module constants, run simulate() for
     every shot, restore the originals on exit.
 
-    Relies on ``ballistics._cd`` and ``ballistics._cl`` resolving the
-    constants at call time from the module's global namespace.
+    Relies on ``ballistics._cd`` and ``ballistics._cl`` resolving
+    ``CD_POLY``/``CL_POLY`` at call time from the module's global namespace.
     """
-    saved = (bl.CD_BASE, bl.CD_SPIN_COEFF, bl.CL_SATURATION, bl.CL_HALF_SP)
-    bl.CD_BASE, bl.CD_SPIN_COEFF, bl.CL_SATURATION, bl.CL_HALF_SP = coeffs
+    saved = (bl.CD_POLY, bl.CL_POLY)
+    bl.CD_POLY, bl.CL_POLY = tuple(coeffs[:3]), tuple(coeffs[3:])
     try:
         results = []
         for s in shots:
             traj = simulate(_build_conditions(s), air_density=TM_FLAT_AIR_DENSITY)
             results.append(traj.carry_yards)
         return results
     finally:
-        bl.CD_BASE, bl.CD_SPIN_COEFF, bl.CL_SATURATION, bl.CL_HALF_SP = saved
+        bl.CD_POLY, bl.CL_POLY = saved
 
 
 def make_loss(shots: List[TMShot], measured: np.ndarray):
-    """Closure that the optimizer can call with a 4-vector."""
+    """Closure that the optimizer can call with a 6-vector."""
     def _loss(x: np.ndarray) -> float:
         preds = simulate_with_coeffs(shots, tuple(x))
         return float(np.sqrt(np.mean((np.asarray(preds) - measured) ** 2)))
@@ -154,7 +150,7 @@ def _loss(x: np.ndarray) -> float:
 def evaluate(
     shots: List[TMShot],
     measured: np.ndarray,
-    coeffs: Tuple[float, float, float, float],
+    coeffs: Coeffs,
 ) -> FitResult:
     preds = simulate_with_coeffs(shots, coeffs)
     rmse = float(np.sqrt(np.mean((np.asarray(preds) - measured) ** 2)))
@@ -258,7 +254,7 @@ def write_scatter(
     preds_default: List[float],
     preds_fit: List[float],
     out_path: Path,
-    fit_coeffs: Tuple[float, float, float, float],
+    fit_coeffs: Coeffs,
 ) -> None:
     import matplotlib
 
@@ -318,10 +314,9 @@ def write_scatter(
         ax.legend(loc="best", fontsize=8)
 
     fig.suptitle(
-        f"Default vs fit (fit: CD_BASE={fit_coeffs[0]:.4f}, "
-        f"CD_SPIN_COEFF={fit_coeffs[1]:.4f}, "
-        f"CL_SATURATION={fit_coeffs[2]:.4f}, "
-        f"CL_HALF_SP={fit_coeffs[3]:.4f})",
+        "Default vs fit (fit: "
+        + ", ".join(f"{n}={v:.4f}" for n, v in zip(PARAM_NAMES, fit_coeffs))
+        + ")",
         fontsize=11,
     )
     fig.tight_layout()
@@ -464,10 +459,12 @@ def _logged_loss(x):
 
     fit_lines = [
         "Optimal coefficients (TM-inputs, RMSE objective):",
-        f"  CD_BASE        = {best_x[0]:.5f}   (default {DEFAULT_COEFFS[0]:.5f})",
-        f"  CD_SPIN_COEFF  = {best_x[1]:.5f}   (default {DEFAULT_COEFFS[1]:.5f})",
-        f"  CL_SATURATION  = {best_x[2]:.5f}   (default {DEFAULT_COEFFS[2]:.5f})",
-        f"  CL_HALF_SP     = {best_x[3]:.5f}   (default {DEFAULT_COEFFS[3]:.5f})",
+        *(
+            f"  {name:<5} = {fit:>9.5f}   (default {default:>9.5f})"
+            for name, fit, default in zip(PARAM_NAMES, best_x, DEFAULT_COEFFS)
+        ),
+        f"  CD_POLY = ({best_x[0]:.4f}, {best_x[1]:.4f}, {best_x[2]:.4f})",
+        f"  CL_POLY = ({best_x[3]:.4f}, {best_x[4]:.4f}, {best_x[5]:.4f})",
         "",
         f"Baseline RMSE: {baseline.rmse:.3f} yd",
         f"Fit RMSE:      {fit.rmse:.3f} yd",
@@ -533,8 +530,7 @@ def _logged_loss(x):
             loso_lines.append(
                 f"{ho_disp:24s}  {len(test_shots):>6d}  {sub_fit_rmse:>9.3f}  "
                 f"{test_rmse:>10.3f}  {test_bias:>+10.3f}  "
-                f"({sub_best[0]:.3f}, {sub_best[1]:.3f}, "
-                f"{sub_best[2]:.3f}, {sub_best[3]:.3f})"
+                "(" + ", ".join(f"{v:.3f}" for v in sub_best) + ")"
             )
         fit_lines.extend(loso_lines)
     elif args.loso:

@@ -74,6 +74,10 @@
 # → 1.1839 kg/m³. We round to 1.184 for the input.
 TM_FLAT_AIR_DENSITY = 1.184
 
+# TrackMan reports apex ("Max Height - Height") in feet while its distance
+# columns are yards. The model returns apex in yards, so convert to compare.
+YD_TO_FT = 3.0
+
 # Club-name → ClubType map. The comparison CSV uses normalized names like
 # "driver", "7-iron", "pw"; the raw TM CSV uses "7 Iron", "Driver", "PW".
 _CLUB_MAP: Dict[str, ClubType] = {
@@ -151,6 +155,10 @@ class TMShot:
     carry_yards: Optional[float]
     timestamp: str
     session: str = ""  # session label (typically derived from the source filename)
+    # TrackMan reports apex as "Max Height - Height" in FEET, unlike the
+    # yard-denominated distance columns. Kept in feet end to end so the
+    # numbers stay comparable to the raw export.
+    apex_feet: Optional[float] = None
 
 
 def _default_session_label(path: Path) -> str:
@@ -201,6 +209,7 @@ def load_trackman(path: Path, session: Optional[str] = None) -> List[TMShot]:
                 carry_yards=_to_float(row.get("Carry Flat - Length")),
                 timestamp=row.get("Date", "") or "",
                 session=session_label,
+                apex_feet=_to_float(row.get("Max Height - Height")),
             ))
     return shots
 
@@ -284,6 +293,11 @@ class ValidationRow:
     measured_carry_yards: float
     model_carry_yards: float
     delta_yards: float  # model - measured
+    # Apex, in FEET, and only where the source carries it (the paired
+    # comparison CSV does not). None means "not measured", never zero.
+    measured_apex_feet: Optional[float] = None
+    model_apex_feet: Optional[float] = None
+    delta_apex_feet: Optional[float] = None
 
 
 def _has_required_inputs(*vals) -> bool:
@@ -312,6 +326,8 @@ def validate_tm_inputs(
             spin_source="measured",
         )
         traj = simulate(conditions, air_density=air_density)
+        # apex_yards is a height in yards; TrackMan's column is feet.
+        model_apex_feet = traj.apex_yards * YD_TO_FT
         session_tag = s.session or "default"
         out.append(ValidationRow(
             source="tm",
@@ -327,6 +343,11 @@ def validate_tm_inputs(
             measured_carry_yards=s.carry_yards,
             model_carry_yards=traj.carry_yards,
             delta_yards=traj.carry_yards - s.carry_yards,
+            measured_apex_feet=s.apex_feet,
+            model_apex_feet=model_apex_feet,
+            delta_apex_feet=(
+                None if s.apex_feet is None else model_apex_feet - s.apex_feet
+            ),
         ))
     return out
 

@@ -38,31 +38,28 @@
 BALL_AREA_M2 = math.pi * BALL_RADIUS_M ** 2
 AIR_DENSITY_STD = 1.225  # kg/m³ at sea level, 15 °C ISA
 
-# Cd = CD_BASE + CD_SPIN_COEFF * Sp
-#   Linear rise with spin parameter Sp = r·ω/v.
-# Cl = CL_SATURATION * Sp / (CL_HALF_SP + Sp)
-#   Hill-type saturating form: Cl → CL_SATURATION as Sp → ∞,
-#   reaches CL_SATURATION/2 at Sp = CL_HALF_SP.
-# These are simple parametric forms consistent with Bearman & Harvey (1976)
-# and Kensrud & Smith (2018) for dimpled balls past the drag crisis
-# (Re ~ 5e4–2e5), which covers the full range of realistic golf shots.
+# Cd = CD_POLY[0] + CD_POLY[1]*Sp + CD_POLY[2]*Sp^2
+# Cl = CL_POLY[0] + CL_POLY[1]*Sp + CL_POLY[2]*Sp^2, clamped at >= 0
+#   Second-order polynomials in the spin parameter Sp = r*omega/v, the form
+#   Ferguson, McNally & McPhee (2022, ISEA 14, doi:10.5703/1288284317493)
+#   fitted to N=1040 shots (GCQuad launch conditions, FlightScope X3
+#   carry/apex/offline, Pro V1, Sp 0.02-0.75). The coefficients below are
+#   their published values, not re-fitted here.
+#   Lift peaks near Sp ~ 0.52 and turns over, as in Bearman & Harvey (1976)
+#   and Lyu, Kensrud & Smith (2018). The previous Hill-type saturating form
+#   was monotonic by construction, so no choice of its two constants could
+#   represent that peak, and it sat ~4 yd low on apex across every club.
+#   Beyond Sp = SP_FIT_MAX both curves are held at their end value rather
+#   than extrapolating a parabola into a region the fit never saw.
 #
-# Fitted with scripts/analysis/sweep_ballistic_coeffs.py against the committed
-# TrackMan capture (session_logs/OpenFlight-Test.Normalized.csv, 24 shots,
-# differential evolution + Nelder-Mead, rho=1.184 to match TrackMan "Flat").
-# Overall carry RMSE against TrackMan: 24.52 -> 3.97 yd. The previous
-# CL_HALF_SP of 0.15 sat well above the low-spin driver regime (driver
-# Sp ~ 0.05-0.08), so for drivers the lift curve never left its low-lift
-# regime and driver carry ran ~37 yd short. Irons and wedges (Sp ~ 0.21-0.63)
-# were already past CL_HALF_SP and ran long instead, which is why the error
-# flipped sign by club. Resulting Cl is ~0.13-0.16 for drivers, rising to
-# ~0.22-0.23 for irons and wedges; the iron/wedge values sit inside the
-# 0.18-0.25 band the cited sources report, while the driver values remain
-# below it.
-CD_BASE = 0.19071
-CD_SPIN_COEFF = 0.31588
-CL_SATURATION = 0.25544
-CL_HALF_SP = 0.04758
+# Validation, TrackMan measured launch conditions in, simulated flight out,
+# Hill form (as re-fit in #229) -> this form:
+#   committed 24-shot capture:            carry RMSE 3.97 -> 2.90 yd
+#   independent 5-session 279-shot set:   carry RMSE 3.48 -> 1.98 yd,
+#                                         apex RMSE 4.41 -> 0.97 yd
+CD_POLY = (0.1304, 0.9287, -0.8259)
+CL_POLY = (0.0504, 1.2031, -1.1490)
+SP_FIT_MAX = 0.75
 
 # Exponential spin decay: ω(t) = ω₀·exp(-rate·t).
 # ~4%/s per Kiratidis & Leinweber (2018); small but matters over ~6 s flights.
@@ -166,12 +163,18 @@ def resolve_launch(shot: Shot) -> Optional[LaunchConditions]:
     )
 
 
+def _poly(coeffs: tuple, sp: float) -> float:
+    return coeffs[0] + coeffs[1] * sp + coeffs[2] * sp * sp
+
+
 def _cd(sp: float) -> float:
-    return CD_BASE + CD_SPIN_COEFF * sp
+    return _poly(CD_POLY, min(sp, SP_FIT_MAX))
 
 
 def _cl(sp: float) -> float:
-    return CL_SATURATION * sp / (CL_HALF_SP + sp) if sp > 0 else 0.0
+    if sp <= 0:
+        return 0.0
+    return max(0.0, _poly(CL_POLY, min(sp, SP_FIT_MAX)))
 
 
 def _derivatives(