Merge pull request #58 from rahil-makadia/dev

binary PCK maps are HERE🎉; improved parallel postprocessing; analytic B-plane partials

.github/workflows/cpp_tests.yml

@@ -23,6 +23,7 @@ jobs:
       - name: Run tests
         run: |
           cd tests/cpp/prop
+          ./pck_map.out
+          ./spk_map.out
           ./apophis.out
           ./didymos.out
-          ./spk_map.out

docs/doxygen/doxygen.config

@@ -61,13 +61,13 @@ PROJECT_BRIEF          = "Gauss-Radau Small-body Simulator"
 # pixels and the maximum width should not exceed 200 pixels. Doxygen will copy
 # the logo to the output directory.
 
-PROJECT_LOGO           = /Users/rahil/grad/projects/keyholes/grss/docs/source/_static/grss-doxygen.png
+PROJECT_LOGO           = ../source/_static/grss-doxygen.png
 
 # With the PROJECT_ICON tag one can specify an icon that is included in the tabs
 # when the HTML document is shown. Doxygen will copy the logo to the output
 # directory.
 
-PROJECT_ICON           = /Users/rahil/grad/projects/keyholes/grss/docs/source/_static/grss-square.svg
+PROJECT_ICON           = ../source/_static/grss-square.svg
 
 # The OUTPUT_DIRECTORY tag is used to specify the (relative or absolute) path
 # into which the generated documentation will be written. If a relative path is

docs/source/api_cpp.rst

@@ -10,6 +10,7 @@ C++ API
    cppsummary/gr15
    cppsummary/interpolate
    cppsummary/parallel
+   cppsummary/pck
    cppsummary/simulation
    cppsummary/spk
    cppsummary/stm
@@ -39,19 +40,22 @@ C++ API
    <tr class="row-even"><td><p><a class="reference internal" href="cppsummary/parallel.html" title="parallel.h"><code class="xref py py-obj docutils literal notranslate"><span class="pre">parallel.h</span></code></a></p></td>
    <td><p>GRSS C++ parallel propagation submodule</p></td>
    </tr>
-   <tr class="row-odd"><td><p><a class="reference internal" href="cppsummary/simulation.html" title="simulation.h"><code class="xref py py-obj docutils literal notranslate"><span class="pre">simulation.h</span></code></a></p></td>
+   <tr class="row-odd"><td><p><a class="reference internal" href="cppsummary/pck.html" title="pck.h"><code class="xref py py-obj docutils literal notranslate"><span class="pre">pck.h</span></code></a></p></td>
+   <td><p>GRSS C++ PCK file submodule</p></td>
+   </tr>
+   <tr class="row-even"><td><p><a class="reference internal" href="cppsummary/simulation.html" title="simulation.h"><code class="xref py py-obj docutils literal notranslate"><span class="pre">simulation.h</span></code></a></p></td>
    <td><p>GRSS C++ propagator simulation submodule</p></td>
    </tr>
-   <tr class="row-even"><td><p><a class="reference internal" href="cppsummary/spk.html" title="spk.h"><code class="xref py py-obj docutils literal notranslate"><span class="pre">spk.h</span></code></a></p></td>
+   <tr class="row-odd"><td><p><a class="reference internal" href="cppsummary/spk.html" title="spk.h"><code class="xref py py-obj docutils literal notranslate"><span class="pre">spk.h</span></code></a></p></td>
    <td><p>GRSS C++ SPK file submodule</p></td>
    </tr>
-   <tr class="row-odd"><td><p><a class="reference internal" href="cppsummary/stm.html" title="stm.h"><code class="xref py py-obj docutils literal notranslate"><span class="pre">stm.h</span></code></a></p></td>
+   <tr class="row-even"><td><p><a class="reference internal" href="cppsummary/stm.html" title="stm.h"><code class="xref py py-obj docutils literal notranslate"><span class="pre">stm.h</span></code></a></p></td>
    <td><p>GRSS C++ state transition matrix submodule</p></td>
    </tr>
