updated method for adding events; continuous events analytic estimation

grss/fit/fit_simulation.py

@@ -1030,7 +1030,7 @@ def _x_dict_to_events(self, x_dict):
             fixed_event = tuple(self.fixed_propsim_params['events'][i])
             fixed_event_impulsive = len(fixed_event) == 5
             fixed_event_continuous = len(fixed_event) == 6
-            if fixed_event_continuous and fixed_event[5] == 0.0:
+            if fixed_event_continuous and fixed_event[4] == 0.0:
                 raise ValueError("Continuous event time constant cannot be 0.")
             assert fixed_event_impulsive or fixed_event_continuous
             event = [None]*5
@@ -1039,12 +1039,13 @@ def _x_dict_to_events(self, x_dict):
                 event[1] = x_dict[f"dvx{i}"] if f"dvx{i}" in x_dict.keys() else fixed_event[1]
                 event[2] = x_dict[f"dvy{i}"] if f"dvy{i}" in x_dict.keys() else fixed_event[2]
                 event[3] = x_dict[f"dvz{i}"] if f"dvz{i}" in x_dict.keys() else fixed_event[3]
+                event[4] = x_dict[f"mult{i}"] if f"mult{i}" in x_dict.keys() else fixed_event[4]
             elif fixed_event_continuous:
                 event[1] = x_dict[f"ax{i}"] if f"ax{i}" in x_dict.keys() else fixed_event[1]
                 event[2] = x_dict[f"ay{i}"] if f"ay{i}" in x_dict.keys() else fixed_event[2]
                 event[3] = x_dict[f"az{i}"] if f"az{i}" in x_dict.keys() else fixed_event[3]
-                event.append(x_dict[f"dt{i}"] if f"dt{i}" in x_dict.keys() else fixed_event[5])
-            event[4] = x_dict[f"mult{i}"] if f"mult{i}" in x_dict.keys() else fixed_event[4]
+                event[4] = x_dict[f"dt{i}"] if f"dt{i}" in x_dict.keys() else fixed_event[4]
+                event.append(fixed_event[5])
             events.append(tuple(event))
         return events
 
