B-plane partials upgraded, JOSS updates

.github/workflows/joss.yml

@@ -20,7 +20,7 @@ jobs:
           journal: joss
           paper-path: ./joss/joss_paper.md
       - name: Upload
-        uses: actions/upload-artifact@v1
+        uses: actions/upload-artifact@master
         with:
           name: paper
           path: ./joss/paper.pdf

environment.yml

@@ -8,4 +8,4 @@ dependencies:
   - jupyter
   - ipykernel
   - pip:
-    -r requirements.txt
+    - -r requirements.txt

grss/fit/fit_optical.py

@@ -952,7 +952,8 @@ def eliminate_obs(obs_df, max_obs_per_night, verbose):
 
 def get_optical_obs(body_id, optical_obs_file=None, t_min_tdb=None,
                     t_max_tdb=None, debias_lowres=True, deweight=True,
-                    eliminate=False, num_obs_per_night=5, verbose=False):
+                    eliminate=False, num_obs_per_night=5, verbose=False,
+                    accept_weights=False):
     """
     Assemble the optical observations for a given body into an array for an orbit fit.
 
@@ -977,6 +978,8 @@ def get_optical_obs(body_id, optical_obs_file=None, t_min_tdb=None,
         Effective/Maximum number of effective observations per night, by default 5
     verbose : bool, optional
         Flag to print out information about the observations, by default False
+    accept_weights : bool, optional
+        Flag to accept the weights from the input file, by default False
 
     Returns
     -------
@@ -1004,7 +1007,8 @@ def get_optical_obs(body_id, optical_obs_file=None, t_min_tdb=None,
                 "Setting biases to zero.")
         opt_idx = obs_df.query("mode != 'RAD'").index
         obs_df.loc[opt_idx, ['biasRA', 'biasDec']] = 0.0
-    obs_df = apply_weighting_scheme(obs_df, verbose)
+    if not accept_weights:
+        obs_df = apply_weighting_scheme(obs_df, verbose)
     if deweight:
         obs_df = deweight_obs(obs_df, num_obs_per_night, verbose)
     elif eliminate:

grss/fit/fit_simulation.py

@@ -563,7 +563,11 @@ def __init__(self, x_init, obs_df, cov_init=None, n_iter_max=10,
         self.fixed_propsim_params['events'] = events if events is not None else []
         self.reject_outliers = True
         self.reject_criteria = [3.0, 2.8]
-        self.num_rejected = 0
+        # number of rejected obs is the count of ['D', 'd'] in selAst
+        sel_ast = self.obs['selAst'].values
+        num_auto_rejected = np.sum(sel_ast == 'D')
+        num_force_rejected = np.sum(sel_ast == 'd')
+        self.num_rejected = num_auto_rejected + num_force_rejected
         self.converged = False
         return None
 
@@ -1525,7 +1529,6 @@ def _get_rms_and_reject_outliers(self, partials, residuals, start_rejecting):
                                 "Default values are chi_reject=3.0 and chi_recover=2.8 "
                                 "(Implemented as FitSimulation.reject_criteria=[3.0, 2.8])")
         full_cov = self.covariance
-        self.num_rejected = 0
         rms_u = 0
         chi_sq = 0
         j = 0
@@ -1556,9 +1559,10 @@ def _get_rms_and_reject_outliers(self, partials, residuals, start_rejecting):
                 # outlier rejection, only reject RA/Dec measurements
                 if size == 2:
                     if residual_chi_squared > chi_reject**2 and sel_ast[i] not in {'a', 'd'}:
+                        if sel_ast[i] == 'A':
+                            self.num_rejected += 1
                         sel_ast[i] = 'D'
-                        self.num_rejected += 1
-                    elif sel_ast[i] == 'D' and residual_chi_squared < chi_recover**2:
+                    elif residual_chi_squared < chi_recover**2 and sel_ast[i] == 'D':
                         sel_ast[i] = 'A'
                         self.num_rejected -= 1
             if sel_ast[i] not in {'D', 'd'}:
@@ -1614,7 +1618,7 @@ def _check_convergence(self, delta_x):
         """
         # check for convergence based on weighted rms
         if self.n_iter > 1:
-            del_rms_convergence = 1e-3
+            del_rms_convergence = 1e-4
             curr_rms = self.iters[-1].weighted_rms
             prev_rms = self.iters[-2].weighted_rms
             del_rms = abs(prev_rms - curr_rms)#/prev_rms
