Add Patera's contour integral for shaped hard bodies

README.md

@@ -81,13 +81,15 @@ projection is what makes a two-dimensional formulation valid: under linear relat
 the secondary crosses that plane in a straight line, so the three-dimensional question
 becomes a two-dimensional one about a region.
 
-The probability is the mass of a bivariate Gaussian inside that region, and it is evaluated
-two independent ways. Foster and Estes (1992) integrate in polar coordinates with adaptive
-quadrature. Alfano (2005a) performs the inner integral analytically with the error function
-and applies Simpson's rule to what is left. The two share nothing but the reduction to
-principal axes, so agreement between them checks both. Chan's series (2008) and a Monte Carlo
-estimator that samples the three-dimensional relative position provide two further checks of
-a different kind.
+The probability is the mass of a bivariate Gaussian inside that region, and it is evaluated three
+independent ways. Foster and Estes (1992) integrate in polar coordinates with adaptive quadrature.
+Alfano (2005a) performs the inner integral analytically with the error function and applies
+Simpson's rule to what is left. Patera (2001) applies Green's theorem and integrates around the
+boundary of the region instead of over its interior, so the region enters only as the curve that
+bounds it and an ellipse costs a different curve rather than a different derivation. The three
+share nothing but the reduction to principal axes, so agreement between them checks all three.
+Chan's series (2008) and a Monte Carlo estimator that samples the three-dimensional relative
+position provide two further checks of a different kind.
 
 ```python
 from conjunction_screening import generate_catalog, run_screening
@@ -285,7 +287,7 @@ uv run ruff format --check .
 uv run mypy
 ```
 
-209 tests cover 95.75 percent of the 1696 statements in the package. Continuous integration
+222 tests cover 95.90 percent of the 1754 statements in the package. Continuous integration
 runs that same command with `--cov-fail-under=93` on Ubuntu and on Windows, which is the
 measured figure rounded down and given two points of headroom, so that a platform difference
 in which branch a filter takes cannot fail a build on its own.
@@ -303,14 +305,14 @@ equal-period circles crossing a quarter of a revolution out of phase. The Lipsch
 the orbit path filter depends on is checked against a finely sampled numerical derivative,
 because a filter whose bound is not a bound could discard a real conjunction.
 
-Other invariants covered: the state transition matrix is symplectic; the time of closest
-approach has zero relative range rate; miss distance is symmetric under swapping the two
-objects; the encounter plane projection preserves the magnitude of a perpendicular relative
-position; the covariance stays symmetric and positive semi-definite through every stage;
-Foster and Alfano agree; Chan is exact for a circular covariance and departs monotonically as
-the aspect ratio grows; the combined hard body contains the Minkowski sum of the two bodies in
-every direction; the dilution curve rises then falls; and a Monte Carlo estimate agrees with
-every analytic value.
+Other invariants covered: the state transition matrix is symplectic; the time of closest approach
+has zero relative range rate; miss distance is symmetric under swapping the two objects; the
+encounter plane projection preserves the magnitude of a perpendicular relative position; the
+covariance stays symmetric and positive semi-definite through every stage; Foster, Alfano, and
+Patera agree on a disc, and Foster and Patera still agree on an outline that is not one; Chan is
+exact for a circular covariance and departs monotonically as the aspect ratio grows; the combined
+hard body contains the Minkowski sum of the two bodies in every direction; the dilution curve
+rises then falls; and a Monte Carlo estimate agrees with every analytic value.
 
 Two rules govern the tolerances. Only values from a converged solve are pinned, and the
 regression module asserts that every pinned event converged, because the state of a
