simulation.h

Functions

std::vector<std::vector<real>> reconstruct_stm(const std::vector<real> &stm)

Reconstruct the STM matrix from the flattened STM vector.

Parameters:

stm – [in] The flattened STM vector.

Returns:

The reconstructed STM matrix.

void get_baseBodyFrame(const int &spiceId, const real &tMjdTDB, std::string &baseBodyFrame)

Get the name of body-fixed frame for a given SPICE ID.

Parameters:

• spiceId – [in] SPICE ID of the body.

• tMjdTDB – [in] Time of observation in Modified Julian Date (TDB).

• baseBodyFrame – [out] Name of the body-fixed frame.

void get_observer_state(const real &tObsMjd, const std::vector<real> &observerInfo, PropSimulation *propSim, const bool tObsInUTC, std::vector<real> &observerState)

Get the observer state for a given time.

Parameters:

• tObsMjd – [in] Observation time in Modified Julian Date.

• observerInfo – [in] Observer information (base body SPICE ID, lon [rad], lat [rad], alt [m])

• propSim – [in] PropSimulation object for querying Ephemeris of base body.

• tObsInUTC – [in] Flag to indicate if the observation time is in UTC (true) or TDB (false).

• observerState – [out] Output observer state (position [DU], velocity [DU/TU]).

struct Constants

#include <simulation.h>

Structure to hold constants used in a PropSimulation.

Param du2m:

Distance unit conversion factor (default: AU to meters).

Param tu2s:

Time unit conversion factor (default: days to seconds).

Param duptu2mps:

Distance unit per time unit conversion factor (default: AU/day to m/s).

Param G:

Gravitational constant (default: kg AU^3/day^2).

Param clight:

Speed of light (default: AU/day).

Param j2000Jd:

Julian date of J2000 epoch.

Param JdMinusMjd:

Julian date minus modified Julian date.

Public Members

real du2m = 149597870700.0L

real tu2s = 86400.0

real duptu2mps = du2m / tu2s

real G = 6.6743e-11L / (du2m * du2m * du2m) * tu2s * tu2s

real clight = 299792458.0L / du2m * tu2s

real j2000Jd = 2451545.0

real JdMinusMjd = 2400000.5

struct IntegrationParameters

#include <simulation.h>

Structure to hold integration parameters used in a PropSimulation.

Param nInteg:

Number of integrated bodies.

Param nSpice:

Number of SPICE bodies (queried from PropSimulation Ephemeris object).

Param nTotal:

Total number of bodies (nInteg + nSpice).

Param n2Derivs:

Number of second derivatives (usually 3 for each body unless STM is propagated).

Param t0:

Initial time.

Param tf:

Final time.

Param dt0:

Initial timestep.

Param dtMin:

Minimum timestep.

Param dtChangeFactor:

Maximum factor by which to change timestep.

Param adaptiveTimestep:

Flag to use adaptive timestep.

Param timestepCounter:

Timestep counter.

Param tolPC:

Tolerance for Gauss-Radau predictor-corrector loop.

Param tolInteg:

Tolerance for integrator.

Public Members

size_t nInteg

size_t nSpice

size_t nTotal

size_t n2Derivs

real t0

real tf

real dt0

real dtMin

real dtChangeFactor

bool adaptiveTimestep

size_t timestepCounter

real tolPC

real tolInteg

class Body

#include <simulation.h>

Parent class for all bodies in a PropSimulation.

Param t0:

Initial time.

Param mass:

Mass of the body.

Param radius:

Radius of the body.

Param J2:

J2 coefficient for oblateness.

Param poleRA:

Right ascension of the pole.

Param poleDec:

Declination of the pole.

Param name:

Name of the body.

Param spiceId:

SPICE ID of the body.

Param pos:

Position of the body.

Param vel:

Velocity of the body.

Param acc:

Acceleration of the body.

Param isPPN:

Flag to indicate if the body is a PPN body (for EIH relativity loop 1).

Param isJ2:

Flag to indicate if the body has J2 oblateness.

Param isNongrav:

Flag to indicate if the body has nongravitational forces.

Param isMajor:

Flag to indicate if the body is a major body (for EIH relativity loop 2).

Param caTol:

