Appendix N Pc Calculation Approaches 1

Appendix N. Pc Calculation Approaches – UPDATE 3/20/24

N.1 Introduction

This appendix discusses the selection of the Probability of Collision (Pc) as the risk assessment 

parameter to use for conjunction assessment requirements compliance. While this parameter is 

widely used in the conjunction assessment industry, issues related to its calculation exist that 

merit extended discussion. The most frequently used analytical techniques to calculate the Pc are 

well established and computationally efficient but include assumptions that restrict their use and 

make the calculations vulnerable to error for a small fraction of conjunctions. Numerical 

techniques also exist, such as Monte Carlo Pc estimation. Predictably, these make fewer 

assumptions and are more widely applicable, but they are much more computationally 

demanding. There are also issues related to the regularization and interpretation of the input data 

to the Pc calculation, some of which are resolved by techniques that are now becoming standard 

practices. In calculating Pc estimates, it is necessary to examine and prepare the input datasets 

carefully and then to choose the calculation approach that is appropriate to the situation. The 

purpose of this appendix is to amplify the Pc calculation-related recommendations by providing 

an extended technical explanation of the data preparation and Pc calculation issues so that 

implementers and users of these calculations can proceed with a better understanding of the 

different options and resultant fidelities of Pc calculation. 

To accomplish this goal, this appendix will address the following technical areas:

• A step-by-step description of the conjunction plane two-dimensional Pc calculation, 

which is the most established and widely used analytical Pc computation technique, and 

its enabling simplifications and assumptions;

• Examination and repair/expansion of input data to the calculation, focusing mostly on the 

orbital state covariance matrices for the two objects in conjunction;

• Discussion of the use of Monte Carlo techniques, which is the numerical method used for 

high-fidelity Pc computation;

• A test to determine whether an analytical or numerical technique should be used for a 

particular conjunction;

• Approaches to choosing for the Pc calculation the hard-body radius, which gives a 

statement of the combined sizes of the primary and secondary objects; and

• Discussion of alternative analytic Pc calculation methods, specifically the two- and threedimensional “Nc” (as opposed to “Pc”) estimation methods that address conjunctions 

affected by curvilinear trajectories and non-Gaussian distribution effects.

• Discussion of correcting Pc estimates for covariance cross correlation effects, which 

affects some conjunctions involving high-drag satellites.

N.2 Conjunction Plane Analytic Pc Calculation

National Aeronautics and Space Administration
Appendix N Pc Calculation Approaches 2

The conjunction plane method of Pc calculation, which is by far the most widely used approach 

in the conjunction assessment industry, was developed for the Space Shuttle Program and first 

described in the literature in 1992 (Foster and Estes). There have been a number of important 

treatments since that time ⎯ e.g., Akella and Alfriend (2000), Patera (2001), Alfano (2005a), 

Chan (2008), Garcia-Pelayo (2016), and Elrod (2019) ⎯ but all rely on the same basic 

methodology: applying reasonable assumptions to enable analytical approximations. While the

particulars vary, these approaches all share the same concept of calculating the Pc estimate by 

integrating over a two-dimensional region on a conjunction encounter plane.

Conjunction plane Pc analysis begins with the states and covariances for the primary and 

secondary objects’ orbital states at TCA. An important set of questions should be addressed 

concerning whether these data, especially the covariances, are truly representative of the 

expected states and uncertainties at TCA or whether the propagation process has distorted them. 

These questions will be addressed in a subsequent section, once a more focused context for them 

has been established through the present discussion of the general procedure. 

The first step is to recognize that Pc calculations depend on the relative positions and 

uncertainties of the two objects, so moving to a framework that views the problem this way is 

helpful. Subtraction of the two objects’ positions produces a relative position vector, the 

magnitude of which is the miss distance between the two objects (at TCA). Similarly, because 

interest is in the relative rather than absolute position uncertainty, it is possible to combine the 

two objects’ covariances to form a joint (relative) covariance and allocate that, if desired, to one 

“end” of the relative position vector (by convention the end for the secondary object), as shown 

in Figure N-1. There are, of course, questions regarding whether the two covariances are 

independent and can be combined by simple addition; this question is addressed in a later 

section, but it is typically acceptable to presume independence and combine the covariances in 

this manner. 

Figure N-1 Relative Positions and Uncertainties of Two Objects

The second step relates to modeling the combined sizes of the primary and secondary objects. 

The general approach is to place a circumscribing sphere about the primary (whose size is well 

known by the O/O because it is their satellite) and then do the same thing for the secondary but 

often via an estimation technique as the secondary is usually a debris object for which there is no 

definitive size information. If these two spheres approach and begin to overlap one another at 

any point during the encounter, then a potential collision has been identified. Again, keeping in 

mind that a relative framework is useful here, the two objects’ size spheres can be combined into 

one super-sphere (also called the “collision sphere”) and placed at one end, by convention the 
Appendix N Pc Calculation Approaches 3

primary object end, of the relative miss vector, as shown in Figure N-2. If the miss vector should 

shrink to be smaller than the radius of the collision sphere (also called the hard-body radius), that 

would be the equivalent of the two original spheres encroaching on each other at TCA. 

Figure N-2 Combining the Sizes of the Primary and Secondary Objects

The third step is to envision the situation at TCA in which all the uncertainty is assigned to the 

secondary object’s end of the relative miss vector held in a fixed position in the mind. The 

primary object end of the relative miss vector is moving along through TCA and bringing with it 

the sphere representing both objects’ combined size. Even though the satellites follow curved 

trajectories and the covariance evolves and changes at each instant, if the encounter is presumed 

to take place extremely quickly—and this is in most conjunctions a good assumption because the 

satellites’ relative velocity usually exceeds 10 km/sec—then two assumptions can be made: the 

trajectories are essentially rectilinear during the encounter period, and the covariances (and thus 

the joint covariance) can be considered static. This means that the encounter can be envisioned as 

in Figure N-3:

Figure N-3 Geometry at TCA

The passing of the primary object by the secondary can be seen as following a “soda straw” 

straight trajectory whose cylindrical radius is the same as that for the hard-body radius and 

whose placement is one miss distance away from the position of the secondary at TCA. Since the 

joint covariance shown as the ellipsoid above represents the uncertainty of one “end” of the 
Appendix N Pc Calculation Approaches 4

miss-distance vector (as shown by the central dot in the above diagram), this dot can be 

presumed to be potentially in any place within the ellipsoid, meaning that in some portion of 

those realizations, it will fall within the soda straw. When this is the case, a collision is 

presumed. The task is to determine the likelihood that the dot will, in the actual realization of this 

conjunction, fall within the cylindrical path (“soda straw”) swept out by the motion of the 

primary object. This probability can be decomposed to be rendered as the product of the 

individual probabilities that each component of the secondary object position (the dot) will fall 

within the soda straw pathway, i.e., if an orthogonal x-y-z coordinate system is presumed, this 

overall probability can be generated as the product of the probability that the x-component of the 

dot’s position will fall within the straw, the y-component of the dot’s position will fall within the 

straw, and the z-component of the dot’s position will fall within the straw. If this coordinate 

system is aligned so that one of the axes (e.g., the z axis) aligns with the direction of the straw, 

because one is assuming rectilinear motion, the soda straw can be presumed to be unbending and 

infinite in length, and as such, it will contain all of the z-component probability of the dot’s 

falling within the straw. As such, the z-component probability in this arrangement will be unity 

and will be what is called “marginalized out,” meaning that the overall probability can be fully 

represented by the probability remaining with the x- and y-components. The entire situation can 

thus be reduced from three to two dimensions and analyzed as a phenomenon that occurs on a

plane that is orthogonal to the soda-straw direction, which is the direction of the relative velocity 

vector. This procedure defines the “conjunction plane,” which can be viewed in two equivalent 

representations as discussed by Chan (2008, see Figure 5.1), and as shown in Figure N-4.

Figure N-4 Two Equivalent Representations of the Conjunction Plane

In the first representation, shown in the left panel of Figure N-4, the center of the secondary 

object is placed at the origin, and the local (x, y) coordinate system configured so that the miss 

vector (which extends a distance of rmiss from the secondary to the primary object) lies along the 

horizontal x axis. The “soda straw” is coming out of the page and represented as the hard-bodyradius circle. With this planar reduction, the Pc is now the frequency with which the missdistance vector will fall within the hard-body-radius circle; this is equivalent to the amount of 

joint covariance probability density (illustrated here using a 1 ellipse) falling within that circle. 

Since the projected joint covariance represents a bivariate Gaussian distribution, the amount of 
Appendix N Pc Calculation Approaches 5

its probability density falling within the hard-body-radius circle is given by the following twodimensional integral (Foster and Estes, 1992)

Pc =

1

√det(2πC)

∫∫ exp (−

r

TC−1r

2

) dx dy

A

(N-1)

in which r is the 21 miss vector, C is the 22 covariance matrix, and A is the hard-body-radius

circular area. There are several different approaches to evaluating this integral. Foster and Estes 

(1992) applied a two-dimensional quadrature technique; this works well, and MATLAB®’s1

adaptive quad2d integrator is quite equal to the task, although perhaps not the most 

computationally efficient of all possible choices. Chan (2008) uses equivalent-area transforms to 

produce a single-dimensional integral, which has a series solution. Garcia-Pelayo et al. (2016) 

also derives a multi-term series approximation. However, even though these series 

approximations are relatively computationally efficient, experience indicates that they 

occasionally produce inaccuracies. Elrod (2019) formulates the integral in terms of 

complementary error functions (to improve accuracy for conjunctions with small Pc values) and 

uses Gauss-Chebyshev quadrature with nodal symmetry (which improves the efficiency of the 

numerical integration considerably). The Foster and Elrod approaches (the former being the 

established standard, and the latter being extremely fast, especially with a vectorized MATLAB 

implementation) are included as the Pc2D_Foster and PcElrod functions in the Pc Computation 

Software Development Kit (SDK) in the NASA CARA software repository. (See Section 7, 

Contact Information in this document for the specific URL.)

In the second conjunction plane representation, shown in the right panel of Figure N-4, the 

secondary object location again is placed at the origin, but in this case the local (x, y) coordinate 

system is configured so that the principal axis of the covariance ellipse lies along the horizontal 

axis. The center of the primary object lies at (xmiss, ymiss), and the axes are oriented so that this 

miss position lies within the 1st quadrant, so that xmiss ≥ 0 and ymiss ≥ 0. Using this representation,

the Pc is given by a one-dimensional integral involving error functions (Alfano, 2005a) 

Pc =

1

√8πσx

∫ {erf(y+) − erf(y−)} exp [−

(u + xmiss)

2

2σx

2

] du

R

−R

(N-2)

with R indicating the combined hard-body radius, u = x

′ − xmiss the integration variable, and

y± = (ymiss ± √R2 − u

2)/(√2σy). The integrand factor in the curly brackets represents the 

difference of two error functions, ∆ = −erf(y−). In many software environments (e.g., 

MATLAB), this factor often can be computed significantly more accurately as a difference of 

complementary error functions, i.e., ∆ = erfc(y−) − erfc(y+), especially when calculating very 

small Pc values (Elrod, 2019). See Abramowitz and Stegun (1970) and Press et al. (1992) for 

details on computing erf(-) and erfc(-) functions.

For most conjunctions, Gauss-Chebyshev quadrature provides an efficient and accurate means to 

calculate the one-dimensional integral in equation (N-2), using an approach similar to that 

1 MATLAB is a registered trademark of The MathWorks, Inc.
Appendix N Pc Calculation Approaches 6

described by Elrod (2019). However, for conjunctions involving relatively large hard-body radii

(or, equivalently, small covariance ellipse σx dimensions), this method can potentially become

inaccurate. Specifically, Gauss-Chebyshev quadrature becomes increasingly inaccurate as the 

hard-body radius grows beyond the limit Rlim = 4σx − min(R). In such cases, which rarely 

occur in practice, an adaptive numerical integrator can be used to calculate an accurate estimate, 

e.g., MATLAB’s integral function. (The function PcConjPlaneCircle of the NASA CARA SDK 

repository implements a vectorized algorithm that automatically determines which of these two 

integration methods should be used to compute Pc values accurately for all input hard-body radii

and covariance values, with a computation speed comparable to that of the PcElrod function for 

most conjunctions. Also, for cases that have hard-body radii well below the limit given above, 

testing indicates that the SDK functions Pc_Foster, PcElrod, and PcConjPlaneCircle all output

nearly identical numerical Pc values, i.e., that typically agree to five digits of precision or more.)

The second conjunction plane representation shown in the right panel of Figure N-4 also 

provides an extremely efficient means to calculate an upper limit estimate for the Pc value. This 

upper bound corresponds to the two-dimensional (2-D) integral over the square that 

circumscribes the hard-body radius circle (as shown in Figure N-4), which has the following 

analytical solution

Pc < Psq =

[erf(X+) − erf(X−)] [erf(Y+) − erf(Y−)]

4

(N-3)

with

X± =

xmiss ± R

√2σx

and Y± =

ymiss ± R

√2σy

(N-4)

Again, in many cases the erf(-) differences in equation (N-3) can be computed more accurately

using erfc(-) differences. Notably, the circumscribing square probability estimate does not 

require any numerical integration at all, but instead only requires the computation of four erf(-) 

or erfc(-) functions, usually making it relatively easy to program into software and significantly

