Merge pull request #79 from rahil-makadia/dev

Initialization internet check; Earth GDL; impulsive/continuous event partials

grss/fit/fit_ades.py

@@ -32,6 +32,7 @@
 }
 ades_grss_columns = {
     'obsTimeMJD': 'float', 'obsTimeMJDTDB': 'float', 'cosDec': 'float',
+    'resChi': 'float',
 }
 ades_column_types = ades_keep_columns | ades_add_columns | ades_grss_columns
 

grss/fit/fit_optical.py

@@ -479,7 +479,8 @@ def apply_debiasing_scheme(obs_df, lowres, verbose):
             unbias_counter += len(group)
             continue
         if cat_code not in biased_catalogs:
-            print(f"\tUnknown star catalog: {cat}")
+            if verbose:
+                print(f"\tUnknown star catalog: {cat}")
             obs_df.loc[group.index, 'biasRA'] = 0.0
             obs_df.loc[group.index, 'biasDec'] = 0.0
             continue

grss/fit/fit_radar.py

@@ -71,7 +71,6 @@ def add_radar_obs(obs_df, t_min_tdb=None, t_max_tdb=None, verbose=False):
         print(f"Read in {num_obs} radar observations from JPL radar API.")
     data = raw_data['data']
     time_range_count = 0
-    pre1972_count = 0
     for i in range(num_obs):
         obs = data[num_obs-i-1]
         date = Time(obs[1], format='iso', scale='utc')
@@ -105,14 +104,6 @@ def add_radar_obs(obs_df, t_min_tdb=None, t_max_tdb=None, verbose=False):
         obs_df.loc[idx,'frq'] = freq
         obs_df.loc[idx,'sigDelay'] = obs_sigma if delay else np.nan
         obs_df.loc[idx,'sigDoppler'] = obs_sigma if doppler else np.nan
-        if date.utc.mjd < 41317:
-            pre1972_count += 1
-            obs_df.loc[idx,'selAst'] = 'd'
-    if pre1972_count > 0:
-        print(f"\tWARNING: {pre1972_count} radar observations were taken before 1972.\n"
-                "\t\tThese residuals cannot be reliably computed because high-accuracy Earth "
-                "orientation kernels are not available for this time period.\n"
-                "\t\tForce deleting them...")
     if verbose:
         print(f"\tFiltered to {num_obs-time_range_count} observations that satisfy the "
                 "time range constraints.")

grss/kernels/get_kernels.py

@@ -37,26 +37,26 @@
 
 # get the earth orientation binary spice kernels and their comments if they are not already present
 # latest earth pck
-latest_earth_pck_exists = os.path.exists(f'{script_dir}/earth_latest_high_prec.bpc')
+latest_earth_pck_exists = os.path.exists(f'{script_dir}/earth_latest.bpc')
 latest_earth_pck_old = (latest_earth_pck_exists and
-                        time.time() - os.path.getmtime(f'{script_dir}/earth_latest_high_prec.bpc')
+                        time.time() - os.path.getmtime(f'{script_dir}/earth_latest.bpc')
                         > 86400)
 if not latest_earth_pck_exists or latest_earth_pck_old:
     print('Downloading latest Earth binary PCK...')