Distance tolerance for close approach detection (default: 0.1 AU).

Subclassed by IntegBody, SpiceBody

Public Functions

void set_J2(real J2, real poleRA, real poleDec)

Set the J2 for the body.

Parameters:

• J2 – [in] Oblateness coefficient.

• poleRA – [in] Right ascension of the pole.

• poleDec – [in] Declination of the pole.

void set_harmonics(real poleRA, real poleDec, int nZon, int nTes, std::vector<real> J, std::vector<std::vector<real>> C, std::vector<std::vector<real>> S)

Set the harmonic coefficients for the body.

Parameters:

• poleRA – [in] Right ascension of the pole.

• poleDec – [in] Declination of the pole.

• nZon – [in] Degree of the zonal coefficients.

• nTes – [in] Order of the tesseral coefficients.

• J – [in] Zonal coefficient vector.

• C – [in] Sectoral coefficient array.

• S – [in] Tesseral coefficient array.

Public Members

real t0

real mass

real radius

real J2 = 0.0L

real poleRA = 0.0L

real poleDec = 90.0L

size_t nZon = 0

size_t nTes = 0

std::vector<real> J = {}

std::vector<std::vector<real>> C = {}

std::vector<std::vector<real>> S = {}

std::string name

int spiceId

real pos[3]

real vel[3]

real acc[3]

bool isPPN = false

bool isJ2 = false

bool isHarmonic = false

bool isNongrav = false

bool isMajor = false

real caTol = 0.1

class SpiceBody : public Body

#include <simulation.h>

Class for SPICE bodies in a PropSimulation.

Param isSpice:

Flag to indicate if the body is a SPICE body.

Public Functions

SpiceBody(std::string name, int spiceId, real t0, real mass, real radius)

Construct a new SpiceBody object.

Parameters:

• name – [in] Name of the body.

• spiceId – [in] SPICE ID of the body.

• t0 – [in] Initial time [MJD TDB].

• mass – [in] Mass of the body [kg].

• radius – [in] Radius of the body [m].

Public Members

bool isSpice = true

struct NongravParameters

#include <simulation.h>

Structure to hold parameters for asteroi/comet nongravitational forces.

Param a1:

Parameter for radial nongravitational forces.

Param a2:

Parameter for transverse nongravitational forces.

Param a3:

Parameter for normal nongravitational forces.

Param a1Est:

Flag to indicate if a1 is estimated (for orbit determination).

Param a2Est:

Flag to indicate if a2 is estimated (for orbit determination).

Param a3Est:

Flag to indicate if a3 is estimated (for orbit determination).

Param alpha:

Normalizing factor for nongravitational forces.

Param k:

Exponent for nongravitational forces.

Param m:

Exponent for nongravitational forces.

Param n:

Exponent for nongravitational forces.

Param r0_au:

Normalizing distance for nongravitational forces.

Public Members

real a1 = 0.0L

real a2 = 0.0L

real a3 = 0.0L

bool a1Est = false

bool a2Est = false

bool a3Est = false

real alpha = 0.1112620426L

real k = 4.6142L

real m = 2.15L

real n = 5.093L

real r0_au = 2.808L

class IntegBody : public Body

#include <simulation.h>

Class for integrated bodies in a PropSimulation.

Param spiceId:

SPICE ID of the body.

Param logCA:

Flag to indicate if close approaches are logged (default: true).

Param isCometary:

Flag to indicate if the body is initialized with a cometary state (as opposed to Cartesian state).

Param initState:

Initial state of the body.

Param isInteg:

Flag to indicate if the body is integrated.

Param isThrusting:

Flag to indicate if the body is thrusting.

Param ngParams:

Nongravitational parameters for the body.

Param n2Derivs:

Number of second derivatives for the body.

Param propStm:

Flag to indicate if the state transition matrix is propagated.

Param stm:

State transition matrix.

Param dCartdState:

Derivatives of initial Cartesian state with respect to initial state.

Public Functions

IntegBody(std::string name, real t0, real mass, real radius, std::vector<real> cometaryState, NongravParameters ngParams)

Construct a new IntegBody object.

Parameters:

• name – [in] Name of the body.

• t0 – [in] Initial time.

• mass – [in] Mass of the body [kg].

