simulation.h

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#