1

DERIVING EVENT THRESHOLDS AND COLLISION PROBABILITY 

FOR AUTOMATED CONJUNCTION ASSESSMENT AT MARS AND 

THE MOON 

Zahi B. Tarzi,* Brian T. Young†, and David S. Berry‡

It is well known that conjunction assessment is necessary in the Earth orbiting 

environment to prevent spacecraft collisions and reduce debris in orbit. Although 

there are many fewer spacecraft in orbit around Mars and the Moon, the number 

of missions to these bodies is increasing rapidly and the consequences of creating 

debris riskier due to the inability to track debris objects in those environments. 

The Multimission Automated Deepspace Conjunction Assessment Process 

(MADCAP) is used at the Jet Propulsion Laboratory (JPL) to perform conjunction 

assessment at Mars, the Moon, and Sun/Earth libration points. Previous papers 

have described this process and its development. Conjunction assessment of 

spacecraft trajectories produces many close approaches which must be screened 

in order to inform users when a specific event may be of concern. Probability of 

Collision (PC) is the main attribute used for screening conjunctions in the Earth 

orbiting environment. Trajectory uncertainty information in the form of covariance data is required to calculate the PC. However, covariance data is not as easily 

obtained for objects in the Lunar and Martian environments. Instead, various attributes of the close approach must be calculated and compared against thresholds 

based on typical expected orbit uncertainty. These thresholds are best informed 

by the individual spacecraft navigation teams who are most knowledgeable about 

the spacecraft’s orbital uncertainties. This paper will explore the development of 

the thresholds used in MADCAP and provide a guideline for navigation teams to 

derive them from their orbit determination analysis. A method for the derivation 

of close approach attribute thresholds from covariance data, when it is available, 

is also described. Finally, collision probability calculations are examined, with 

exploration of methods for approximating probabilities when covariance data are 

limited. 

INTRODUCTION 

The ongoing problem with space debris in the Earth orbital regime has been well documented. 

There are currently no known debris fields at the Moon and Mars, but the creation of such debris 

fields would be highly undesirable due to the infeasibility of tracking such objects from Earth. 

There are a growing number of missions in the Martian and Lunar orbital environments, with many 

 

*

 Navigation Engineer, Mission Design and Navigation Section, Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive Pasadena, CA 91109. 

†

 Navigation Engineer, Mission Design and Navigation Section, Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive Pasadena, CA 91109. 

‡

 Program Area Manager, Mission Design and Navigation Section, Jet Propulsion Laboratory, California Institute of 

Technology, 4800 Oak Grove Drive Pasadena, CA 91109. 

(Preprint) AAS 22-042
 2

more planned in the near future. The creation of debris fields in these environments would greatly 

jeopardize future operations in those orbital regimes for both robotic and human missions. Preventing such debris fields is imperative to safely continue spacecraft operations. 

NASA currently expends resources to monitor the Earth's orbital debris environment via its 

Conjunction Assessment and Risk Analysis (CARA) program located at the Goddard Space Flight 

Center (GSFC). 1, 2, 3 The Multimission Automated Deepspace Conjunction Assessment Process 

(MADCAP) has been used at the Jet Propulsion Laboratory to perform conjunction assessment at 

Mars and the Moon since 2011, and Sun/Earth libration points since 2020. The MADCAP process 

has been described in previous publications. 4, 5, 6 The general process is reviewed here, before 

entering a more detailed explanation of how thresholds are used to categorize conjunction events. 

The derivation of these thresholds from orbit uncertainty data delivered by spacecraft navigation 

teams is investigated. Updates in covariance processing capabilities and Probability of Collision 

(PC) calculations are also covered. Areas for future work are also identified, with some enhancements already in progress. 

MADCAP GENERAL DESCRIPTION 

The fundamental design concepts of MADCAP have not materially changed since their conception in 2011. A parameter file is setup for each orbital environment to be analyzed, which allows 

conjunction analyses in any orbital environment without modifying the underlying software. The 

parameters fall into a few general classes: environment (central body, coordinate system); bodies 

within the environment (active spacecraft, inactive spacecraft, natural bodies); thresholds used to 

classify conjunction events and control report generation; options for detailed reports and plots; 

and email lists for report participants. The main parameters that establish the orbital environment 

are the specification of the central body and a list of at least two spacecraft (or other bodies including natural satellites or debris). 

MADCAP requires information about the trajectory of the orbiting bodies in order to find close 

conjunctions, and information about the uncertainty of those trajectories in order to evaluate the 

risk of collision. In Earth orbit, there are many ways to compute the orbits of spacecraft, including 

an onboard GPS receiver providing continual updates to the position, and active and passive tracking using radar, optical, and radiometric methods. Passive radar methods can also be used for 

inactive spacecraft, and orbital elements for all non-classified spacecraft are available through TwoLine Element (TLE) files provided by the Department of Defense. At lunar distances and beyond, 

however, the options are fewer, and ephemerides of these spacecraft are obtained primarily from 

the project’s ground-based navigation teams computing their orbits using active radiometric tracking. The number of antennas that can perform this tracking are limited; the primary ones used by 

NASA are those of the Deep Space Network (DSN). At the time of MADCAP’s conception, all 

spacecraft in these environments used DSN tracking data. Thus, the navigation teams for these 

spacecraft regularly upload updated trajectory ephemeris files to the DSN Service Preparation Subsystem (SPS) portal; a central repository convenient for the DSN and other users to easily access 

the orbital information. The SPS portal also allows for autonomous downloads of these files, which 

makes it a convenient place for MADCAP to collect these files from in order to conduct conjunction 

assessment. In cases where SPS is not used, trajectory information can be loaded manually or 

through specialized automated interfaces that can be worked on a case-by-case basis. 

Ephemeris data is stored on SPS in two formats: SPICE (Spacecraft Planet Instrument C-matrix 

Event) SPK (Spacecraft and Planet Kernel) 10 and CCSDS (Consultative Committee for Space Data 

Systems) OEM (Orbit Ephemeris Message) 11. The SPK is a standardized format which allows any 

user with the publicly available SPICE software to access the trajectory information. The OEM 
 3

stores ephemeris information in a standardized text format that is both machine and human readable. A third format, which is occasionally used internally at JPL if ephemeris information is provided independently and not through the SPS portal, is the MONTE (Mission Design and Operations Navigation Toolkit Environment) BOA12 (Binary Object Archive) format. MONTE is JPL's 