-    os.system((f'curl --silent --show-error -o {script_dir}/earth_latest_high_prec.bpc '
+    os.system((f'curl --silent --show-error -o {script_dir}/earth_latest.bpc '
                 f'{NAIF_SITE}/pck/earth_latest_high_prec.bpc'))
 # historical earth pck
-historical_earth_pck_exists = os.path.exists(f'{script_dir}/earth_720101_230601.bpc')
+historical_earth_pck_exists = os.path.exists(f'{script_dir}/earth_historic.bpc')
 if not historical_earth_pck_exists:
-    print('Downloading historical Earth binary PCK...')
-    os.system((f'curl --silent --show-error -o {script_dir}/earth_720101_230601.bpc '
-                f'{NAIF_SITE}/pck/earth_720101_230601.bpc'))
+    print('Downloading historic Earth binary PCK...')
+    os.system((f'curl --silent --show-error -o {script_dir}/earth_historic.bpc '
+                f'{NAIF_SITE}/pck/earth_620120_240827.bpc'))
 # predicted earth pck
-predicted_earth_pck_exists = os.path.exists(f'{script_dir}/earth_200101_990825_predict.bpc')
+predicted_earth_pck_exists = os.path.exists(f'{script_dir}/earth_predict.bpc')
 if not predicted_earth_pck_exists:
     print('Downloading predicted Earth binary PCK...')
-    os.system((f'curl --silent --show-error -o {script_dir}/earth_200101_990825_predict.bpc '
-                f'{NAIF_SITE}/pck/earth_200101_990825_predict.bpc'))
+    os.system((f'curl --silent --show-error -o {script_dir}/earth_predict.bpc '
+                f'{NAIF_SITE}/pck/earth_200101_990827_predict.bpc'))
 # generic frame kernels
 generic_frame_kernel_exists = os.path.exists(f'{script_dir}/pck00011.tpc')
 if not generic_frame_kernel_exists:

grss/utils.py

@@ -2,6 +2,7 @@
 import os
 import time
 import gzip
+import requests
 
 __all__ = [ 'default_kernel_path',
             'grss_path',
@@ -11,6 +12,14 @@
 grss_path = os.path.dirname(os.path.abspath(__file__))
 default_kernel_path = f'{grss_path}/kernels/'
 
+def _connected_to_internet(url='http://www.google.com/', timeout=5):
+    try:
+        _ = requests.head(url, timeout=timeout)
+        return True
+    except requests.ConnectionError:
+        print("No internet connection available.")
+    return False
+
 def _download_ades_files(file):
     """
     Download a text file of the observatory code info from the Minor Planet Center.
@@ -66,6 +75,10 @@ def initialize():
     None : NoneType
         None
     """
+    internet = _connected_to_internet()
+    if not internet:
+        print("Please connect to the internet to download/update the necessary data files.")
+        return None
     # run the get_debiasing_data.py script
     os.system(f'python {grss_path}/debias/get_debiasing_data.py')
     # run the get_kernels.py script

grss/version.txt

@@ -1 +1 @@
-4.2.0
+4.3.2

include/gr15.h

@@ -79,16 +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.
+ */
+void event_timestep_check(PropSimulation *propSim, real &dt);
 
 /**
  * @brief 15th-order Gauss-Radau integrator for the PropSimulation.
- * 
- * @param[inout] propSim PropSimulation object for the integration.
  */
 void gr15(PropSimulation *propSim);
 

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,57 @@ 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 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;
+    bool isContinuous = false;
+    bool eventEst = false;
     size_t bodyIndex;
-};
+    size_t xIntegIndex;
+    bool hasStarted = false;
 
-/**
- * @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:
-    std::vector<real> deltaV = {0.0L, 0.0L, 0.0L};
-    real multiplier = 1.0L;
+    // for impulsive events
+    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 = 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.
-     * 
-     * @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(PropSimulation *propSim, const real &t, std::vector<real> &xInteg);
+};
+
+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 +484,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;
@@ -502,10 +525,9 @@ class PropSimulation {
      */
     void remove_body(std::string name);
     /**
-     * @brief Add an impulse event to the simulation.
+     * @brief Add an event to the simulation.
      */
-    void add_event(IntegBody body, real tEvent, std::vector<real> deltaV,
-                   real multiplier = 1.0L);
+    void add_event(Event event);
     /**
      * @brief Set the values of the PropSimulation Constants object.
      */
@@ -524,7 +546,7 @@ class PropSimulation {
         bool convergedLightTime = false,
         std::vector<std::vector<real>> observerInfo =
             std::vector<std::vector<real>>(),
-        bool adaptiveTimestep = true, real dt0 = 0.0L,
+        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);
     /**

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

joss/joss_paper.bib

@@ -83,6 +83,34 @@ @article{Milani1999
 author = {Andrea Milani and Giovanni B. Valsecchi},
 }
 
+@Article{Farnocchia2019,
+  author    = {Farnocchia, Davide and Eggl, Siegfried and Chodas, Paul W. and Giorgini, Jon D. and Chesley, Steven R.},
+  journal   = {Celestial Mechanics and Dynamical Astronomy},
+  title     = {Planetary encounter analysis on the B-plane: a comprehensive formulation},
+  year      = {2019},
+  issn      = {1572-9478},
+  month     = aug,
+  number    = {8},
+  volume    = {131},
+  doi       = {10.1007/s10569-019-9914-4},
+  publisher = {Springer Science and Business Media LLC},
+}
+
+@article{Park2021,
+doi = {10.3847/1538-3881/abd414},
+url = {https://dx.doi.org/10.3847/1538-3881/abd414},
+year = {2021},
+month = {feb},
+publisher = {The American Astronomical Society},
+volume = {161},
+number = {3},
+pages = {105},
+author = {Ryan S. Park and William M. Folkner and James G. Williams and Dale H. Boggs},
+title = {The JPL Planetary and Lunar Ephemerides DE440 and DE441},
+journal = {The Astronomical Journal},
+abstract = {The planetary and lunar ephemerides called DE440 and DE441 have been generated by fitting numerically integrated orbits to ground-based and space-based observations. Compared to the previous general-purpose ephemerides DE430, seven years of new data have been added to compute DE440 and DE441, with improved dynamical models and data calibration. The orbit of Jupiter has improved substantially by fitting to the Juno radio range and Very Long Baseline Array (VLBA) data of the Juno spacecraft. The orbit of Saturn has been improved by radio range and VLBA data of the Cassini spacecraft, with improved estimation of the spacecraft orbit. The orbit of Pluto has been improved from use of stellar occultation data reduced against the Gaia star catalog. The ephemerides DE440 and DE441 are fit to the same data set, but DE441 assumes no damping between the lunar liquid core and the solid mantle, which avoids a divergence when integrated backward in time. Therefore, DE441 is less accurate than DE440 for the current century, but covers a much longer duration of years −13,200 to +17,191, compared to DE440 covering years 1550–2650.}
+}
+
 @article{Holman2023,
 doi = {10.3847/PSJ/acc9a9},
 year = {2023},