-   <tr class="row-even"><td><p><a class="reference internal" href="cppsummary/timeconvert.html" title="timeconvert.h"><code class="xref py py-obj docutils literal notranslate"><span class="pre">timeconvert.h</span></code></a></p></td>
+   <tr class="row-odd"><td><p><a class="reference internal" href="cppsummary/timeconvert.html" title="timeconvert.h"><code class="xref py py-obj docutils literal notranslate"><span class="pre">timeconvert.h</span></code></a></p></td>
    <td><p>GRSS C++ time conversion submodule</p></td>
    </tr>
-   <tr class="row-odd"><td><p><a class="reference internal" href="cppsummary/utilities.html" title="utilities.h"><code class="xref py py-obj docutils literal notranslate"><span class="pre">utilities.h</span></code></a></p></td>
+   <tr class="row-even"><td><p><a class="reference internal" href="cppsummary/utilities.html" title="utilities.h"><code class="xref py py-obj docutils literal notranslate"><span class="pre">utilities.h</span></code></a></p></td>
    <td><p>GRSS C++ utilities submodule</p></td>
    </tr>
    </tbody>

docs/source/cppsummary/pck.rst

@@ -0,0 +1,6 @@
+pck.h
+=====
+
+.. doxygenfile:: pck.h
+   :project: GRSS
+

docs/source/examples.rst

@@ -16,3 +16,5 @@ Examples
    tests/python/fit/eggl.ipynb
    tests/python/fit/farnocchia.ipynb
    tests/python/fit/shantanunaidu.ipynb
+
+Note: The OD examples show a statistically significant difference between GRSS and JPL SBDB orbit solutions. This stems from the fact that the two handle astrometry from the Gaia Focused Product Release (FPR) differently. These differences will be reconciled in the future.

grss/fit/fit_simulation.py

@@ -1019,7 +1019,7 @@ def _x_dict_to_state(self, x_dict):
             arg_peri = x_dict['w']
             inc = x_dict['i']
             cometary_elements = [ecc, peri_dist, time_peri,
-                                    omega*np.pi/180.0, arg_peri*np.pi/180.0, inc*np.pi/180.0]
+                                    omega, arg_peri, inc]
             state = cometary_elements
         else:
             raise ValueError("fit_cartesian or fit_cometary must be True")
@@ -1147,11 +1147,13 @@ def _get_perturbed_state(self, key):
             elif key in {'vx', 'vy', 'vz'}:
                 fd_pert = 1e-10
         elif self.fit_cometary:
-            fd_pert = 1e-7
+            fd_pert = 1e-9
             if key == 'e':
                 fd_pert = 7.5e-9
             elif key == 'q':
                 fd_pert = 1e-9
+            elif key == 'tp':
+                fd_pert = 1e-7
         if key in {'a1', 'a2', 'a3'}:
             fd_pert = 1e-12
         if key[:4] == 'mult':
@@ -1329,19 +1331,8 @@ def _get_computed_obs(self, prop_sim_past, prop_sim_future, integ_body_idx):
         return computed_obs
 
     def _get_analytic_stm(self, t_eval, prop_sim):