signature astrodynamics computing platform, used for space mission design and in-flight navigation, with ephemeris information stored in the binary BOA file format.

In addition to the ephemerides of the spacecraft, it is important for calculating collision risks to 

include the uncertainty associated with the position estimates stored in the ephemeris files. This is 

normally computed by the spacecraft’s navigation team as part of the orbit determination (OD) 

process. Ideally, MADCAP would use this uncertainty estimate as it represents the best estimate 

of the position errors, represented by a 3x3 position covariance matrix which can be evaluated at 

any time along the orbit. However, because the DSN does not utilize this information, navigation 

teams do not include it in the files submitted to SPS. Furthermore, there is no mechanism in the 

SPK or BOA formats to provide the covariance – only the OEM format allows a user to include 

this. When the formal covariance data is not available, an approximation of the uncertainty can be 

computed using time-varying polynomials of orbit radial and timing uncertainties relative to the 

trajectory creation time. These are also specified by spacecraft navigation teams, with assistance 

from the MADCAP team, and used to determine the uncertainties of the state at a given conjunction. The formal covariance information is preferred because it is a more accurate representation 

of the true errors. However, to date, covariance information has not been consistently delivered, 

so polynomials are generally being used to evaluate collision risk. Future work is planned to automate the delivery of trajectory and covariance files to MADCAP outside of SPS in order to handle 

the increasing number spacecraft arriving in the Martian and Lunar environments which are not 

tracked by the DSN. 

 MADCAP Process 

MADCAP is automatically initiated and downloads the latest ephemerides from the DSN SPS 

portal that were prepared by the navigation teams for tracking purposes. Two basic types of files 

are downloaded from SPS: "predicts grade" ephemerides and "scheduling grade" ephemerides. The 

predicts grade ephemerides represent the navigation team's best estimate of the spacecraft trajectory; these ephemerides are used in the generation of DSN pointing and frequency predicts. The 

scheduling grade ephemerides may represent a lower fidelity predict used for scheduling antenna 

time or a "reference trajectory". The reference ephemeris files are usually longer duration files that 

represent a reference, baseline, or nominal trajectory, and often include some future planned maneuvers. The reference ephemeris files are typically updated less frequently than predicts grade 

ephemerides. 

The SPS metadata is parsed to find the most recent predicts grade SPK and OEM files available. 

Files are checked against those downloaded in the last run to determine if they are new. If the most 

recent OEM and SPK files are based on the same submission, then the OEM file will be downloaded and later checked for covariance data (position covariance matrices at specified times). The 

SPK file is used for conjunction assessment analysis while the OEM file is used to obtain covariance data to calculate thresholds. If no matching predicts grade OEM file is found, MADCAP will 

then search for a matching “ProjectUse” grade OEM file to use. The “ProjectUse” grade files are 

not used by the DSN for tracking purposes, but are intended for projects to be able to provide 

supplemental information . Alternatively, a manually specified trajectory file can be used for conjunction assessment instead of downloading the latest file from SPS. One additional ephemeris is 

allowed to be specified per spacecraft. This can be either the latest scheduling grade SPK file on 
 4

SPS, or a local file. Covariance data are not considered for additional ephemerides. A flow chart 

showing these possible trajectory and covariance file sources is shown in Figure 1. 

Figure 1. MADCAP Trajectory and Covariance File Sources for a Given Spacecraft 

Given the ephemerides of two orbiting objects, the next step is to find the close approaches 

between them. MADCAP performs pairwise comparisons of the ephemerides of the bodies listed 

in the input parameter file. Analysis for all combinations of two objects is performed (e.g., for 

spacecraft A, B and C, results are calculated for A-B, A-C, and B-C). Comparisons are made over 

the duration of the overlapping time period of the two ephemeris files analyzed (with a reduced 

overlap during times when the objects are nearly coplanar or very long duration). The pairwise 

analysis will be performed once using the nominal ephemeris file specified for each body and then 

again for all possible combinations of the additional files specified (1 additional per body, thus a 

total of 4 analyses per pair is possible). For each analysis, a search is conducted for local minimum 

relative distances between the two bodies analyzed; each relative minimum is considered a “Close 

Approach Event”. Times of the events and various orbit attributes are printed to an output table, 

plotted, and used to classify events into two main categories based on impact risk: ‘Red’ and ‘All’. 

These are described in detail in the next section. 

For the case of two non-coplanar orbits, an impact can only occur where the two orbit planes 

intersect. Therefore, the impact risk can be characterized by two uncertainties: that of the time at 

which the primary spacecraft crosses the orbit plane of the secondary spacecraft, and that of the 

radius from the central point of the orbit at the time of that orbit plane crossing. With a known time 

and radius of both spacecraft at that point, the feasibility of an impact can be determined based on 

whether the difference in timing and the difference in radius is within a certain number of standard 

deviations of the uncertainty, typically three. Given this, there are three main close approach event 

attributes used to analyze the events: close approach distance (CAD), orbit crossing distance 

(OXD), and orbit crossing timing (OXT). The CAD is the relative distance between the two bodies 
 5

at the time of the closest approach (t஼஺). The OXD is the minimum distance between the orbits of 

the two bodies. For non-coplanar orbits, this occurs at the two points where the orbits cross each 

other. Since the orbits are slowly changing over time, the orbit crossings that are within one orbital 

period of the t஼஺ are evaluated, and the one with the minimum OXD is reported for the close approach event. The OXD value can be represented by the following equation: 

OXD = rଵ

(tை௑ଵ) − rଶ(tை௑ଶ) (1)

Where r(tை௑) represents the orbital radius of each body evaluated at the time it is at the orbit 

crossing location. The “sign” of this number is meant to convey information about which orbit is 

higher at the point at which the orbits cross. So, for a “Primary-Secondary” body pairing, if the 

orbit crossing distance is positive, then the primary orbit is above the secondary orbit at the crossing 

time. If it’s negative, then the primary orbit is below the secondary orbit. A visual representation 

of CAD and OXD is shown in Figure 2, where r⃑ை௑ represents the body’s orbit radial direction at 

the time of the crossing. The third attribute, OXT, is defined as the difference between the times 

that the two bodies are at the crossing. This can be represented as: 

OXT = tை௑ଵ − tை௑ଶ (2)

If the value is positive, then the primary arrives later than secondary; if it is negative then the reverse 

is true. 

Figure 2. Main Close Approach Event Attributes 