more computationally efficient. These considerations could be important in some circumstances 

(e.g., for computations performed on orbiting satellites), prompting the use of Psq as an

approximation for Pcitself. For instance, CDMs generated by the USSPACECOM ASW

processing system often reports Psq as the estimate of conjunction’s collision probability. (An 

optional, non-default mode of the function PcConjPlaneCircle in the NASA CARA SDK 

repository implements an efficient vectorized algorithm that calculates Psq values).

It is perhaps helpful at this point to review the four assumptions employed by the conjunction 

plane Pc calculation methods given by equations (N-1) through (N-4), because alternative Pc 

estimation techniques will be needed when these assumptions do not inhere:

1. Statistical Independence: The two objects’ uncertainty distributions are statistically 

independent so that the joint covariance can be obtained by simple addition of the two 

covariances. This assumption is largely true but can break down for objects that share global 

atmospheric density forecast error in a manner that influences the conjunction. This issue will 

be discussed as an isolated topic in a subsequent section of this appendix.
Appendix N Pc Calculation Approaches 7

2. Rectilinear Motion: The conjunction circumstances are such that it is reasonable to presume 

rectilinear motion and static covariances during the encounter. These conditions inhere for 

most conjunctions between Earth-orbiting satellites. The objects’ relative velocity at TCA is 

in some places used as an indication of the reasonability of these assumptions, but this 

parameter alone is not sufficient to identify situations in which the conjunction plane Pc

approximation will miscarry.

3. Gaussian Position Distributions and Negligible Velocity Uncertainties: The position vector 

state errors for each satellite at TCA follow Gaussian distributions, and the velocity vector 

state errors are negligibly small. When combined, these assumptions lead to the conjunction 

plane representations shown in Figure N-4. Conjunctions that do not satisfy these 

assumptions are addressed by alternate analytical and Monte Carlo methods, discussed in the 

sections below.

4. Temporally Isolated Event: The conjunction presents a single, well-defined event so that the 

collision likelihood can be ascertained by examining that single instance. Objects that stay in 

close proximity for extended periods accumulate risk throughout long interactions, rather 

than quickly accumulating risk at or near a well-defined TCA. A different approach is also 

required for their evaluation, which is discussed in the sections below.

As mentioned previously, the issue of statistical independence (i.e., assumption 1 above) will be 

discussed in a subsequent section of this appendix. To address conjunctions that do not satisfy 

any of the other assumptions, two alternative (but more computationally intensive) analytical 

methods are available: the “three-dimensional Nc” method, which research indicates can be 

applied to conjunctions that violate assumptions 2, 3 and/or 4 above, and the “two-dimensional 

Nc” method, applicable to temporally isolated conjunctions that violate assumptions 2 and/or 3.

Notably, the three-dimensional Nc method could, in principal, be applied to all conjunctions to

estimate Pc values. However, this is not justified because most LEO satellite conjunction Pc 

values can be estimated accurately using the much more efficient conjunction plane methods 

described above. Additionally, most of the remaining conjunctions can be estimated accurately 

using the two-dimensional Nc method, which, although much slower than the conjunction plane 

methods, is still significantly faster than the full three-dimensional Nc method. (The function 

PcConjPlaneUsageViolation soon to be posted on the NASA CARA Software Development Kit

(SDK) repository provides an algorithm that determines if a given conjunction violates one or 

more of conjunction plane Pc estimation assumptions, and if so, which of the other available Pc 

estimation methods are appropriate.)

N.3 Three-Dimensional Nc Method Analytic Pc Calculations

The relatively infrequent conjunctions that do not satisfy the conjunction plane method

assumptions discussed in the previous section must be addressed with a different methodology,

and in response to this need, several authors have formulated semi-analytical approaches 

relaxing some or all of these assumptions. Coppola (2012) proposed a method for single 

encounters designed to account for non-linear orbital motion and velocity uncertainties, resulting 

in an approximation for the probability rate, Ṗ

c

(t), calculated using integration over the surface 

of a unit sphere. When combined with a one-dimensional time integration, this yields a “three-
Appendix N Pc Calculation Approaches 8

dimensional Pc” approximation. Chan (2015) contested Coppola’s formulation, arguing that a 

proper approach must employ a set of random variables associated with a time-invariant 

Probability Density Function (PDF). NASA CARA implemented the three-dimensional Pc

method in software (Hall et al. 2017a) and subsequently discovered that, for some conjunctions, 

it can produce Pc estimates that differ significantly from high-fidelity Monte Carlo Pc

calculations, even though all the required three-dimensional Pc assumptions are satisfied.

Shelton and Junkins (2019) provided a key insight into why the original three-dimensional Pc

approximation fails in certain situations. Their analysis indicates that accurate Pc approximations 

require that the state uncertainty PDFs of the two satellites be estimated accurately in the volume 

of space where they overlap the most. The original Coppola (2012) three-dimensional Pc

formulation did not incorporate this concept, but Hall (2021) reformulated the method to do so 

explicitly. For single encounters, the reformulated approach approximates the curvilinear motion 

of each satellite using a first-order Taylor series expansion, not centered on the mean orbital 

state, but centered instead on a state that coincides with the maximum overlap of the PDFs for

the two satellites. The analysis demonstrates that such “peak overlap” states can be determined 

using an iterative calculation that converges quickly. The formulation derives an expression for 

“Nc” ⎯ the statistically expected number of collisions ⎯ which equals Pc for single, temporally 

isolated conjunctions, but that may exceed Pc for multi-encounter interactions. The resulting 

“three-dimensional Nc” method entails a total of three numerical integrations, one over time and 

two over the surface of a sphere. The outermost, time integration expression has the form

Nc

(τa, τb) = ∫ Ṅ

c

(t) dt

τb

τa

(N-5)

This expression estimates Nc

(τa, τb), the number of collisions statistically expected to occur at 

some time during the risk assessment interval τa ≤ t < τb, which can represent either a short 

duration closely bracketing a single close-approach encounter, or an extended duration spanning 

multiple encounters. The expected collision number is closely related to the collision probability. 

In fact, they are equal for single, isolated encounters between well-tracked satellites. The

collision rate for such a temporally isolated conjunction is expressed as a two-dimensional 

integral over the unit sphere

Ṅ

c

(t) = R

2 ∫ ∫ [νt(r̂) MVN(Rr̂, r̆t

,Ãt)] sin(θ) dθ dφ

π

0

2π

0

(N-6)

In this equation, the radial unit vector is given by r̂ = [cos(φ) sin(θ), sin(φ) sin(θ), cos(θ)]

T

, 

so the surface of the unit sphere is spanned by the azimuthal angle 0 ≤ φ < 2π and the axial 

angle 0 ≤ θ ≤ π. The leading factor of R

2

represents the square of the combined hard-body 

radius, meaning that the expression actually represents an integration over the surface area of the 

collision sphere. As explained in Hall (2021), the integrand function in the square brackets is the 

product of two factors. The first in an average velocity term, νt(r̂), which is a function of time

(as indicated by the t subscript) and calculated using the TCA states and covariances of the 

primary and secondary objects. The second factor, MVN(Rr̂, r̆t,Ãt), represents a multi-variate 
Appendix N Pc Calculation Approaches 9

normal (MVN) function which depends on r̂ and R, as well as a 31 relative position/velocity 

state vector, r̆t, and an associated 33 Ãt covariance matrix, both of which are also calculated 

from the TCA states and covariances. Lebedev quadrature (Lebedev and Laikov, 1999) provides 

an efficient method for numerical integration over the unit sphere. (The function Pc3D_Hall of 

the NASA CARA SDK repository computes conjunction Pc estimates, along with associated Nc 

and Nc rate estimates, using the Hall (2021) three-dimensional Nc method as summarized in 

equations (N-5) and (N-6) above.)

Figure N-5 Pc Rate and Cumulative Pc for a Non-Rectilinear Conjunction

Figure N-5 shows the time-dependent collision rate (left) and the cumulative Pc (right) 

calculated using the three-dimensional Nc method, as applied to an archived CARA conjunction 

that has a high relative velocity (~13 km/s), but that fails to satisfy the rectilinear motion and 

Gaussian PDF assumptions. This specific event does not satisfy these assumptions because it 

involves an object in a highly eccentric orbit with a conjunction occurring in the relatively tightly 

curved part of the trajectory near perigee, a phenomenology discussed by Hall (2018). In this 

case, the three-dimensional Nc method calculates a Pc value which accurately matches the 

Monte Carlo from TCA method Pc estimate, but that is a factor of fifteen larger than the 

erroneous conjunction plane Pc estimate, calculated in this case using the two-dimensional 

Foster and Estes (1992) method. It is of interest that the peak point of risk accumulation occurs 

nine seconds before the TCA, with essentially all the risk accumulated by 8 seconds before TCA. 

Such a result can seem counterintuitive at first, for one would initially expect that the point of 

highest risk would always be at the TCA. However, what is operative is the alignment between 

the geometric miss distance and the covariance. If at TCA very little of the position uncertainty 

lies along the relative miss vector, then that miss vector is a strong statement of the actual miss; 

and if the vector is somewhat larger than the hard-body radius (HBR), then the collision risk at 

that point is quite low. If, however, somewhat earlier or later than TCA the combined covariance 

(which, it must be remembered, is constantly changing position) does more substantially align 

along the miss vector, then more possible instantiations of the true miss are likely to be smaller 

than the HBR (even with the nominal miss vector larger than the expected miss at TCA), so the 

risk at that point is actually higher. In the more extreme cases, such as that shown in the above 

example (and even more strongly in the one below), most of the risk accumulation can be 
Appendix N Pc Calculation Approaches 10

relatively far from TCA, meaning that methods that examine the situation only at TCA can 

misrepresent this risk.

Figure N-6 Pc Rate and Cumulative Pc for a Low Velocity Conjunction

Figure N-6 shows the collision rate (left) and the cumulative Pc (right) for low relative velocity 

conjunction (~3 m/s), that fails to satisfy the rectilinear motion, temporally isolated, and 

Gaussian PDF assumptions. In this case, the three-dimensional Nc algorithm estimates a Pc 

value which matches the Monte Carlo from TCA method Pc estimate, but that is a factor of 330 

larger than the erroneous conjunction plane Pc estimate.

N.4 Two-dimensional Nc Method Analytic Pc Calculations

As formulated in equations (N-5) and (N-6), the three-dimensional Nc method requires the 

computation of a time-series of many unit-sphere integrations, each calculated using Lebedev 

quadrature. However, as described by Hall et al. (2023), for temporally isolated conjunctions 

(i.e., those that are not too “extended” or “blended” in time as described by Hall 2021), the 

integration over time can be approximated analytically, ultimately yielding a single unit-sphere 

integral. This relatively efficient “two-dimensional Nc” (or 2D-Nc) method yields accurate Pc 

estimates for high velocity, temporally isolated conjunctions that are affected significantly by

curvilinear trajectory effects (such as the example shown in Figure N-5). However, 2D-Nc is not 

applicable to some low velocity events that are too temporally extended or blended (as shown in 

Figure N-6); these kind of conjunctions still require the use of the three-dimensional Nc method 

or Monte Carlo estimation. The function Pc2D_Hall of the NASA CARA SDK repository 

computes Pc estimates for temporally isolated conjunctions using the two-dimensional Nc 

method.

Hall et al. (2023) describe a multistep algorithm that uses the semi-analytical 2D-Pc, 2D-Nc, and 

3D-Nc methods along with the from-TCA and from-epoch Monte Carlo methods, in order to 

calculate accurate Pc values regardless whether they are affected by curvilinear trajectory effects 

or by low-velocity encounters. The multistep algorithm sequentially applies these five 

increasingly accurate Pc estimation methods, but only as required to ensure computational 

efficiency, by conservatively evaluating potential usage violations for each.
Appendix N Pc Calculation Approaches 11

N.5 Comparison of Pc Estimates for Temporally Isolated Conjunctions

Figure N-7 shows a comparison of collision probabilities for 63,603 temporally isolated 

conjunctions extracted from the NASA CARA database for the period 2017-05-01 and 

2019-08-15 and for events with 2D Pc > 10-7(Hall, 2021). The vertical axes on all three panels 

plot Monte Carlo (MC) estimates for the collision probability⎯specifically, Monte Carlo from

TCA method Pc estimates, which are also referred to in this figure as two-body Monte Carlo 

method Pc estimates (i.e., TBMC-Pc estimates, as described in more detail in section N.13). The

error bars in Figure N-7 show 95% confidence intervals estimated using the Clopper-Pearson 

method (1934); several of the conjunctions had zero hits registered in the Monte Carlo

simulations, which are represented in Figure N-7 using downward-pointing triangles and a 

single-sided error bar. The horizontal axes plot the corresponding semi-analytical 

approximations: two-dimensional Pc on the left graph, three-dimensional Nc in the center, and 

two-dimensional Nc on the right. The colored points on each plot indicate the results of a 

binomial proportion statistical test that evaluates the agreement between the estimates. 

Specifically, black points in Figure N-7 indicate analytical Pc estimates that agree reasonably 

well with the Monte Carlo estimates as they do not violate a null-hypothesis that the two are 

equal at a p-value  10-3significance level. However, those highlighted in yellow do violate the 

hypothesis at this significance level, and those in red at a more stringent level of p-value  10-6. 

