From 4db523325c02ebe93b881f9fabdd6c88a0ab25b7 Mon Sep 17 00:00:00 2001 From: cvrt-jh Date: Mon, 14 Sep 2026 16:54:24 +0200 Subject: [PATCH 1/2] fix(ballistics): replace Hill lift form with Ferguson quadratic Cd/Cl The saturating Hill form is monotonic, so it cannot represent the measured lift peak near Sp ~ 0.5 and left every club ~4 yd low at apex. Switch Cd and Cl to the second-order polynomials Ferguson, McNally & McPhee (2022) fitted to 1040 shots, using their published coefficients unchanged, and hold both curves at Sp = 0.75 so in-flight spin parameters beyond the fitted range do not extrapolate into negative lift or drag. Carry RMSE on the committed capture 3.97 -> 2.90 yd; on an independent 279-shot set 3.48 -> 1.98 yd carry and 4.41 -> 0.97 yd apex. Regression budgets ratcheted, sweep tool moved to six parameters, docs updated. Fixes #230 --- docs/how-it-works/ballistics.md | 20 +++--- docs/reference/constants.md | 7 +- scripts/analysis/sweep_ballistic_coeffs.py | 70 +++++++++----------- src/openflight/ballistics.py | 55 +++++++-------- tests/test_ballistics_trackman_regression.py | 27 +++++--- 5 files changed, 92 insertions(+), 87 deletions(-) diff --git a/docs/how-it-works/ballistics.md b/docs/how-it-works/ballistics.md index 708f519cb..f57b455fd 100644 --- a/docs/how-it-works/ballistics.md +++ b/docs/how-it-works/ballistics.md @@ -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 | diff --git a/docs/reference/constants.md b/docs/reference/constants.md index 0e2702e2a..9156d52bd 100644 --- a/docs/reference/constants.md +++ b/docs/reference/constants.md @@ -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 | diff --git a/scripts/analysis/sweep_ballistic_coeffs.py b/scripts/analysis/sweep_ballistic_coeffs.py index 35514a5e1..d5e9260f4 100644 --- a/scripts/analysis/sweep_ballistic_coeffs.py +++ b/scripts/analysis/sweep_ballistic_coeffs.py @@ -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,16 +119,16 @@ 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: @@ -140,11 +136,11 @@ def simulate_with_coeffs( 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: diff --git a/src/openflight/ballistics.py b/src/openflight/ballistics.py index 8f64c6b88..c7ba343bc 100644 --- a/src/openflight/ballistics.py +++ b/src/openflight/ballistics.py @@ -36,31 +36,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. @@ -184,12 +181,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( diff --git a/tests/test_ballistics_trackman_regression.py b/tests/test_ballistics_trackman_regression.py index 4d301d81c..fe3c18bfd 100644 --- a/tests/test_ballistics_trackman_regression.py +++ b/tests/test_ballistics_trackman_regression.py @@ -36,24 +36,29 @@ TRACKMAN_CSV = _REPO_ROOT / "session_logs" / "OpenFlight-Test.Normalized.csv" -# Per-club RMSE ceilings in yards. Ratcheted down after the aero coefficients -# were re-fit against this same capture (see ballistics.py CD_/CL_ constants). +# Per-club RMSE ceilings in yards. Ratcheted twice: after the aero +# coefficients were re-fit against this capture (#229), and again when the +# Cd/Cl functional form moved to the Ferguson quadratics (#230). # -# before re-fit after re-fit -# club rmse bias rmse bias -# 7-iron 11.64 +11.00 2.96 +1.52 -# driver 38.60 -37.18 1.75 -0.45 -# pitching wedge 13.56 +11.02 6.26 -2.96 -# OVERALL 24.52 -5.05 3.97 -0.45 +# pre-#229 #229 (Hill) quadratic +# club rmse bias rmse bias rmse bias +# 7-iron 11.64 +11.00 2.96 +1.52 1.33 +1.30 +# driver 38.60 -37.18 1.75 -0.45 2.38 -1.33 +# pitching wedge 13.56 +11.02 6.26 -2.96 4.48 +1.83 +# OVERALL 24.52 -5.05 3.97 -0.45 2.90 +0.58 +# +# Driver moved up within its ceiling: the quadratic coefficients are an +# external published fit, not tuned to this capture, and the driver ceiling +# was left where it was rather than raised. # # Ceilings sit slightly above the measured values to absorb float/platform # drift. Lower them again if the model improves; never raise them. RMSE_BUDGET_YARDS = { "driver": 3.0, - "7-iron": 4.0, - "pitching wedge": 7.5, + "7-iron": 2.0, + "pitching wedge": 5.5, } -OVERALL_RMSE_BUDGET_YARDS = 5.0 +OVERALL_RMSE_BUDGET_YARDS = 3.5 # The reference capture is a fixed, committed file: if the shot count changes, # the fixture changed and every budget above needs re-deriving. From 2308b02b16c3f00fc9eda4eccf2600a16cc251b4 Mon Sep 17 00:00:00 2001 From: cvrt-jh Date: Wed, 16 Sep 2026 12:35:15 +0200 Subject: [PATCH 2/2] test(ballistics): cover Cd/Cl clamp and add apex error budget The Sp clamp added with the Ferguson quadratics had no test: the full suite passed with both clamps removed, even though the raw parabolas give Cd = -0.35 at Sp 1.5 (negative drag accelerates the ball) and Cl < 0 above Sp 1.09. Add coefficient-level guards for zero-spin lift, non-negative lift and positive drag above the fitted range, that both curves are held flat past SP_FIT_MAX, and that the fitted region below it is untouched. Reverting the clamp fails three of them. Also add an apex error budget against the committed TrackMan capture, which ships apex as "Max Height - Height" in feet. Apex is out-of-sample evidence since nothing was ever fitted against it: overall RMSE is 15.05 ft on the Hill form against 2.56 ft on the quadratics, with the 7-iron at 18.15 -> 1.53 ft and the pitching wedge at 18.00 -> 2.81 ft. Restoring the Hill form fails all four apex budgets. Apex is plumbed through the validator's own loader rather than parsed separately, so the manual tool and the test cannot disagree about the fixture. Budgets are seeded above measured to absorb float drift, matching the carry budgets, and a coverage test pins apex presence so a changed fixture cannot turn them into vacuous passes over an empty list. --- scripts/analysis/validate_ballistics.py | 21 +++++ tests/test_ballistics.py | 85 ++++++++++++++++++++ tests/test_ballistics_trackman_regression.py | 64 +++++++++++++++ 3 files changed, 170 insertions(+) diff --git a/scripts/analysis/validate_ballistics.py b/scripts/analysis/validate_ballistics.py index d10a68180..2161cb84a 100644 --- a/scripts/analysis/validate_ballistics.py +++ b/scripts/analysis/validate_ballistics.py @@ -73,6 +73,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] = { @@ -150,6 +154,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: @@ -200,6 +208,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 @@ -283,6 +292,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: @@ -311,6 +325,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", @@ -326,6 +342,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 diff --git a/tests/test_ballistics.py b/tests/test_ballistics.py index 5cc6df705..71585779d 100644 --- a/tests/test_ballistics.py +++ b/tests/test_ballistics.py @@ -6,9 +6,17 @@ import pytest from openflight.ballistics import ( + BALL_RADIUS_M, + CD_POLY, + CL_POLY, CLUB_TYPICAL_SPIN_RPM, + MPH_TO_MPS, + SP_FIT_MAX, SPIN_DECAY_RATE, LaunchConditions, + _cd, + _cl, + _poly, resolve_launch, simulate, ) @@ -177,3 +185,80 @@ def test_zero_launch_angle_does_not_crash(self): def test_total_distance_includes_rollout(self): traj = simulate(_driver()) assert traj.total_yards > traj.carry_yards + + +class TestAeroCoefficientSafeguards: + """Guards on the Cd/Cl quadratics outside their fitted range. + + The Ferguson coefficients are fitted over Sp <= SP_FIT_MAX. Extrapolated + past it both parabolas turn over and cross zero (Cl above Sp ~1.09, Cd + above ~1.25), which would mean lift pulling the ball down and drag + *accelerating* it. Sp is recomputed every integration step from + r*omega/v, and v decays faster than spin, so a normal lob wedge reaches + Sp ~1.5 near apex - inside the negative-drag region. Holding both curves + at their end value is the guard; these tests are what make its removal + fail. + """ + + def test_zero_spin_produces_no_lift(self): + assert _cl(0.0) == 0.0 + + def test_lift_never_negative_above_fit_range(self): + # Raw Cl(1.5) is about -0.73; clamped it must not pull the ball down. + for sp in (SP_FIT_MAX, 1.0, 1.5, 3.0): + assert _cl(sp) >= 0.0, f"Cl({sp}) = {_cl(sp):.4f} is negative" + + def test_drag_stays_positive_above_fit_range(self): + # Raw Cd crosses zero near Sp 1.25; negative drag is unphysical. + for sp in (SP_FIT_MAX, 1.3, 1.5, 3.0): + assert _cd(sp) > 0.0, f"Cd({sp}) = {_cd(sp):.4f} is not positive" + + def test_both_curves_held_flat_above_fit_range(self): + # Beyond the fitted range the value is pinned to the endpoint, so the + # curves are constant there rather than continuing to fall. + for sp in (1.0, 1.5, 3.0): + assert _cd(sp) == pytest.approx(_cd(SP_FIT_MAX)) + assert _cl(sp) == pytest.approx(_cl(SP_FIT_MAX)) + + def test_curves_unclamped_inside_fit_range(self): + # The guard must not disturb the fitted region it sits above. + for sp in (0.1, 0.25, 0.5): + assert _cd(sp) == pytest.approx(_poly(CD_POLY, sp)) + assert _cl(sp) == pytest.approx(_poly(CL_POLY, sp)) + + def test_lob_wedge_reaches_clamped_region_in_flight(self): + """A real lob wedge drives Sp past the fitted range in flight. + + This is why the clamp is not hypothetical: 55 mph / 40 deg / + 10000 rpm peaks near Sp 1.5, where the raw parabola gives + Cd = -0.35 and Cl = -0.73. The clamp is what holds that at the + Sp 0.75 endpoint instead. + + Note this asserts reachability, not the clamp's presence: removing + the clamp *lowers* peak Sp (less drag keeps the ball faster), so it + cannot serve as a red control. The coefficient tests above do that. + """ + cond = LaunchConditions( + ball_speed_mph=55.0, + launch_angle_v=40.0, + launch_angle_h=0.0, + spin_rpm=10000, + spin_axis_deg=0.0, + spin_source="measured", + ) + traj = simulate(cond) + + # Recorded points are downsampled, so this is a lower bound on the + # true peak - which only makes the assertion stricter. + peak_sp = max( + BALL_RADIUS_M * (p.spin_rpm * 2 * math.pi / 60.0) / (p.speed_mph * MPH_TO_MPS) + for p in traj.points + if p.speed_mph > 1e-6 + ) + assert peak_sp > SP_FIT_MAX, ( + f"lob wedge peak Sp {peak_sp:.2f} no longer exceeds {SP_FIT_MAX}; " + f"the clamp would be untested by any realistic shot" + ) + assert 30.0 < traj.carry_yards < 90.0, ( + f"lob wedge carry {traj.carry_yards:.1f} yd is not plausible" + ) diff --git a/tests/test_ballistics_trackman_regression.py b/tests/test_ballistics_trackman_regression.py index fe3c18bfd..2aec1b0f6 100644 --- a/tests/test_ballistics_trackman_regression.py +++ b/tests/test_ballistics_trackman_regression.py @@ -60,6 +60,25 @@ } OVERALL_RMSE_BUDGET_YARDS = 3.5 +# Apex ceilings in FEET (TrackMan's "Max Height - Height" column is feet). +# Apex is the strongest evidence the quadratic form is right: nothing was ever +# fitted against it, so the improvement is out-of-sample. +# +# club Hill (#229) quadratic +# 7-iron 18.15 1.53 +# driver 5.08 3.17 +# pitching wedge 18.00 2.81 +# OVERALL 15.05 2.56 +# +# Seeded above measured to absorb float/platform drift, matching the carry +# budgets above. Lower them when the model improves; never raise them. +APEX_RMSE_BUDGET_FEET = { + "driver": 4.0, + "7-iron": 2.0, + "pitching wedge": 3.5, +} +OVERALL_APEX_RMSE_BUDGET_FEET = 3.0 + # The reference capture is a fixed, committed file: if the shot count changes, # the fixture changed and every budget above needs re-deriving. EXPECTED_SHOT_COUNT = 24 @@ -107,3 +126,48 @@ def test_per_club_carry_rmse_within_budget(validation_rows, club): f"(bias {stats['mean']:+.2f}, max |delta| {stats['max_abs']:.2f}, " f"n={stats['n']})." ) + + +def _apex_deltas(rows, club=None): + """Apex errors in feet, skipping rows whose source carried no apex.""" + return [ + r.delta_apex_feet + for r in rows + if r.delta_apex_feet is not None and (club is None or r.club == club) + ] + + +def test_reference_capture_has_apex_measurements(validation_rows): + """Pin apex coverage: a fixture silently losing the column would turn + every apex budget below into a vacuous pass over an empty list.""" + deltas = _apex_deltas(validation_rows) + assert len(deltas) == EXPECTED_SHOT_COUNT, ( + f"{len(deltas)} of {EXPECTED_SHOT_COUNT} reference shots carry an apex " + f"measurement. The fixture's 'Max Height - Height' column changed." + ) + + +def test_overall_apex_rmse_within_budget(validation_rows): + """Model apex vs TrackMan apex, fed TrackMan's own launch conditions.""" + stats = _stats(_apex_deltas(validation_rows)) + assert stats["rmse"] <= OVERALL_APEX_RMSE_BUDGET_FEET, ( + f"Overall apex RMSE {stats['rmse']:.2f} ft exceeds budget " + f"{OVERALL_APEX_RMSE_BUDGET_FEET} ft (bias {stats['mean']:+.2f}, " + f"max |delta| {stats['max_abs']:.2f}, n={stats['n']}). " + f"The ballistic model's flight height regressed." + ) + + +@pytest.mark.parametrize("club", sorted(APEX_RMSE_BUDGET_FEET)) +def test_per_club_apex_rmse_within_budget(validation_rows, club): + """Per-club apex budgets: trajectory shape is club-dependent, and the + quadratic's biggest win (irons and wedges) must not silently erode.""" + deltas = _apex_deltas(validation_rows, club) + assert deltas, f"No reference shots with apex for club {club!r}" + stats = _stats(deltas) + budget = APEX_RMSE_BUDGET_FEET[club] + assert stats["rmse"] <= budget, ( + f"{club}: apex RMSE {stats['rmse']:.2f} ft exceeds budget {budget} ft " + f"(bias {stats['mean']:+.2f}, max |delta| {stats['max_abs']:.2f}, " + f"n={stats['n']})." + )