In the case the two spacecraft are in coplanar or nearly coplanar (hereafter referred to as coplanar) obits, an updated algorithm was developed to calculate minimum orbit distances. This algorithm uses a brute force method to compare distances at different points along both orbits to find 

the minimum distance between them. It has an increased run time, but provides more accurate 

results when orbits are coplanar. Thus, a check has been inserted to determine when the orbits 

being compared are coplanar (angular momentum vectors within a parameter specified angle). For 

each close approach event, if the orbits are determined to be coplanar, the updated algorithm is used 

to calculate minimum orbit distances and timing which are reported as orbit crossing distances and 
 6

timing in both the detailed tables and summary reports. To avoid long run times for bodies which 

are often coplanar, the coplanar algorithm is only used for events within 60 days from the current 

time. Another way MADCAP avoids long run times is to not analyze natural bodies against each 

other (e.g. Phobos v Deimos is not analyzed unless specifically requested). 

A summary report is generated for each MADCAP run which highlights important results from 

the conjunction analysis. The report is generated in HTML (Hyper Text Markup Language) format 

and sent out via email to a list of interested parties. Text messages which highlight conjunctions 

of concern can also be sent to specified users. More detailed plots and tables of the close approach 

events are also generated and can be sent out via email for each body pairing. An example summary 

report is included in the Appendix. 

THRESHOLDS 

The close approach attributes mentioned in the previous section are checked against thresholds 

to place events into two categories: ‘Red’ and ‘All’. The ‘Red’ event list is comprised of significant, near-term events which may require further investigation, while the ‘All’ event list is composed of any future events which may be of interest. These thresholds are specified by the individual spacecraft navigation teams based on the criteria explained in the following sections. 

 ‘Red’ Events 

The category of ‘Red’ Events is intended to represent significant, near-term predicted conjunction events. This category includes data for all of the events that meet the ‘Red’ thresholds and 

occur within 14 days of the analysis time (defined as the time that a MADCAP run was initiated). 

Only pairings that involve at least one active spacecraft using the primary specified ephemeris files 

(predicts grade for those on SPS) are included in the ‘Red’ Event category. Pairings involving 

inactive spacecraft are not categorized as ‘Red’. 

As mentioned above, the three types of threshold categories are orbit crossing distance (OXD), 

orbit crossing timing (OXT), and close approach distance (CAD) thresholds. A conjunction event 

is considered ‘Red’ when the OXD and OXT for that event are less than the ‘Red’ thresholds. The 

CAD threshold is not considered for ‘Red’ Events, but the CAD value is still listed for reference in 

the summary report. The ‘Red’ thresholds correspond to the 3-sigma uncertainties as follows: 

OXD-radial position uncertainty, OXT-downtrack timing uncertainty. When covariance files become available, they will be used to calculate these values based on an interpolation and mapping 

of the position covariance matrices which bracket the event in time (see the next section for detailed 

description of covariance processing). The threshold values for both bodies in a pairing are calculated for each event, and the Root Sum Square (RSS) of the two values is then used as the event 

threshold. In the absence of formal covariance files, ‘Red’ thresholds are based on a quadratic fit 

of the 3-sigma uncertainty values (as provided by the spacecraft navigation team) as a function of 

time to the event. Figure 3 shows a flowchart of this process, where the threshold is represented 

by τ. 

The polynomial equations used to calculate ‘Red’ threshold values for bodies without covariance data specified are described as follows: 

RED OXD Threshold = 3σ௥ = OXD0 + OXD1 ∗ t + OXD2 ∗ t

ଶ

 [km] (3)

RED OXT Threshold = 3σ௧ = OXT0 + OXT1 ∗ t + OXT2 ∗ t

ଶ

 [km] (4)

 Where t = Close Approach Time − Ephemeris Submit Time (in days) (5)
 7

Figure 3. Threshold Sources 

The six coefficients above must be specified for each body (except for inactive spacecraft which 

are not considered for ‘Red’ Events). Bodies using constant thresholds simply have zero values 

for OXD1, OXD2, OXT1, and OXT2. The time in days, t, stated above can be any positive real 

number (not limited to integer value). In order to obtain the values of these coefficients, spacecraft 

navigation teams must fit a second order polynomial to their 3σ worst case nominal radial (for 

OXD) and timing (for OXT) uncertainties over the 14 days after last available tracking data. 

“Worst case nominal” is meant to reflect the largest 3σ uncertainties experienced during nominal 

spacecraft operations (Ex: certain attitude configuration or pointing mode, not safe mode). Alternatively, spacecraft navigation teams can deliver their worst case nominal radial and timing uncertainty data to MADCAP analysts to produce a curve fit. In certain cases, the polynomial fit may 

be altered from “best-fit” of the data to a curve that better reflects the intentions of its use (Ex: 

avoid negative thresholds, better fit in first 10 days vs last 4 days). In such cases, the spacecraft 

navigation teams will be notified of this alteration and must approve of the final fit. Figure 4 shows 

four examples of worst case nominal orbit uncertainty data that was fit to a second order polynomial 

to obtain MADCAP thresholds. 

These fits are an estimate of the uncertainty of the predicted trajectory based on the available 

data when the prediction was made. The uncertainty of the prediction is based on the amount of 

time that has passed since the last available tracking data used for the prediction, regardless of when 

the analysis was run. The ephemeris submit time to SPS provides a satisfactory general approximation of the data cutoff time since the actual data cutoff is not available in the ephemeris file. If 

a manually specified ephemeris file is used instead of one submitted to SPS, the analysis time will 

be used in place of the submit time. 

The threshold coefficients for a body are colored blue and noted with an “*” in the summary 

report if they have changed since the last MADCAP run. The polynomial coefficients used are 

listed in a table in the report, and the threshold source for each ‘Red’ event is also reported in body 
 8

pair order (Ex: P-P, P-C, etc. where P=polynomial and C=covariance). The coefficients are also 

colored blue if they have changed since the last run. 

Figure 4. 2nd Order Polynomial Fits of 3σ Uncertainties of Select Mars Orbiters 

‘All’ Events 

The ‘All’ Events category is intended to represent any future events which may be of interest to 

the recipients of the summary report. It includes data for all of the events identified in the analysis 

that have orbit crossing and close approach distances less than the constant ‘All’ OXD and ‘All’ 

CAD thresholds for all pairings of active spacecraft. The OXT threshold is not used for the ‘All’ 

events, but the timing values are still listed in the summary report for reference. For any particular 

pairing of bodies, the analysis is performed for the overlap of the two ephemeris files used, from 