-        stm_state_full = np.array(prop_sim.interpolate(t_eval))
-        stm = stm_state_full[6:42].reshape((6,6))
-        if self.fit_cometary:
-            stm[:, 3:6] /= 180.0/np.pi # covert partial w.r.t. rad -> partial w.r.t. deg
-        if len(stm_state_full) > 42:
-            param_block = stm_state_full[42:].reshape((6, -1), order='F')
-            stm = np.hstack((stm, param_block))
-            num_params = param_block.shape[1]
-            if num_params > 0:
-                bottom_block = np.zeros((num_params, 6+num_params))
-                bottom_block[:,6:] = np.eye(num_params)
-                stm = np.vstack((stm, bottom_block))
-        return stm
+        stm = prop_sim.interpolate(t_eval)[6:]
+        return libgrss.reconstruct_stm(stm)
 
     def _get_analytic_partials(self, prop_sim_past, prop_sim_future):
         """
@@ -1763,10 +1754,19 @@ def print_summary(self, iter_idx=-1):
         final_sol = data.x_nom
         with np.errstate(divide='ignore'):
             for i, key in enumerate(init_sol.keys()):
-                print(f"{key}\t\t\t{init_sol[key]:.11e}\t\t{init_variance[i]:.11e}",
-                        f"\t\t{final_sol[key]:.11e}\t\t{final_variance[i]:.11e}",
-                        f"\t\t{final_sol[key]-init_sol[key]:+.11e}"
-                        f"\t\t{(final_sol[key]-init_sol[key])/init_variance[i]:+.3f}")
+                init_val = init_sol[key]
+                init_unc = init_variance[i]
+                final_val = final_sol[key]
+                final_unc = final_variance[i]
+                if self.fit_cometary and key in ['om', 'w', 'i']:
+                    init_val *= 180/np.pi
+                    init_unc *= 180/np.pi
+                    final_val *= 180/np.pi
+                    final_unc *= 180/np.pi
+                print(f"{key}\t\t\t{init_val:.11e}\t\t{init_unc:.11e}",
+                        f"\t\t{final_val:.11e}\t\t{final_unc:.11e}",
+                        f"\t\t{final_val-init_val:+.11e}"
+                        f"\t\t{(final_val-init_val)/init_unc:+.3f}")
         return None
 
     def plot_summary(self, auto_close=False):
@@ -1938,7 +1938,7 @@ def _generate_simulated_obs(ref_sol, ref_cov, ref_ng_info, events, modified_obs_
         arg_peri = ref_sol['w']
         inc = ref_sol['i']
         cometary_elements = [ecc, peri_dist, time_peri,
-                                omega*np.pi/180.0, arg_peri*np.pi/180.0, inc*np.pi/180.0]
+                                omega, arg_peri, inc]
         target_body = libgrss.IntegBody("body_simulated_obs", ref_sol['t'],
                                         ref_sol['mass'], ref_sol['radius'], cometary_elements,
                                         ref_cov, nongrav_params)

grss/fit/fit_utils.py

@@ -415,8 +415,8 @@ def get_sbdb_elems(tdes, cov_elems=True):
         elements[key] = full_elements_dict[key]
         if key == 'tp':
             elements[key] = full_elements_dict[key] - 2400000.5
-        # if key in {'om', 'w', 'i'}:
-        #     elements[key] = full_elements_dict[key]*np.pi/180
+        if key in {'om', 'w', 'i'}:
+            elements[key] = full_elements_dict[key]*np.pi/180
     return elements
 
 def get_sbdb_info(tdes):
@@ -445,6 +445,9 @@ def get_sbdb_info(tdes):
     # covariance matrix for cometary orbital elements
     cov_keys = raw_data['orbit']['covariance']['labels']
     cov_mat = (np.array(raw_data['orbit']['covariance']['data'])).astype(float)
+    # convert covariance matrix angle blocks from degrees to radians
+    cov_mat[3:6, :] *= np.pi/180
+    cov_mat[:, 3:6] *= np.pi/180
     # nongravitatinoal constants for target body
     nongrav_data = raw_data['orbit']['model_pars']
     hdr = [param['name'] for param in nongrav_data]

grss/kernels/planets_big16_de431_1950_2350.tm

@@ -11,8 +11,7 @@
 
     PATH_SYMBOLS = ( 'GRSS' )
 
-    KERNELS_TO_LOAD = ( '$GRSS/latest_leapseconds.tls',
-                        '$GRSS/pck00011.tpc',
+    KERNELS_TO_LOAD = ( '$GRSS/pck00011.tpc',
                         '$GRSS/earth_200101_990825_predict.bpc',
                         '$GRSS/earth_720101_230601.bpc',
                         '$GRSS/earth_latest_high_prec.bpc' )

grss/kernels/planets_big16_de441_1950_2350.tm

@@ -11,8 +11,7 @@
 
     PATH_SYMBOLS = ( 'GRSS' )
 
-    KERNELS_TO_LOAD = ( '$GRSS/latest_leapseconds.tls',
-                        '$GRSS/pck00011.tpc',
+    KERNELS_TO_LOAD = ( '$GRSS/pck00011.tpc',
                         '$GRSS/earth_200101_990825_predict.bpc',
                         '$GRSS/earth_720101_230601.bpc',
                         '$GRSS/earth_latest_high_prec.bpc' )