Function p_hit_bivariate 

pub fn p_hit_bivariate(
    var_drop: f64,
    var_wind: f64,
    cov: f64,
    target: TargetGeometryV1,
) -> f64

P(impact falls inside target) for a bivariate normal impact distribution centred at the origin – the nominal (zero-mean) trajectory solution – with drop variance var_drop, windage variance var_wind, and drop/windage covariance cov. target is always centred on that same origin; see TargetGeometryV1’s doc for why there is no separate aim-point offset.

The math (pinned, MBA-1347 spec section 6.2)

A correlated covariance does NOT let the rectangle probability separate into a product of two normal_cdf differences (uncorrelated_rectangle_matches_the_separable_closed_form in this module’s tests exists specifically to demonstrate the separable form is only valid at zero correlation, and a_strongly_correlated_case_differs_materially_from_the_wrong_separable_approximation shows how far a real, moderately-correlated case departs from it). Instead this integrates over the drop axis and, at each drop value u, applies the CONDITIONAL normal distribution of windage given that drop:

P = integral of phi(u) * [Phi(beta(u)) - Phi(alpha(u))] du,

where phi is the drop marginal’s density, Phi is normal_cdf, and alpha(u)/beta(u) are the target’s windage bounds at drop u (constant for a rectangle; the circle’s chord for a circle), expressed in units of the conditional windage standard deviation and offset by the conditional mean rho * (sigma_windage / sigma_drop) * u.

Two degenerate covariances, handled explicitly (not merely “does not panic”)

• Zero total variance (var_drop <= 0.0 && var_wind <= 0.0): a deterministic impact exactly at the origin, which is always the target’s own centre here – inside by definition, so this returns 1.0 outright without touching the quadrature below.
• Zero drop variance alone (var_drop <= 0.0, var_wind > 0.0): drop is deterministic at 0 but windage is not. The quadrature above integrates OVER the drop axis, which cannot represent a Dirac delta; instead of letting sd -> 0 silently zero out every quadrature node’s density (which the naive translation of the formula above does, and which would wrongly report 0.0 regardless of target size – caught by drop_deterministic_windage_random_matches_closed_form_not_hardcoded_zero in this module’s tests), this evaluates the windage marginal directly at drop = 0.

A degenerate target

Checked FIRST, before either degenerate-covariance branch above: a target with no positive area (width_m <= 0.0 || height_m <= 0.0 for a rectangle; radius_m <= 0.0 for a circle – after the same .max(0.0) treatment negative dimensions get elsewhere) can never be hit, returning 0.0 unconditionally, REGARDLESS of the covariance. Without this guard, a zero-size target combined with zero total variance would fall into the “deterministic impact” branch above and report 1.0 – technically defensible under a boundary-inclusive convention (a point impact exactly at a zero-size target’s own centre), but indistinguishable from a caller’s degenerate-target bug silently reading as total confidence, the single most misleading number this report could produce. error_budget_with_target reaches this with an empty sources list (zero total variance from having nothing to accumulate) more easily than it might seem, so this is checked unconditionally rather than only when sources happens to be empty.

A single declared source (one nonzero-sigma axis) produces a RANK-1 covariance – cov exactly +-sigma_drop * sigma_windage before the correlation clamp below – which is the routine case error_budget_with_target hits every time it prices “if this source alone were perfected” against a row with exactly one OTHER remaining source (see SourceContributionV1::p_hit_gain_if_perfect’s doc). rho is clamped to [-0.999_999, 0.999_999] so the conditional variance below is never exactly zero, avoiding a division by zero without needing a third special case for perfect correlation.

Why the quadrature is PANELLED, not a single 20-node call over the whole domain

The spec pins “fixed-order Gauss-Legendre quadrature over a truncated +/-6 sigma domain,” but a single 20-node rule spread across the WHOLE [-6 sigma_drop, 6 sigma_drop] interval is not merely imprecise, it is badly wrong for realistic inputs, for two DIFFERENT reasons, both found by comparing that naive translation against an independent fine (4000-point-per-smooth- piece composite Simpson) reference over a broad sweep of target sizes and correlations, not guessed:

First: the target boundary is a genuine discontinuity in the naive formulation. Written as “integrate over the whole +/-6 sigma domain, contributing zero outside the target’s own drop extent,” the integrand jumps from a generic nonzero value to zero exactly at the target’s edge. Gauss-Legendre quadrature assumes smoothness across its whole panel; a hidden jump degrades it to first-order accuracy. Measured on this module’s own pinned zero-correlation rectangle test (drop/windage sigma 0.10 m/0.20 m, target 0.30 m x 0.40 m): the naive single-panel translation is wrong by 8.6e-3 against the closed form (the test requires < 1e-6) – for a SMALL circular target relative to sigma (radius one-fifth of sigma: radius 0.02 with sigma 0.1 in both axes) the naive version places every one of its 20 nodes outside the target entirely and returns exactly 0.0 against a true value near 0.020.