the present time through the end of the overlap, or for a specified maximum number of days (as 

decided by the appropriate input parameters). The notes on sign convention and epoch time in the 

above ‘Red’ Events section apply to the ‘All’ Events as well. The ‘All’ Events thresholds are 

provided by the spacecraft teams in order to cover any upcoming events they would like to see 

listed in the summary report. The ‘Red’ Thresholds are assumed to be tighter than the ‘All’ thresholds, such that only the events in the ‘All’ Events category are screened for inclusion in the ‘Red’ 

Events category. 

Data for active spacecraft and natural bodies are included in this category. Data for inactive 

spacecraft are not included, but events from analysis using ‘additional’ trajectories are encompassed. However, the ‘additional’ trajectory events are only displayed if the event comes after the 
 9

end of the time span covered by the ‘main’ trajectory file (see Figure 1 for source of ‘main’ vs 

‘additional’ files). This is done to give precedence to the ‘main’ file and not confuse users with 

multiple listings for the same event. 

DERIVING THRESHOLDS FROM COVARIANCE DATA 

The polynomials described in the previous section are approximations taken from computations 

of expected uncertainties, and are not updated frequently, so necessarily represent a worst case 

nominal scenario. This causes events to be flagged that may be much safer than this approach 

would imply. For that reason, there has been increased interest in the delivery of covariances generated by the navigation team’s OD filter as part of the standard delivery that would allow these 

uncertainties to be computed in a more robust way. 

 Characterizing the Desired Uncertainties 

The first component to consider is the “timing” error, more explicitly defined as the uncertainty 

of the time of the event at which the primary body crosses the orbit plane of the secondary body. 

Typically, because it is assumed that variation in the anomaly of the orbit is the most significant 

error in any orbit prediction, this has been approximated by projecting the position covariance into 

the velocity direction, to determine the downtrack position error, and dividing by the nominal speed 

to convert to timing. However, this computation is an approximation, since off-axis position errors 

and velocity errors also affect the geometry of this crossing, though in practice the effect of this is 

very small. It is important to clarify that the position covariance that is projected into the velocity 

direction is the covariance at the epoch of the nominal plane crossing; that is, it is not constrained 

to the plane crossing event, which would by definition reduce the uncertainty drastically by eliminating the velocity-direction position errors from the uncertainty. 

The second component is the “radial” error. More explicitly, there is interest in the uncertainty 

of the radius of the primary body at the time it crosses the orbit plane of the secondary body. This 

is calculated by projecting the position covariance at the event of the plane crossing onto the nominal radius vector. The computation from the constrained covariance is important, because the radius of the orbit will shift with the downtrack variation of the orbit according to the flight path 

angle, so that there is a high degree of correlation between the two components. Only when the 

downtrack portion of the covariance is constrained to the plane of the other orbit does the radial 

component correctly reflect the value that potentially poses a danger to the two bodies. 

Throughout the rest of this document, the term “constrained” covariance will refer to the covariance of the state at the time-varying plane-crossing event, while the term “unconstrained” covariance will refer to the variation of the state at the epoch of the nominal plane crossing. 

Proposed Algorithm 

In the previous section, a tension stands out. The timing error is naturally defined by the unconstrained covariance, while the radial error is naturally defined by the constrained. Previous efforts 

attempted to project a single covariance, either constrained or unconstrained, onto the position and 

velocity vectors to determine the radial and timing uncertainties, but this method was destined for 

failure, either drastically increasing the radial errors or drastically reducing the timing errors. The 

navigation teams will not necessarily know the trajectory of the other conjunction object in advance, so constraining covariances to the plane-crossing event is not practical. Fixed-anomaly constraint mappings would be possible, though asking teams to deliver two sets of covariances does 

not seem practical either, and unconstrained covariances are easier to both generate and explain. 

Therefore, it is recommended for projects to deliver unconstrained covariances of the position and 

velocity in EME2000 at fixed times, which can be converted to radial and timing errors using the 
 10

formulae described here.7 EME2000 is the commonly used Earth Centered Inertial (ECI) coordinate frame defined with the Earth's Mean Equator and Mean Equinox (MEME) at 12:00 Terrestrial 

Time on 1 January 2000. The x-axis is pointed toward the mean equinox, the z-axis pointed toward 

the Earth's rotation axis (at that time), and y-axis is the cross product of x and z. 

Timing Error. In general, an “event” is the time at which some function f(xሬ⃑,t) = 0 where xሬ⃑ are 

the constant orbit parameters defining a trajectory, such as the state at a specific time. Given a 

nominal trajectory and the epoch that meets this constraint, defined by xሬ⃑෡ and t̂, the event function 

can be linearized about that point as 

f൫xሬ⃑෡ + xሬ⃑,t̂+ ∆t൯ = f൫xሬ⃑෡,t̂൯ +

డ௙

డxሬ⃑

∆xሬ⃑ +

డ௙

డ௧ ∆t = 0 +

డ௙

డxሬ⃑

∆xሬ⃑ +

డ௙

డ௧ ∆t (6)

Solving this linearized form for when it equals zero yields the partial derivative of the event time 

relative to the parameters: 

డ௧

డxሬ⃑ =

∆௧

∆xሬ⃑ = −

డ௙⁄డxሬ⃑

డ௙ డ௧ ⁄

(7)

Now, we can define the function for the plane crossing by recognizing that the plane crossing event 

is defined by the epoch at which rሬ⃑ଵ, the position of the primary body, is perpendicular to hሬ⃑

ଶ, the 

angular momentum vector of secondary body, which is assumed to be constant for the duration of 

interest for a given conjunction (so that we avoid considering its time derivatives). Defining the 

function and taking the derivative gives the partial derivative of the event time relative to the state 

at the nominal epoch of the plane-crossing 

fை௑ = rሬ⃑ଵ ∙ hሬ⃑

ଶ

(8)

డ௧ೀ೉

డrሬ⃑భ

= −

డ௙ೀ೉⁄డrሬ⃑భ

డ௙ೀ೉⁄డ௧ = −

hሬ⃑మ

vሬ⃑భ∙hሬ⃑మ

(9)

where vሬ⃑ଵ represents the velocity vector of the primary body. Applying this transform to the position covariance gives the uncertainty of the timing event: 

σ௧ଵ =

ටhሬ⃑

మ

೅Pభhሬ⃑మ

vሬ⃑భ∙

ሬሬhሬ⃑మ

