earth GDL; continuous events

grss/version.txt

@@ -1 +1 @@
-4.2.0
+4.2.1

include/gr15.h

@@ -79,11 +79,28 @@ void refine_b(std::vector<real> &b, std::vector<real> &e, const real &dtRatio,
               const size_t &dim);
 
 /**
- * @brief Check if any events need to be applied after the current timestep.
+ * @brief Check if any impulsive events need to be applied after the current timestep.
  */
-void check_and_apply_events(PropSimulation *propSim, const real &t,
-                            real &tNextEvent, size_t &nextEventIdx,
-                            std::vector<real> &xInteg);
+void check_and_apply_impulsive_events(PropSimulation *propSim, const real &t,
+                                      std::vector<real> &xInteg);
+
+/**
+ * @brief Check if any continuous events need to be applied after the current timestep.
+ */
+void check_continuous_events(PropSimulation *propSim, const real &t);
+
+/**
+ * @brief Top level function to check for events after the current timestep.
+ */
+void check_events(PropSimulation *propSim, const real &t, std::vector<real> &xInteg);
+
+/**
+ * @brief Check whether timestep is too large that it would skip an event.
+ * 
+ * @param[in] propSim PropSimulation object for the integration.
+ * @param[in] dt Current timestep.
+ */
+void event_timestep_check(PropSimulation *propSim, real &dt);
 
 /**
  * @brief 15th-order Gauss-Radau integrator for the PropSimulation.

include/grss.h

@@ -20,7 +20,8 @@ void PropSimulation::integrate() {
 
     // flip vectors if integration is backwards in time
     if (this->integParams.t0 > this->integParams.tf) {
-        std::reverse(this->events.begin(), this->events.end());
+        std::reverse(this->eventMngr.impulsiveEvents.begin(), this->eventMngr.impulsiveEvents.end());
+        std::reverse(this->eventMngr.continuousEvents.begin(), this->eventMngr.continuousEvents.end());
         std::reverse(this->xObserver.begin(), this->xObserver.end());
         std::reverse(this->observerInfo.begin(), this->observerInfo.end());
         std::reverse(this->tEval.begin(), this->tEval.end());

include/simulation.h

@@ -222,34 +222,50 @@ class IntegBody : public Body {
  * @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 threshold Threshold for the continuous impulse to reach in time tEvent+dt.
+ * @param isContinuous Flag to indicate if the event is continuous.
+ * @param isHappening Flag to indicate if the continuous event is happening.
  */
 class Event {
    private:
    public:
     real t;
     std::string bodyName;
     size_t bodyIndex;
-};
+    size_t xIntegIndex;
 
-/**
- * @brief Class for impulse events in a PropSimulation.
- * 
- * @param deltaV Delta-V for the impulse.
- * @param multiplier Multiplier for the Delta-V.
- */
-class ImpulseEvent : public Event {
-   private:
-   public:
+    // for impulsive events
     std::vector<real> deltaV = {0.0L, 0.0L, 0.0L};
     real multiplier = 1.0L;
+
+    // for continuous ejecta events
+    real dt = 1.0L;
+    real threshold = 0.999L;
+    bool isContinuous = false;
+    bool isHappening = false;
+    real c = 3.8L;
     /**
      * @brief Apply the impulse event to the body.
-     * 
-     * @param[in] t Time of the event.
-     * @param[inout] xInteg State of the body.
-     * @param[in] propDir Direction of propagation.
      */
-    void apply(const real &t, std::vector<real> &xInteg, const real &propDir);
+    void apply_impulsive(const real &t, std::vector<real> &xInteg, const real &propDir);
+};
+
+class EventManager {
+   private:
+   public:
+    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;
 };
 
 /** 
@@ -461,7 +477,7 @@ class PropSimulation {
     IntegrationParameters integParams;
     std::vector<SpiceBody> spiceBodies;
     std::vector<IntegBody> integBodies;
-    std::vector<ImpulseEvent> events;
+    EventManager eventMngr;
     std::vector<CloseApproachParameters> caParams;
     std::vector<ImpactParameters> impactParams;
     real t;
@@ -506,6 +522,11 @@ class PropSimulation {
      */
     void add_event(IntegBody body, real tEvent, std::vector<real> deltaV,
                    real multiplier = 1.0L);
+    /**
+     * @brief Add an ejecta event to the simulation.
+     */
+    void add_event(IntegBody body, real tEvent, std::vector<real> deltaV,
+                   real multiplier = 1.0L, real dt = 1.0L, real threshold = 0.999L);
     /**
      * @brief Set the values of the PropSimulation Constants object.
      */

src/force.cpp

