Merge pull request #55 from rahil-makadia/dev

better spk memory handling; switch wget to curl
rahil-makadia · Mar 27, 2024 · 7772035 · 7772035
2 parents c7ecc78 + 96c1dbc
commit 7772035
Show file tree

Hide file tree

Showing 18 changed files with 320 additions and 240 deletions.
diff --git a/.gitignore b/.gitignore
@@ -183,6 +183,10 @@ cython_debug/
 *.app
 
 ################################### Individual user exceptions
+# temporary files
+temp*
+temp
+
 # Sphinx documentation
 docs/_build/
 docs/source/_autosummary/

diff --git a/README.md b/README.md
@@ -80,13 +80,7 @@ Once the GRSS library has been installed, it can be imported into a Python scrip
    import grss
 ```
 
-The first time the library is imported, it will download some data files such as NAIF SPICE kernels and the data needed to debias optical astrometry. This should should take a few minutes, and if the download was completed, the following message will be printed for each file:
-
-```console
-   YYYY-MM-DD HH:MM:SS URL:url-of-downloaded-file [filesize] -> path/to/downloaded/file [1]
-```
-
-Once these files are available to the library, you are ready to use GRSS to its full potential!
+The first time the library is imported, it will download some data files such as NAIF SPICE kernels and the data needed to debias optical astrometry. This should should take a few minutes. Once these files are available to the library, you are ready to use GRSS to its full potential!
 
 Check out the [examples](https://rahil-makadia.github.io/grss/examples.html) on the [GRSS website](https://rahil-makadia.github.io/grss/) to get started.
 

diff --git a/docs/source/start.rst b/docs/source/start.rst
@@ -68,10 +68,4 @@ Once the GRSS library has been installed, it can be imported into a Python scrip
 
    import grss
 
-The first time the library is imported, it will download some data files such as NAIF SPICE kernels and the data needed to debias optical astrometry. This should should take a few minutes, and if the download was completed, the following message will be printed for each file:
-
-.. code-block:: console
-
-   YYYY-MM-DD HH:MM:SS URL:url-of-downloaded-file [filesize] -> path/to/downloaded/file [1]
-
-Once these files are available to the library, you are ready to use GRSS to its full potential!
+The first time the library is imported, it will download some data files such as NAIF SPICE kernels and the data needed to debias optical astrometry. This should should take a few minutes. Once these files are available to the library, you are ready to use GRSS to its full potential!
diff --git a/extern/get_cspice.py b/extern/get_cspice.py
@@ -154,7 +154,8 @@ def _download(self):
             HTTPS GET call to the NAIF server to download the required CSPICE
             distribution package.
         """
-        cmd = f"wget --no-clobber {self._rcspice}"
+        # cmd = f"wget --no-clobber {self._rcspice}"
+        cmd = f"curl -O {self._rcspice} --silent --show-error"
         subprocess.run(cmd,shell=True,check=True)
         return
 

diff --git a/grss/debias/get_debiasing_data.py b/grss/debias/get_debiasing_data.py
@@ -8,10 +8,14 @@
 # if lowres_data.tgz or lowres_data/ already exist, do not download
 if not os.path.exists(f'{cwd}/lowres_data.tgz') and not os.path.exists(f'{cwd}/lowres_data/'):
     FTP_FILE = 'debias_2018.tgz'
-    os.system(f'wget --no-verbose --no-clobber {FTP_SITE}/{FTP_FILE} -O {cwd}/lowres_data.tgz')
+    cmd = f'curl --silent --show-error -o {cwd}/lowres_data.tgz {FTP_SITE}/{FTP_FILE}'
+    print('Downloading low-resolution debiasing data')
+    os.system(cmd)
 if not os.path.exists(f'{cwd}/hires_data.tgz') and not os.path.exists(f'{cwd}/hires_data/'):
     FTP_FILE = 'debias_hires2018.tgz'
-    os.system(f'wget --no-verbose --no-clobber {FTP_SITE}/{FTP_FILE} -O {cwd}/hires_data.tgz')
+    cmd = f'curl --silent --show-error -o {cwd}/hires_data.tgz {FTP_SITE}/{FTP_FILE}'
+    print('Downloading high-resolution debiasing data')
+    os.system(cmd)
 
 # extract the data files
 # create the directories if they don't exist, delete the files if they do

diff --git a/grss/kernels/get_kernels.py b/grss/kernels/get_kernels.py
@@ -6,55 +6,83 @@
 # get the path to the directory containing this script
 script_dir = os.path.dirname(os.path.realpath(__file__))
 
-GRSS_SITE = 'https://github.com/rahil-makadia/grss/raw/dev/grss/kernels'
+GRSS_SITE = 'https://raw.githubusercontent.com/rahil-makadia/grss/main/grss/kernels'
 NAIF_SITE = 'https://naif.jpl.nasa.gov/pub/naif/generic_kernels'
 SSD_SITE = 'https://ssd.jpl.nasa.gov/ftp'
 # get the custom spice kernels if they are not already present
 # de431 planets + big16 1950-2350
-os.system((f'wget --no-verbose --no-clobber {GRSS_SITE}/planets_big16_de431_1950_2350.bsp '
-            f'-O {script_dir}/planets_big16_de431_1950_2350.bsp'))
-os.system((f'wget --no-verbose --no-clobber {GRSS_SITE}/planets_big16_de431_1950_2350.tm '
-            f'-O {script_dir}/planets_big16_de431_1950_2350.tm'))
+if not os.path.exists(f'{script_dir}/planets_big16_de431_1950_2350.bsp'):
+    print('Downloading combined DE430/431 planets and big16 asteroid SPK kernels...')
+    os.system((f'curl -H "Accept: application/vnd.github+json" '
+                f'-L --silent --show-error -o {script_dir}/planets_big16_de431_1950_2350.bsp '
+                f'{GRSS_SITE}/planets_big16_de431_1950_2350.bsp'))
+if not os.path.exists(f'{script_dir}/planets_big16_de431_1950_2350.tm'):
+    print('Downloading combined DE430/431 planets and big16 asteroid meta-kernel...')
+    os.system((f'curl -H "Accept: application/vnd.github+json" '
+                f'-L --silent --show-error -o {script_dir}/planets_big16_de431_1950_2350.tm '
+                f'{GRSS_SITE}/planets_big16_de431_1950_2350.tm'))
 # de441 planets + big16 1950-2350
-os.system((f'wget --no-verbose --no-clobber {GRSS_SITE}/planets_big16_de441_1950_2350.bsp '
-            f'-O {script_dir}/planets_big16_de441_1950_2350.bsp'))
-os.system((f'wget --no-verbose --no-clobber {GRSS_SITE}/planets_big16_de441_1950_2350.tm '
-            f'-O {script_dir}/planets_big16_de441_1950_2350.tm'))
+if not os.path.exists(f'{script_dir}/planets_big16_de441_1950_2350.bsp'):
+    print('Downloading combined DE440/441 planets and big16 asteroid kernels...')
+    os.system((f'curl -H "Accept: application/vnd.github+json" '
+                f'-L --silent --show-error -o {script_dir}/planets_big16_de441_1950_2350.bsp '
+                f'{GRSS_SITE}/planets_big16_de441_1950_2350.bsp'))
+if not os.path.exists(f'{script_dir}/planets_big16_de441_1950_2350.tm'):
+    print('Downloading combined DE440/441 planets and big16 asteroid meta-kernel...')
+    os.system((f'curl -H "Accept: application/vnd.github+json" '
+                f'-L --silent --show-error -o {script_dir}/planets_big16_de441_1950_2350.tm '
+                f'{GRSS_SITE}/planets_big16_de441_1950_2350.tm'))
 
 # get the generic spice kernels if they are not already present
 # de430 planets + de431 big16
-os.system((f'wget --no-verbose --no-clobber {NAIF_SITE}/spk/planets/de430.bsp '
-            f'-O {script_dir}/de430.bsp'))
-os.system((f'wget --no-verbose --no-clobber {SSD_SITE}/xfr/sb431-n16s.bsp '
-            f'-O {script_dir}/sb431-n16s.bsp'))
+if (not os.path.exists(f'{script_dir}/de430.bsp') or
+    not os.path.exists(f'{script_dir}/sb431-n16s.bsp')):
+    print('Downloading generic DE430/431 planets and big16 asteroid kernels...')
+    os.system((f'curl --silent --show-error -o {script_dir}/de430.bsp '
+                f'{NAIF_SITE}/spk/planets/de430.bsp'))
+    os.system((f'curl --silent --show-error -o {script_dir}/sb431-n16s.bsp '
+                f'{SSD_SITE}/xfr/sb431-n16s.bsp'))
 # de440 planets + de441 big16
-os.system((f'wget --no-verbose --no-clobber {NAIF_SITE}/spk/planets/de440.bsp '
-            f'-O {script_dir}/de440.bsp'))
-os.system((f'wget --no-verbose --no-clobber {SSD_SITE}/xfr/sb441-n16s.bsp '
-            f'-O {script_dir}/sb441-n16s.bsp'))
+if (not os.path.exists(f'{script_dir}/de440.bsp') or
+    not os.path.exists(f'{script_dir}/sb441-n16s.bsp')):
+    print('Downloading generic DE440/441 planets and big16 asteroid kernels...')
+    os.system((f'curl --silent --show-error -o {script_dir}/de440.bsp '
+                f'{NAIF_SITE}/spk/planets/de440.bsp'))
+    os.system((f'curl --silent --show-error -o {script_dir}/sb441-n16s.bsp '
+                f'{SSD_SITE}/xfr/sb441-n16s.bsp'))
 
 # get the latest spice leap second kernel
-os.system((f'wget --no-verbose --no-clobber {NAIF_SITE}/lsk/latest_leapseconds.tls '
-            f'-O {script_dir}/latest_leapseconds.tls'))
+if not os.path.exists(f'{script_dir}/latest_leapseconds.tls'):
+    print('Downloading latest leap second kernel...')
+    os.system((f'curl --silent --show-error -o {script_dir}/latest_leapseconds.tls '
+                f'{NAIF_SITE}/lsk/latest_leapseconds.tls'))
 # get the earth orientation binary spice kernels and their comments if they are not already present
 # latest earth pck
-os.system((f'wget --no-verbose --no-clobber {NAIF_SITE}/pck/earth_latest_high_prec.cmt '
-            f'-O {script_dir}/earth_latest_high_prec.cmt'))
-os.system((f'wget --no-verbose --no-clobber {NAIF_SITE}/pck/earth_latest_high_prec.bpc '
-            f'-O {script_dir}/earth_latest_high_prec.bpc'))
+if not os.path.exists(f'{script_dir}/earth_latest_high_prec.bpc'):
+    print('Downloading latest Earth binary PCK...')
+    os.system((f'curl --silent --show-error -o {script_dir}/earth_latest_high_prec.cmt '
+                f'{NAIF_SITE}/pck/earth_latest_high_prec.cmt'))
+    os.system((f'curl --silent --show-error -o {script_dir}/earth_latest_high_prec.bpc '
+                f'{NAIF_SITE}/pck/earth_latest_high_prec.bpc'))
 # historical earth pck
-os.system((f'wget --no-verbose --no-clobber {NAIF_SITE}/pck/earth_720101_230601.bpc '
-            f'-O {script_dir}/earth_720101_230601.bpc'))
-os.system((f'wget --no-verbose --no-clobber {NAIF_SITE}/pck/earth_720101_230601.cmt '
-            f'-O {script_dir}/earth_720101_230601.cmt'))
+if not os.path.exists(f'{script_dir}/earth_720101_230601.bpc'):
+    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'))
+    os.system((f'curl --silent --show-error -o {script_dir}/earth_720101_230601.cmt '
+                f'{NAIF_SITE}/pck/earth_720101_230601.cmt'))
 # predicted earth pck
-os.system((f'wget --no-verbose --no-clobber {NAIF_SITE}/pck/earth_200101_990825_predict.bpc '
-            f'-O {script_dir}/earth_200101_990825_predict.bpc'))
-os.system((f'wget --no-verbose --no-clobber {NAIF_SITE}/pck/earth_200101_990825_predict.cmt '
-            f'-O {script_dir}/earth_200101_990825_predict.cmt'))
+if not os.path.exists(f'{script_dir}/earth_200101_990825_predict.bpc'):
+    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_200101_990825_predict.cmt '
+                f'{NAIF_SITE}/pck/earth_200101_990825_predict.cmt'))
 # generic frame kernels
-os.system((f'wget --no-verbose --no-clobber {NAIF_SITE}/pck/pck00011.tpc '
-            f'-O {script_dir}/pck00011.tpc'))
+if not os.path.exists(f'{script_dir}/pck00011.tpc'):
+    print('Downloading generic frame kernels...')
+    os.system((f'curl --silent --show-error -o {script_dir}/pck00011.tpc '
+                f'{NAIF_SITE}/pck/pck00011.tpc'))
 
 # run this code if the no-tm-overwrite flag argument is not present
 if len(sys.argv) > 1:

diff --git a/grss/prop/prop_utils.py b/grss/prop/prop_utils.py
@@ -181,6 +181,7 @@ def plot_ca_summary(prop_sim, flyby_body, central_body='Earth',
         if body.name == flyby_body:
             break
         j += 2*body.n2Derivs
+    prop_sim.map_ephemeris()
     for idx, time in enumerate(t_stack):
         central_body_state = prop_sim.get_spiceBody_state(time, central_body)
         ca_body_state = x_stack[idx]
@@ -190,6 +191,7 @@ def plot_ca_summary(prop_sim, flyby_body, central_body='Earth',
         rel_state[idx,3] = ca_body_state[j+3] - central_body_state[3]
         rel_state[idx,4] = ca_body_state[j+4] - central_body_state[4]
         rel_state[idx,5] = ca_body_state[j+5] - central_body_state[5]
+    prop_sim.unmap_ephemeris()
     rel_dist = np.linalg.norm(rel_state[:,:3], axis=1)
     rel_radial_vel = np.sum(rel_state[:,:3]*rel_state[:,3:6], axis=1) / rel_dist
     scale_factor, units = get_scale_factor(central_body)

diff --git a/grss/version.txt b/grss/version.txt
@@ -1 +1 @@
-3.8.1
+3.8.2
diff --git a/include/force.h b/include/force.h
@@ -31,8 +31,8 @@ void force_ppn_eih(const PropSimulation *propSim, std::vector<real> &accInteg,
 /**
  * @brief Compute the acceleration of the system due to the J2 zonal harmonic.
  */
-void force_J2(const PropSimulation *propSim, std::vector<real> &accInteg,
-              std::vector<STMParameters> &allSTMs);
+void force_J2(const real &t, const PropSimulation *propSim,
+              std::vector<real> &accInteg, std::vector<STMParameters> &allSTMs);
 
 /**
  * @brief Compute the acceleration of the system due to the nongravitational forces.

diff --git a/include/grss.h b/include/grss.h
@@ -17,6 +17,7 @@ void PropSimulation::integrate() {
     if (!this->parallelMode) furnsh_c(this->DEkernelPath.c_str());
     gr15(this);
     if (!this->parallelMode) unload_c(this->DEkernelPath.c_str());
+    this->unmap_ephemeris();
 
     // flip vectors if integration is backwards in time
     if (this->integParams.t0 > this->integParams.tf) {

diff --git a/include/simulation.h b/include/simulation.h
@@ -6,6 +6,12 @@
 // forward declaration for main PropSimulation class
 class PropSimulation;
 
+/**
+ * @brief Get the name of body-fixed frame for a given SPICE ID.
+ */
+void get_baseBodyFrame(const int &spiceId, const real &tMjdTDB,
+                       ConstSpiceChar *&baseBodyFrame);
+
 /**
  * @brief Get the observer state for a given time.
  */
@@ -406,6 +412,7 @@ class PropSimulation {
    public:
     std::string name;
     std::string DEkernelPath;
+    Ephemeris ephem;
     /**
      * @brief Construct a new PropSimulation object.
      */
@@ -415,7 +422,14 @@ class PropSimulation {
      * @brief Construct a new PropSimulation object from a reference simulation.
      */
     PropSimulation(std::string name, const PropSimulation &simRef);
-    Ephemeris ephem;
+    /**
+     * @brief Memory map the ephemeris files.
+     */
+    void map_ephemeris();
+    /**
+     * @brief Unmap the ephemeris files.
+     */
+    void unmap_ephemeris();
     /**
      * @brief Get the state of a SpiceBody in the PropSimulation at a given time.
      */

diff --git a/include/spk.h b/include/spk.h
@@ -56,7 +56,7 @@ struct CacheItem {
  */
 struct EphemerisCache {
     double t;
-    struct CacheItem items[SPK_CACHE_ITEM_SIZE];
+    CacheItem items[SPK_CACHE_ITEM_SIZE];
 };
 
 /**
@@ -70,21 +70,6 @@ struct EphemerisCache {
  */
 #define RECORD_LEN 1024
 
-/**
- * @brief Structure to hold all the data for the SPK files in a PropSimulation.
- *
- * @param mb Main body ephemeris data.
- * @param sb Small body ephemeris data.
- * @param nextIdxToWrite Next index to write to in the cache.
- * @param cache Cache of recently queried SPK data.
- */
-struct Ephemeris {
-    struct SpkInfo *mb;
-    struct SpkInfo *sb;
-    size_t nextIdxToWrite = -1;
-    EphemerisCache cache[SPK_CACHE_SIZE];
-};
-
 /**
  * @brief Structure to hold the data for a single body in an SPK file.
  *
@@ -117,24 +102,51 @@ struct SpkTarget {
  * @param map Memory map of the SPK file.
  * @param len Length of the memory map.
  */
-struct SpkInfo {
-    struct SpkTarget *targets;
+struct DafInfo {
+    SpkTarget* targets;
     int num;
     int allocatedNum;
     void *map;
     size_t len;
 };
 
+/**
+ * @brief Structure to hold all the data for the SPK files in a PropSimulation.
+ *
+ * @param mb Main body ephemeris data.
+ * @param sb Small body ephemeris data.
+ * @param nextIdxToWrite Next index to write to in the cache.
+ * @param cache Cache of recently queried SPK data.
+ */
+struct Ephemeris {
+    std::string mbPath;
+    std::string sbPath;
+    DafInfo* mb = nullptr;
+    DafInfo* sb = nullptr;
+    size_t nextIdxToWrite = -1;
+    EphemerisCache cache[SPK_CACHE_SIZE];
+};
+
+/**
+ * @brief Free the memory allocated for an DafInfo structure.
+ */
+void daf_free(DafInfo *pl);
+
+/**
+ * @brief Initialise a DAF file.
+ */
+DafInfo* daf_init(const std::string &path, const std::string &type);
+
 /**
  * @brief Initialise an SPK file.
  */
-SpkInfo* spk_init(const std::string &path);
+DafInfo* spk_init(const std::string &path);
 
 /**
  * @brief Compute pos, vel, and acc for a given body
- * at a given time using an SpkInfo structure.
+ * at a given time using an DafInfo structure.
  */
-void spk_calc(SpkInfo *pl, double epoch, int spiceId, double *out_x,
+void spk_calc(DafInfo *pl, double epoch, int spiceId, double *out_x,
              double *out_y, double *out_z, double *out_vx, double *out_vy,
              double *out_vz, double *out_ax, double *out_ay, double *out_az);
 

diff --git a/src/approach.cpp b/src/approach.cpp
@@ -633,49 +633,7 @@ void CloseApproachParameters::print_summary(int prec){
  */
 void ImpactParameters::get_impact_parameters(PropSimulation *propSim){
     ConstSpiceChar *baseBodyFrame;
-    switch (this->centralBodySpiceId) {
-        case 10:
-            baseBodyFrame = "IAU_SUN";
-            break;
-        case 1:
-        case 199:
-            baseBodyFrame = "IAU_MERCURY";
-            break;
-        case 2:
-        case 299:
-            baseBodyFrame = "IAU_VENUS";
-            break;
-        case 399:
-            baseBodyFrame = "ITRF93";
-            // High precision frame is not defined before 1972 JAN 01 00:00:42.183 TDB
-            if (this->t < 41317.0004882291666666666L) {
-                baseBodyFrame = "IAU_EARTH";
-            }
-            break;
-        case 499:
-            baseBodyFrame = "IAU_MARS";
-            break;
-        case 599:
-            baseBodyFrame = "IAU_JUPITER";
-            break;
-        case 699:
-            baseBodyFrame = "IAU_SATURN";
-            break;
-        case 799:
-            baseBodyFrame = "IAU_URANUS";
-            break;
-        case 899:
-            baseBodyFrame = "IAU_NEPTUNE";
-            break;
-        case 999:
-            baseBodyFrame = "IAU_PLUTO";
-            break;
-        default:
-            std::cout << "get_impact_parameters: Given impacted body: "
-                      << this->centralBody << std::endl;
-            throw std::invalid_argument("Given base body not supported");
-            break;
-    }
+    get_baseBodyFrame(this->centralBodySpiceId, this->t, baseBodyFrame);
     real tET;
     mjd_to_et(this->t, tET);
     SpiceDouble rotMatSpice[6][6];