(10)

where Pଵ represents the unconstrained position covariance of the primary body. The timing error 

can be understood as the unconstrained covariance projected onto the normal vector of the orbital 

plane, scaled by the velocity projected onto that same normal vector. When the largest principal 

axis of the covariance is also aligned with the velocity vector, this is equivalent to our previouslydefined approximation. Note that this can easily be computed from the unconstrained covariance 

and the nominal trajectories of the two spacecraft. 

Radial Error. The radial error can be understood as the constrained position covariance projected onto the position vector direction. In order to compute the transition matrix from the unconstrained state to a constrained state: 

డxሬ⃑೎

డqሬ⃑ =

డxሬ⃑ೠ

డqሬ⃑ +

డxሬ⃑ೠ

డ௧

డ௧೐ೡ೐೙೟

డqሬ⃑

(11)
 11

where xሬ⃑௖ represents the constrained state, xሬ⃑௨ represents the unconstrained state. Assuming that qሬ⃑

is the unconstrained state at the epoch of the node crossing xሬ⃑௨, we can get a transition matrix from 

the unconstrained to the constrained state: 

డxሬ⃑೎

డxሬ⃑ೠ

=

డxሬ⃑ೠ

డxሬ⃑ೠ

+

డxሬ⃑ೠ

డ௧

డ௧ೀ೉

డxሬ⃑ೠ

= I଺௫଺ + ൤

vሬ⃑ଵ

aሬ⃑ଵ

൨൤

ିhሬ⃑మ

೅

vሬ⃑భ∙hሬ⃑మ

0

்

൨ (12)

where I଺௑଺ represents the six by six identity matrix. The acceleration, aሬ⃑ଵ, is multiplied by zero 

and has no impact. From this transformation, the radial uncertainty can be computed by projecting 

the constrained position covariance onto the radial direction, for a final formulation (eliminating 

zero terms) where: 

σ௥ଵ = ටФሬ⃑Pଵ

ሬФ்⃑ (13)

 where Фሬ⃑ =

rሬ⃑భ

೅

௥భ

൬Iଷ௫ଷ −

vሬ⃑భhሬ⃑మ

೅

vሬ⃑భ∙hሬ⃑మ

൰ (14)

Correlations. While it is tempting, and mathematically straightforward, to compute the correlations of these values and potentially compute the Mahalanobis norm to use as a threshold, that is 

not recommended at this time. Because the radial and timing errors become highly correlated after 

a few orbits, the semiminor axis of this joint covariance is very small. JPL navigation computes 

these covariances using a linearized state transition matrix propagation, which leads to small errors 

in these values, as can be demonstrated with simple two-body propagation Monte Carlo simulations. With small correlations this would be fine, but as the correlations grow, small linearization 

errors lead to large areas that do not overlap between the 3σ ellipses represented by linearized and 

true covariances. Therefore, while we trust the ‘corners’ of the bounding rectangle of the covariance 

in the radial/timing space to be accurate to within a few percent, the width and orientation of the 

ellipse within that rectangle are extremely sensitive and cannot be accurately used for computation. 

Results. The radial and timing errors were computed using these equations for a Mars orbiter over 

a three-week period of time with results shown in Figure 5. The computed errors are compared to 

those obtained via the polynomial fits of uncertainty described in the previous section showing that 

the polynomial fit used was slightly conservative in this case. 

Figure 5. Covariance derived timing and radial errors compared to polynomials for a Mars orbiter 

Uncertainty in Secondary Orbit Plane. The timing and radial error equations above were derived 

without considering the uncertainty in the angular momentum vector of the secondary spacecraft, 
 12

hሬ⃑

ଶ. The uncertainty in angular momentum for orbiters is typically very small, but can be accounted 

for if the information is available. However, covariance data is currently still rarely reported to 

MADCAP and information on orbit normal uncertainty has not historically been requested of missions by the MADCAP team. Testing based on data from a current Mars orbiter has shown that 

induced uncertainties in the primary spacecraft’s radial error by the secondary body’s normal error 

can be up to a few hundred meters in certain cases (induced timing uncertainty is much smaller, 

sub-second level). Though not negligible, these errors represent a small fraction of the overall 

uncertainty and are currently unaccounted for by MADCAP. Future methods for inclusion could 

involve obtaining the angular momentum uncertainty of the secondary body from covariance data 

in OEM files, requesting a constant worst case normal uncertainty from spacecraft navigation 

teams, or assuming a constant worst case normal uncertainty for all bodies analyzed by MADCAP. 

These methods are still under consideration for future MADCAP enhancements. 

Interpolation of Covariances. In order to obtain the unconstrained covariances at the close approach time from data delivered at arbitrary intervals, MADCAP must interpolate the covariance 

data delivered via OEM file. A method based on the two-body state transition matrix (STM) mapping is used to perform this interpolation. The covariances at epochs which bracket the close approach time are linearly mapped to this time using the two-body STM: 

P(t௘) = Ψ௘௜P(t௜)Ψ௘௜

்

 (15)

Where Ψei represents the STM calculated to transform the states from the time being 

interpolated from, t௜, to the event time, t௘: 

Ψ௘௜ =

ௗxሬ⃑(௧೐)

ௗxሬ⃑(௧೔)

 

(16)

The weighted average of these mappings is then computed based on the relative proximity of 

the entry times to the close approach time: 

P(t௘) = ωP(t௘)௔ + [1 − ω]P(t௘)௕ (17)

Where 

ω = (t௘ − t௔)⁄(t௕ − t௔) (18)

And t௔ and t௕ represent the interpolated from times. 

This method was compared to three other candidate interpolation methods for two orbiters in 

the Mars environment: linear interpolation of EME2000 cartesian covariances, linear interpolation 

of covariances in a velocity-aligned frame (y-axis pointed in the spacecraft velocity direction, zaxis pointed in the angular momentum direction, and x-axis defined by the cross-product of y and 

z), and linear interpolation of EME2000 covariances of equinoctial conic elements. Comparisons 

were made to truth values at known plane-crossing events, interpolated from covariances generated 

at time steps from 1 minute to 30 minutes. The maximum error in computed 1σ uncertainties for 

conjunctions involving the analyzed orbiters are shown in Table 1 for covariances reported at these 

intervals. The STM mapping interpolation performs best among the tested methods and should 

hold valid across a range of orbits and missions. The two-body gravitational acceleration used in 

this method is dominant over other forces by orders of magnitude, and the interpolation smooths 