@@ -1727,15 +1731,17 @@ def filter_lsq(self, verbose=True):
             self._check_convergence(delta_x)
             if self.converged:
                 if self.reject_outliers and start_rejecting:
-                    print(f"Converged after rejecting outliers. Rejected {self.num_rejected} out",
-                            f"of {len(self.optical_idx)} optical observations.")
+                    if verbose:
+                        print(f"Converged after rejecting outliers. Rejected {self.num_rejected}",
+                                f"out of {len(self.optical_idx)} optical observations.")
                     break
                 msg = "Converged without rejecting outliers."
                 if self.reject_outliers:
                     msg += " Starting outlier rejection now..."
                     start_rejecting = True
                     self.converged = False
-                print(msg)
+                if verbose:
+                    print(msg)
                 if not self.reject_outliers:
                     break
         if self.n_iter == self.n_iter_max and not self.converged:
@@ -1791,7 +1797,7 @@ def print_summary(self, iter_idx=-1, out_stream=sys.stdout):
                 str_to_print += f"{key}\t\t\t{init_val:.11e}\t\t{init_unc:.11e}"
                 str_to_print += f"\t\t{final_val:.11e}\t\t{final_unc:.11e}"
                 str_to_print += f"\t\t{final_val-init_val:+.11e}"
-                str_to_print += f"\t\t{(final_val-init_val)/init_unc:+.3f}\n"
+                str_to_print += f"\t\t{(final_val-init_val)/final_unc:+.3f}\n"
         print(str_to_print, file=out_stream)
         return None
 

grss/prop/prop_utils.py