@@ -37,6 +37,12 @@ static void force_nongrav(const PropSimulation *propSim, std::vector<real> &accI
  */
 static void force_thruster(const PropSimulation *propSim, std::vector<real> &accInteg);
 
+/**
+ * @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);
+
 /**
  * @param[in] t Time [TDB MJD]
  * @param[in] xInteg State vector
@@ -115,6 +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);
     #ifdef PRINT_FORCES
     forceFile.open("cpp.11", std::ios::app);
     forceFile << std::setw(10) << "total_acc" << std::setw(25) << accInteg[0]
@@ -636,3 +643,32 @@ static void force_thruster(const PropSimulation *propSim,
     forceFile.close();
     #endif
 }
+
+/**
+ * @param[in] propSim PropSimulation object.
+ * @param[inout] accInteg State derivative vector.
+ */
+static void force_continuous_event(const real &t, const PropSimulation *propSim,
+                                   std::vector<real> &accInteg) {
+    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){
+                const Event event = propSim->eventMngr.continuousEvents[i];
+                const size_t starti = event.xIntegIndex / 2;
+                const real tPastEvent = t - event.t;
+                // f = np.sin(np.pi*x/(2*time_for_val))
+                // df = np.pi/(2*time_for_val)*np.cos(np.pi*x/(2*time_for_val))
+                const real preFac = PI/(2*event.dt);
+                real accFac = preFac*cos(preFac*tPastEvent);
+                if (accFac > preFac || accFac < 0.0){
+                    accFac = 0.0;
+                }
+                accInteg[starti + 0] += accFac * event.deltaV[0] * propDir;
+                accInteg[starti + 1] += accFac * event.deltaV[1] * propDir;
+                accInteg[starti + 2] += accFac * event.deltaV[2] * propDir;
+            }
+        }
+    }
+}

src/gr15.cpp