over small differences associated with higher fidelity models. If this approximation proves insufficient for certain time periods or missions, the performance can be improved by reducing the time 
 13

step. MADCAP does not require covariance data to be delivered in fixed time steps, so covariance 

data can be reported at finer spacing for different parts of an orbit if needed. 

Table 1. Maximum Interpolation Error in 1σ Uncertainty 

Interpolation 

Method 

Error 

Type 

1-min 

Interval 

2-min 

Interval 

5-min 

Interval 

10-min 

Interval 

30-min 

Interval 

EME2000 

Cartesian 

Radial (km) 1.3e+01 2.8e+01 6.9e+01 1.4e+02 4.1e+02 

Timing (s) 5.1e-02 2.1e-01 1.3e+00 5.2e+00 4.4e+01 

VelocityAligned 

Radial (km) 1.0e-04 1.0e-04 2.0e-03 6.0e-03 5.9e-02 

Timing (s) 7.0e-03 3.0e-02 1.9e-01 7.6e-01 7.4e+00 

EME2000 

Conic 

Radial (km) 3.0e-03 1.1e-02 7.2e-02 2.9e-01 1.1e+00 

Timing (s) 7.0e-03 1.8e-02 1.1e-01 3.9e-01 1.3e+00 

EME2000 

STM Map 

Radial (km) 2.7e-06 2.7e-06 5.1e-05 2.8e-06 1.8e-05 

Timing (s) 1.5e-07 1.9e-07 1.1e-04 2.3e-06 2.0e-05 

The above equations were derived with the assumption of plane-crossing events and do not hold 

in the case of coplanar orbiting bodies. A separate method was derived for use in such cases, the 

details of which are described in Reference 7. The performance of this method compares reasonably well with errors obtained via the plane-crossing method described above for non-coplanar 

events analyzed for a specific set of Mars orbiters. Future analysis is planned to ensure acceptable 

performance across a range of orbital environments and uncertainties. 

PROBABILITY OF COLLISION 

In the past, collision probability calculations were not performed for a few reasons. The main 

factor in not calculating PC is that covariance data has not been provided for most bodies in the 

Martian and Lunar environments. This is slowly changing as missions will begin to deliver covariance information for use in MADCAP. However, this covariance data is self-reported and thus of 

unknown methodology. It is possibly less accurate than the radar-based uncertainty information 

available for Earth orbiters. For example, a mission reporting overly conservative covariance data 

could lead to artificially low PC values. Nevertheless, mission teams are the most knowledgeable 

of their own uncertainties and will be more likely to accurately report information when the usage 

is explained. A final factor in the previous absence of reporting PC values is to avoid collision 

avoidance decisions being made based on one number without considering the larger context of the 

conjunction being analyzed. However, this can be prevented by clarifying to users that the PC is 

to be used as an additional conjunction event metric, and not the sole determining number to base 

decisions on. 

MADCAP now has the capability to report PC values when the appropriate covariance data is 

available from at least one of the bodies involved in a conjunction pair. In this section, we describe 

how PC is computed for the cases where covariance data is available for both spacecraft or for only 

one spacecraft (PC will not be computed if neither spacecraft has covariance information). The 
 14

flowchart in Figure 7 shows the process followed. We will first describe how covariance information is obtained from polynomials in the case that the formal covariance matrix is not available. 

Constructing Covariance from Polynomials 

When covariance data is not available, a pseudo-covariance matrix can be constructed from the 

‘Red’ threshold polynomials previously described and used for PC calculations. This can provide 

a more realistic probability than that for the “worst case” which can be artificially large, especially 

for large miss distances. However, it may conflict with the primary purpose of the polynomials: to 

be used as thresholds. Missions tend to choose conservative polynomials to reflect their worst case 

nominal uncertainties so that any conjunction events of concern will be flagged ‘Red’. Yet, artificially conservative covariances lead to artificially low PC values. To avoid inaccurate covariance 

data produced via polynomials, the default MADCAP setting is to not convert the polynomials into 

pseudo-covariance data since not reporting PC is considered better than reporting inaccurate information. The conversion of polynomials into pseudo-covariances can be turned on for each body 

independently. This allows the feature to only be used when the corresponding polynomials are 

known to provide more accurate PC information than the alternative. 

The pseudo-covariance matrix is constructed by using Equations (19-21) to compute the position 

uncertainties in the body velocity aligned frame. This frame is composed of a y-axis defined in the 

body velocity direction, z-axis in the angular momentum vector direction, and x-axis defined by 

the cross product of y and z (x - crosstrack, y - downtrack, z - normal). The uncertainty in the 

normal direction is assumed to be reasonably approximated by that in the crosstrack direction. This 

may not hold for all missions, in which case the pseudo-covariance matrix should not be used. 

Future MADCAP enhancements could allow for direct use of polynomial orbit-normal uncertainties provided by the spacecraft navigation teams. 

3σ஼ோைௌௌ்ோ஺஼௄ = OXD0 + OXD1 ∗ t + OXD2 ∗ t

ଶ

 (19)

3σ஽ைௐே்ோ஺஼௄ = [OXT0 + OXT1 ∗ t + OXT2 ∗ t

ଶ

]v̿ (20)

σேைோெ஺௅ = σ஼ோைௌௌ்ோ஺஼௄ (21)

The position errors can then be used to create a diagonal pseudo-covariance matrix as follows: 

P = ቎

σ஼௥௢௦௦௧௥௔௖௞

ଶ 0 0

0 σ஽௢௪௡௧௥௔௖௞

ଶ 0

0 0 σே௢௥௠௔௟

ଶ

቏ (22)

 This pseudo-covariance matrix can then be rotated into the Collison Frame and added to 

the position covariance matrix from the secondary body. The Collison Frame is defined by x-axis 

pointed in the relative miss direction (body1 to body2), y-axis pointed in the relative velocity direction, and z-axis defined by the cross-product of x and y. The x-z plane is known as the “collision 

plane”. 
 15

Figure 6. Collision Frame 

Calculating Probability of Collision 

In the case where covariance data for both spacecraft are available (either the formal or pseudocovariance constructed as described above), the covariance matrices and hard body radii from both 

bodies are combined through addition. Depending on the source of the covariance data, the matrices can either be added together in EME2000 frame and then rotated into the Collison Frame (see 

Figure 6), or rotated into the Collision Frame and then combined. Once in the Collision Frame, the 

problem is reduced into the 2 dimensions of the collision plane. This allows for the use of Foster’s 