Overall, the two-dimensional conjunction plane Pc comparison plotted in the left graph contains 

254 yellow and 436 red points, which both significantly exceed the number of disagreements 

expected from purely statistical variations, even though together they represent a small fraction 

(~1%) of the original conjunctions. These disagreements represent conjunctions that violate one 

or more of the assumptions required for conjunction plane Pc estimation. The center graph

clearly shows that the three-dimensional Nc method matches the Monte Carlo Pc estimates 

better, producing only 66 yellow and zero red points. Finally, the two-dimensional Nc method

plotted on the right produces very nearly the same results as the three-dimensional Nc method

but requires much less computation time. The three comparisons shown in Figure N-7 indicate

that, for temporally isolated conjunctions, the two- and three-dimensional Nc methods are 

consistent with one another and match Monte Carlo Pc estimates significantly better than the 

conjunction plane Pc estimation method.
Appendix N Pc Calculation Approaches 12

Figure N-7 Comparison of Monte Carlo Collision Probabilities with the

Two-Dimensional Pc Method (left), the Three-Dimensional Nc

Method (center), and the Two-Dimensional Nc Method (right) for a Large 

Set of CARA Conjunctions

N.6 Pc Calculations for Multi-encounter Interactions

An additional feature of the two- and three-dimensional Nc methods is their ability to address 

explicitly the amalgamated risk of repeating conjunctions. While most conjunctions are 

temporally isolated, there are two conjunction types that exhibit different behaviors. The first 

type is produced by objects in orbits in synodic alignment that generate a temporal sequence of

conjunctions once every revolution or multiple of a revolution (e.g., two nearly circular orbits at 

different inclinations that produce conjunctions at one or both nodal crossing points). The second 

type involves objects that orbit close to each other for extended periods, generating extended 

interactions with multiple close approaches (e.g., two satellites in nearly the same orbit, but with 

only a slight difference in inclination and/or eccentricity). In both cases, if each of the multiple 

encounters is considered separately, situations can arise in which each encounter in the series 

falls below a Pc mitigation threshold, but the combined risk of all of the encounters exceeds that 

threshold. The two- and three-dimensional Nc methods account for such multi-encounter 

interactions by providing estimates for the total expected number of collisions, and upper and 

lower bounds for the probability of collision (Hall, 2021; Hall et al., 2023). The total statistically 

expected number of collisions for a multi-encounter interaction is given by a summation of Nc 

values for each conjunction in the sequence

Nc

(τa, τb) = ∑Nc,k

K

k=1

(N-7)

with the index k = 1...K indicating the close approaches applicable to the risk assessment 

interval τa ≤ t < τb, and Nc,k representing the associated expected number of collisions for 

each. The upper and lower Pc bounds for the combined interaction are given by

max(Nc,k) ≤ Pc

(τa, τb) ≤ 1 −∏(1 − Nc,k)

K

k=1

(N-8)

which also implies that Pc

(τa, τb) ≤ Nc(τa, τb). Figure N-8 shows an example of a multiencounter interaction, in which a pair of satellites experience four conjunction events over about 

a five-hour period. Each of the individual conjunctions produces a Pc value in the upper yellow 

region (between 10

-5

and 10

-4

), as plotted in the bottom panel. The solid line in the top panel

shows the upper limit of the cumulative Pc for the interaction; notably, after the third 

conjunction, the cumulative Pc exceeds 10

-4 — a value frequently selected as a risk mitigation 

threshold. So, while these events would not necessarily prompt a mitigation action if examined 

individually, when considered collectively they do appear to represent a situation of sufficiently 

high collision likelihood to warrant mitigation. In such a case, it is advisable to run a Monte 

Carlo investigation (discussed in a subsequent section of this appendix) to verify that the upper-
Appendix N Pc Calculation Approaches 13

limit Pc value generated by the method is in fact representative of the actual cumulative risk. The 

utility of the analytic calculation, however, should be clear: if the upper-bound calculation for 

repeating events is found to be below the mitigation threshold, then there is no need to marshal 

computationally intensive methods (such as Monte Carlo), for it has already been demonstrated

no mitigation action is warranted.

Figure N-8 Cumulative Three-Dimensional Nc Risk Over 

Four Repeating Conjunctions

N.7 Input Covariance Data Considerations

Most calculations are only as good as the input data that drive them, and Pc calculation is no 

exception. Appendix P discusses orbit determination quality and propagation issues for 

individual objects and addresses the question of the circumstances under which the state estimate 

and covariance might be considered sufficiently poor so as not to constitute a basis for 

conjunction risk mitigation actions. The purpose of this section is to address the routine 

improvement and basic error checking extended to covariances as part of the Pc calculation. 

These activities fall into three basic groups: correction for known problems in propagation, 

covariance formation and stability issues, and correlation between primary and secondary 

covariances. Each of these topics will be treated in turn.
Appendix N Pc Calculation Approaches 14

N.7.1 Propagation Error Compensation

Historically, accuracy analysis of state estimate and uncertainty estimates has focused on 

products that emerge directly from the orbit determination fit. These direct orbit-determination

products include best-estimate states and covariances of the satellite at an epoch that usually 

coincides with the acquisition time of the last tracking metric observation incorporated into the 

analysis. While these direct orbit-determination products are of course foundational, it is 

important to remember that most SSA activities, and conjunction assessment in particular, are 

not usually conducted with these direct products but rather with predictions propagated from the 

orbit-determination epoch solution, often over a non-trivial duration that spans many orbital 

revolutions into the future. The batch orbit-determination analysis method used by the DOD 

produces a formation covariance that represents the expected at-epoch state uncertainty based on 

the number, quantity, and temporal spacing of the incorporated metric observations; when the 

state is propagated forward, a parallel process can also be used to propagate the covariance 

forward in time. The same dynamical models used for the orbit-determination analysis as well as 

the state propagation itself are used to perform this covariance propagation, although in a 

linearized way. This means that the propagated covariance will be sized (mostly) appropriately 

for both the propagation duration and the final prediction point in the orbit.

Despite the use of appropriate models to propagate the covariance forward in time, a number of 

additional sources of error manifest themselves during the propagation interval yet are not part of 

the dynamical model used during the fit; these errors are therefore neither incorporated into the 

orbit-determination-produced covariance nor added as part of the regular propagation process. 

Because of the prevalence of such outside-of-model errors, techniques have been developed to 

account for them, the most familiar of which is the addition of process noise during propagation. 

Originally developed to account for acceleration errors that, due to model inadequacy, were to 

some degree known, this method propagates a noise matrix alongside the propagated covariance 

and combines both matrices as part of the overall process. The result is a covariance that is larger 

than it would have been otherwise to account for these (characterized) errors in the force 

model(s). A second approach, which is the one used by DOD in the propagation of their orbit 

prediction products, applies parameters to the covariance before propagation to guide the 

propagation process in producing a more realistic result. Because this is the approach reflected in 

the CDM covariances that conjunction assessment practitioners receive from the DOD, it will be 

described here in some detail. 

Orbit determination makes a distinction between “solved-for” parameters that are actually 

estimated during an orbit-determination activity, and “consider” parameters that are not “solved 

for” but represent a priori information that is “considered” as part of the orbit-determination

process. In the present case, the use of the term “consider parameter” is somewhat non-nominal 

in that it is referring not to additions or alterations made during the fit but to modifications to the 

fit-produced covariance so that when it is propagated, it may give a more realistic representation 

of the expected state errors. For DOD covariances, two different consider parameters are applied 

to compensate for two distinct expected errors during propagation: atmospheric density forecast 

error and satellite frontal area uncertainty.
Appendix N Pc Calculation Approaches 15

Atmospheric drag is a significant force that affects satellite orbits with perigee heights less than 

1000 km, and the calculation of the expected atmospheric drag at any particular moment requires 

an estimate of the neutral atmospheric density that that satellite will encounter at that moment. 

Because the atmospheric models that generate this estimate are driven by space weather indices, 

the ability to predict these indices accurately affects the fidelity of the predicted atmospheric 

density and thus the atmospheric drag. Unfortunately, it is difficult to predict future space 

weather indices well, primarily because they are affected by activities on the part of the sun’s 

surface that has not yet rotated into view from the Earth. This particular issue was studied with 

the Jacchia-Bowman High Accuracy Satellite Drag Model (HASDM) 2009, which is the 

atmospheric density model used by DOD, by comparing predicted densities (using predicted 

space weather indices) to actual densities and producing polynomial fits of the relative density 

error as a function of satellite perigee height and the Ap (major magnetic storms list) and Dst 

(disturbance storm-time) space weather indices. Figure N-9 shows the behavior of these 

polynomial fits divided into four different classes of Ap/Dst activity; y-axis values are omitted 

here to allow full releasability of the figure:

Figure N-9 Behavior of Relative Density Error by Perigee Height and

Solar Activity

These fits produce a variance term that can be added to the ballistic coefficient variance in the 

covariance: because in the drag equation the ballistic coefficient and the atmospheric density 

estimate are multiplicatively coupled, changing one of these parameters has the same effect as 

changing the other. When the covariance is propagated, this augmented term will appropriately 

expand the other covariance terms.

The amount of drag acceleration a satellite encounters is also governed by the frontal area that 

the satellite presents to the atmosphere; this makes intuitive sense (amount of resistance is a 

function of area presented to the resisting fluid) and is reflected in the ballistic coefficient (B) 

equation
Appendix N Pc Calculation Approaches 16

M

A

B = CD

(N-9)

in which CD is the (dimensionless) drag coefficient, which indicates the degree to which the 

satellite surface is susceptible to drag resistance; M is the satellite mass; and A is the satellite 

frontal area. As A increases, B increases as well, increasing the overall drag acceleration value. 

Stabilized satellites should manifest a stable B term, but rotating satellites, because their frontal 

area term is continuously changing, can exhibit a range of B terms. Three outcomes are possible 

depending on the rapidity of the rotation: the effect can be washed out during the fit (because the 

rotation is so rapid that an average value is quite representative), the effect can be not relevant 

during the fit (because the rotation is so slow that it does not affect near-term look-back and 

look-forward), or the effect can be such that the rotation does affect results fit-to-fit. It is this last 

case for which compensation is helpful. A history of regularized B histories for individual 

satellites is examined and a relative error and variance calculated, and this variance is added to 

the drag variance in the covariance as a corrective term whose influence will be realized in 

propagation.

There is some additional subtlety regarding the exact way these consider parameters are applied. 

A typical propagation consists of two conceptual stages: the first stage is the propagation forward 

from the epoch time of the orbit determination to the present time, which can make use of 

measured and thus highly reliable space weather indices; and the second stage is from the present 

time to the desired final propagation time, which has to use predicted space weather indices and 

the errors that these introduce. The two consider parameters are thus applied at different times. 

Because satellite rotation and its resultant uncertainty will occur for the entire interval from 

epoch time to the propagation end point, that consider parameter is applied at epoch. 

Atmospheric density forecast error, however, is encountered only forward in time from the 

present moment, so it is added only for that portion of the propagation. Figure N-10 outlines the 

two-phase application of these consider parameters:

DCP = Dynamic Consider Parameter

Figure N-10 Two-phase Application of Consider Parameters

If the CDMs generated by the DOD are used for conjunction assessment applications, the good 

news is that all the consider parameter activities described above are already performed—the 

propagated covariances that the CDM contains have had these two consider parameters applied 

during the covariance propagation executed at the DOD work center. If one is working with 
Appendix N Pc Calculation Approaches 17

epoch state estimates, which is sometimes necessary when performing Pc calculations with 

Monte Carlo techniques, then manual application of the consider parameters may be necessary. 

This issue is discussed at greater length in the section that addresses Monte Carlo Pc

calculations.

N.7.2 Defective Covariances

There are several ways in which a covariance contained in a CDM can be defective. 

• Null Covariances. All-zero, or null, covariances are sometimes observed, usually arising 

from a conjunction assessment screening for which the O/O-provided ephemeris does not 

contain covariance data. In such a case, it is possible to compute the Pc either using only 

the one covariance that the CDM message contains or by applying a special technique 

that determines the maximum Pc possible presuming that the null covariance could take 

on any possible value (developed and described in Frisbee 2015). 

• Default Covariances. Default covariances are diagonal covariances that contain a value 

of ten earth radii squared for each of the three position variances. The presence of this 

covariance indicates that a precision, special-perturbation solution for the object was not 

possible; the state estimate provided arose from a general-perturbation solution, and an 

orbit-determination-produced covariance was not available. Such a result is not a 

precision solution and does not constitute a basis for conjunction risk mitigation actions.

• Non-Positive-Semidefinite Covariances. Another defective covariance type found in 

CDMs, now quite rare due to improvements to the DOD operational system, is a 

covariance that fails to be positive semidefinite. A positive semidefinite matrix is one that 

contains no negative eigenvalues. Because the covariance represents a hyperellipsoid of 

actual state error information, it must have a set of eigenvalues all greater than or equal to 

zero for error information to consist of real rather than imaginary quantities. The orbitdetermination mechanism that generates the covariance should always produce at least a 

positive semidefinite matrix, for the linear algebra function involves the product of a 

square matrix and its transpose, and one can prove that this procedure always produces a 

positive semidefinite result. Due to either numerical precision limits or interpolation, the 

provided matrix is sometimes not positive semidefinite. If the 2 × 2 projection of the 