@@ -1070,20 +1071,42 @@ def _check_and_add_events(self, prop_sim_past, prop_sim_future, integ_body, even
         prop_sim_future : libgrss.PropSimulation object
             PropSimulation object for the future.
         """
-        for event in events:
-            t_event = event[0]
-            dvx_or_ax = event[1]
-            dvy_or_ay = event[2]
-            dvz_or_az = event[3]
-            multiplier = event[4]
-            event_args = (integ_body, t_event, [dvx_or_ax, dvy_or_ay, dvz_or_az], multiplier)
-            if len(event) == 6:
-                dt = event[5]
-                event_args += (dt,)
+        est_keys = self.x_nom.keys()
+        for i, event_list in enumerate(events):
+            t_event = event_list[0]
+            dvx_or_ax = event_list[1]
+            dvy_or_ay = event_list[2]
+            dvz_or_az = event_list[3]
+            multiplier_or_dt = event_list[4]
+            continuous_flag = event_list[5] if len(event_list) == 6 else False
+            event = libgrss.Event()
+            event.t = t_event
+            event.bodyName = integ_body.name
+            event.isContinuous = continuous_flag
+            if continuous_flag:
+                event.expAccel0 = [dvx_or_ax, dvy_or_ay, dvz_or_az]
+                event.tau = multiplier_or_dt
+                event.expAccel0Est = (
+                    f"ax{i}" in est_keys and
+                    f"ay{i}" in est_keys and
+                    f"az{i}" in est_keys
+                )
+                event.tauEst = f"dt{i}" in est_keys
+                event.eventEst = event.expAccel0Est or event.tauEst
+            else:
+                event.deltaV = [dvx_or_ax, dvy_or_ay, dvz_or_az]
+                event.multiplier = multiplier_or_dt
+                event.deltaVEst = (
+                    f"dvx{i}" in est_keys and
+                    f"dvy{i}" in est_keys and 
+                    f"dvz{i}" in est_keys
+                )
+                event.multiplierEst = f"mult{i}" in est_keys
+                event.eventEst = event.deltaVEst or event.multiplierEst
             if t_event < self.t_sol:
-                prop_sim_past.add_event(*event_args)
+                prop_sim_past.add_event(event)
             else:
-                prop_sim_future.add_event(*event_args)
+                prop_sim_future.add_event(event)
         return prop_sim_past, prop_sim_future
 
     def _get_perturbed_state(self, key):

grss/version.txt

@@ -1 +1 @@
-4.3.1
+4.3.2

include/simulation.h

@@ -234,19 +234,27 @@ class Event {
    public:
     real t;
     std::string bodyName;
+    bool isContinuous = false;
+    bool eventEst = false;
     size_t bodyIndex;
     size_t xIntegIndex;
+    bool hasStarted = false;
 
     // for impulsive events
-    std::vector<real> deltaV = {0.0L, 0.0L, 0.0L};
-    real multiplier = 1.0L;
+    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;
 
     // for continuous ejecta events
-    std::vector<real> expAccel0 = {0.0L, 0.0L, 0.0L};
-    real tau = 1.0L;
-    bool isContinuous = false;
-    bool isHappening = 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;
+    /**
+     * @brief Empty constructor for the Event class.
+     */
+    Event() {};
     /**
      * @brief Apply the impulse event to the body.
      */
@@ -517,15 +525,9 @@ class PropSimulation {
      */
     void remove_body(std::string name);
     /**
-     * @brief Add an impulse event to the simulation.
-     */
-    void add_event(IntegBody body, real tEvent, std::vector<real> deltaV,
-                   real multiplier);
-    /**
-     * @brief Add an ejecta event to the simulation.
+     * @brief Add an event to the simulation.
      */
-    void add_event(IntegBody body, real tEvent, std::vector<real> expAccel0,
-                   real multiplier, real tau);
+    void add_event(Event event);
     /**
      * @brief Set the values of the PropSimulation Constants object.
      */

include/stm.h

@@ -71,4 +71,10 @@ void stm_nongrav(STMParameters &stmParams, const real &g,
                  const real &dy, const real &dz, const real &dvx,
                  const real &dvy, const real &dvz, real *rVec, real *nVec);
 
+/**
+ * @brief Compute the derivatives of a continuous event with respect to position, velocity, and parameters.
+ */
+void stm_continuous_event(STMParameters &stmParams,
+                          const PropSimulation *propSim, const size_t &eventIdx,
+                          const real &tPastEvent, const real &postFac);
 #endif

src/force.cpp

@@ -41,7 +41,7 @@ static void force_thruster(const PropSimulation *propSim, std::vector<real> &acc
  * @brief Compute the acceleration of the system due to a continuous event.
  */
 static void force_continuous_event(const real &t, const PropSimulation *propSim,
-                                   std::vector<real> &accInteg);
+                                   std::vector<real> &accInteg, STMParameters* allSTMs);
 
 /**
  * @param[in] t Time [TDB MJD]
@@ -121,7 +121,7 @@ void get_state_der(PropSimulation *propSim, const real &t,
     force_J2(propSim, accInteg, allSTMs);
     force_nongrav(propSim, accInteg, allSTMs);
     force_thruster(propSim, accInteg);
-    force_continuous_event(t, propSim, accInteg);
+    force_continuous_event(t, propSim, accInteg, allSTMs);
     #ifdef PRINT_FORCES
     forceFile.open("cpp.11", std::ios::app);
     forceFile << std::setw(10) << "total_acc" << std::setw(25) << accInteg[0]
@@ -648,27 +648,31 @@ static void force_thruster(const PropSimulation *propSim,
  * @param[in] t Time [TDB MJD]
  * @param[in] propSim PropSimulation object.
  * @param[inout] accInteg State derivative vector.
+ * @param[in] allSTMs STMParameters vector for IntegBodies in the simulation.
  */
 static void force_continuous_event(const real &t, const PropSimulation *propSim,
-                                   std::vector<real> &accInteg) {
+                                   std::vector<real> &accInteg, STMParameters* allSTMs) {
     if (!propSim->eventMngr.allConEventDone){
         const bool forwardProp = propSim->integParams.tf > propSim->integParams.t0;
         const real propDir = forwardProp ? 1.0 : -1.0;
         for (size_t i = 0; i < propSim->eventMngr.continuousEvents.size(); i++){
-            if (propSim->eventMngr.continuousEvents[i].isHappening){
+            if (propSim->eventMngr.continuousEvents[i].hasStarted){
                 const size_t starti =
                     propSim->eventMngr.continuousEvents[i].xIntegIndex / 2;
                 const real tPastEvent = t - propSim->eventMngr.continuousEvents[i].t;
-                const real postFac = exp(-tPastEvent/propSim->eventMngr.continuousEvents[i].tau);
+                const real postFac = exp(-tPastEvent/propSim->eventMngr.continuousEvents[i].tau) * propDir;
                 accInteg[starti + 0] += propSim->eventMngr.continuousEvents[i].expAccel0[0] *
-                    postFac *
-                    propSim->eventMngr.continuousEvents[i].multiplier * propDir;
+                    postFac;
                 accInteg[starti + 1] += propSim->eventMngr.continuousEvents[i].expAccel0[1] *
-                    postFac *
-                    propSim->eventMngr.continuousEvents[i].multiplier * propDir;
+                    postFac;
                 accInteg[starti + 2] += propSim->eventMngr.continuousEvents[i].expAccel0[2] *
-                    postFac *
-                    propSim->eventMngr.continuousEvents[i].multiplier * propDir;
+                    postFac;
+                const size_t bodyIdx = propSim->eventMngr.continuousEvents[i].bodyIndex;
+                if (propSim->integBodies[bodyIdx].propStm &&
+                    propSim->eventMngr.continuousEvents[i].eventEst) {
+                    stm_continuous_event(allSTMs[bodyIdx], propSim, i,
+                                            tPastEvent, postFac);
+                }
             }
         }
     }

src/gr15.cpp

@@ -262,9 +262,9 @@ void check_continuous_events(PropSimulation *propSim, const real &t) {
                 eventEnded = t < propSim->eventMngr.continuousEvents[i].t;
             }
             if (eventStarted && !eventEnded) {
-                propSim->eventMngr.continuousEvents[i].isHappening = true;
+                propSim->eventMngr.continuousEvents[i].hasStarted = true;
             } else {
-                propSim->eventMngr.continuousEvents[i].isHappening = false;
+                propSim->eventMngr.continuousEvents[i].hasStarted = false;
             }
             if (!eventEnded) {
                 allDone = false;
@@ -311,7 +311,7 @@ void event_timestep_check(PropSimulation *propSim, real &dt) {
                 tNextEvent = fmax(tNextEvent, tNextConEvent);
             }
             // break up the timestep if early in a continuous event
-            if (propSim->eventMngr.continuousEvents[i].isHappening) {
+            if (propSim->eventMngr.continuousEvents[i].hasStarted) {
                 const real timeConstant = propSim->eventMngr.continuousEvents[i].tau;
                 const real timeConstantFac = 5.0L;
                 const real numSegments = 100.0L;