Merge pull request #83 from rahil-makadia/dev

update pip install pipeline; impact location comparison routines

Author: rahil-makadia

.github/workflows/joss.yml

@@ -18,7 +18,7 @@ jobs:
         uses: openjournals/openjournals-draft-action@master
         with:
           journal: joss
-          paper-path: ./joss/joss_paper.md
+          paper-path: ./joss/paper.md
       - name: Upload
         uses: actions/upload-artifact@master
         with:

.github/workflows/python_tests.yml

@@ -23,6 +23,9 @@ jobs:
           python-version: 3.11
       - name: Install GRSS (including dependencies)
         run: |
+          cd ./extern/
+          python3 get_cspice.py
+          cd ..
           python3 -m pip install --upgrade pip
           pip3 install "pybind11[global]>=2.10.0"
           pip3 install .

build_cpp.sh

@@ -1,4 +1,6 @@
 #!/bin/bash
+# make error stop the script
+set -e
 
 clean=0
 docs=0

grss/prop/prop_utils.py

@@ -12,6 +12,7 @@
             'plot_ca_summary',
             'plot_bplane',
             'plot_earth_impact',
+            'compare_earth_impact',
 ]
 
 earth_obliq = 84381.448/3600.0*np.pi/180.0
@@ -904,3 +905,95 @@ def plot_earth_impact(impact_list, print_ellipse_params=False, sigma_points=None
                     y=0.93)
     plt.show()
     return None
+
+def compare_earth_impact(coord1, coord2, sig1, sig2, labels, zoom_size=5.0e3, save_name=None):
+    """
+    Show two impact circumstances and compare them.
+
+    Parameters
+    ----------
+    coord1 : np.ndarray
+        [time, lon, lat] of the first impact
+    coord2 : np.ndarray
+        [time, lon, lat] of the second impact
+    sig1 : np.ndarray
+        [time, lon, lat] uncertainty of the first impact
+    sig2 : np.ndarray
+        [time, lon, lat] uncertainty of the second impact
+    labels : list
+        List of labels for the two impacts
+    zoom_size : float, optional
+        Size of the zoomed in plot, in km, by default 5.0e3.
+    save_name : str, optional
+        Name to save the plots as, by default None.
+
+    Returns
+    -------
+    None : NoneType
+        None
+    """
+    _, lon1, lat1 = coord1
+    _, lon2, lat2 = coord2
+    fig, ax = plt.subplots(1, 3, figsize=(10, 4), dpi=150, gridspec_kw={'width_ratios': [4, 3, 1]})
+    size = zoom_size*1e3
+    m = Basemap(width=size,height=size,
+                rsphere=(6378137.00,6356752.3142),
+                resolution='i',area_thresh=size/100,projection='lcc',
+                lat_0=lat1,lon_0=lon1,ax=ax[0])
+    m.bluemarble()
+    m.drawcountries()
+    lat_step = lon_step = 3 if zoom_size <= 1e3 else 10
+    if zoom_size <= 500:
+        lat_step = lon_step = 1
+    if zoom_size <= 100:
+        lat_step = lon_step = 0.5
+    if zoom_size > 5e3:
+        lon_step = 20
+        lat_step = 20
+    # draw parallels and meridians.
+    parallels = np.arange(-90, 91, lat_step)
+    # labels = [left,right,top,bottom]
+    m.drawparallels(parallels,labels=[True,False,False,False], color='gray')
+    meridians = np.arange(0, 361, lon_step)
+    m.drawmeridians(meridians,labels=[False,False,False,True], color='gray')
+    x_lon1,y_lat1 = m(lon1,lat1)
+    m.plot(x_lon1,y_lat1,'b.')
+    x_lon2,y_lat2 = m(lon2,lat2)
+    m.plot(x_lon2,y_lat2,'r.')
+    # draw line between two points
+    m.plot([x_lon1,x_lon2], [y_lat1,y_lat2], 'g--')
+
+    # plot error bars for grss solution, and plot GRSS solution
+    n_sig = 1
+    deg2m = np.pi/180*6378137.0
+    sig1[1:] *= deg2m
+    sig2[1:] *= deg2m
+    diff = coord2 - coord1
+    # diff[0] *= 86400*1e3
+    diff[1:] *= deg2m
+    diff[1] *= np.cos(lat1*np.pi/180)
+
+    # add two subplots
+    ax[1].set_xlabel('Longitude Difference [m]')
+    ax[1].set_ylabel('Latitude Difference [m]')
+    ax[1].errorbar(0.0, 0.0, xerr=sig1[1]*n_sig, yerr=sig1[2]*n_sig,
+                    fmt='s', capsize=6, markeredgewidth=1.5, lw=1.5, label=labels[0])
+    ax[1].errorbar(diff[1], diff[2], xerr=sig2[1]*n_sig, yerr=sig2[2]*n_sig,
+                    fmt='o', capsize=6, markeredgewidth=1.5, lw=1.5, label=labels[1])
+    ax[1].legend(ncol=1)
+
+    ax[2].set_ylabel('Time Difference [ms]')
+    ax[2].errorbar(0.0, 0.0, yerr=sig1[0]*n_sig,
+                    fmt='s', capsize=6, markeredgewidth=1.5, lw=1.5, label=labels[0])
+    ax[2].errorbar(0.0, diff[0], yerr=sig2[0]*n_sig,
+                    fmt='o', capsize=6, markeredgewidth=1.5, lw=1.5, label=labels[1])
+    ax[2].set_xticks([])
+    ax[2].yaxis.tick_right()
+    ax[2].yaxis.set_label_position('right')
+
+    fig.tight_layout()
+    if save_name is not None:
+        plt.savefig(save_name, bbox_inches='tight')
+    plt.suptitle(f'Impact Circumstance Comparison ({n_sig}$\\sigma$)', y=1.05)
+    plt.show()
+    return None