covariance matrix into the conjunction plane is not positive semidefinite, the twodimensional Pc calculation is not possible. If the full 6 × 6 or 8 × 8 matrix is not positive 

semidefinite, then Monte Carlo sampling on the entire matrix is not possible either. 

As such, some attention must be paid to this issue of positive semidefinite matrix conditioning. A 

recent paper on this subject (Hall 2017b) examined the issue in some detail and compared 

different matrix repair algorithms to minimally adjust the covariance to make it positive 

semidefinite compliant; it found that most repair approaches yield equivalent answers in terms of 

the resultant calculated Pc. An “eigenvalue clipping” procedure was developed in which any 

negative eigenvalues (which are almost always small) are set to a small positive or zero value, as 

required. 
Appendix N Pc Calculation Approaches 18

The CARA operational implementation of this method proceeds parsimoniously, namely by 

directing such repair only to the level needed to perform a calculation. For example, a covariance 

used for the two-dimensional Pc calculation would neither be tested for positive semidefinite 

compliance in its full 8 x 8 form nor its position-only 3 x 3 form; instead, the 2 x 2 projection of 

the joint covariance into the conjunction plane would be tested and repaired only when necessary 

to enable the two-dimensional Pc calculation. To do otherwise is to make repairs and potentially 

alter the result when this is not strictly necessary. 

The Pc Computation CARA SDK includes the source code for identifying the need for and 

making the covariance matrix repairs described above; it is available in the NASA CARA 

software repository. (See Section 7, Contact Information in this document for the specific URL.)

N.7.3 Covariance Correlation

For nearly all the broader conjunction assessment conduct during the past decade, practitioners 

operated with the presumption that the primary and secondary objects’ covariances could be 

considered uncorrelated. Not only was this the “right” answer in that it greatly simplified the Pc

computation because the joint covariance could be formed by simple addition of the two 

covariances, but there was also a good intuitive justification for the presumption. Because the 

focus had been on the two objects’ orbit-determination fits, which are based on separate sets of 

observations, there was no expectation that there would exist any significant correlation between 

the two objects’ covariances. The principal source of potentially correlated error was presumed 

to be uncharacterized but correlated errors in space sensor observations used by both primary and 

secondary objects. Because most primaries receive many observations from many different 

sensors, it was seen as unlikely that this particular source would introduce much correlation. 

Correlation between covariances was thus expected to be small, and conjunction assessment

operators proceeded as though it were.

With the initiative several years ago to include outside-of-fit prediction error characterization 

into DOD satellite covariances (see the above section on Propagation Error Compensation), the 

issue of covariance cross-correlation began to be rethought. The principal outside-of-fit 

prediction error is global atmospheric density forecast error due to inadequate space weather 

index prediction. Because this is a global error, it is likely to be shared among large classes of 

objects, some of which might constitute both the primary and secondary objects in a conjunction. 

As discussed previously, this global density forecast error has been parameterized by satellite 

perigee height and predicted geomagnetic index, so the degree of such error, both identified 

separately and injected into each DOD-furnished covariance by means of a consider parameter, 

is known for each satellite. It is possible to determine the degree of shared error from this source 

and account for it when forming the joint covariance.

A study by Casali et al. (2018) provides a full development of this theory and presents results for 

an evaluation set of conjunction data. Essentially, one has to determine the global density 

forecast error relevant to each satellite and the degree to which the associated drag error induced 

by this density forecast error will manifest itself in position error relevant to the particular 

conjunction. The governing equation is the following:

Pm = Ps + Pp - s/g p/g Gs G

T

p - s/g p/g Gp G

T

s (N-10)
Appendix N Pc Calculation Approaches 19

in which Pm is the decorrelated joint covariance at TCA, Ps is the secondary covariance at TCA, 

Pp is the primary covariance at TCA, σs/g is the density forecast error germane to the secondary 

satellite, σp/g is the density forecast error germane to the primary satellite, Gs is the sensitivity 

vector mapping drag acceleration error to secondary satellite position error at TCA, and Gp is the

sensitivity vector mapping drag acceleration error to primary satellite position error at TCA. One 

could wonder how the conjunction assessment practitioner would come upon all of the needed 

data to effect the proposed compensation. As described in detail in the next section, a recent 

enhancement to the DOD operational system has placed all of this information in the CDM itself, 

allowing the direct calculation of the decorrelated joint covariance. The CARA Pc Calculation 

SDK, available on the NASA CARA software repository, also contains both a math specification 

outlining this calculation and source code to perform it. (See Section 7, Contact Information in 

this document for the URL.) Hall (2021) describes how the density forecast errors and sensitivity 

vectors can be used to estimate decorrelated joint covariances for the two- and three-dimensional 

Nc methods.

A heuristic probing of the situation reveals that, conjunction by conjunction, different levels of 

Pc changes are possible due to cross-correlation remediation. Orbits that are insensitive to 

atmospheric drag are little affected by density forecast error and will have Pc estimates that, as 

expected, also change little. Head-on conjunctions are expected to be left mostly unaffected as 

well, for the components of the error that govern the Pc are not subject to density forecast error 

perturbation. Crossing events are perhaps the most susceptible to cross-correlation effects, 

especially if the drag level of both satellites is similar.

The plot in Figure N-11 profiles 250 conjunctions in which the primary and secondary satellites 

are of non-negligible drag (i.e., Energy Dissipation Rate (EDR) values greater than 0.0006 

watts/kg; see Hejduk 2008 for an explanation of energy dissipation rate) and plots the ratio of the 

Pc calculated with the decorrelated joint covariance to that of the Pc calculated with the 

unmodified joint covariance. One can see that for somewhat more than half of the events, the 

ratio hovers near unity, meaning that the conjunction is little affected by the compensation. For 

about one-third of the cases, the decrease in Pc is notable, in many instances more than an order 

of magnitude. For the remaining cases, there is an increase in Pc from a factor of 1.5 to 5. 
Appendix N Pc Calculation Approaches 20

Figure N-11 Profiles of 250 Conjunctions with Primary and Secondary 

Satellites of Non-negligible Drag2

While the Pc value changes by less than a factor of 1.5 for most of the conjunctions, a sufficient 

number are affected more substantially therefore justifying the integration of this additional 

consideration into the Pc computation, especially because it is a straightforward calculation from 

data provided directly in the CDM.

N.7.4 Conjunction Data Message DCP Uncertainty and Sensitivity 

Vector Explanation

The content of Conjunction Data Message files produced by the United States Space Force’s 

ASW software has been changed to enable conjunction risk assessments that account for 

covariance correlation. Specifically, ASW release versions 19.2 and after insert additional 

information within CDMs to enable users to include the atmospheric drag covariance crosscorrelation effect when estimating a conjunction’s probability of collision, as calculated from 

DCP uncertainty values and sensitivity vectors (as described in the previous section). This 

section provides an overview of these CDM modifications, and a brief description of their usage.

CDM Mean Position/Velocity State Vectors

CDMs can contain the mean position/velocity state vectors of the primary and secondary 

satellites at TCA in a couple of reference frames, as described in more detail in the CDM 

specification documents by Consultative Committee for Space Data Systems (specifically, 

CCSDS 502.0-B-1 and 502.0-B-2, listed in the references section). For Pc computation, an 

inertial reference frame must be used, sometimes requiring these original CDM states to be 

converted into the inertial reference frame. This section denotes the resulting mean inertial 

position and velocity vectors at TCA for the primary and secondary objects as (rp̅, v̅p) and 

( r̅s, v̅s), respectively.

2 From Casali et al. 2018.
Appendix N Pc Calculation Approaches 21

CDM Position/Velocity State Covariance Matrices

CDMs specify position/velocity state covariance matrices for the primary and secondary 

satellites using the radial-transverse-normal (RTN) coordinate frame. A CDM file specifies this 

symmetric 66 covariance for each object using keyword values for its 21 non-redundant matrix 

elements as follows

C =

[

 

 

 

 

 

 

CR,R CT,R CN,R CṘ

,R CṪ,R CṄ,R

CT,R CT,T CN,T CṘ

,T CṪ,T CṄ,T

CN,R CN,T CN,N CṘ

,N CṪ,N CṄ,N

CṘ

,R CṘ,T CṘ,N CṘ,Ṙ CṪ,Ṙ CṄ,Ṙ

CṪ

,R CṪ,T CṪ,N CṪ,Ṙ CṪ,Ṫ CṄ,Ṫ

CṄ

,R CṄ,T CṄ,N CṄ,Ṙ CṄ,Ṫ CṄ,Ṅ]

 

 

 

 

 

 

 (N-11)

The RTN-frame covariance C can be transformed into a 66 inertial frame position/velocity 

covariance P by applying the following equation

P = M C MT (N-12)

with the 66 transformation matrix M having the form

M = [

M 03x3

03x3 M

] (N-13)

with 03x3 representing a 33 matrix of zeros. The 33 orthonormal matrix M rotates vectors 

from the pseudo-inertial (non-rotating, instantaneously frozen) RTN frame into the inertial frame

M = [R̂ T̂ N̂] (N-14)

with column vectors given by the three RTN unit vectors, calculable from the object’s inertial 

mean position and velocity vectors as follows

R̂ = r̅⁄|r̅| and N̂ = (r̅× v̅)⁄|r̅× v̅| and T̂ = N̂ × R̂ (N-15)

The 66 inertial frame covariance in eq. (N-12) can each be decomposed into three 33 submatrices 

P = [

A B

T

B C

] (N-16)

with A representing the marginalized covariance of the position vector, C the marginalized 

covariance of the velocity vector, and B position-velocity cross correlations.

Eqs. (N-11) to (N-16) can be used to calculate inertial frame position/velocity state covariance 

matrices at TCA for the primary and secondary objects involved in a conjunction, Pp and Ps, 

respectively, as well as the marginalized position covariance matrices, Ap and As. (Note: 

because RTN is an object-specific frame of reference, these calculations must employ different 

rotation matrices, Mp and Ms, respectively.)

Uncorrelated and Correlated Relative Position Covariance Matrices
Appendix N Pc Calculation Approaches 22

Collision probability estimation using the 2D-Pc method requires the conjunction’s inertial 

relative position miss-vector, rm̅ = r̅s−rp̅, along with the associated miss-vector covariance 

matrix, Am. If the primary and secondary position vectors are statistically independent (i.e., 

uncorrelated), then the relative position miss-vector covariance is given by their sum

Am = Ap + As (N-17)

This approach provides a viable approximation for 2D-Pc estimation for conjunctions in which 

Ap and As have sufficiently weak correlation. However, as described in the previous section, 

Casali et al. (2018) have demonstrated that some conjunctions have stronger covariance 

correlations, due to shared atmospheric density forecast components arising from the ASW

global density portion of the Dynamic Consider Parameter. In these cases, the miss-vector 

covariance can be corrected by removing the cross-correlated components as follows

Am = Ap + As − σp/gσs/g[GpGs

T + GsGpT

] (N-18)

with σp/g and σs/g denoting the atmospheric density 1-sigma relative uncertainties for the 

primary and secondary, respectively. The vectors Gs and Gp represent the sensitivity of the 

miss-vector covariance on global density relative uncertainties. (Note: see eq. (11) of Casali et 

al. (2018) and the related discussion for a derivation of eq. (N-18), and a more detailed 

explanation of its components. Also, this section uses the symbol “A” for 33 position 

covariances, instead of the symbol “P” which this section uses to represent full 66 

position/velocity covariances. This change of notation means that eq. (N-18) is simply a 

restatement of eq. (N-10), the latter of which uses the original Casali et al. (2018) notation.)

CDM files produced by ASW versions 19.2 and later provide the sigma values and sensitivity 

vectors required to calculate corrected relative position covariance matrices using eq. (N-18). 

These data have been added within “comment” lines in the primary and secondary object 

portions of the CDM, with the format

COMMENT DCP Density Forecast Uncertainty=2.45030416E-01

COMMENT DCP Sensitivity Vector RTN Pos=-1.67441647E+01 3.68889831E+02 1.63797508E-01 [m]

COMMENT DCP Sensitivity Vector RTN Vel=-3.98670591E-01 1.09452965E-02 -4.83454839E-04 [m/sec]

The first occurrence of parameters with this specific format appears in the primary object section 

of the CDM and provides σp/g (the DCP density forecast 1-sigma uncertainty), Gp

RTN (the 31 

DCP position sensitivity vector, expressed in the primary’s RTN frame), and Hp

RTN (the 31 

DCP velocity RTN-frame sensitivity vector). The second occurrence provides the corresponding 

DCPs for the secondary, i.e., σs/g, Gs

RTN and HsRTN. The RTN frame sensitivity vectors provided 

in the CDM can be converted to inertial frame vectors using the transformation matrix defined in 

eq. (N-14) separately for each object:

Gp = MpGp

RTN and Gs = MsGsRTN and Hp = MpHpRTN and Hs = MsHsRTN (N-19)

Uncorrelated and Correlated Relative Position/Velocity Covariance Matrices
Appendix N Pc Calculation Approaches 23

For statistically independent primary and secondary states, the relative position/velocity missstate covariance is given by the sum of the covariances for the two objects

Pm = Pp + Ps (N-20)

This approach provides a viable approximation for conjunctions in which the primary and 

secondary position/velocity states have sufficiently weak correlation. In other cases, the missvector covariance can be corrected by removing the cross-correlated components as follows