• radius – [in] Radius of the body [m].

• cometaryState – [in] Initial heliocentric ecliptic cometary state of the body.

• ngParams – [in] Nongravitational parameters for the body.

IntegBody(std::string name, real t0, real mass, real radius, std::vector<real> pos, std::vector<real> vel, NongravParameters ngParams)

Construct a new IntegBody object.

Parameters:

• name – [in] Name of the body.

• t0 – [in] Initial time.

• mass – [in] Mass of the body [kg].

• radius – [in] Radius of the body [m].

• pos – [in] Initial barycentric position of the body [AU].

• vel – [in] Initial barycentric velocity of the body [AU/day].

• ngParams – [in] Nongravitational parameters for the body.

void prepare_stm()

Prepare the state transition matrix for propagation.

Public Members

int spiceId = -99999

bool logCA = true

bool isCometary = false

std::vector<real> initState

std::vector<real> initCart

bool isInteg = true

bool isThrusting = false

NongravParameters ngParams

size_t n2Derivs = 3

bool propStm = false

std::vector<real> stm

std::vector<std::vector<real>> dCartdState

class Event

#include <simulation.h>

Class for events in a PropSimulation.

Param t:

Time of the event.

Param bodyName:

Name of the IntegBody for the event

Param bodyIndex:

Index of the IntegBody in the PropSimulation.

Param xIntegIndex:

Starting index of the body’s state in the flattened state vector.

Param deltaV:

Delta-V for the event.

Param multiplier:

Multiplier for the Delta-V.

Param dt:

Time duration for the continuous event.

Param isContinuous:

Flag to indicate if the event is continuous.

Param isHappening:

Flag to indicate if the continuous event is happening.

Public Functions

inline Event()

Empty constructor for the Event class.

void apply_impulsive(PropSimulation *propSim, const real &t, std::vector<real> &xInteg)

Apply the impulse event to the body.

Parameters:

• propSim – [in] PropSimulation object.

• t – [in] Time of the event.

• xInteg – [inout] State of the body.

Public Members

real t

std::string bodyName

bool isContinuous = false

bool eventEst = false

size_t bodyIndex

size_t xIntegIndex

bool hasStarted = false

std::vector<real> deltaV = std::vector<real>(3, std::numeric_limits<real>::quiet_NaN())

real multiplier = std::numeric_limits<real>::quiet_NaN()

bool deltaVEst = false

bool multiplierEst = false

std::vector<real> expAccel0 = std::vector<real>(3, std::numeric_limits<real>::quiet_NaN())

real tau = std::numeric_limits<real>::quiet_NaN()

bool expAccel0Est = false

bool tauEst = false

class EventManager

#include <simulation.h>

Public Members

std::vector<Event> impulsiveEvents = {}

std::vector<Event> continuousEvents = {}

size_t nextImpEventIdx = 0

size_t nextConEventIdx = 0

real tNextImpEvent = std::numeric_limits<real>::quiet_NaN()

real tNextConEvent = std::numeric_limits<real>::quiet_NaN()

size_t nImpEvents = 0

size_t nConEvents = 0

bool allConEventDone = true

struct BPlaneParameters

#include <simulation.h>

Structure to hold parameters for a close approach B-plane.

Param x:

X-coordinate of the B-plane.

Param y:

Y-coordinate of the B-plane.

Param z:

Z-coordinate of the B-plane (usually 0).

Param dx:

Derivatives of x with respect to the state at the close approach.

Param dy:

Derivatives of y with respect to the state at the close approach.

Public Members

real x

real y

real z

std::vector<real> dx = std::vector<real>(6, std::numeric_limits<real>::quiet_NaN())

std::vector<real> dy = std::vector<real>(6, std::numeric_limits<real>::quiet_NaN())

class CloseApproachParameters

#include <simulation.h>

Class for close approach parameters in a PropSimulation.

Param t:

Time of the close approach or impact (from child class).

Param xRel:

Relative state at the close approach or impact (from child class).

Param tCA:

Time of the close approach.

Param xRelCA:

Relative state at the close approach.

Param tMap:

Time of the B-plane map.

Param xRelMap:

Relative state at the B-plane map time.

Param dist:

Distance at the close approach.

Param vel:

