added integrator upgrades and close approach calcs

.github/workflows/cpp_tests.yml

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

.github/workflows/joss.yml

@@ -0,0 +1,25 @@
+name: JOSS build
+on:
+  push:
+    branches:
+      - main
+      - dev
+
+jobs:
+  paper:
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v3
+        with:
+          submodules: true
+      - name: Build draft PDF
+        uses: openjournals/openjournals-draft-action@master
+        with:
+          journal: joss
+          paper-path: ./joss/joss_paper.md
+      - name: Upload
+        uses: actions/upload-artifact@v1
+        with:
+          name: paper
+          path: ./joss/paper.pdf

.github/workflows/python_tests.yml

@@ -26,7 +26,8 @@ jobs:
         run: |
           python3 -m pip install jupyter ipykernel
           python3 -m ipykernel install --user
-          jupyter nbconvert --to script tests/python/prop/*.ipynb
-          python3 tests/python/prop/*.py
-          jupyter nbconvert --to script tests/python/fit/*.ipynb
-          python3 tests/python/fit/*.py
+          cd tests/python
+          jupyter nbconvert --to script fit/*.ipynb
+          jupyter nbconvert --to script prop/*.ipynb
+          python fit/*.py
+          python prop/*.py

.gitignore

@@ -190,6 +190,9 @@ cython_debug/
 # gperftools profiler output
 *.prof
 
+# pytest settings are exempt from ignoring
+!pytest.ini
+
 # VS code setup directory
 .vscode/
 

CMakeLists.txt

@@ -10,7 +10,8 @@ project(grss VERSION ${ver})
 # TODO: add gprof profiling, iee standards for floats
 set(CMAKE_CXX_STANDARD 11)
 set(CMAKE_CXX_STANDARD_REQUIRED True)
-add_compile_options(-std=c++11 -DLONGDOUBLE -O2 -fPIC -Werror -Wall -Wextra -pedantic -Wno-unused-but-set-variable -Wno-unused-result)
+add_compile_options(-std=c++11 -O3 -fPIC) # operational flags
+# add_compile_options(-std=c++11 -DLONGDOUBLE -O2 -fPIC -Werror -Wall -Wextra -pedantic -Wno-unused-but-set-variable -Wno-unused-result)
 # if not clang, add -Wno-maybe-uninitialized
 if(NOT CMAKE_CXX_COMPILER_ID MATCHES "Clang")
     add_compile_options(-Wno-maybe-uninitialized)

docs/source/api.rst

@@ -5,10 +5,6 @@ API
    :template: custom-module.rst
    :recursive:
 
-   grss.fit_optical
-   grss.fit_radar
-   grss.fit_simulation
-   grss.fit_utils
-   grss.prop_simulation
-   grss.prop_utils
+   grss.fit
+   grss.prop
    grss.utils

grss/__init__.py

@@ -1,12 +1,9 @@
 """GRSS: the Gauss-Radau Small-body Simulator"""
-from . import prop
 from . import fit
+from . import prop
 from . import utils
 
-__all__ = ['prop', 'fit', 'utils']
-
 # get version from version.txt
-with open(f'{utils.grss_path}/version.txt', 'r', encoding='utf-8') as ver_file:
-    __version__ = ver_file.read().strip()
+__version__ = open(f'{utils.grss_path}/version.txt', 'r', encoding='utf-8').read().strip()
 
 utils.initialize()

grss/fit.py → grss/fit/__init__.py

@@ -1,6 +1,5 @@
-"""GRSS orbit determination code"""
-#pylint: disable=wildcard-import,unused-wildcard-import
-from .fit_utils import *
+"""GRSS orbit determination subpackage"""
 from .fit_optical import *
 from .fit_radar import *
 from .fit_simulation import *
+from .fit_utils import *

grss/fit_optical.py → grss/fit/fit_optical.py

@@ -1,5 +1,4 @@
 """Optical observation handling for the GRSS orbit determination code"""
-import os
 from astropy.time import Time
 from astroquery.mpc import MPC
 from astroquery.gaia import Gaia
@@ -9,6 +8,7 @@
 import spiceypy as spice
 
 from .fit_utils import get_ra_from_hms, get_dec_from_dms, radec2icrf, icrf2radec, rec2lat
+from .. import utils
 
 __all__ = [ 'get_ades_optical_obs_array',
             'get_gaia_optical_obs_array',
@@ -238,7 +238,8 @@ def _get_gaia_query_results(body_id, release='gaiadr3', verbose=False):
         print(f"Found {len(res)} observations from Gaia DR3.")
     return res
 
-def get_gaia_optical_obs_array(body_id, de_kernel_path, t_min_tdb=None, t_max_tdb=None, verbose=False):
+def get_gaia_optical_obs_array(body_id, de_kernel_path, t_min_tdb=None,
+                                t_max_tdb=None, verbose=False):
     """
     Assemble the optical observations for a given body from Gaia DR3.
 
@@ -282,12 +283,24 @@ def get_gaia_optical_obs_array(body_id, de_kernel_path, t_min_tdb=None, t_max_td
         cosdec = np.cos(data['dec']*np.pi/180)
         obs_array[i, 0] = obs_time.utc.mjd
         obs_array[i, 1] = data['ra']*3600
-        obs_array[i, 1]-= data['ra_error_systematic']/cosdec/1000
         obs_array[i, 2] = data['dec']*3600
-        obs_array[i, 2]-= data['dec_error_systematic']/1000
-        obs_array[i, 3] = data['ra_error_random']/cosdec/1000
-        obs_array[i, 4] = data['dec_error_random']/1000
-        obs_array[i, 5] = data['ra_dec_correlation_random']
+        ra_sig_sys = data['ra_error_systematic']/cosdec/1000
+        ra_sig_rand = data['ra_error_random']/cosdec/1000
+        dec_sig_sys = data['dec_error_systematic']/1000
+        dec_sig_rand = data['dec_error_random']/1000
+        corr_sys = data['ra_dec_correlation_systematic']
+        corr_rand = data['ra_dec_correlation_random']
+        cov_sys = np.array([[ra_sig_sys**2, corr_sys*ra_sig_sys*dec_sig_sys],
+                            [corr_sys*ra_sig_sys*dec_sig_sys, dec_sig_sys**2]])
+        cov_rand = np.array([[ra_sig_rand**2, corr_rand*ra_sig_rand*dec_sig_rand],
+                            [corr_rand*ra_sig_rand*dec_sig_rand, dec_sig_rand**2]])
+        cov = cov_sys + cov_rand
+        ra_sig = cov[0,0]**0.5
+        dec_sig = cov[1,1]**0.5
+        corr = cov[0,1]/ra_sig/dec_sig
+        obs_array[i, 3] = ra_sig
+        obs_array[i, 4] = dec_sig
+        obs_array[i, 5] = corr
         pos_x = data['x_gaia_geocentric']*au2km
         pos_y = data['y_gaia_geocentric']*au2km
         pos_z = data['z_gaia_geocentric']*au2km
@@ -391,9 +404,8 @@ def apply_debiasing_scheme(obs_array_optical, star_catalog_codes, observer_codes
     columns = []
     for cat in biased_catalogs:
         columns.extend([f"{cat}_ra", f"{cat}_dec", f"{cat}_pm_ra", f'{cat}_pm_dec'])
-    cwd = os.path.dirname(os.path.realpath(__file__))
-    debias_lowres_path = os.path.join(cwd,'debias/lowres_data')
-    debias_hires_path = os.path.join(cwd,'debias/hires_data')
+    debias_lowres_path = utils.grss_path+'/debias/lowres_data'
+    debias_hires_path = utils.grss_path+'/debias/hires_data'
     if lowres:
         biasfile = f'{debias_lowres_path}/bias.dat'
         nside = 64

grss/fit_simulation.py → grss/fit/fit_simulation.py

@@ -1,10 +1,11 @@
 """Simulation classes for the Python GRSS orbit determination code"""
-#pylint: disable=no-member, too-many-lines
+# pylint: disable=no-name-in-module, too-many-lines
 import numpy as np
 import matplotlib.pyplot as plt
 import spiceypy as spice
+from astropy.time import Time
 
-from . import prop
+from .. import prop
 from .fit_utils import get_observer_info, get_radec
 
 __all__ = [ 'FitSimulation',
@@ -288,7 +289,7 @@ def plot_residuals(self, t_arr, ra_residuals, dec_residuals, ra_cosdec_residuals
         ax1.plot(t_arr[is_rejected], dec_residuals[is_rejected], 'ro',
                     markersize=2*markersize, markerfacecolor='none')
         ax1.legend()
-        ax1.set_xlabel('MJD [UTC]')
+        ax1.set_xlabel('Time [UTC]')
         ax1.set_ylabel('Residuals, O-C [arcsec]')
         ax1.grid(True, which='both', axis='both', alpha=0.2)
         if show_logarithmic:
@@ -314,7 +315,7 @@ def plot_residuals(self, t_arr, ra_residuals, dec_residuals, ra_cosdec_residuals
         ax3.plot(t_arr, doppler_residuals, '.', mfc='C3', mec='C3', label='Doppler',
                     markersize=radar_scale*markersize)
         ax3.legend()
-        ax3.set_xlabel('MJD [UTC]')
+        ax3.set_xlabel('Time [UTC]')
         ax3.set_ylabel(r'Residuals, O-C [$\mu$s, Hz]')
         ax3.grid(True, which='both', axis='both', alpha=0.2)
         if show_logarithmic:
@@ -411,28 +412,31 @@ def plot_chi(self, t_arr, ra_chi, dec_chi, delay_chi, doppler_chi,
         plt.suptitle(msg, y=0.95)
         if plot_chi_squared:
             plt.subplot(1,2,1)
-        plt.plot(t_arr, ra_chi, '.', markersize=markersize, label='RA')
-        plt.plot(t_arr, dec_chi, '.', markersize=markersize, label='Dec')
-        plt.plot(t_arr[is_rejected], ra_chi[is_rejected], 'ro',
-                    markersize=2*markersize, markerfacecolor='none')
-        plt.plot(t_arr[is_rejected], dec_chi[is_rejected], 'ro',
-                    markersize=2*markersize, markerfacecolor='none')
-        plt.plot(t_arr, delay_chi, '.', mfc='C2', mec='C2',
-                    markersize=radar_scale*markersize, label='Delay')
-        plt.plot(t_arr, doppler_chi, '.', mfc='C3', mec='C3',
-                    markersize=radar_scale*markersize, label='Doppler')
-        plt.axhline(-sigma_limit, c='khaki', linestyle='--', alpha=1.0,
+        if not np.all(np.isnan(ra_chi)) and not np.all(np.isnan(dec_chi)):
+            plt.plot(t_arr, ra_chi, '.', markersize=markersize, label='RA')
+            plt.plot(t_arr, dec_chi, '.', markersize=markersize, label='Dec')
+            plt.plot(t_arr[is_rejected], ra_chi[is_rejected], 'ro',
+                        markersize=2*markersize, markerfacecolor='none')#, label='Rejected Obs.')
+            plt.plot(t_arr[is_rejected], dec_chi[is_rejected], 'ro',
+                        markersize=2*markersize, markerfacecolor='none')
+        if not np.all(np.isnan(delay_chi)):
+            plt.plot(t_arr, delay_chi, '.', mfc='C2', mec='C2',
+                        markersize=radar_scale*markersize, label='Delay')
+        if not np.all(np.isnan(doppler_chi)):
+            plt.plot(t_arr, doppler_chi, '.', mfc='C3', mec='C3',
+                        markersize=radar_scale*markersize, label='Doppler')
+        plt.axhline(-sigma_limit, c='red', linestyle='--', alpha=0.5,
                         label=fr'$\pm{sigma_limit:.0f}\sigma$')
-        plt.axhline(sigma_limit, c='khaki', linestyle='--', alpha=1.0)
-        plt.axhline(-2*sigma_limit, c='red', linestyle='--', alpha=0.5,
-                        label=fr'$\pm{2*sigma_limit:.0f}\sigma$')
-        plt.axhline(2*sigma_limit, c='red', linestyle='--', alpha=0.5)
-        plt.legend(ncol=3)
-        plt.xlabel('MJD [UTC]')
-        plt.ylabel(r'$\chi$, (O-C)/$\sigma$ $[\cdot]$')
+        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.grid(True, which='both', axis='both', alpha=0.2)
         if show_logarithmic:
             plt.yscale('log')
+        else:
+            lim = np.max(np.abs(plt.ylim()))
+            plt.ylim(-lim, lim)
         if plot_chi_squared:
             plt.subplot(1,2,2)
             plt.plot(t_arr, ra_chi_squared, '.', markersize=markersize, label='RA')
@@ -446,7 +450,7 @@ def plot_chi(self, t_arr, ra_chi, dec_chi, delay_chi, doppler_chi,
             plt.plot(t_arr, doppler_chi_squared, '.', mfc='C3', mec='C3',
                         markersize=radar_scale*markersize, label='Doppler')
             plt.legend(ncol=2)
-            plt.xlabel('MJD [UTC]')
+            plt.xlabel('Time [UTC]')
             plt.ylabel(r'$\chi^2$, (O-C)$^2/\sigma^2$ $[\cdot]$')
             plt.grid(True, which='both', axis='both', alpha=0.2)
             # if show_logarithmic: plt.yscale('log')
@@ -493,6 +497,7 @@ def plot_iteration_summary(self, show_logarithmic=False, title=None, plot_chi_sq
         plt.rcParams['font.size'] = 12
         plt.rcParams['axes.labelsize'] = 12
         t_arr = self.obs_array[:, 0]
+        t_arr = Time(t_arr, format='mjd', scale='utc').utc.datetime
         if show_logarithmic:
             ra_residuals = np.abs(self.all_info['ra_res'])
             dec_residuals = np.abs(self.all_info['dec_res'])

grss/prop/__init__.py

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

grss/utils.py

@@ -1,4 +1,4 @@
-"""General utilities for the GRSS package"""
+"""GRSS utilities subpackage"""
 import os
 import sys
 import datetime
@@ -52,7 +52,7 @@ def _download_obs_codes_file():
     req = request("GET", url, timeout=30)
     text = req.content.replace(b'<pre>',b'').replace(b'</pre>',b'')
 
-    with open(f'{grss_path}/obs_codes.txt', "wb") as file:
+    with open(f'{grss_path}/fit/obs_codes.txt', "wb") as file:
         now = datetime.datetime.now()
         file.write((f'{now.year}, {now.month}, {now.day}, '
                     f'{now.hour}, {now.minute}, {now.second}\n'.encode()))
@@ -67,8 +67,8 @@ def get_obs_codes_file():
     week ago.
     """
     # check if file exists
-    if os.path.isfile(f'{grss_path}/obs_codes.txt'):
-        with open(f'{grss_path}/obs_codes.txt', 'r', encoding='utf-8') as file:
+    if os.path.isfile(f'{grss_path}/fit/obs_codes.txt'):
+        with open(f'{grss_path}/fit/obs_codes.txt', 'r', encoding='utf-8') as file:
             now = datetime.datetime.now()
             data = file.read()
             last_updated = data.split('\n')[0]
@@ -91,12 +91,10 @@ def initialize():
     None : NoneType
         None
     """
-    # get the file directory
-    cwd = os.path.dirname(os.path.abspath(__file__))
     # run the get_debiasing_data.py script
-    os.system(f'python {cwd}/debias/get_debiasing_data.py')
+    os.system(f'python {grss_path}/debias/get_debiasing_data.py')
     # run the get_kernels.py script
-    os.system(f'python {cwd}/kernels/get_kernels.py')
+    os.system(f'python {grss_path}/kernels/get_kernels.py')
     # get the observatory codes file
     get_obs_codes_file()
     return None

grss/version.txt

@@ -1 +1 @@
-0.9.8
+0.9.9

include/force.h

@@ -4,11 +4,12 @@
 
 std::vector<real> get_state_der(const real &t, const std::vector<real> &xInteg,
                                 propSimulation *propSim);
-void force_newton(const propSimulation *propSim, real *accInteg);
-void force_ppn_simple(const propSimulation *propSim, real *accInteg);
-void force_ppn_eih(const propSimulation *propSim, real *accInteg);
-void force_J2(const propSimulation *propSim, real *accInteg);
-void force_nongrav(const propSimulation *propSim, real *accInteg);
-void force_thruster(const propSimulation *propSim, real *accInteg);
+void force_newton(const propSimulation *propSim, std::vector<real> &accInteg);
+void force_ppn_simple(const propSimulation *propSim,
+                      std::vector<real> &accInteg);
+void force_ppn_eih(const propSimulation *propSim, std::vector<real> &accInteg);
+void force_J2(const propSimulation *propSim, std::vector<real> &accInteg);
+void force_nongrav(const propSimulation *propSim, std::vector<real> &accInteg);
+void force_thruster(const propSimulation *propSim, std::vector<real> &accInteg);
 
 #endif

include/gr15.h

@@ -37,9 +37,8 @@ real get_initial_timestep(const real &t, const std::vector<real> &xInteg0,
 void compute_g_and_b(const std::vector<std::vector<real> > &AccIntegArr,
                      const size_t &hIdx, real *g,
                      std::vector<std::vector<real> > &b, const size_t &dim);
-void refine_b(std::vector<std::vector<real> > &b,
-              real *e, const real &dtRatio,
-              const size_t &dim, const size_t &timestepCounter);
+void refine_b(std::vector<std::vector<real> > &b, real *e, const real &dtRatio,
+              const size_t &dim);
 void check_and_apply_events(propSimulation *propSim, const real &t,
                             real &tNextEvent, size_t &nextEventIdx,
                             std::vector<real> &xInteg);

include/simulation.h

@@ -89,6 +89,7 @@ class propSimulation {
 
     // memory-mapped ephemeris
     Ephemeris ephem;
+    std::vector<real> get_spiceBody_state(const real t, const std::string &bodyName);
 
     // constants
     Constants consts;