Pm = Pp + Ps − σp/gσs/g[ΓpΓs

T + ΓsΓpT

] (N-21)

with σp/g and σs/g again denoting the atmospheric density 1-sigma relative uncertainties for the 

primary and secondary objects, respectively. The 61 vector Γp = [Gp

T HpT

]

T

represents the 

sensitivity of the primary object’s miss-state covariance on global density relative uncertainties. 

The secondary object’s 61 sensitivity vector, Γs, is defined similarly. (Again, see Casali et al. 

(2018) for more detail on these DCP uncertainties and sensitivity vectors.) The RTN frame 

sensitivity vectors provided in the CDM can be converted to inertial frame vectors using the 

transformation matrix defined in eq. (N-13) for each object

Γp = MpΓp

RTN and Γs = MsΓsRTN (N-22)

Pc Estimates with and without Covariance Correlation Correction

Conjunction 2D-Pc values calculated using the miss-vector covariance in eq. (N-18) represent 

collision probabilities corrected for global atmospheric cross-correlation effects. These can 

differ from the uncorrected 2D-Pc values calculated using the covariance in eq. (N-17). As 

shown in Figure N-11, analysis of archived conjunctions indicates that this correction usually 

does not change Pc values appreciably, except in a minority of conjunctions that have both

elevated drag energy dissipation rates and an appropriate combination of orbital geometries. 

Among this minority, however, the corrections can potentially elevate Pc values by a factor of 

1.5 or more (as shown on the right side of the plot in Figure N-11), meaning that accurate and 

conservative risk assessments for these cases rely on applying the covariance cross-correlation 

corrections made possible by the ASW system’s recent CDM modifications.

For low-velocity or multi-conjunction interactions, the statistically expected number of collisions 

(i.e., the “3D-Nc” value) can be calculated using a miss-state covariance matrix of the same form 

as given in eq. (N-21), which then can be used to estimate the Pc value for the interaction, as 

explained in more detail by Hall (2021).

N.8 Monte Carlo Pc Calculation Techniques

Analytic approaches to Pc calculation are more computationally efficient than Monte Carlo 

methods, especially the conjunction plane two-dimensional Pc method. However, as discussed 

previously, analytic methods require certain enabling assumptions that are not necessarily valid 

for all conjunctions. Monte Carlo approaches require fewer enabling assumptions, but they are 

not typically employed as the first method of computation. Instead, Monte Carlo methods are 

usually reserved for cases that are suspected of violating the enabling assumptions of the analytic 
Appendix N Pc Calculation Approaches 24

methods. As described by Hall et al. (2018) and Hall (2021), there are two strains of the Monte 

Carlo method that are regularly employed:

• “Monte Carlo from epoch,” in which the orbit-determination-epoch mean states and 

covariances are used to generate sampled states, and potentially long (i.e., multiday) 

high-fidelity propagations are required for each Monte Carlo trial.

• “Monte Carlo from TCA,” in which the mean states and covariances predicted at TCA 

are used to generate sampled states, and only relatively short propagations to each trial’s 

new TCA are required.

Each of these two approaches will be discussed in turn.

N.8.1 Monte Carlo from Epoch

The principal appeal of calculating the Pc using the Monte Carlo from epoch approach is that it 

requires almost no simplifying or restrictive assumptions, making it is as close to a “gold 

standard” for Pc estimation as can be devised. The input information includes the states and 

covariances for the primary and secondary objects at their respective orbit-determination epoch 

times, the combined hard-body radius of the two objects, and an ensemble of environmental 

datasets required for the high-fidelity propagations (such as predicted space weather indices and 

atmospheric density debiasing data; see Hall et al. 2018). The most elaborate instantiation of this 

technique uses the full eight-dimensional orbital state vectors (each containing six coordinate or 

element variables, plus drag and solar radiation pressure variables) along with the associated 8x8 

covariance matrices for the primary and secondary satellites. 

For each Monte Carlo trial, a state perturbation is obtained by performing a random draw from 

the distribution using the covariance matrix to generate the associated Gaussian distribution for 

each variable. These perturbations are then used to alter the initial states by adding each to the 

appropriate epoch mean state estimates. To do this for the primary satellite, for example, one 

would generate random perturbations for the eight variables representing the primary satellite’s 

state based on their distribution as defined by the covariance, and then create a new, sampled 

state vector by adding the (signed) perturbations to the mean epoch state vector. This same 

procedure would also be performed for the secondary satellite. This sampling process generates 

epoch states for both the primary and secondary objects that represent statistically reasonable 

alternatives for their actual states. These two sampled epoch states are then propagated forward, 

a TCA identified, and a check performed to determine whether the miss distance at TCA is 

smaller than the hard-body radius; if it is, the trial results in a simulated collision and a “hit” is 

registered; if not, it is considered a “miss.” This sampling/ propagation/ hit-registration 

procedure is then repeated for a large number of Monte Carlo trials, and the final number of hits 

divided by the total number of trials constitutes the Monte Carlo Pc estimate. There are 

algorithms that can be applied to estimate the confidence interval on the Monte Carlo Pc after a 

given number of trials due to event counting uncertainties (e.g., MATLAB’s binofit function).

This procedure seems straightforward enough, and in many respects it is. But there are subtleties 

that require attention, especially if the technique is deployed for LEO conjunctions:
Appendix N Pc Calculation Approaches 25

• For the result to be valid, the same force models and force model settings must be used 

for the Monte Carlo propagations as were used when generating the original orbitdetermination solution. While it often is not difficult to apply the same general force 

model settings, there does need to be overall compatibility between the orbitdetermination engine and the Monte Carlo propagator, such as the same set of 

geopotential Legendre polynomial coefficients (not just the same geopotential order) and, 

most critically for LEO, the same atmospheric density model. While it may be possible to 

deviate somewhat from this level of compatibility and still obtain reasonably accurate 

outcomes, the “gold standard” propriety of the result is lost.

• Correlation between the primary and secondary covariances, as described in the previous 

section, should be considered. This correlation can be modeled during the random draw 

process by forcing correlation in the primary and secondary objects’ drag perturbations.

• Monte Carlo from epoch often requires extremely long computation run times. The run 

time is a function of the actual Pc, since this will determine how often hits are likely to 

occur and the number of trials required to obtain a result with a desired confidence level. 

Without marshalling extensive high-performance computing, the Pc levels that can be 

explored with this method have to remain relatively large (e.g., 10-5or above). 

Computation times for smaller Pc events can be prohibitively long. For example, 

validating a Pc estimate of ~10

-7

at the 95% confidence level for a conjunction with TCA 

five days from the orbit determination epoch time would require an estimated two years 

of execution time on a 20-CPU, reasonably modern server (Hall et al. 2018). For this 

reason, it is typically necessary to reserve the Monte Carlo from epoch method for larger 

Pc events only, which means that one must trust analytic methods to identify a candidate 

subset of conjunctions for follow-up Monte Carlo analysis.

• High-fidelity orbital state propagations require the most processing during typical Monte 

Carlo from epoch computations, so when applying this approach, there is a temptation to 

“reuse” propagations to gain computational efficiency. Suppose that ten perturbations 

were performed for the primary satellite and ten propagated ephemerides were generated, 

and the same was done for the secondary object as well. Ephemeris #1 for the primary 

could be compared to ephemeris #1 for the secondary to determine whether a hit 

occurred, ephemeris #2 for the primary to ephemeris #2 for the secondary, ephemeris #3 

for the primary to ephemeris #3 for the secondary, etc., and have as a result ten 

comparisons/trials. To go further, one could compare ephemeris #1 for the primary to 

ephemeris #2 for the secondary, ephemeris #1 for the primary to ephemeris #3 for the 

secondary, etc., and realize 100 comparisons/trials from merely 20 total propagations. 

Such a procedure certainty seems advantageous, given that processing time is the limiting 

factor to the deployment of Monte Carlo from epoch, and following such a procedure will 

produce a result that converges on the true Pc. The drawback is that such reuse of 

samples violates the conditions for the proper application of the formulae to calculate Pc 

uncertainty confidence intervals. Monte Carlo results without associated reliable 

confidence intervals are not operationally useful because it is never known how close one 

is to the true Pc value. Schilling et al. (2016) discuss this issue and confirm it to be a 
Appendix N Pc Calculation Approaches 26

problem, and they recommend some estimation techniques that allow (large) confidence 

intervals to be constructed around Monte Carlo products that avail themselves of sample 

reuse. The CARA implementation has avoided any such sample reuse to ensure that the 

“gold standard” status of the results not be in question and to produce more narrow 

confidence intervals.

Due to all the above considerations, but especially the run-time requirements, Monte Carlo from 

epoch is usually reserved only for those cases that require it for accurate Pc estimation.

A study effort discussed in greater detail in the next section determined that the Monte Carlo 

from epoch method appears to be needed only when the two objects stay in such close proximity 

that they experience a sequence of multiple, comparable-risk close approaches during a multiday 

risk assessment projection interval. For closely spaced co-orbiting objects, these conjunctions 

may also become effectively “blended” in time with one another such that collision probability 

accumulates appreciably even during the mid-point times between the encounters rather than just 

during short bursts very near the close approaches (Baars et al. 2019). In such cases, two 

impediments arise to estimating accurate Pc values using methods formulated for temporally 

isolated close approaches. First, there is no clear, single TCA at which to evaluate the collision 

likelihood. While one could in fact find the unique point of closest approach between the two 

nominal trajectories for the entire sequence and perform a two-dimensional Pc calculation for 

that encounter, there is no guarantee that another encounter in the sequence may actually possess 

a higher Pc due to different covariance orientations even though it has a larger nominal miss 

distance. Second, calculating a single encounter Pc at each of the close approaches and then 

amalgamating these using the following formula (derived from DeMorgan’s Law of 

Complements)

PTot = 1 − ∏ (1 − Pci)

n

i=1

(N-23)

potentially overestimates the overall collision likelihood because it presumes that the individual 

events in the sequence are statistically independent, but in fact they may not be, especially if 

blended in time. This total probability estimate matches that given in equation (N-8) and 

provides an upper bound for the amalgamated risk over the time interval of interest. For 

maximum efficiency, Monte Carlo from epoch would then be run operationally only in cases in 

which this upper bound is above a mitigation threshold and there is interest in determining 

whether the higher-fidelity Monte Carlo calculation would reduce this value to one much closer 

to or below the threshold. Monte Carlo from epoch can also be run if there is any question about 

the overall rectitude of the Pc calculation. As stated earlier, lower-Pc conjunctions may present 

intractable Monte Carlo execution times, but if one wishes only to ensure that the Monte Carlo 

Pc falls below the mitigation threshold (rather than establish a high-fidelity Pc value), this can 

usually be accomplished with far fewer Monte Carlo trials.

Optimal application of the Monte Carlo from epoch Pc estimation method does not entail Pc

evaluation over a relatively short time interval bracketing a single conjunction’s TCA, but rather 

over a more extended interval that spans several close approach encounters. For example, the 

collision likelihood between two objects would not be evaluated at the nominal TCA for a single 

conjunction, but, perhaps, over a risk assessment interval projecting forward seven days. This 
Appendix N Pc Calculation Approaches 27

multiday interval would not only include the original nominal TCA but also a sequence of other 

encounters between the primary and secondary as well. In this case, each Monte Carlo trial 

would be propagated forward seven days and a hit registered at the earliest time that the hardbody radius is violated, if such a time exists (Hall et al. 2018 and Hall 2021). Temporal risk plots 

can be produced using the sequence of hits registered during all the trials, an example of which is 

shown in Figure N-12 (which shows the same conjunction sequence as Figure N-8 from the 

earlier discussion of two- and three-dimensional Nc calculation methods). The pink shaded area

in Figure N-12 shows the Monte Carlo Pc estimation confidence region, and the pink line shows 

the best estimate Monte Carlo result. As can be seen, the black upper line, which is the upper 

bound estimate from the three-dimensional Nc function, is within the confidence interval of the 

Monte Carlo results and thus is a reasonable actual realization of the repeating conjunctions’ 

cumulative risk.

Figure N-12 Three-Dimensional Nc Temporal Risk Plot with 

Monte Carlo from Epoch Result Overlay (in Pink)

Because it is complicated to set up the execution environment for the Monte Carlo from epoch

calculation, and because “gold standard” results require assembling extensive environmental data 

and software settings identical to the original DOD orbit-determination solutions, it is envisioned 

that the ability to run this strain of Monte Carlo estimation will remain with NASA CARA. 

However, a more computationally efficient mode of Monte Carlo estimation, which is 

serviceable for several different applications and is easier to obtain and employ operationally, is 

described in the next section.
Appendix N Pc Calculation Approaches 28

N.8.2 Monte Carlo from TCA

A much more computationally efficient variant on Monte Carlo from epoch, which has been 

used by conjunction assessment practitioners for some time, is Monte Carlo conducted from 

TCA, or two-body Monte Carlo (TBMC)-Pc estimation. As the definition implies, the Monte 

Carlo simulation begins with the primary and secondary objects’ equinoctial element states 

propagated to TCA. Perturbation and sampling of both states is conducted much as described 

earlier for Monte Carlo from epoch, and each sampled primary and secondary state is propagated 

both forward and backward in time to find the pair’s TCA and determine whether the 

corresponding miss distance is smaller than the hard-body radius (backward propagations are 

required to register hits that occur at times before TCA). The simplification arises from the fact 