Relative velocity at the close approach.

Param vInf:

Hyperbolic excess velocity of the close approach.

Param flybyBody:

Name of the flyby body.

Param flybyBodyIdx:

Index of the flyby body in the PropSimulation.

Param centralBody:

Name of the central body.

Param centralBodyIdx:

Index of the central body in the PropSimulation.

Param centralBodySpiceId:

SPICE ID of the central body.

Param impact:

Flag to indicate if the close approach is an impact.

Param tPeri:

Time of periapsis.

Param tLin:

Linearized time of intersection on the B-plane.

Param bVec:

B-vector of the B-plane.

Param bMag:

Magnitude of the B-vector.

Param gravFocusFactor:

Gravitational focus factor to account for nonlinear effects from central body.

Param kizner:

Kizner B-plane parameters.

Param opik:

Öpik B-plane parameters.

Param scaled:

Scaled Kizner B-plane parameters.

Param mtp:

Modified Target Plane parameters.

Param dtLin:

Partial derivatives of the (linearized intersection minus map) time with respect to the state at the close approach.

Param dt:

Partial derivatives of the time of closest approach with respect to the state at the close approach.

Subclassed by ImpactParameters

Public Functions

void get_ca_parameters(PropSimulation *propSim, const real &tMap)

Get the close approach parameters.

Parameters:

• propSim – [in] PropSimulation object.

• tMap – [in] Time of the B-plane map.

void print_summary(int prec = 8)

Print a summary of the close approach parameters.

Parameters:

prec – [in] Precision of the output.

Public Members

real t

std::vector<real> xRel

real tCA

std::vector<real> xRelCA

real tMap

std::vector<real> xRelMap

real dist

real vel

real vInf

std::string flybyBody

int flybyBodyIdx

std::string centralBody

int centralBodyIdx

int centralBodySpiceId

bool impact

real tPeri

real tLin

std::vector<real> bVec = std::vector<real>(3, 0.0L)

real bMag

real gravFocusFactor

BPlaneParameters kizner

BPlaneParameters opik

BPlaneParameters scaled

BPlaneParameters mtp

std::vector<real> dtLin = std::vector<real>(6, 0.0L)

std::vector<real> dt = std::vector<real>(6, 0.0L)

class ImpactParameters : public CloseApproachParameters

#include <simulation.h>

Class for impact parameters in a PropSimulation.

Param xRelBodyFixed:

Relative state at the impact time in body-fixed frame.

Param lon:

Longitude of the impact [rad].

Param lat:

Latitude of the impact [rad].

Param alt:

Altitude of the impact [km].

Public Functions

void get_impact_parameters(PropSimulation *propSim)

Get the impact parameters.

Parameters:

propSim – [in] PropSimulation object.

void print_summary(int prec = 8)

Print a summary of the impact parameters.

Parameters:

prec – [in] Precision of the output.

Public Members

std::vector<real> xRelBodyFixed = std::vector<real>(6, std::numeric_limits<real>::quiet_NaN())

real lon

real lat

real alt

struct InterpolationParameters

#include <simulation.h>

Structure to hold interpolation parameters for a PropSimulation.

Param interpIdx:

Index of the interpolation parameters.

Param t0:

Initial time of the interpolation window.

Param dt0:

Timespan of the interpolation window.

Param b0:

Interpolation coefficients for the interpolation window.

Param xInteg0:

States at the initial time of the interpolation window.

Param accInteg0:

Accelerations at the initial time of the interpolation window.

Param tStack:

Stack of integrator epochs.

Param xIntegStack:

Stack of states at integrator epochs.

Param bStack:

Stack of interpolation coefficients at integrator epochs.

Param accIntegStack:

Stack of accelerations at integrator epochs.

Public Members

size_t interpIdx = 0

real t0

real dt0

std::vector<real> b0

std::vector<real> xInteg0

std::vector<real> accInteg0

std::vector<real> tStack

std::vector<std::vector<real>> xIntegStack

std::vector<std::vector<real>> bStack

std::vector<std::vector<real>> accIntegStack

class PropSimulation

#include <simulation.h>

Class for a PropSimulation.

Param name:

Name of the simulation.

Param DEkernelPath:

Path to DE kernels.