@@ -365,6 +367,11 @@ Methods:
   Space Center, August 1992. Stable record:
   [Stanford SearchWorks 13354320](https://searchworks.stanford.edu/view/13354320). Source of
   the polar quadrature formulation of the two-dimensional probability of collision.
+- Patera, R. P. "General Method for Calculating Satellite Collision Probability." Journal of
+  Guidance, Control, and Dynamics, Vol. 24, No. 4, 2001, pp. 716 to 722.
+  DOI [10.2514/2.4771](https://doi.org/10.2514/2.4771). Source of the reduction of the area
+  integral to a contour integral around the boundary of the hard body outline, which is the
+  analytic method that does not assume the outline is a circle.
 - Alfano, S. "A Numerical Implementation of Spherical Object Collision Probability." The
   Journal of the Astronautical Sciences, Vol. 53, No. 1, 2005, pp. 103 to 109.
   DOI [10.1007/BF03546397](https://doi.org/10.1007/BF03546397). Source of the reduction of

docs/design-notes.md

@@ -111,11 +111,36 @@ each term is one minus a partial sum of a Poisson series, which cancels catastro
 the ratio of hard body radius to covariance scale is small, and that ratio is small in every
 real conjunction.
 
+Patera, "General Method for Calculating Satellite Collision Probability", Journal of Guidance,
+Control, and Dynamics 24(4), 2001, DOI 10.2514/2.4771, converts the area integral into a
+contour integral around the boundary of the region with Green's theorem. In the plane scaled
+by the two principal standard deviations the density is isotropic and has a vector potential
+in closed form, `(1 - exp(-r^2 / 2)) / r^2` times the perpendicular of the radius vector
+measured from the miss point, so the probability becomes an integral in one variable along a
+closed curve. That integrand is periodic and analytic in the variable, which is the regime the
+trapezoidal rule converges geometrically in, so nothing adaptive is needed: doubling from 128
+nodes meets the 1e-11 tolerance at 256 on every case in the test suite. Patera agrees with
+Foster to 2.5e-15 relative over the eight comparison cases, the level Alfano reaches, and
+unlike Alfano it keeps that agreement once the region stops being a disc.
+
+One numerical point decides whether the contour method is usable at screening depths. Its
+kernel splits as `1 / r^2` minus `exp(-r^2 / 2) / r^2`, and the first part integrates to the
+winding number of the outline about the miss point, which is one when the miss vector lies
+inside the body and zero otherwise. Left inside the quadrature it cancels away: on a
+probability of 1e-234 the terms it contributes are of order one and the sum returns 5e-21,
+which is noise. Subtracted in closed form, what is left is the decaying part alone and the
+value is right. The subtracted form divides by zero when the miss point lies exactly on the
+outline, where the other form is well behaved, so both are kept and the one whose integrand
+has the smaller peak is used. That comparison needs no threshold: the two forms differ by a
+term that is known exactly, so the smaller integrand is by construction the one whose sum
+cancels less.
+
 The region being integrated over is a disc only when both objects are spheres. A hard body may
 be given as an ellipsoid instead, in which case the region is the shadow it casts along the
-relative velocity and the polar quadrature takes a radial limit that varies with the angle.
-Foster and Monte Carlo accept that region; Alfano and Chan reject it, for reasons recorded
-under closed limitations below.
+relative velocity: the polar quadrature takes a radial limit that varies with the angle, and
+the contour integral runs around an ellipse in place of a circle. Foster, Patera, and Monte
+Carlo accept that region; Alfano and Chan reject it, for reasons recorded under closed
+limitations below.
 
 The Monte Carlo estimator samples the three-dimensional relative position from the combined
 covariance and projects each draw onto the plane normal to the relative velocity, rather than
@@ -250,8 +275,9 @@ direction of approach. That is the whole content of the limitation.
 
 The probability integral then runs over that ellipse. Foster's polar quadrature needs one
 change: the upper limit of the radial integral becomes `R(theta)`, and the angle moves from
-the inner to the outer variable. Monte Carlo needs one change: the hit test becomes a
-quadratic form in the plane rather than a distance in three dimensions.
+the inner to the outer variable. Patera's contour integral needs one change: the curve it runs
+around becomes the ellipse. Monte Carlo needs one change: the hit test becomes a quadratic
+form in the plane rather than a distance in three dimensions.
 
 What it cost. Four things, none of them hidden.
 
@@ -260,10 +286,12 @@ expands a non-central chi-square tail about a circular region; in both derivatio
 enters before the density does, so neither extends by changing a limit. They raise on a
 non-circular cross section rather than substituting a disc of the same area, because a
 silently substituted disc would produce a number that looks like a cross check of the Foster
-result and is not one. The consequence is that the strongest validation available, two
-independent quadratures agreeing to 4.3e-16, exists only for spheres. For a non-spherical body
-the cross check is Foster against Monte Carlo, which is four binomial standard errors wide
-rather than at machine precision.
+result and is not one. Two of the five methods are therefore unavailable for a shaped body,
+and that is the cost. It is not a loss of cross validation: Patera's contour integral does not
+use the circle either, so it follows an elliptical outline by changing the curve rather than
+the formulation, and it agrees with Foster to 1.0e-15 relative on the elliptical cases in the
+test suite, the level the two disc methods reach on a disc. Monte Carlo stays the check that
+also covers the covariance square root and the projection, four binomial standard errors wide.
 
 The combined body is an outer approximation, so the probability it produces is an
 overestimate. The size of that is measured rather than asserted: for a 30 m by 3 m by 3 m body

src/conjunction_screening/__init__.py

@@ -21,6 +21,7 @@
     CloseApproachSettings,
     FosterMethod,
     MonteCarloMethod,
+    PateraMethod,
     ProbabilityMethod,
     ProbabilityResult,
     find_close_approaches,
@@ -51,6 +52,7 @@
     "KeplerianElements",
     "MonteCarloMethod",
     "OrbitState",
+    "PateraMethod",
     "ProbabilityMethod",
     "ProbabilityResult",
     "ScreeningConfig",

src/conjunction_screening/algorithm/__init__.py

@@ -39,10 +39,12 @@
     CHAN,
     FOSTER,
     MONTE_CARLO,
+    PATERA,
     AlfanoMethod,
     ChanMethod,
     FosterMethod,
     MonteCarloMethod,
+    PateraMethod,
     ProbabilityMethod,
     ProbabilityResult,
 )
@@ -61,6 +63,7 @@
     "CHAN",
     "FOSTER",
     "MONTE_CARLO",
+    "PATERA",
     "PATH_FILTER",
     "PERIGEE_APOGEE_FILTER",
     "TIME_FILTER",
@@ -74,6 +77,7 @@
     "FosterMethod",
     "MaximumProbability",
     "MonteCarloMethod",
+    "PateraMethod",
     "PathSeparation",
     "ProbabilityMethod",
     "ProbabilityResult",