method to calculate the probability of collision (described in detail in the referenced paper). 8

If covariance data is only available for one of the bodies, then the Frisbee Method is used to 

calculate the worst case 2D collision probability. 9 An upper bound on the PC can be obtained by 

assuming the 1σ error ellipse for the missing body is a degenerate ellipse (straight line) with semimajor axis equal to the miss distance. When linearly combined with the known covariance in the 

2D collision frame (x-z), this results in the worst case PC: 

P஼ைெ஻ூோ஽ = Pଵ + ൤

x଴

ଶ 0

0 0

൨ (23)

Where x଴ represents the miss distance. This problem is not well formed for large miss distances, 

(x଴ ≫ σ௫ଵ). MADCAP checks for this condition and does not calculate a PC if the miss distance 

is too large. The current value used for “large miss distance” in MADCAP is (x଴ > 10,000σ௫ଵ). 

This value was arrived at heuristically with future work necessary to determine more specifically 

when the assumption breaks down. 

The probability of collision is printed in the summary report for ‘Red’ events as an additional 

conjunction parameter. The probability value is followed by a descriptor pair of letters which explain what method was used to calculate the probability. If covariance data is available for both 

bodies, then Foster’s Method is used to calculate a 2D collision probability value, and “C-C” will 

be displayed. If covariance data is only available for one of the bodies, then the Frisbee Method is 

used to calculate the worst case 2D collision probability, and “C-N” or N-C” will be displayed. 

The order corresponds to the order the pair is specified in the first column, with “C” indicating that 

covariance data was available and used, while “N” indicates that no data was available. If no covariance data is available for either body, then the probability of collision cannot be calculated and 

“No Data” is displayed. If a polynomial derived covariance matrix is used, the letter “P” will be 

displayed in place of “C”. This process is shown graphically in Figure 7. 
 16

Figure 7. Flow Charts for a) Covariance Origin and b) PC Calculation 

CONCLUSION 

This paper has presented an update to the automated process used at NASA/JPL to conduct 

conjunction assessment at Mars, the Moon, and Sun-Earth Libration points. The process of deriving thresholds and collision probabilities from orbit uncertainties to screen conjunctions in these 

environments has been thoroughly described. The MADCAP screening and reporting process can 

only work as well as the inputs used to generate these thresholds and probabilities. The methods 

described here for threshold formation will enable current and future spacecraft navigation teams 

to provide the relevant data required for effective conjunction screening. This will ensure a safer 

orbital environment for all spacecraft operating in these increasingly crowded environments. 
 17

ACKNOWLEDGMENTS 

The authors would like to thank Dr. Shyam Bhaskaran for his comments and suggestions regarding the structure and organization of this paper. This work was carried out at the Jet Propulsion 

Laboratory, California Institute of Technology, Pasadena, California, under contract to the National 

Aeronautics and Space Administration. United States Government sponsorship acknowledged. All 

rights reserved. Reference herein to any specific commercial product, process, or service by trade 

name, trade-mark, manufacturer, or otherwise, does not constitute or imply its endorsement by the 

United States Government or the Jet Propulsion Laboratory, California Institute of Technology. © 

2022. California Institute of Technology. Government sponsorship acknowledged 

APPENDIX: EXAMPLE MOON SUMMARY REPORT 

An example summary report is displayed below (the ‘All’ table has been abridged for length): 

Analysis Time: 2021-12-31 16:47:32 UTC

RED Threshold Updates: 0

ALL Threshold Updates: 0

Ephemeris Updates: 4

Conjunction Assessment Bodies and Types

Body Name Type 

1 LRO Active 

1r LRO Active/Reference

2 ARTEMIS-P1 Active 

2r ARTEMIS-P1 Active/Reference

3 ARTEMIS-P2 Active 

3r ARTEMIS-P2 Active/Reference

4 CH2O Active 

5 Ouna Inactive 

6 CH1 Inactive 

Red (Conjunction Data < 'Red' Thresholds and Event < 14 days from 

Analysis Time)

Bodies OXD value/limit 

(km)

OXT value/limit 

(sec)

CAD value/limit 

(km)

Collision 

Probability

CA Epoch 

(UTC-SCET)

1-4 -0.4 3.2 4P 3.2 8.3 1P 163.4 - - no data 2022-01-13 

07:31:48
 18

1-4 -0.4 3.2 4P 3.2 8.3 1P 4.1 - - no data 2022-01-13 

08:31:26

All (Conjunction Data < 'All' Thresholds for <= 100 days)

Bodies OXD (km) OXT (sec) CAD (km) CA Epoch (UTC-SCET)

1-2 -152.5 -212.8 305.8 2022-01-03 10:52:35

1-4 -9.7 367.1 472.8 2022-01-03 23:57:37

1-4 -9.4 307.8 485.5 2022-01-04 00:55:55

1-4 -9.4 307.8 396.8 2022-01-04 01:55:18

1-4 -9.2 248.3 409.1 2022-01-04 02:53:36

Notes

OXD means "Orbit Crossing Distance". OXT means "Orbit Crossing Timing". CAD 

means "Close Approach Distance". 

Data for active spacecraft and natural bodies are displayed in the tables 

above. Data for inactive spacecraft are not displayed, but they are available 

in the conjunction metric tables and plots, which have been stored in the output directory listed below. Data for reference trajectories are not considered 

for Red events, but are considered in the All section for events not covered by 

the predicts file. Reference trajectories use the same thresholds as the nominal trajectories. 

The note after the probability value listed for RED events refers to the source 

of covariance data used for the probability calculation, with the following 

definitions: 

C - Covariance data provided by the mission in OEM format. 

P - A covariance converted from the RED threshold polynomial coefficients 

listed below. 

N - No covariance data, worst case covariance assumed for this body. 

No Data - Probability could not be calculated due to a lack of covariance data 

from both bodies. 

For more information, please see the point of contact listed below. 

Analysis time: 2021-12-31 16:47:32 UTC 

Active spacecraft: LRO, ARTEMIS-P1, ARTEMIS-P2, CH2O 

Natural bodies: None 

Inactive spacecraft: Ouna, CH1 

Output directory: /nav/home/jplmdnav/MADCAP/Moon/archive

Point of contact: MADCAP_Moon@jpl.nasa.gov

MADCAP build: 3.1.2 
 19

Red Thresholds -- Polynomial Coefficients

Body Name OXD0 

(km)

OXD1 

(km/t)