that, since one is beginning from TCA, the propagations required will be short. This means that 

an efficient two-body motion propagation scheme usually provides an accurate trajectory 

approximation, and this, combined with the very short propagation times, vastly improves the 

computational efficiency of the calculation—by a factor of 10,000 to 100,000 according to the 

study by Hall et al. (2018). This specific method of Monte Carlo from TCA is also referred to as 

“two-body Monte Carlo” Pc (TBMC-Pc) estimation. To use the TBMC-Pc method, the 

conjunction duration needs to be short so that one may safely presume a single, unblended event 

that does not require the Monte Carlo from epoch method. As a second condition, one must have 

confidence that both objects’ states and covariances propagated to TCA are good representations 

of the states and state uncertainties at that point. Usually, there is reasonable confidence in the 

mean state estimates themselves, but the covariances are a different matter: a number of studies 

(e.g., DeMars et al. 2014) have indicated that propagated covariances represented in Cartesian 

space fail to represent the actual uncertainty distributions, due both to the potential failure of the 

linearized dynamics to remain representative over long propagation intervals and, more 

importantly, a mismatch between elongated in-track uncertainty volumes and the forced 

representation of these uncertainty volumes as Cartesian ellipsoids. The latter problem is 

illustrated in Figure N-13. The actual in-track error volume should follow the curved trajectory 

of the orbit, but the Cartesian covariance is limited to the rectilinear representation shown: as the 

elongation grows in the in-track direction (which occurs for longer propagations), the mismatch 

between the two representations also increases.
Appendix N Pc Calculation Approaches 29

Figure N-13 Mismatch between Elongated In-track Covariances and 

Forced Cartesian Ellipsoidal Representation

To address the latter problem, the results from a study conducted by Sabol et al. (2010) are both 

important and extremely helpful. This study addressed directly the question of the optimal choice 

of orbital state representation for covariances, finding that it is not the specific state 

representation in which the propagation is executed but rather the one in which the propagated 

covariance is rendered that ultimately governs the realism of the uncertainty distribution and 

associated error volume. Specifically, if the covariance is rendered and used in a curvilinear state 

representation, such as equinoctial elements, then it tends to represent the error volume much 

more accurately than if it is transformed and used in Cartesian coordinates. The surprising result 

is that a non-representative Cartesian covariance transformed into an equinoctial covariance 

becomes a representative covariance. Furthermore, taking random samples using the equinoctial 

state representation and performing the non-linear conversion of each sample to Cartesian 

coordinates generates a point cloud in the Cartesian frame that also approximates the true error 

volume much more accurately. 

This latter procedure allows the Monte Carlo from TCA method to employ more realistic 

uncertainty volumes, at least with respect to orbital state representation-related mismatches. The 

detailed procedure is the following:

1. Convert both objects’ states and covariances at TCA to equinoctial elements.

2. Generate a set of perturbations for each object based on the equinoctial covariances.

3. Combine these with the mean equinoctial states to generate sampled equinoctial states for 

the primary and secondary.

4. Convert these sampled states from equinoctial elements to Cartesian coordinates using 

the non-linear transformation.

r
Appendix N Pc Calculation Approaches 30

5. Propagate the Cartesian states for both the primary and secondary using two-body 

equations of motion to find the new TCA, which may precede or follow the nominal 

TCA.

6. Determine whether the new miss distance is less than the hard-body radius, and if so, 

register a hit at the time that the hard-body radius sphere is violated.

7. Repeat steps 5-6 until the entire set of Monte Carlo sampling trials has been processed.

8. The Pc is the number of hits divided by the total number of trials, and the confidence 

interval can be calculated from an appropriate formula.

This approach seems reasonable enough; but it would be presumptuous to assert, without further 

study, that it is truly robust, especially since the question of the durability of the linearized 

dynamics typically used to propagate covariances was not directly addressed. As it is, additional 

study efforts have been performed to verify that it is indeed sufficiently representative for 

conjunction assessment applications, and they are described below.

The first of these study efforts was performed as part of the previously cited analysis by Hall et 

al. (2018). A set of 373 high-Pc conjunctions was selected and evaluated with both Monte Carlo 

from epoch and Monte Carlo from TCA, and the comparative results are shown in Figure N-14. 

The top window is a scatter plot of the Pc calculated by Monte Carlo from epoch versus that 

from Monte Carlo from TCA. The intersection of each “plus” sign gives the scatter-plot point, 

and the length of the plus-symbol tails indicates the uncertainty of the calculation. One can see 

that the agreement is strong because all the points are close to the dashed y-x line that would 

indicate perfect equality. The bottom window plots the base-ten logarithm of the ratio of Pc 

values estimated using the Monte Carlo from epoch method to those estimated with the Monte 

Carlo from TCA method. The largest deviations are about 0.2 of an order of magnitude in Pc, 

which is considered to be below operational significance. A separate statistical test for similarity 

of results produced p-values all less than 10-3, indicating that one should reject a hypothesis that 

these results arise from different distributions. Good agreement is thus observed between the 

abbreviated Monte Carlo from TCA and Monte Carlo from epoch results.
Appendix N Pc Calculation Approaches 31

Figure N-14 Comparative Results of 373 High-Pc Conjunctions3

The second study effort involved a comparison between the Monte Carlo from TCA Pc 

estimation method and the two-dimensional Pc method (Hall 2019b), which was later extended 

to compare the from TCA method to two- and three-dimensional Nc methods (Hall 2021). In 

these studies, 63,603 temporally isolated conjunction events with two-dimensional Pc values 

greater than 10-7 were subjected to Pc calculation by the two-dimensional Pc method, as well as 

the two- and three-dimensional Nc methods, and all of these were compared to the Monte Carlo 

from TCA method. The comparative results are shown in Figure N-15 (which refers to Pc values 

estimated using the Monte Carlo from TCA method as TBMC-Pc values). These plots are similar 

to those already shown in Figure N-7, except in this case the colored diamonds represent the 

worst kind of Pc estimation failure, i.e., conjunctions in which the analytical calculations

significantly underestimate the Monte Carlo from TCA method Pc values. The left panel shows 

that, although the two-dimensional Pc method performs reasonably well for the vast majority of 

temporally isolated conjunctions, it underestimates from TCA Pc values by a factor of 1.5 or 

more in 0.258% of the investigated cases. (See the top legend of the left graph.) For those cases 

that showed large disparities, the subset that had true Pc values in the tractable range for Monte 

Carlo from epoch were validated with this methodology; and in each case the Monte Carlo from 

epoch reruns matched the output from the Monte Carlo from TCA. To the degree that nonrepresentative covariances may be responsible for two-dimensional Pc failures (due to coordinate 

frame mismatches), Monte Carlo from TCA certainly appears to be able to recover the true Pc. 

3 From Hall et al. 2018.
Appendix N Pc Calculation Approaches 32

As an aside, some of these differences between the conjunction plane two-dimensional Pc

estimates and the Monte Carlo method estimates are considerable; there are several cases in 

which the two-dimensional Pc estimate understates the Monte Carlo from TCA estimate by more 

than an order of magnitude, as shown in the red in the left panel of Figure N-15. The center and 

left panels of Figure N-15 show that the three- and two-dimensional Nc methods do not similarly 

underestimate the from TCA Pc method values, except in a few statistically insignificant cases.

Figure N-15 Comparative Results of 63,603 Conjunction Events4

Originally, Monte Carlo from TCA was advanced as a robust and computationally tractable way 

to ensure reliable Pc calculations in the face of occasional miscarriage of the two-dimensional Pc

algorithm. However, Figure N-15 demonstrates that the development and testing of the two- and

three-dimensional Nc analytic methods have led NASA CARA to recommend that these methods

supplant routine use of Monte Carlo from TCA. One may still resort to this more 

computationally efficient Monte Carlo method if desired, but testing indicates that the Nc

calculation methods outperform Monte Carlo from TCA (in that it accurately matches Pc 

estimates for isolated conjunctions and also provides a cumulative risk upper bound for repeating 

conjunctions) and is at least an order of magnitude more computationally efficient, especially for 

events with small Pc values. Monte Carlo from TCA is a capability available for download from 

the NASA CARA software repository,

5 and O/Os who are currently using it or a similar 

implementation are not urged to take it out of service, but as a Pc calculation approach, it does 

not seem to offer any enduring advantage over the two- and three-dimensional Nc methods.

Finally, a dedicated study on covariance Gaussianity was recently conducted (Lechtenberg 

2019b). The same set of 44,000 conjunctions used by Hall (2019b) was analyzed for multivariate 

normality of the position covariance in Cartesian coordinates. The methodology was to convert 

the covariance to equinoctial elements, generate a set of random position samples from this 

covariance, convert the sample set to Cartesian coordinates, and apply a multivariate normality 

4 From Hall 2019b and Hall 2021.

5 Specifically, the MCWorkbench SDK function estimates TBMC-Pc conjunction values using CDM test files as 

input. See Section 7, Contact Information in this document for the CARA software repository URL.
Appendix N Pc Calculation Approaches 33

test (here the Henze-Zirkler test) to assess compliance to this distribution. At a 5% significance 

level, only 60% of the cases could be considered to conform to a multivariate Gaussian 

distribution in Cartesian coordinates. In principle, such conjunctions would be considered 

suspect cases whose two-dimensional Pc results would be considered doubtful. As it is, since 

such a large fraction of the investigated cases show good agreement between the Monte Carlo at 

TCA and the two-dimensional Pc, clearly the Gaussianity of the covariances in the Cartesian 

framework does not matter appreciably for the Pc calculation. This result corroborates that of an 

earlier study by Ghrist and Plakalovic (2012), which reported similar findings.

How can this be, given that the covariance is integral to the two-dimensional Pc calculation? It is 

important to remember that the analyses above are restricted to high-Pc events. For the 

probability of collision to be high, significant overlap needs to exist between the primary and 

secondary object covariances; this means that the central parts of the covariances, which are 

where most of the probability density lies, have to overlap substantially. Such a requirement 

makes the behavior of the tails of the covariance much less important, and it is in the tails that 

non-Gaussian behavior is most strongly manifested. Even though a good number of conjunctions 

fail tests for covariance Gaussianity in Cartesian coordinates, for high-Pc events, this result does 

not appear to affect the rectitude of the calculated Pc. A difference probably is observable for 

lower-Pc conjunctions, but because these conjunctions are not operationally significant, it is not 

operationally important to identify or characterize this phenomenon. 

There is the lingering question of the eventual failure of the linearized dynamics in propagating 

covariances since DOD covariances are propagated through a pre- and post-multiplication by a 

linearized state transition matrix. It is agreed that, given sufficiently long propagations, such an 

eventuality should arise. However, the propagation duration required for this problem to manifest 

itself substantially is believed to be much longer than encountered for most conjunctions—on the 

order of weeks. This is why more attention is paid to this phenomenon in other areas of SSA 

such as uncorrelated track processing for which propagations of 30 days or more may be 

required, whereas it is rare for propagations longer than ten days to take place in conjunction 

assessment. Of course, fresh tracking data are always appreciated as they shorten the propagation 

interval and lend additional confidence to the solution.

N.9 Choosing an Appropriate Hard-Body Radius

As discussed in the first section of this appendix, the hard-body radius represents the combined 

size of the primary and secondary objects; the word “radius” is used because this combined size 

is typically envisioned as a sphere, and the radius of a sphere is a convenient linear 

representation of its size. The hard-body radius is needed for the Pc calculation because it 

represents the circle/sphere within which a simultaneous penetration by both the primary and 

secondary objects’ trajectories will constitute a presumed collision. It is not just a required input 

to the Pc calculation, however, it is also one of the governing parameters of the calculation: the 

Pc value represents an integral over the area of the hard-body radius circle (or surface area of the 

hard-body radius sphere) and thus in many circumstances, varies roughly in proportion to the 

square of the hard-body radius, so an increase of the hard-body radius by a factor of three 

increases the calculated Pc by a factor of nine, or nearly one order of magnitude. Because of this 

sensitivity, it is important not to overstate the hard-body radius simply for “conservatism”
Appendix N Pc Calculation Approaches 34

because the effect on the calculated Pc can be considerable. The best overall strategy for 

applying conservatism, should one wish that, is to apply it at the end of the process by lowering 

the Pc threshold at which mitigation actions should occur. Injecting conservatism into different 

places throughout the Pc calculation makes it difficult to determine how much conservatism has 

actually been introduced, whereas addressing this desire through a modification of the Pc

threshold for mitigation allows it to be understood precisely. 

Mashiku and Hejduk (2019) recently completed a study of different hard-body radius calculation 

and setting strategies, and a streamlined summary of the possibilities examined is provided 

below:

1. Large a priori value. For some years it was standard practice simply to choose a hardbody radius value that was notably larger than the expected actual size of the combined 

primary and secondary object; 20 meters is a value that was typically used. Perhaps in 

the early days of conjunction assessment, this was an acceptable initial screening 

strategy to identify potential serious events, and the hard-body radius would then be 

reduced when analysis began in earnest. As the space catalog has grown in size, and 

especially with the recent growth through the deployment of the Space Fence radar, this 

particular strategy merely creates additional false alarms that needlessly burden the 

process. Nearly all O/Os have moved away from this hard-body radius strategy.