@@ -327,8 +327,7 @@ def partials_to_ellipse(ca, au2units, n_std, print_ellipse_params, units, analyt
     units : str
         units of the plot
     analytic_info : tuple
-        Tuple of analytic B-plane information (initial covariance
-        and state conversion partial derivatives)
+        Tuple of analytic B-plane information (just initial covariance for now)
 
     Returns
     -------
@@ -341,47 +340,31 @@ def partials_to_ellipse(ca, au2units, n_std, print_ellipse_params, units, analyt
         init_cov = np.array(analytic_info[0])
     else:
         init_cov = analytic_info[0]
-    if not isinstance(analytic_info[1], np.ndarray):
-        part_cart_part_state = np.array(analytic_info[1])
-    else:
-        part_cart_part_state = analytic_info[1]
-    part_state_part_cart = np.linalg.inv(part_cart_part_state)
 
-    stm_flat = ca.xRel[6:]
-    stm_tca_t0 = np.array(libgrss.reconstruct_stm(stm_flat)) @ part_state_part_cart
     stm_flat = ca.xRelMap[6:]
-    stm_tmap_t0 = np.array(libgrss.reconstruct_stm(stm_flat)) @ part_state_part_cart
-    if ca.t == ca.tMap:
-        stm_tca_tmap = np.eye(stm_tca_t0.shape[0])
-    else:
-        stm_tca_tmap = stm_tca_t0 @ np.linalg.inv(stm_tmap_t0)
-
-    extra_zeros = [0]*(stm_tca_tmap.shape[0]-6)
-    part_tlin_minus_t = ca.dTLinMinusT + extra_zeros
-    part_tlin = part_tlin_minus_t @ stm_tca_tmap
-    part_t = ca.dt + extra_zeros
-    # part_tlin += part_t
-    part_kizner = np.zeros((3,stm_tca_tmap.shape[0]))
-    part_kizner[:2,:] = np.array([ca.kizner.dx+extra_zeros,ca.kizner.dy+extra_zeros])@stm_tca_tmap
-    part_kizner[2,:] = part_tlin
-    part_scaled = np.zeros((3,stm_tca_tmap.shape[0]))
-    part_scaled[:2,:] = np.array([ca.scaled.dx+extra_zeros,ca.scaled.dy+extra_zeros])@stm_tca_tmap
-    part_scaled[2,:] = part_tlin
-    part_opik = np.zeros((3,stm_tca_tmap.shape[0]))
-    part_opik[:2,:] = np.array([ca.opik.dx+extra_zeros,ca.opik.dy+extra_zeros])@stm_tca_tmap
-    part_opik[2,:] = part_tlin
-    part_mtp = np.zeros((3,stm_tca_tmap.shape[0]))
-    part_mtp[:2,:] = np.array([ca.mtp.dx+extra_zeros,ca.mtp.dy+extra_zeros])@stm_tca_tmap
-    part_mtp[2,:] = part_tlin
-
-    init_cart_cov = part_cart_part_state @ init_cov @ part_cart_part_state.T
-    cov_tmap = stm_tmap_t0 @ init_cart_cov @ stm_tmap_t0.T
+    stm_tmap_t0 = np.array(libgrss.reconstruct_stm(stm_flat))# @ part_state_part_cart
+
+    extra_zeros = [0]*(stm_tmap_t0.shape[0]-6)
+    part_kizner = np.zeros((3,stm_tmap_t0.shape[0]))
+    part_kizner[:2,:] = np.array([ca.kizner.dx+extra_zeros,ca.kizner.dy+extra_zeros])
+    part_kizner[2,:] = ca.dtLin + extra_zeros
+    part_scaled = np.zeros((3,stm_tmap_t0.shape[0]))
+    part_scaled[:2,:] = np.array([ca.scaled.dx+extra_zeros,ca.scaled.dy+extra_zeros])
+    part_scaled[2,:] = ca.dtLin + extra_zeros
+    part_opik = np.zeros((3,stm_tmap_t0.shape[0]))
+    part_opik[:2,:] = np.array([ca.opik.dx+extra_zeros,ca.opik.dy+extra_zeros])
+    part_opik[2,:] = ca.dtLin + extra_zeros
+    part_mtp = np.zeros((3,stm_tmap_t0.shape[0]))
+    part_mtp[:2,:] = np.array([ca.mtp.dx+extra_zeros,ca.mtp.dy+extra_zeros])
+    part_mtp[2,:] = ca.dt + extra_zeros
+
+    cov_tmap = stm_tmap_t0 @ init_cov @ stm_tmap_t0.T
     cov_kizner = part_kizner @ cov_tmap @ part_kizner.T
     cov_scaled = part_scaled @ cov_tmap @ part_scaled.T
     cov_opik = part_opik @ cov_tmap @ part_opik.T
     cov_mtp = part_mtp @ cov_tmap @ part_mtp.T
 
-    t_dev = (part_t @ cov_tmap @ part_t)**0.5
+    t_dev = cov_mtp[2,2]**0.5
 
     cov_kizner = cov_kizner[:2,:2]*au2units**2
     cov_opik = cov_opik[:2,:2]*au2units**2
@@ -865,9 +848,10 @@ def plot_earth_impact(impact_list, print_ellipse_params=False, sigma_points=None
         size = zoom_size*1e3
         m = Basemap(width=size,height=size,
                     rsphere=(6378137.00,6356752.3142),
-                    resolution='l',area_thresh=size/100,projection='lcc',
+                    resolution='i',area_thresh=size/100,projection='lcc',
                     lat_0=mean_lat,lon_0=mean_lon)
-        m.shadedrelief()
+        m.bluemarble()
+        m.drawcountries()
         lat_step = lon_step = 3 if zoom_size <= 1e3 else 10
         if zoom_size <= 500:
             lat_step = lon_step = 1
@@ -893,10 +877,12 @@ def plot_earth_impact(impact_list, print_ellipse_params=False, sigma_points=None
     x_lon,y_lat = m(lon,lat)
     m.plot(x_lon,y_lat,'r.')
     lat_half = 'N' if mean_lat > 0 else 'S'
+    if mean_lat < 0:
+        mean_lat = np.abs(mean_lat)
+    if save_name is not None:
+        plt.savefig(save_name, dpi=500, bbox_inches='tight')
     plt.suptitle('Impact Location at 100km altitude w.r.t Earth: '
                     f'{mean_lat:0.2f}$^o${lat_half}, {mean_lon:0.2f}$^o$E',
                     y=0.93)
-    if save_name is not None:
-        plt.savefig(save_name, dpi=500, bbox_inches='tight')
     plt.show()
     return None

grss/version.txt

@@ -1 +1 @@
-4.1.4
+4.2.0

include/simulation.h

@@ -296,7 +296,7 @@ struct BPlaneParameters {
  * @param opik Öpik B-plane parameters.
  * @param scaled Scaled Kizner B-plane parameters.
  * @param mtp Modified Target Plane parameters.
- * @param dTLinMinusT Partial derivatives of the (linearized intersection minus map) time with respect to the state at the close approach.
+ * @param dtLin Partial derivatives of the (linearized intersection minus map) time with respect to the state at the close approach.
  * @param dt Partial derivatives of the time of closest approach with respect to the state at the close approach.
  */
 class CloseApproachParameters {
@@ -326,7 +326,7 @@ class CloseApproachParameters {
     BPlaneParameters opik;
     BPlaneParameters scaled;
     BPlaneParameters mtp;
-    std::vector<real> dTLinMinusT = std::vector<real>(6, 0.0L);
+    std::vector<real> dtLin = std::vector<real>(6, 0.0L);
     std::vector<real> dt = std::vector<real>(6, 0.0L);
     /**
      * @brief Get the close approach parameters.

include/utilities.h

@@ -181,6 +181,12 @@ void vabs_max(const real *v, const size_t &dim, real &max);
 void mat_vec_mul(const std::vector<std::vector<real>> &A,
                  const std::vector<real> &v, std::vector<real> &Av);
 
+/**
+ * @brief Multiply a vector by a matrix.
+ */
+void vec_mat_mul(const std::vector<real> &v, real **A,
+                 const size_t &dim, std::vector<real> &vA);
+
 /**
  * @brief Multiply two matrices.
  */

joss/joss_paper.bib

@@ -1,3 +1,18 @@
+@article{Makadia2022,
+doi = {10.3847/PSJ/ac7de7},
+url = {https://dx.doi.org/10.3847/PSJ/ac7de7},
+year = {2022},
+month = {aug},
+publisher = {The American Astronomical Society},
+volume = {3},
+number = {8},
+pages = {184},
+author = {Rahil Makadia and Sabina D. Raducan and Eugene G. Fahnestock and Siegfried Eggl},
+title = {Heliocentric Effects of the DART Mission on the (65803) Didymos Binary Asteroid System},
+journal = {The Planetary Science Journal},
+abstract = {The Double Asteroid Redirect Test (DART) is NASA’s first kinetic impact–based asteroid deflection mission. The DART spacecraft will act as a projectile during a hypervelocity impact on Dimorphos, the secondary asteroid in the (65803) Didymos binary system, and alter its mutual orbital period. The initial momentum transfer between the DART spacecraft and Dimorphos is enhanced by the ejecta flung off the surface of Dimorphos. This exchange is characterized within the system by the momentum enhancement parameter, β, and on a heliocentric level by its counterpart, β ⊙. The relationship between β and the physical characteristics of Dimorphos is discussed here. A nominal set of Dimorphos physical parameters from the design reference asteroid and impact circumstances from the design reference mission are used to initialize the ejecta particles for dynamical propagation. The results of this propagation are translated into a gradual momentum transfer onto the Didymos system barycenter. A high-quality solar system propagator is then used to produce precise estimates of the post-DART encounters between Didymos and Earth by generating updated close approach maps. Results show that even for an unexpectedly high β ⊙, a collision between the Didymos system and Earth is practically excluded in the foreseeable future. A small but significant difference is found in modeling the overall momentum transfer when individual ejecta particles escape the Didymos system, as opposed to imparting the ejecta momentum as a single impulse at impact. This difference has implications for future asteroid deflection campaigns, especially when it is necessary to steer asteroids away from gravitational keyholes.}
+}
+
 @article{Makadia2024,
 author    = {Rahil Makadia and Steven R. Chesley and Davide Farnocchia and Shantanu P. Naidu and Damya Souami and Paolo Tanga and Kleomenis Tsiganis and Masatoshi Hirabayashi and Siegfried Eggl},
 journal   = {The Planetary Science Journal},
@@ -56,13 +71,6 @@ @inproceedings{Chodas1999
 url = {https://hdl.handle.net/2014/16816}
 }
 
-@misc{pybind11,
-author = {Wenzel Jakob and Jason Rhinelander and Dean Moldovan},
-year = {2017},
-url = {https://github.com/pybind/pybind11},
-title = {pybind11 -- Seamless operability between C++11 and Python}
-}
-
 @article{Milani1999,
 title = {The Asteroid Identification Problem: II. Target Plane Confidence Boundaries},
 journal = {Icarus},
@@ -71,7 +79,7 @@ @article{Milani1999
 pages = {408-423},
 year = {1999},
 issn = {0019-1035},
-doi = {https://doi.org/10.1006/icar.1999.6135},
+doi = {10.1006/icar.1999.6135},
 author = {Andrea Milani and Giovanni B. Valsecchi},
 }