added b-plane plotting routines

grss/fit/fit_simulation.py

@@ -430,7 +430,7 @@ def plot_chi(self, t_arr, ra_chi, dec_chi, delay_chi, doppler_chi,
         plt.axhline(sigma_limit, c='red', linestyle='--', alpha=0.5)
         plt.legend(ncol=2)
         plt.xlabel('Time [UTC]')
-        plt.ylabel(r'Weighted Residuals, (O-C)/$\sigma$ $[\cdot]$')
+        plt.ylabel(r'Weighted Residuals, (O-C) $[\sigma]$')
         plt.grid(True, which='both', axis='both', alpha=0.2)
         if show_logarithmic:
             plt.yscale('log')
@@ -451,7 +451,7 @@ def plot_chi(self, t_arr, ra_chi, dec_chi, delay_chi, doppler_chi,
                         markersize=radar_scale*markersize, label='Doppler')
             plt.legend(ncol=2)
             plt.xlabel('Time [UTC]')
-            plt.ylabel(r'$\chi^2$, (O-C)$^2/\sigma^2$ $[\cdot]$')
+            plt.ylabel(r'$\chi^2$')
             plt.grid(True, which='both', axis='both', alpha=0.2)
             # if show_logarithmic: plt.yscale('log')
             plt.yscale('log')

grss/prop/__init__.py

@@ -1,3 +1,4 @@
 """GRSS orbit propagation subpackage"""
 from .prop_simulation import *
+from .prop_unscented import *
 from .prop_utils import *

grss/prop/prop_unscented.py

@@ -0,0 +1,36 @@
+"""Unscented transformations for the GRSS orbit propagation code"""
+import numpy as np
+
+def generate_sigma_points(sol, cov):
+    """
+    Generate sigma points for propagating uncertainties using the unscented
+    transformation.
+
+    Parameters
+    ----------
+    sol : dict
+        solution dictionary from the SBDB API
+    cov : numpy.ndarray
+        covariance matrix of the solution
+
+    Returns
+    -------
+    sigma_points : list
+        list of dictionaries containing each of the sigma points
+    """
+    assert len(sol)-1 == cov.shape[0]
+    assert cov.shape[0] == cov.shape[1]
+    dim = cov.shape[0]
+    sqrt_cov = np.linalg.cholesky(cov)
+    sol_array = np.array([sol[key] for key in sol if key != 't'])
+    sigma_points = [list(sol_array)]
+    fac = dim**0.5
+    for i in range(dim):
+        plus = sol_array + fac*sqrt_cov.T[i]
+        minus = sol_array - fac*sqrt_cov.T[i]
+        sigma_points.append(list(plus))
+        sigma_points.append(list(minus))
+    # convert sigma points to dict
+    for i in range(2*dim+1):
+        sigma_points[i] = {key: sigma_points[i][j-1] for j, key in enumerate(sol) if key != 't'}
+    return sigma_points