OXD2 

(km/t^2)

OXT0 

(sec)

OXT1 

(sec/t)

OXT2 

(sec/t^2)

1 LRO 0.1500 0.0125 0.0005 1.8750 0.2671 0.0184

2 ARTEMISP1 0.1370 -0.0024 0.0089 0.0100 2.8154 0.0045

3 ARTEMISP2 0.1370 -0.0024 0.0089 0.0100 2.8154 0.0045

4 CH2O 0.0000 0.2509 0.0000 0.0000 0.1490 0.0005

Red OX Distance Threshold = OXD0 + (OXD1 * t) + (OXD2 * t^2) [km] 

Red OX Timing Threshold = OXT0 + (OXT1 * t) + (OXT2 * t^2) [sec] 

where t = CA Epoch - Ephemeris File Submit Time (in days) 

Red thresholds are based on 3-sigma values. Thresholds listed as "P" are based 

on a quadratic fit of the 3-sigma values as a function of time to the event. 

The polynomial coefficients used are listed in the table above. Lines for coefficients which have been updated since the last run are colored blue, and each 

line's body is marked with an "*". Thresholds listed as "C" are based on 3-

sigma covariance data provided by the mission. 

All Thresholds -- Constants

Body Name OXD (km) CAD (km)

1 LRO 1 40

2 ARTEMIS-P1 500 500

3 ARTEMIS-P2 500 500

4 CH2O 500 500

All OX Distance Threshold = OXD 

All CA Distance Threshold = CAD 

All thresholds are always constants. The constants used are listed in the table 

above. Lines for constants which have been updated since the last run are colored blue, and each line's body is marked with an "*". 

Ephemerides

Body Ephemeris Submitted Begin End 

1* 28day_20211231_01.bsp 2021-12-31 

12:34:41 UTC 

31-DEC-2021 

00:00:00 UTC 

28-JAN-2022 

00:00:00 UTC 

1r 558day_20211215_01.bsp 2021-12-15 

13:39:10 UTC 

15-DEC-2021 

00:00:00 UTC 

26-JUN-2023 

00:00:00 UTC 
 20

2* 192.THEMIS_B.SHORT_TERM.2021_365

.oem.bsp_V0.1 

2021-12-31 

10:41:59 UTC 

31-DEC-2021 

00:00:00 UTC 

30-JAN-2022 

00:00:00 UTC 

2r 192.THEMIS_B.LONG_TERM.2021_104.

oem.bsp_V0.1 

2021-04-14 

17:57:45 UTC 

14-APR-2021 

00:00:00 UTC 

13-APR-2027 

00:00:00 UTC 

3* 193.THEMIS_C.SHORT_TERM.2021_365

.oem.bsp_V0.1 

2021-12-31 

10:43:08 UTC 

31-DEC-2021 

00:00:00 UTC 

30-JAN-2022 

00:00:00 UTC 

3r 193.THEMIS_C.LONG_TERM.2021_104.

oem.bsp_V0.1 

2021-04-14 

18:01:40 UTC 

14-APR-2021 

00:00:00 UTC 

13-APR-2027 

00:00:00 UTC 

4* ISRO-CH2-2021-12-29-OD893-365-

v1.xsp.bsp 

2021-12-31 

12:08:16 UTC 

29-DEC-2021 

22:00:00 UTC 

19-JAN-2022 

00:00:00 UTC 

5 ouna_191201_230101_150608_SMM07100

31456-jpl-ekl.bsp 

Analysis 

Time 

01-DEC-2019 

00:59:23 UTC 

31-DEC-2022 

23:58:23 UTC 

6 traj_ch1_010920-011022.bsp Analysis 

Time 

31-AUG-2020 

23:58:50 UTC 

30-SEP-2022 

23:58:50 UTC 

Ephemeris files for the bodies analyzed are listed in the table above. Lines for 

files which have been updated since the last run are colored blue, and each line's 

body is marked

1

 NASA Conjunction Assessment Risk Analysis, URL: https://satellitesafety.gsfc.nasa.gov/CARA.html 

2

 Newman, L.K. "Conjunction Assessment Risk Analysis (CARA) Program Overview." https://foiaelibrary.gsfc.nasa.gov/_assets/doclibBidder/tech_docs/CARA%20overview%20briefing.pdf, 7 August 2013 

3

 Newman, Lauri Kraft, "The NASA Robotic Conjunction Assessment Process: Overview and Operational Experiences", 

Proceedings of the 59th International Astronautical Congress 2008, IAC-08-A.6.2.6 

4

 Berry, D. S., Guinn, J. R., Tarzi, Z. B., & Demcak, S. W., “Automated Spacecraft Conjunction at Mars and the Moon” 

Paper 1262983, 12th International Conference on Space Operations, Stockholm, Sweden, June 2012 

5

 Tarzi, Z., Berry, D., Roncoli, R., “An Updated Process for Automated Deepspace Conjunction Assessment” Paper AAS 

14-373, 25th International Symposium on Space Flight Dynamics, Munich, Germany, October 2015 

6

 Berry, D., Tarzi, Z., Roncoli, R., Wilson, R., “Automated Spacecraft Conjunction Assessment at Mars and the MoonA Five Year Update” Paper 2810036, 14th International Conference on Space Operations, Marseille, France, May 2018 

7

 Young, B., “Application of Filter-Derived Covariance to MADCAP”, JPL IOM 392C-20-007, July 2020 

8

 Foster, J.L. and Estes, H.S., “A Parametric Analysis of Orbital Debris Collision Probability and Maneuver Rate for 

Space Vehicles,” NASA JSC 25898, August 1992 

9

 Frisbee, J., “An Upper Bound on High Speed Satellite Collision Probability When Only One Object Has Position Uncertainty Information”, Paper AAS 15-717, AAS/AIAA Astrodynamics Specialist Conference, August 09, 2015 - August 

13, 2015, Vail, CO 

10 SPK Required Reading, https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/FORTRAN/req/spk.html 

11 Orbit Data Messages. Issue 2. Recommendation for Space Data System Standards (Blue Book), CCSDS 502.0-B-2. 

Washington, D.C.: CCSDS, November 2009, https://public.ccsds.org/Pubs/502x0b2c1e2.pdf 

12 Evans, S., Monte: The Next Generation of Mission Design & Navigation Software, The 6th International Conference 

on Astrodynamics Tools and Techniques (ICATT) proceedings 2016, Darmstadt, Germany. 

REFERENCES