@@ -225,30 +225,104 @@ void refine_b(std::vector<real> &b, std::vector<real> &e, const real &dtRatio, c
 /**
  * @param[in] propSim PropSimulation object for the integration.
  * @param[in] t Current time.
- * @param[out] tNextEvent Time of the next event.
- * @param[out] nextEventIdx Index of the next event.
  * @param[in] xInteg Vector of integrated states after the current timestep.
  */
-void check_and_apply_events(PropSimulation *propSim, const real &t,
-                            real &tNextEvent, size_t &nextEventIdx,
-                            std::vector<real> &xInteg) {
-    while (nextEventIdx < propSim->events.size() && t == tNextEvent) {
+void check_and_apply_impulsive_events(PropSimulation *propSim, const real &t,
+                                      std::vector<real> &xInteg) {
+    size_t *nextImpEventIdx = &propSim->eventMngr.nextImpEventIdx;
+    real *tNextImpEvent = &propSim->eventMngr.tNextImpEvent;
+    while (*nextImpEventIdx < propSim->eventMngr.nImpEvents && t == *tNextImpEvent) {
         // apply events for the state just reached by the integrator
-        real propDir;
-        if (propSim->integParams.t0 < propSim->integParams.tf) {
-            propDir = 1.0L;
+        bool forwardProp = propSim->integParams.tf > propSim->integParams.t0;
+        real propDir = forwardProp ? 1.0L : -1.0L;
+        propSim->eventMngr.impulsiveEvents[*nextImpEventIdx].apply_impulsive(t, xInteg, propDir);
+        // update next event index and time
+        (*nextImpEventIdx)++;
+        if (*nextImpEventIdx < propSim->eventMngr.nImpEvents) {
+            *tNextImpEvent = propSim->eventMngr.impulsiveEvents[*nextImpEventIdx].t;
         } else {
-            propDir = -1.0L;
+            *tNextImpEvent = propSim->integParams.tf;
         }
-        propSim->events[nextEventIdx].apply(t, xInteg, propDir);
-        // update next event index and time
-        nextEventIdx++;
-        if (nextEventIdx < propSim->events.size()) {
-            tNextEvent = propSim->events[nextEventIdx].t;
+    }
+}
+
+/**
+ * @param[in] propSim PropSimulation object for the integration.
+ * @param[in] t Current time.
+ */
+void check_continuous_events(PropSimulation *propSim, const real &t) {
+    if (!propSim->eventMngr.allConEventDone){
+        bool allDone = true;
+        bool forwardProp = propSim->integParams.tf > propSim->integParams.t0;
+        for (size_t i = 0; i < propSim->eventMngr.nConEvents; i++) {
+            // const bool eventStarted = t >= propSim->eventMngr.continuousEvents[i].t;
+            // const bool eventEnded = t > propSim->eventMngr.continuousEvents[i].t + propSim->eventMngr.continuousEvents[i].dt;
+            bool eventStarted, eventEnded;
+            if (forwardProp) {
+                eventStarted = t >= propSim->eventMngr.continuousEvents[i].t;
+                eventEnded = t >= (propSim->eventMngr.continuousEvents[i].t + propSim->eventMngr.continuousEvents[i].dt);
+            } else {
+                eventStarted = t <= (propSim->eventMngr.continuousEvents[i].t + propSim->eventMngr.continuousEvents[i].dt);
+                eventEnded = t < propSim->eventMngr.continuousEvents[i].t;
+            }
+            if (eventStarted && !eventEnded) {
+                propSim->eventMngr.continuousEvents[i].isHappening = true;
+            } else {
+                propSim->eventMngr.continuousEvents[i].isHappening = false;
+            }
+            if (!eventEnded) {
+                allDone = false;
+            }
+        }
+        if (allDone) {
+            propSim->eventMngr.allConEventDone = true;
+            propSim->eventMngr.tNextConEvent = propSim->integParams.tf;
+        }
+    }
+}
+
+/**
+ * @param[in] propSim PropSimulation object for the integration.
+ * @param[in] t Current time.
+ * @param[in] xInteg Vector of integrated states after the current timestep.
+ */
+void check_events(PropSimulation *propSim, const real &t, std::vector<real> &xInteg){
+    check_and_apply_impulsive_events(propSim, t, xInteg);
+    check_continuous_events(propSim, t);
+}
+
+/**
+ * @param[in] propSim PropSimulation object for the integration.
+ */
+void event_timestep_check(PropSimulation *propSim, real &dt) {
+    real tNextEvent = propSim->integParams.tf;
+    const bool forwardProp = propSim->integParams.tf > propSim->integParams.t0;
+    if (propSim->eventMngr.nImpEvents > 0) {
+        const real tNextImpEvent = propSim->eventMngr.tNextImpEvent;
+        if (forwardProp) {
+            tNextEvent = fmin(tNextEvent, tNextImpEvent);
         } else {
-            tNextEvent = propSim->integParams.tf;
+            tNextEvent = fmax(tNextEvent, tNextImpEvent);
+        }
+    }
+    if (!propSim->eventMngr.allConEventDone) {
+        for (size_t i = 0; i < propSim->eventMngr.nConEvents; i++) {
+            if (propSim->eventMngr.continuousEvents[i].isHappening) {
+                // if any continuous event is happening, set dt to 1/10th of the event duration
+                dt = propSim->eventMngr.continuousEvents[i].dt/10.0L;
+                dt = forwardProp ? dt : -dt;
+            }
+            if (forwardProp && propSim->t < propSim->eventMngr.continuousEvents[i].t) {
+                tNextEvent = fmin(tNextEvent, propSim->eventMngr.continuousEvents[i].t);
+            } else if (!forwardProp && propSim->t > propSim->eventMngr.continuousEvents[i].t+propSim->eventMngr.continuousEvents[i].dt) {
+                tNextEvent = fmax(tNextEvent, propSim->eventMngr.continuousEvents[i].t+propSim->eventMngr.continuousEvents[i].dt);
+            }
         }
     }
+    if ((forwardProp && propSim->t + dt > tNextEvent) ||
+        (!forwardProp && propSim->t + dt < tNextEvent)) {
+        dt = tNextEvent - propSim->t;
+    }
 }
 
 /**
@@ -342,21 +416,8 @@ void gr15(PropSimulation *propSim) {
     std::vector<real> g(7 * dim, 0.0);
     std::vector<real> e(7 * dim, 0.0);
     real dtReq;
-    real tNextEvent = propSim->integParams.tf;
-    size_t nextEventIdx = 0;
-    if (t == propSim->integParams.t0) {
-        nextEventIdx = 0;
-    }
-    if (propSim->events.size() != 0) {
-        tNextEvent = propSim->events[0].t;
-    }
-    check_and_apply_events(propSim, t, tNextEvent, nextEventIdx, xInteg);
-    if ((propSim->integParams.tf > propSim->integParams.t0 &&
-         t + dt > tNextEvent) ||
-        (propSim->integParams.tf < propSim->integParams.t0 &&
-         t + dt < tNextEvent)) {
-        dt = tNextEvent - t;
-    }
+    check_events(propSim, t, xInteg);
+    event_timestep_check(propSim, dt);
     propSim->interpParams.tStack.push_back(t);
     propSim->interpParams.xIntegStack.push_back(xInteg);
     std::vector<real> accInteg0(dim, 0.0);
@@ -430,11 +491,12 @@ void gr15(PropSimulation *propSim) {
             approx_xInteg(xInteg0, accInteg0, dt, 1.0, b, dim,
                         propSim->integBodies, xInteg, xIntegCompCoeffs);
             t += dt;
+            check_events(propSim, t, xInteg);
             get_state_der(propSim, t, xInteg, accInteg0);
-            check_and_apply_events(propSim, t, tNextEvent, nextEventIdx,
-                                    xInteg);
             propSim->interpParams.tStack.push_back(t);
             propSim->interpParams.xIntegStack.push_back(xInteg);
+            propSim->t = t;
+            propSim->xInteg = xInteg;
             propSim->integParams.timestepCounter++;
             refine_b(b, e, dtReq/dt, dim);
             check_ca_or_impact(propSim, t-dt, xInteg0, t, xInteg);
@@ -445,14 +507,7 @@ void gr15(PropSimulation *propSim) {
                 keepStepping = 0;
             }
             dt = dtReq;
-            if ((propSim->integParams.tf > propSim->integParams.t0 &&
-                 t + dt > tNextEvent) ||
-                (propSim->integParams.tf < propSim->integParams.t0 &&
-                 t + dt < tNextEvent)) {
-                dt = tNextEvent - t;
-            }
-            propSim->t = t;
-            propSim->xInteg = xInteg;
+            event_timestep_check(propSim, dt);
             oneStepDone = 1;
         }
     }