grss/version.txt

@@ -1 +1 @@
-4.3.4
+4.3.5

include/simulation.h

@@ -388,8 +388,6 @@ class ImpactParameters : public CloseApproachParameters {
     real lon;
     real lat;
     real alt;
-    std::vector<real> dlon = std::vector<real>(6, std::numeric_limits<real>::quiet_NaN());
-    std::vector<real> dlat = std::vector<real>(6, std::numeric_limits<real>::quiet_NaN());
     /**
      * @brief Get the impact parameters.
      */

joss/joss_paper.bib joss/paper.bibjoss/joss_paper.bib renamed to joss/paper.bib

@@ -26,6 +26,15 @@ @article{Makadia2024
 publisher = {The American Astronomical Society},
 }
 
+@article{Makadia2024b,
+author    = {Rahil Makadia and Davide Farnocchia and Steven R. Chesley and Siegfried Eggl},
+journal   = {The Planetary Science Journal},
+title     = {{Gauss-Radau Small-body Simulator (GRSS): An Open-Source Library for Planetary Defense}},
+year      = {2024},
+volume    = {Submitted},
+publisher = {The American Astronomical Society},
+}
+
 @article{Rein2014,
 author = {Rein, Hanno and Spiegel, David S.},
 title = "{ias15: a fast, adaptive, high-order integrator for gravitational dynamics, accurate to machine precision over a billion orbits}",

joss/joss_paper.md joss/paper.mdjoss/joss_paper.md renamed to joss/paper.md

@@ -12,12 +12,12 @@ authors:
     corresponding: true
     orcid: 0000-0001-9265-2230
     affiliation: 1
-  - name: Steven R. Chesley
-    orcid: 0000-0003-3240-6497
-    affiliation: 2
   - name: Davide Farnocchia
     orcid: 0000-0003-0774-884X
     affiliation: 2
+  - name: Steven R. Chesley
+    orcid: 0000-0003-3240-6497
+    affiliation: 2
   - name: Siegfried Eggl
     orcid: 0000-0002-1398-6302
     affiliation: 1
@@ -27,7 +27,7 @@ affiliations:
   - name: Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA
     index: 2
 date: XX August 2024
-bibliography: joss_paper.bib
+bibliography: paper.bib
 aas-doi: 10.3847/PSJ/xxxxx
 aas-journal: Planetary Science Journal
 ---
@@ -40,7 +40,7 @@ Modeling the motion of small solar system bodies is of utmost importance when lo
 
 In this paper, we present ``GRSS``, the Gauss-Radau Small-body Simulator, an open-source library for orbit determination and propagation of small bodies in the solar system. ``GRSS`` is an open-source software library with a C++11 foundation and a Python binding. The propagator is based on the ``IAS15`` algorithm, a 15<sup>th</sup> order integrator based on Gauss-Radau quadrature [@Rein2014]. Only the particles of interest are integrated within ``GRSS`` to reduce computational cost. The states for the planets and 16 largest main-belt asteroids are computed using memory-mapped JPL digital ephemeris kernels [@Park2021] as done in the ``ASSIST`` orbit propagator library [@Holman2023]. In addition to the propagator, the C++ portion of the library also has the ability to predict impacts and calculate close encounter circumstances using various formulations of the B-plane [@Kizner1961; @Opik1976; @Chodas1999; @Milani1999; @Farnocchia2019].
 
-The C++ functionality is exposed to Python through a binding generated using ``pybind11``[^1]. The Python layer of ``GRSS`` uses the propagator as the foundation to compute the orbits of small bodies from a given set of optical and/or radar astrometry from the Minor Planet Center[^2], the JPL Small Body Radar Astrometry database[^3], and the Gaia Focused Product Release solar system observations database[^4]. Additionally, the orbit determination modules also have the ability to fit especially demanding measurements such as stellar occultations. These capabilities of the ``GRSS`` library have already been used to study the the heliocentric changes in the orbit of the (65803) Didymos binary asteroid system as a result of the DART impact [@Makadia2022; @Makadia2024].
+The C++ functionality is exposed to Python through a binding generated using ``pybind11``[^1]. The Python layer of ``GRSS`` uses the propagator as the foundation to compute the orbits of small bodies from a given set of optical and/or radar astrometry from the Minor Planet Center[^2], the JPL Small Body Radar Astrometry database[^3], and the Gaia Focused Product Release solar system observations database[^4]. Additionally, the orbit determination modules also have the ability to fit especially demanding measurements such as stellar occultations. These capabilities of the ``GRSS`` library have already been used to study the the heliocentric changes in the orbit of the (65803) Didymos binary asteroid system as a result of the DART impact [@Makadia2022; @Makadia2024] and for analyzing the impact locations of two asteroids, 2024 BX<sub>1</sub> and 2024 RW<sub>1</sub> [@Makadia2024b].
 
 ``GRSS`` will continue to be developed in the future, with anticipated contributions including the ability to perform mission studies for future asteroid deflections. ``GRSS`` is publicly available to the community through the Python Package Index and the source code is available on GitHub[^5]. Therefore, ``GRSS`` is a reliable and efficient tool that the community has access to for studying the dynamics of small bodies in the solar system.
 

setup.py

@@ -31,14 +31,17 @@ class CMakeBuild(build_ext):
     """
 
     def build_extension(self, ext: CMakeExtension) -> None:
-        subprocess.run(["./build_cpp.sh"], cwd=ext.sourcedir, check=True)
+        result = subprocess.run(["./build_cpp.sh"], cwd=ext.sourcedir, check=True)
+        # check return code of subprocess
+        if result.returncode != 0:
+            raise subprocess.CalledProcessError(returncode=result.returncode, cmd=result.args)
         binary_created = [f for f in os.listdir(f"{ext.sourcedir}/build/")
                             if f.startswith("libgrss")]
         if not binary_created:
             raise FileNotFoundError("libgrss binary for C++ source code not found "
                                     "in cmake build directory")
         os.system(f"cp {ext.sourcedir}/build/libgrss.* {self.build_lib}/grss/")
-        return
+        return None
 
 setup(
     ext_modules=[CMakeExtension("libgrss")],

src/approach.cpp

@@ -63,46 +63,6 @@ static void rec_to_geodetic(const real &x, const real &y, const real &z,
     }
 }
 
-static void rec_to_geodetic_jac(const real &lon, const real &lat, const real &h,
-                                std::vector< std::vector<real> > &jac) {
-    // calculate derivatives of geodetic coordinates with respect to the
-    // rectangular coordinates by inversion of the derivatives of the
-    // rectangular coordinates with respect to the geodetic coordinates
-    const real a = EARTH_RAD_WGS84;
-    const real slat = sin(lat);
-    const real clat = cos(lat);
-    const real slon = sin(lon);
-    const real clon = cos(lon);
-    const real flat = 1.0 - EARTH_FLAT_WGS84;
-    const real flat2 = flat * flat;
-    const real g = sqrt(clat*clat + flat2*slat*slat);
-    const real g2 = g * g;
-    const real dg_dlat = (-1.0 + flat2) * slat * clat / g;
-    std::vector< std::vector<real> > jacInv = std::vector< std::vector<real> >(3, std::vector<real>(3, 0.0));
-    // partials of rectangular coordinates w.r.t. geodetic longitude
-    jacInv[0][0] = - (h + a/g) * slon * clat;
-    jacInv[1][0] = (h + a/g) * clon * clat;
-    jacInv[2][0] = 0.0;
-    // partials of rectangular coordinates w.r.t. geodetic latitude
-    jacInv[0][1] = (-a*dg_dlat/g2) * clon * clat - (h + a/g) * clon * slat;
-    jacInv[1][1] = (-a*dg_dlat/g2) * slon * clat - (h + a/g) * slon * slat;
-    jacInv[2][1] = (-flat2*a*dg_dlat/g2) * slat + (h + flat2*a/g) * clat;
-    // partials of rectangular coordinates w.r.t. geodetic altitude
-    jacInv[0][2] = clon * clat;
-    jacInv[1][2] = slon * clat;
-    jacInv[2][2] = slat;
-    // invert the Jacobian
-    std::vector< std::vector<real> > jacSmall = std::vector< std::vector<real> >(3, std::vector<real>(3, 0.0));
-    mat3_inv(jacInv, jacSmall);
-    // fill the full 2x6 Jacobian for just longitude and latitude
-    jac.resize(2, std::vector<real>(6, 0.0));
-    for (size_t i = 0; i < 2; i++) {
-        for (size_t j = 0; j < 3; j++) {
-            jac[i][j] = jacSmall[i][j];
-        }
-    }
-}
-
 /**
  * @param[inout] propSim PropSimulation object for the integration.
  * @param[in] tOld Time at the previous integrator epoch.
@@ -1030,13 +990,11 @@ void ImpactParameters::get_impact_parameters(PropSimulation *propSim){
     x = this->xRelBodyFixed[0];
     y = this->xRelBodyFixed[1];
     z = this->xRelBodyFixed[2];
-    std::vector<std::vector<real>> jac(2, std::vector<real>(6, 0.0));
     if (this->centralBodySpiceId == 399){
         rec_to_geodetic(x, y, z, lon, lat, alt);
-        rec_to_geodetic_jac(lon, lat, alt, jac);
     } else {
         const real dist = sqrt(x*x + y*y + z*z);
-        lat = atan2(z, sqrt(x*x + y*y));
+        lat = asin(z/dist);
         lon = atan2(y, x);
         if (lon < 0.0) {
             lon += 2 * PI;
@@ -1048,39 +1006,10 @@ void ImpactParameters::get_impact_parameters(PropSimulation *propSim){
             centralBodyRadius = propSim->spiceBodies[this->centralBodyIdx - propSim->integParams.nInteg].radius;
         }
         alt = dist-centralBodyRadius;
-        // calculate derivatives of latitudinal coordinates with respect to the
-        // rectangular coordinates by inversion of the derivatives of the
-        // rectangular coordinates with respect to the latitudinal coordinates
-        std::vector<std::vector<real>> jacInv(3, std::vector<real>(3, 0.0));
-        // partials of rectangular coordinates with respect to longitude
-        jacInv[0][1] = -dist * sin(lon) * cos(lat);
-        jacInv[1][1] = dist * cos(lon) * cos(lat);
-        jacInv[2][1] = 0.0;
-        // partials of rectangular coordinates with respect to latitude
-        jacInv[0][2] = -dist * cos(lon) * sin(lat);
-        jacInv[1][2] = -dist * sin(lon) * sin(lat);
-        jacInv[2][2] = dist * cos(lat);
-        // partials of rectangular coordinates with respect to radius
-        jacInv[0][3] = cos(lon) * cos(lat);
-        jacInv[1][3] = sin(lon) * cos(lat);
-        jacInv[2][3] = sin(lat);
-        std::vector<std::vector<real>> jacSmall(3, std::vector<real>(3, 0.0));
-        mat3_inv(jacInv, jacSmall);
-        for (size_t k = 0; k < 2; k++) {
-            for (size_t k2 = 0; k2 < 3; k2++) {
-                jac[k][k2] = jacSmall[k][k2];
-            }
-        }
     }
     this->lon = lon;
     this->lat = lat;
     this->alt = alt;
-    std::vector<std::vector<real>> jac_inertial(2, std::vector<real>(6, 0.0));
-    mat_mat_mul(jac, rotMat, jac_inertial);
-    for (size_t k = 0; k < 6; k++) {
-        this->dlon[k] = jac_inertial[0][k];
-        this->dlat[k] = jac_inertial[1][k];
-    }
 }
 
 /**