Param spkEphem:

SPK ephemeris object for the simulation.

Param pckEphem:

PCK ephemeris object for the simulation.

Param unsafePersistentMemoryMap:

Flag to indicate if ephemeris memory maps should not be unmapped (unsafe, default: false).

Param consts:

Constants object for the simulation.

Param integParams:

Integration parameters for the simulation.

Param spiceBodies:

Vector of SpiceBody objects in the simulation.

Param integBodies:

Vector of IntegBody objects in the simulation.

Param events:

Vector of ImpulseEvent objects in the simulation.

Param caParams:

Vector of CloseApproachParameters objects in the simulation.

Param impactParams:

Vector of ImpactParameters objects in the simulation.

Param isPreprocessed:

Flag to indicate if the simulation is preprocessed.

Param t:

Current time of the simulation.

Param xInteg:

Current state of the simulation.

Param interpParams:

Interpolation parameters for the simulation.

Param tEvalUTC:

Flag to indicate if the evaluation times are in UTC.

Param evalApparentState:

Flag to indicate if the apparent state is evaluated.

Param evalMeasurements:

Flag to indicate if measurements are evaluated.

Param convergedLightTime:

Flag to indicate if converged light time needs to be computed.

Param xObserver:

Observer states array.

Param observerInfo:

Observer information array.

Param tEvalMargin:

Margin for interpolation outside integration arc.

Param tEval:

Vector of times at which to evaluate the integrated state.

Param obsType:

Vector of observation type flags (0=none, 1=delay, 2=doppler, 3=gaia).

Param lightTimeEval:

Vector of computed light times.

Param xIntegEval:

Vector of integrated states at evaluation times.

Param opticalObs:

Vector of optical observations.

Param opticalPartials:

Vector of optical observation partials.

Param opticalObsCorr:

Vector of photocenter-barycenter corrections for optical observations.

Param radarObs:

Vector of radar observations.

Param radarPartials:

Vector of radar observation partials.

Public Functions

PropSimulation(std::string name, real t0, const int defaultSpiceBodies, std::string DEkernelPath)

Construct a new PropSimulation object.

Parameters:

• name – [in] Name of the simulation.

• t0 – [in] Initial time.

• defaultSpiceBodies – [in] Default SPICE bodies to load. (0=empty sim with DE440 ephemeris, 430/421=DE430 Sun+planets+Moon+Pluto+big16 asteroids, 440/441=DE440 Sun+planets+Moon+Pluto+big16 asteroids)

• DEkernelPath – [in] Path to SPICE metakernel

PropSimulation(std::string name, const PropSimulation &simRef)

Construct a new PropSimulation object from a reference simulation.

Parameters:

• name – [in] Name of the simulation.

• simRef – [in] Reference simulation to copy.

void map_ephemeris()

Memory map the ephemeris files.

void unmap_ephemeris()

Unmap the ephemeris files.

std::vector<real> get_spiceBody_state(const real t, const std::string &bodyName)

Get the state of a SpiceBody in the PropSimulation at a given time.

Parameters:

• t – [in] Time at which to get the state.

• bodyName – [in] Name of the body.

Returns:

std::vector<real> State of the body at time t.

std::vector<real> interpolate(const real t)

Get the integrated state from the simulation at a given time.

Parameters:

t – [in] Time to interpolate the integrated state from the PropSimulation at.

Returns:

std::vector<real> The interpolated geometric barycentric state at time t.

void add_spice_body(SpiceBody body)

Add a SpiceBody to the simulation.

Parameters:

body – [in] SpiceBody object to add to the simulation.

void add_integ_body(IntegBody body)

Add an IntegBody to the simulation.

Parameters:

body – [in] IntegBody object to add to the simulation.

void remove_body(std::string name)

Remove a body from the simulation.

Parameters:

name – [in] Name of the body to remove.

void add_event(Event event)

Add an event to the simulation.

Parameters:

event – [in] Event object to add to the simulation.

void set_sim_constants(real du2m = 149597870700.0L, real tu2s = 86400.0L, real G = 6.6743e-11L / (149597870700.0L * 149597870700.0L * 149597870700.0L) * 86400.0L * 86400.0L, real clight = 299792458.0L / 149597870700.0L * 86400.0L)