2. Circumscribing sphere. The use of a circumscribing sphere to set the primary object’s 

hard-body radius is perhaps the most commonly used present operational technique, 

which admits of two typical variants: placing the sphere’s center at the center of mass of 

the primary satellite and defining the radius by the length from this center to the 

satellite’s most distant extremity; or allowing the center point to float freely and then 

defining the smallest circumscribing sphere. The overall sphere size then has to be 

increased by the expected size of the secondary object using either an averaged value or 

an estimate for the particular secondary encountered. This size could be obtained either 

from known or published dimensions (for intact spacecraft or rocket bodies) or estimated

from remote sensing data, such as radar cross-section or photometric brightness

measurements. In the latter case, the estimated hard-body radius size of the secondary 

may also have an associated uncertainty estimate, which can also be incorporated into 

the Pc estimation process (as discussed in more detail below).

3. Maximum projected area into any plane. Since the circumscribing sphere described 

above most often ends up being projected into the conjunction plane, it is instructive to 

examine in more detail the implications of such a projection. Clearly the sphere itself 

will project as a circle, but the projection of the three-dimensional spacecraft inside will 

necessarily be smaller in area than the projected circle, and for some spacecraft shapes 

and orientations, it will be substantially smaller in area. In this latter case, the substantial 

“white space” within the projected circular area not occupied by the primary could 

justifiably be excluded in the Pc estimation process, especially for the most common 

debris-encounter scenario when the incoming secondary object is much smaller than the 

primary asset. A straightforward way to address this issue that does require knowledge 

of the satellite’s actual orientation in relation to the conjunction plane is simply to 
Appendix N Pc Calculation Approaches 35

determine in advance the maximum area that the satellite can possibly project into any 

plane (based on a three-dimensional CAD6 model of the satellite) and use a hard-body

radius circle of that equivalent area (which of course must then be enlarged to include 

the estimated size of the secondary) for Pc estimation. It is true that this is a conservative 

formulation in that it uses the maximum possible projected area, but this is often 

substantially smaller than the area of the projected circumscribing sphere. One could 

argue that uncertainty is introduced by using the equivalent circular area rather than a 

contour integral to perform the integration over the actual projected shape, but individual 

exploratory examples show that this difference is usually negligibly small, and in any 

case, the most conservative projection approach should compensate for any differences 

in the shape chosen to represent the hard-body radius area for the Pc calculation.

4. Projection into actual conjunction plane. The most accurate, and at the same time the 

most difficult, approach is to perform a projection of the primary satellite’s shape into 

the actual conjunction plane. Specifically, this requires a three-dimensional CAD model 

of the satellite plus knowledge of its inertial attitude and the orientation of any major 

articulating components at TCA along with a calculation to project the resulting shape 

into the conjunction plane. Once this projection is obtained, its boundaries have to be 

augmented to account for the expected size of the secondary, and the integration of the 

joint covariance probability density can take place over this figure via contour 

integration or over a more convenient shape of equivalent area. Chan (2019) recently 

proposed a method to decompose complex shapes into individual triangles and use an 

equivalent-area method to evaluate the Pc for each triangle; the composite Pc is simply 

the sum of the Pc values for these individual trials of decomposition.

Each successive approach among the four presented brings greater precision to the hard-body

radius determination but at the same time, additional complexity. A reasonable via media would 

appear to be approach 3) above, which keeps the hard-body radius value grounded in reality and 

free from excessive conservatism but avoids the difficulties of gathering and maintaining shape 

and attitude data to enable a detailed projection calculation for each conjunction.

To facilitate Pc estimation, CARA has undertaken an effort to estimate hard-body radii of 

unknown orbiting objects based on radar cross section (RCS) measurements obtained by the 

Space Fence radar system (Baars and Hall 2022; Hall and Baars 2022). Even after accounting for 

radar calibration irregularities and occasional outlier data points, such RCS measurements have 

considerable point-to-point scatter, and the process of converting the RCS data into hard-body 

radius estimates introduces considerable additional uncertainty. The resulting hard-body radius 

uncertainty PDFs are non-Gaussian, although they can be roughly approximated using a lognormal distribution. Notably, RCS-based hard-body radius estimates are typically uncertain by 

factors of two to three, so it is essential that Pc estimates account for these uncertainties. For an 

unknown secondary involved in a conjunction, the RCS-based uncertainty PDF can be used to 

calculate a mean hard-body radius estimate, R̅

2, and an associated variance, σR2

2

. As mentioned 

previously, for many conjunctions, the estimated Pc value varies roughly in proportion to the 

6 Computer-aided Design
Appendix N Pc Calculation Approaches 36

square of the hard-body radius. In these cases, the appropriate hard-body radius to use for Pc 

estimation is given by the effective combined hard-body radius

Reff = √(R1 + R̅

2)

2 + σR2

2

(N-24)

with R1 representing the primary object’s hard-body radius, usually established using methods 2 

or 3 listed above. (The current recommendation for the less often used method 4 is to increase

the projected shape of the primary object outwardly on the conjunction plane by the length R̅

2 +

3σR2to produce a conservative Pc estimate, although CARA continues to research a more 

accurate approach in these cases.) So for most cases, CARA recommends calculating Pc values 

using the combined hard-body radius given in equation (N-24) for all conjunctions involving 

secondary objects with uncertain sizes, partly because of the formula’s relative simplicity, but

primarily because using the effective radius approach reproduces reasonably well the more 

accurate RCS-based Pc estimates calculated using Monte Carlo analyses. (See Hall and Baars 

2022 for details.)

Currently, Space Fence RCS measurements sufficient for size estimation are available for over 

90% of CARA conjunctions involving unknown secondary objects (Baars and Hall 2022). 

Analysis indicates that about 98% of such conjunctions involve secondary objects with mean 

hard-body radius estimates less than 35 cm, and only about 0.3% have R̅

2 > 1.5 m. CARA uses a 

multi-faceted approach to estimate secondary sizes. For known secondary objects (e.g., active or 

retired payloads, rocket-bodies), the approach uses the known hard-body radii with zero

uncertainty, tabulated along with a reference for the source of the size information. For unknown 

objects with sufficient RCS data (which includes a large fraction of tracked LEO orbital debris

objects), the CARA approach tabulates the mean estimated hard-body radius and the 1-sigma 

uncertainty, which can be used in equation (N-24). For reference, this approach also tabulates the

median hard body radius values, as well as the associated 95% and 99% confidence intervals, as

estimated from the actual RCS-based uncertainty PDF. For unknown objects with insufficient 

RCS data, one can assign a reasonably conservative default size estimate of 1.5 m.
Appendix R. List of Works Cited 37

Authored and Other References

1. Abramowitz, M. and Stegun, I. A. 1970, Handbook of Mathematical Functions, New York, 

NY: Dover Publications, Inc.

2. Akella, M. R. and K. T. Alfriend. 2000. “The Probability of Collision Between Space 

Objects.” Journal of Guidance, Control and Dynamics, Vol. 23 No. 5 (September-October 

2000) pp. 769-772.

3. Alfano, S. 2004. “Orbital Covariance Interpolation.” 14th AAS/AIAA Space Flight 

Mechanics Conference, Maui HI, February 2004.

4. Alfano, S. 2005a. “A Numerical Implementation of Spherical Object Collision Probability.”

Journal of the Astronautical Sciences, Vol. 53 No. 1 (January-March 2005), pp. 103-109.

5. Alfano, S. 2005b. “Relating Position Uncertainty to Maximum Conjunction Probability.”

Journal of the Astronautical Sciences, Vol. 53 No. 2 (April-June 2005), pp. 193-205.

6. Baars, L., Hall D., and Casali, S. 2019. “Assessing GEO and LEO Repeating Conjunctions 

using High-Fidelity Brute Force Monte Carlo Simulations.” 2019 AAS Astrodynamics 