Fix: restrict the integration domain to the target’s own drop extent intersected with the +/-6 sigma truncation ([lo, hi] below) – outside that range the contribution is EXACTLY zero, not approximately zero, so there is nothing to lose by not integrating there at all. This alone brings the zero-correlation rectangle case to ~1e-16 (machine precision, since the integrand reduces to a constant windage factor times a plain Gaussian bump, which 20-point Gauss-Legendre integrates essentially exactly).

Second: near-perfect correlation creates a separate, INTERNAL sharp transition the target boundary fix does not touch. As |rho| -> 1, the conditional windage standard deviation sigma_w * sqrt(1 - rho^2) -> 0, so the bracketed [Phi(beta(u)) - Phi(alpha(u))] factor above becomes an increasingly steep (though still, short of exactly rho = +-1, smooth) sigmoid in u, centred wherever the conditional mean crosses the window bound – a location that can fall anywhere inside the domain, not just at its edges. This is not a rare input: EVERY “if this source alone were perfected” comparison in error_budget_with_target evaluates a covariance with exactly one remaining source, which is exactly rank-1 (rho at the +-0.999_999 clamp). Measured on that exact shape (a real single-source covariance, sigma_drop = 18.5, sigma_windage = 6.0, rho clamped to -0.999_999, swept over target heights/widths from 0.05x to 20x each sigma): the boundary-restricted-but-still-single-panel quadrature is wrong by up to 0.28 against the fine reference at the swept extremes (height 20x sigma_drop, width 1x sigma_windage), and by up to 7.4e-2 even restricted to height/width within 0.5x-4x of the natural (sigma_drop, sigma_windage) scale – nowhere near a contrived corner.

Fix: when |rho| is non-negligible, ALSO split the domain at the (closed-form) drop value(s) where the conditional mean crosses the window’s bound – +-half_width * sigma_drop / (rho * sigma_windage) for a rectangle (the window bound is constant), or +-radius / sqrt(1 + (rho * sigma_windage / sigma_drop)^2) for a circle (from solving the circle’s own chord equation) – giving the sigmoid its own smooth sub-panel instead of sharing one with the flat shoulder on either side. This brings the worst case measured over a broad synthetic stress sweep (several (sigma_drop, sigma_windage) magnitudes including the realistic 18.5/6.0 pair above, target heights/widths from 0.05x to 20x sigma, |rho| up to 0.999_999, both shapes) down under 1e-3 (6.1e-4 for rectangles, 2.3e-4 for circles) – and realistic (roughly comparable width/height, moderate correlation) target shapes measured one to two further orders of magnitude better than that worst case.

This is a more careful IMPLEMENTATION of the pinned formula, not a different one: the number of panel BOUNDARIES is bounded at compile time (at most 4: the two domain edges – the target’s own edge only ever contributes 0 extra boundaries since it already bounds [lo, hi] – plus up to two correlation-crossing points), so there are at most 3 panels, each integrated by the exact same 20-node rule named in the spec. The result therefore stays deterministic and its cost stays bounded (at most 3 panels * 20 nodes * 2 normal_cdf calls per node = 120 evaluations of normal_cdf, negligible next to the real trajectory solves error_budget_with_target needs to build the covariance in the first place).

Bounded and monotone

Always clamped to [0.0, 1.0] before returning. The TRUE integral is monotone non-decreasing in target size for a fixed covariance (a bigger rectangle or circle strictly contains a smaller one centred at the same origin, so the region of integration only grows) – but that is a property of the exact mathematical integral, not something the panelled quadrature gets for free: growing the target moves the panel boundaries (the domain-restriction edge, and, at nonzero rho, the correlation-crossing points), which relocates every one of the 20 nodes within the affected panels, so the COMPUTED value is not automatically a monotone functional of target size the way the true integral is. This is verified BY TEST across a range of target sizes and shapes and, since the near-degenerate correlated regime is exactly where panel boundaries move the most, across rho in {0.0, 0.9, 0.999_999} too – see p_hit_is_bounded_and_grows_with_target_size (circle, verbatim from the spec), p_hit_grows_with_target_size_for_a_rectangle_too, and p_hit_is_monotone_in_target_size_across_a_sweep_including_near_rank_one_correlation in this module’s tests – not guaranteed by construction of the floating-point implementation.