Set the values of the PropSimulation Constants object.

Parameters:

• du2m – [in] Distance unit conversion factor (default: AU to meters).

• tu2s – [in] Time unit conversion factor (default: days to seconds).

• G – [in] Gravitational constant (default: kg AU^3/day^2).

• clight – [in] Speed of light (default: AU/day).

void set_integration_parameters(real tf, std::vector<real> tEval = std::vector<real>(), bool tEvalUTC = false, bool evalApparentState = false, bool convergedLightTime = false, std::vector<std::vector<real>> observerInfo = std::vector<std::vector<real>>(), bool adaptiveTimestep = true, real dt0 = 1.0L, real dtMin = 1.0e-4L, real dtChangeFactor = 0.25L, real tolInteg = 1.0e-11L, real tolPC = 1.0e-16L)

Set the values of the PropSimulation IntegrationParameters object.

Parameters:

• tf – [in] Final time.

• tEval – [in] Vector of times at which to evaluate the integrated state (default: empty).

• tEvalUTC – [in] Flag to indicate if the evaluation times are in UTC (default: false).

• evalApparentState – [in] Flag to indicate if the apparent state is evaluated (default: false).

• convergedLightTime – [in] Flag to indicate if converged light time needs to be computed (default: false).

• observerInfo – [in] Observer information array (default: empty).

• adaptiveTimestep – [in] Flag to use adaptive timestep (default: true).

• dt0 – [in] Initial timestep (default: 0.0).

• dtMin – [in] Minimum timestep (default: 0.005).

• dtChangeFactor – [in] Maximum factor by which to change timestep (default: 0.25).

• tolInteg – [in] Tolerance for integrator (default: 1e-11).

• tolPC – [in] Tolerance for Gauss-Radau predictor-corrector loop (default: 1e-16).

std::vector<real> get_sim_constants()

Get the constants used in the simulation.

Returns:

std::vector<real> Vector of simulation constants.

std::vector<real> get_integration_parameters()

Get the integration parameters used in the simulation.

Returns:

std::vector<real> Vector of simulation integration parameters.

void integrate()

Integrate the simulation to the final time.

void extend(real tf, std::vector<real> tEvalNew = std::vector<real>(), std::vector<std::vector<real>> observerInfoNew = std::vector<std::vector<real>>())

Extend the simulation to a new final time.

Parameters:

• tf – [in] New final time.

• tEvalNew – [in] New vector of times at which to evaluate the integrated state.

• observerInfoNew – [in] New vector of observer information.

void save(std::string filename)

Save the simulation to a file.

Parameters:

filename – [in] Name of the file to save the simulation to.

Public Members

std::string name

std::string DEkernelPath

SpkEphemeris spkEphem

PckEphemeris pckEphem

bool unsafePersistentMemoryMap = false

Constants consts

IntegrationParameters integParams

std::vector<SpiceBody> spiceBodies

std::vector<IntegBody> integBodies

EventManager eventMngr

std::vector<CloseApproachParameters> caParams

std::vector<ImpactParameters> impactParams

real t

std::vector<real> xInteg

InterpolationParameters interpParams

bool tEvalUTC = false

bool evalApparentState = false

bool evalMeasurements = false

bool convergedLightTime = false

std::vector<std::vector<real>> xObserver

std::vector<std::vector<real>> observerInfo

real tEvalMargin = 0.0L

std::vector<real> tEval

std::vector<int> obsType

std::vector<std::vector<real>> lightTimeEval

std::vector<std::vector<real>> xIntegEval

std::vector<std::vector<real>> opticalObs

std::vector<std::vector<real>> opticalObsDot

std::vector<std::vector<real>> opticalPartials

std::vector<std::vector<real>> opticalObsCorr

std::vector<std::vector<real>> radarObs

std::vector<std::vector<real>> radarPartials

Private Functions

void prepare_for_evaluation(std::vector<real> &tEval, std::vector<std::vector<real>> &observerInfo)

Prepare the simulation for evaluation.

Parameters:

• tEval – [in] Vector of times at which to evaluate the integrated state.

• observerInfo – [in] Observer information array.

void preprocess()

Preprocess the simulation before integration.

Private Members

bool isPreprocessed = false