Specialist Conference (Paper #19-612), Portland ME, August 2019.

7. Baars, L. and Hall, D. 2022 “Processing Space Fence Radar Cross-Section Data to Produce 

Size and Mass Estimates” 2022 AAS Astrodynamics Specialist Conference (Paper #22-586), 

Charlotte, N.C., August 2022.

8. Balch, M. S., R. Martin, and S. Ferson. 2019. “Satellite Conjunction Assessment and the 

False Confidence Theorem.” Proceedings of the Royal Society A: Mathematical, Physical 

and Engineering Sciences, Vol. 475, Issue 2227, 2019.

(https://doi.org/10.1098/rspa.2018.0565).

9. Butt, Y. et al. 2021. See above: National Geospatial-Intelligence Agency (NGA), Cis-Lunar 

SSA Technical Steering Group. “Proposed Message Set For All Altitude Regimes.” 

(DRAFT). November 2021.

https://www.nasa.gov/sites/default/files/atoms/files/cislunar_ssa_proposed_message_setnov2021_rev2_without_emailspdf.pdf

10. Bassa, C. G., Hainaut, O. R., and Galadi-Enriquezz, D. 2022. “Analytical Simulations of the 

Effect of Satellite Constellations on Optical and Near-infrared Observations.” Astronomy & 

Astrophysics Vol. 657, P. A75, Jan 2022.

11. Carpenter, J. R. and F. L. Markey. 2014. “Wald Sequential Probability Ratio Test for Space 

Object Conjunction Assessment.” Journal of Guidance, Control, and Dynamics Vol. 37 No. 

5 (September-October 2014), pp. 1385-1396.

12. Carpenter 2019. 

13. Carpenter 2020. 
Appendix R. List of Works Cited 38

14. Casali, S., D. Hall, D. Snow, M. Hejduk, L. Johnson, B. Skrehart, and L. Baars. 2018. “Effect 

of Cross-Correlation of Orbital Error on Probability of Collision Determination.” 2018 

AAS/AIAA Astrodynamics Specialist Conference (Paper #18-272), Snowbird UT, August 

2018.

15. Consultative Committee for Space Data Systems (CCSDS) 508.0-B-1. Conjunction Data 

Message (CDM) Blue Book. Issue 1. Recommended Standard for Conjunction Data 

Messages. Washington, D.C.: CCSDS, June 2013. 

https://public.ccsds.org/Pubs/508x0b1e2c2.pdf.

16. CCSDS 502.0-B-2. Orbit Data Messages Blue Book. Issue 2. Recommendation for Space 

Data System Standards7. Washington, D.C.: CCSDS, November 2009. 

https://public.ccsds.org/Pubs/502x0b2c1e2.pdf.

17. Chan, F. C. 2008. Spacecraft Collision Probability. El Segundo, CA: The Aerospace Press.

18. Chan, F. C. 2015. “Formulation of Collision Probability with Time-Dependent Probability 

Distribution Functions,” AAS/AIAA Space Flight Mechanics Meeting (Paper #15-233), 

Williamsburg, VA, January 2015.

19. Chan, F. C. 2019. “Collision Probability for Polygonal Cross Sections.” 2019 AAS 

Astrodynamics Specialist Conference (Paper #19-911), Portland ME, August 2019.

20. Clopper, C. and E. S. Pearson. 1934. “The Use of Confidence or Fiducial Limits Illustrated in 

the Case of the Binomial.” Biometrika, Vol. 26 (1934), pp 404-413.

21. Coppola, V.T. 2012. “Including Velocity Uncertainty in the Probability of Collision between 

Space Objects.” AAS/AIAA Spaceflight Mechanics Meeting, Charleston SC, Paper 12-247, 

February 2012.

22. D’Agostino, R. B. and M. A. Stephens. 1986. Goodness-of-Fit Techniques. New York:

Marcel Dekker, Inc.

23. DeMars, K., Y. Cheng, and M. Jah. 2014. “Collision Probability with Gaussian Mixture 

Orbit Uncertainty.” Journal of Guidance, Control, and Dynamics, Vol. 37 No. 3 (2014), pp. 

979-985.

24. Elkantassi and Davison 2022. 

25. Elrod, C. 2019. “Computational Bayesian Methods Applied to Complex Problems in Bio and 

Astro Statistics.” Doctoral Dissertation, Baylor University, July 2019.

26. Flohrer, T., Lemmens, S., Bastida Virgili, B., Krag, H., Klinkrad, H., Parrilla, E., Sanchez, 

N., Oliveira, J., Pina, F. 2013. “DISCOS – Current Status and Future Developments.” Proc. 

6

th European Conference on Space Debris, (ESA SP-723), Darmstadt, Germany, April 2013.

27. Folta, David, Natasha Bosanac, Ian Elliott, Laurie Mann, Rebecca Mesarch, and Jose 

Rosales. 2002. See above: NASA Technical Publication NASA/TP–20220014814, 

7 There is a version 3 Orbit Data Messages (ODM) that will be published relatively soon (early 2023), but it doesn't 

change the Orbit Ephemeris Message (OEM) in a material way.
Appendix R. List of Works Cited 39

“Astrodynamics Convention and Modeling Reference for Lunar, Cislunar, and Libration 

Point Orbits.” 1 December 2022.

28. Foster, J. L., and H. S. Estes. 1992. “A Parametric Analysis of Orbital Debris Collision 

Probability and Maneuver Rate for Space Vehicles.” NASA/JSC-25898 (August 1992).

29. Frisbee, J. 2015. “An Upper Bound on High-Speed Satellite Collision Probability when only 

One Object has Position Uncertainty Information.” 2015 AAS/AIAA Astrodynamics 

Specialist Conference (Paper #15-717), Vail CO, August 2015.

30. Garcia-Pelayo, R. and Hernando-Ayuso, J. 2016. “Series for the Collision Probability in 

Short-Encounter Model,” Journal of Guidance Control and Dynamics, Vol. 39, No. 8, pp. 

1904-1912, Aug. 2016.

31. Ghrist, R. W., and D. Plakalovic. 2012. “Impact of Non-Gaussian Error Volumes on 

Conjunction Assessment Risk Analysis.” 2012 AIAA/AAS Astrodynamics Specialist 

Conference, Minneapolis MN, August 2012.

32. Grubbs, F.E. 1969. “Procedures for detecting Outlying Observations in Samples.” 

Technometrics 11 (1969), 1-21.

33. Hall, D.T., M. D. Hejduk, and L. C. Johnson. 2017a. “Time Dependence of Collision 

Probabilities During Satellite Conjunctions.” AAS Space Flight Mechanics Meeting (Paper # 

17-271), San Antonio, TX, February 5-9, 2017.

34. Hall, D.T., M. D. Hejduk, and L. C. Johnson. 2017b. “Remediating Non-Positive-Definite 

State Covariances for Collision Probability Estimation.” 2017 AAS/AIAA Astrodynamics 

Specialist Conference, Stevenson WA, August 2017.

35. Hall, D.T., S. J. Casali, L. C. Johnson, B. B. Skrehart, and L. G. Baars. 2018. “High-Fidelity 

Collision Probabilities Estimated using Brute Force Monte Carlo Simulations.” 2018 

AAS/AIAA Astrodynamics Specialist Conference (Paper #18-244), Snowbird UT, August 

2018.

36. Hall, D.T. 2019a. “Determining Appropriate Risk Remediation Thresholds from Empirical 

Conjunction Data using Survival Probability Methods.” 2019 AAS Astrodynamics Specialist 

Conference (Paper #19-631), Portland ME, August 2019.

37. Hall, D.T. 2019b. “Implementation Recommendations and Usage Boundaries for the TwoDimensional Probability of Collision Calculation.” 2019 AAS Astrodynamics Specialist 

Conference (Paper #19-632), Portland ME, August 2019.

38. Hall, D.T. and Matney, M.M. 2000. “Collision Probabilities for Keplerian Orbits.” Space 

Debris, Vol. 2. (2000), pp.161-198.

39. Hall, D.T. 2021. “Expected Collision Rates for Tracked Satellites,” Journal of Spacecraft 

and Rockets, Vol. 58 No. 3 (2021), pp. 715-728.

40. Hall, D.T. 2022. “Semi-Empirical Astronomical Light Pollution Evaluation of Satellite 

Constellations.” Journal of the Astronautical Sciences, .69, pp.1893-1928, 2022. 

doi.org/10.1007/s40295-022-00358-4.
Appendix R. List of Works Cited 40

41. Hall, D. and Baars, L. 2022, “Collision and Fragmentation Probabilities for Satellites with 

Size and Mass Estimates Based on Radar Observations.” In preparation, Nov. 2022.

41.42. Hall, D. T, Baars, L. B., and Casali, S. J. 2023 “A Multistep Probability of Collision 

Computational Algorithm.” 2023 AAS/AIAA Astrodynamics Specialist Conference (Paper 

#23-398), Big Sky MT, August 2023. 

42.43. Hametz, M. E. and B. A. Beaver. 2013. “A Geometric Analysis to Protect Manned Assets 

from Newly Launched Objects – COLA Gap Analysis.” AAS 13-360.

43.44. Hejduk, M.D. 2008. “Space Catalogue Accuracy Modeling Simplifications.” 2008 AAS 

Astrodynamics Specialists Conference, Honolulu, HI, August 2008.

44.45. Hejduk, M.D., S. J. Casali, D. A. Cappellucci, N. L. Ericson, and D. E. Snow. 2013. “A 

Catalogue-Wide Implementation of General Perturbations Orbit Determination Extrapolated 

from Higher-Order Theory Solutions.” 2013 AAS/AIAA Space Flight Mechanics Meeting, 

Kauai, HI. February 2013. 

45.46. Hejduk, M.D., D. Plakalovic, M. E. Hametz, L. K. Newman, J. C. Ollivierre, B. A. 

Beaver, and R. C. Thompson. 2014. “Launch COLA Operations: An Examination of Data 

Products, Procedures, and Thresholds, Revision A.” NASA Special Technical Report. 12 

January 2014.

46.47. Hejduk, M.D., L. K. Newman, R. L. Besser, and D. A. Pachura. 2015. “Predicting Space 

Weather Effects on Close Approach Events.” 2015 AMOS Technical Conference, Kihei HI, 

September 2015.

47.48. Hejduk, M.D., F. Laporte, M. Moury, T. S. Kelso, L. K. Newman, and R. Shepperd.

2017. “Consideration of Collision ‘Consequence’ in Satellite Conjunction Assessment and 

Risk Analysis.” 26th International Symposium on Space Flight Dynamics, Matsuyama, 

Japan, June 2017.

48.49. Hejduk, M.D., D. E. Snow, and L. K. Newman. 2019. “Satellite Conjunction Assessment 

Risk Analysis for ‘Dilution Region’ Events: Issues and Operational Approaches.” Fifth 

Annual Space Traffic Management Conference Progress through Collaboration, Austin TX, 

February 2019.

49.50. Hejduk, M.D. and D. E. Snow. 2019. “Satellite Conjunction ‘Probability,’ ‘Possibility,’ 

and ‘Plausibility’: A Categorization of Competing Conjunction Assessment Risk Analysis 

Paradigms.” 2019 AAS/AIAA Astrodynamics Specialist Conference (Paper #19-652), 

Portland ME, August 2019.

50.51. Hoots, F. R. and R. L. Roehrich. 1980. “Spacetrack Report No.3 – Models for 

Propagation of NORAD Element Sets,” December 1980.

51.52. Jenkin, A.B., J. P. McVey, and G. E. Peterson. 2020. “Solution of long-coast re-entry 

COLA problems with COLA gap methods” Journal of Space Safety Engineering Volume 7, 

Issue 2, June 2020, Pages 105-112.
Appendix R. List of Works Cited 41

52.53. Johnson, N.L., P. H. Krisko, J. C. Liou, and P. D. Anz-Meador. 2001. “NASA’s New 

Breakup Model of EVOLVE 4.0.” Advances in Space Research, Vol. 28 No. 9 (2001), pp. 

1377-1384.

53.54. Kaplan, S. and B. Garrick. 1981. “On the Quantitative Definition of Risk.” Risk Analysis, 

Vol. 1 No. 1 (1981), pp. 11-27.

54.55. Kessler, D. J. 1981. “Derivation of the Collision Probability between Orbiting Objects: 

The Liftetimes of Jupiter’s Outer Moons,” Icarus, Vol. 48 (1981), pp. 39-48.

55.56. Lebedev, V., and Laikov, D. 1999. “A Quadrature Formula for the Sphere of the 131st

Algebraic Order of Accuracy,” Doklady Mathematics, Vol. 59, No. 3 (1999), pp. 477-481.

56.57. Lechtenberg, T. F. and M. D. Hejduk. 2019. “Assessment and Validation of Collision 

‘Consequence’ Method of Assessing Orbital Regime Risk Posed by Potential Satellite 

Conjunctions.” 2019 Hypervelocity Impact Symposium, Destin FL, March 2019, HVIS2019-

061.

57.58. Lechtenberg, T. 2019a. “An Operational Algorithm for evaluating Satellite Collision 

Consequence.” 2019 AAS Astrodynamics Specialist Conference (Paper #19-669), Portland 

ME, August 2019.

58.59. Lechtenberg, T., 2019b. “Multivariate Normality of Cartesian-Framed Covariances:

Evaluation and Operational Significance.” 2019 AAS Astrodynamics Specialist Conference 

(Paper #19-671), Portland ME, August 2019.

59.60. Li, Z., Ziebart, M., Bhattarai, S., Harrison, D., Grey, S. 2018. “Fast Solar Radiation 

Pressure Modelling with Ray Tracing and Multiple Reflections.” Advances in Space 

Research, Vol. 61, No. 9, pp. 2352-2365.

60.61. Mashiku, A.K. and M. D. Hejduk. 2019. “Recommended Methods for Setting Mission 

Conjunction Analysis Hard Body Radii.” 2019 AAS Astrodynamics Specialist Conference 

(Paper #19-702), Portland ME, August 2019.

61.62. Mclean, F., Lemmens, S., Funke, Q., Braun, V. 2017. “DISCOS 3: An Improved Model 

for ESA’s Database and Information System Characterising Objects in Space.” Proc. 7th

European Conference on Space Debris, Darmstadt, Germany, June 2017.

62.63. Patera, R.P. 2001. “General Method for Calculation Satellite Collision Probability. 2001. 

” Journal of Guidance, Control, and Dynamics, Vol. 24, No. 4 (July-August 2001), pp. 716-

722.

63.64. Press, W. H., Teukolsky, S. A., Vettering, W. T., and Flannery, B. B. 1992. Numerical 

Recipes in FORTRAN: The Art of Scientific Computing, 2nd ed., Cambridge University Press, 

New York, NY, 1992.

64.65. Rosenberg, F. 2003. “Radar Characterization of RORSAT Coolant Spheres and their 

Potential Use for Calibration and Validation of Debris Tracking Sensors.” MIT Lincoln 

Laboratory Project Report STK-264, November 2003.
Appendix R. List of Works Cited 42

65.66. Rosner, B. 1975. “On the Detection of Many Outliers.” Technometrics 17 (1975), 221-

227.

66.67. Rosner, B. 1977. “Percentage Points for the RST Many Outlier Procedure.” 

Technometrics 19 (1977), 307-312.

67.68. Sabol, C., T. Sukut, K. Hill, K. Alfriend, B. Wright, Y. Li, and P. Schumacher. 2010. 

“Linearized orbit covariance generation and propagation analysis via simple Monte Carlo 

simulations.” AAS/AIAA Space Flight Mechanics Meeting (paper AAS-10-134), San Diego, 

CA, February 2010.

68.69. Schilling, B., Y. Taleb,, J. Carpenter, M. Balducci, and T. Williams. 2016. “Operational 

Experience with the Wald Sequential Probability Ratio Test for Conjunction Assessment 

from the Magnetospheric MultiScale Mission.” 2016 AIAA/AAS Astrodynamics Specialist 

Conference, Long Beach CA, August 2016.

69.70. Seitzer, P. 2020. (U MI) “Large Constellations of LEO Satellites and Astronomy.” 

Advanced Maui Optical and Space (AMOS) Surveillance Technologies Conference, 

September 15-18, 2020.

70.71. Shelton, C.T. and J. L. Junkins. 2019. “Probability of Collision between Space Objects 

Including Model Uncertainty.” Acta Astronautica, Vol.155 (2019), pp.462-471.

71.72. Sivia, D.S. and Skilling, J. 2006. Data Analysis: A Bayesian Tutorial. 2

nd ed., Oxford, 

UK, Oxford University Press.

72.73. Stokely, C.L., Foster, J.L., Stanbery, E.G., Benbrook, J.R., and Juarez, Q. 2006. 

“Haystack and HAX Radar Measurements of the Orbital Debris Environment.” National 

Aeronautics and Space Administration Report, Orbital Debris Program Office, (JSC-62815), 

Nov. 2006.

73.74. Tarzi, Zahi B., Brian T. Young, David S. Berry. 2022. "Deriving Event Thresholds and 

Collision Probability for Automated Conjunction Assessment at Mars and the Moon", AAS 

22-042, 2022 American Astronomical Association (AAS)/American Institute of Aeronautics 

and Astronautics (AIAA) Astrodynamics Specialist Conference, Charlotte, NC, August, 

2022. https://www.nasa.gov/cara/madcap.

74.75. Tobiska, W.K., D. Knipp, W. J. Burke, D. Bouwer, J. Bailey, D. Odstrcil, M. P. Hagan, 

J. Gannon, and B. R. Bowman. 2013. “The Anemomilos Prediction Methodology for Dst.”

Space Weather (Vol. 11, No. 9 (September 2013) pp 490-508. DOI: 10.1002/swe.20094.

75.76. Walker, C. et al. 2020, “Impact of Satellite Constellations on Optical Astronomy and 

Recommendations Toward Mitigations.” SatCon-1 Workshop Report, NSF NOIRLab, Aug 

2020.

76.77. Walker, C. et al. 2021. “Executive Summary.” Report of the SatCon-2 WorkShop, NSF 

NOIRLab, Jul 2021.

77.78. Vallado 1981. 
Appendix R. List of Works Cited 43

78.79. Vallado, D.A. 2001. Fundamentals of Astrodynamics and Applications. 2

nd ed. El 

Segundo, CA, Microcosm Press.

79.80. Woodburn, J. and S. Tanygin. 2014. “Coordinate Effects on the use of Orbit Error 

Uncertainty.” International Symposium on Space Flight Dynamics 24 (Baltimore, MD), May 

2014.

80.81. Zaidi, W.H. and Hejduk, M.D. 2016. “EOS System Covariance Realism.” 2016 

AAS/AIAA Astrodynamics Specialist Conference, Long Beach CA, September 2016.
Appendix R. List of Works Cited 44

National Aeronautics and Space Administration

Headquarters