Merge pull request #67 from rahil-makadia/dev

computing radius from SBDB info; orbital element updates
rahil-makadia · May 10, 2024 · b44f0c3 · b44f0c3
2 parents e63ab24 + 9ff13f9
commit b44f0c3
Show file tree

Hide file tree

Showing 26 changed files with 870 additions and 840 deletions.
diff --git a/.github/workflows/cpp_tests.yml b/.github/workflows/cpp_tests.yml
@@ -5,6 +5,7 @@ on:
       - main
       - dev
     paths:
+      - '.github/workflows/cpp_tests.yml'
       - 'include/**'
       - 'src/**'
       - 'tests/cpp/**'

diff --git a/.github/workflows/joss.yml b/.github/workflows/joss.yml
@@ -5,6 +5,7 @@ on:
       - main
       - dev
     paths:
+      - '.github/workflows/joss.yml'
       - 'joss/**'
 
 jobs:

diff --git a/.github/workflows/pypi.yml b/.github/workflows/pypi.yml
@@ -24,7 +24,7 @@ jobs:
       - name: Set up Python
         uses: actions/setup-python@v5
         with:
-          python-version: "3.10"
+          python-version: "3.11"
       - name: Set up Xcode
         if: matrix.os == 'macos-14'
         uses: maxim-lobanov/setup-xcode@v1.6.0

diff --git a/.github/workflows/python_tests.yml b/.github/workflows/python_tests.yml
@@ -5,25 +5,30 @@ on:
       - main
       - dev
     paths:
+      - '.github/workflows/python_tests.yml'
       - 'include/**'
       - 'src/**'
       - 'grss/**'
       - 'tests/python/**'
 
 jobs:
   prop-and-fit:
-    runs-on: ubuntu-latest
+    runs-on: macos-14
     steps:
       - name: Checkout repo
         uses: actions/checkout@v4
       - name: Set up Python
         uses: actions/setup-python@v5
         with:
-          python-version: 3.8
+          python-version: 3.11
+      - name: Set up Xcode
+        uses: maxim-lobanov/setup-xcode@v1.6.0
+        with:
+          xcode-version: '15.1.0'
       - name: Install GRSS (including dependencies)
         run: |
           python3 -m pip install --upgrade pip
-          sudo apt-get install doxygen
+          brew install doxygen
           source initialize.sh
           pip3 install .
       - name: Run tests

diff --git a/grss/fit/fit_ades.py b/grss/fit/fit_ades.py
@@ -30,7 +30,7 @@
     'sigDelay': 'float', 'sigDoppler': 'float',
 }
 ades_grss_columns = {
-    'obsTimeMJD': 'float', 'cosDec': 'float',
+    'obsTimeMJD': 'float', 'obsTimeMJDTDB': 'float', 'cosDec': 'float',
 }
 ades_column_types = ades_keep_columns | ades_add_columns | ades_grss_columns
 

diff --git a/grss/fit/fit_optical.py b/grss/fit/fit_optical.py
@@ -136,11 +136,13 @@ def create_optical_obs_df(body_id, optical_obs_file=None, t_min_tdb=None,
         if not os.path.exists(obs_data):
             raise FileNotFoundError(f"File {obs_data} does not exist.")
     # read in the data and add extra columns if not present
-    obs_df = (
-        pd.read_xml(obs_data, dtype=ades_column_types)
-        # compute the obsTimeMJD time from obsTime
-        .assign(obsTimeMJD=lambda x:Time(x['obsTime'].to_list(),format='isot',scale='utc').utc.mjd)
-    )
+    obs_df = pd.read_xml(obs_data, dtype=ades_column_types)
+    obs_times = Time(obs_df['obsTime'].to_list(), format='isot', scale='utc')
+    obs_df['obsTimeMJD'] = obs_times.utc.mjd
+    obs_df['obsTimeMJDTDB'] = obs_times.tdb.mjd
+    if 'deprecated' in obs_df:
+        # drop rows with deprecated discovery observations
+        obs_df.query("deprecated != 'x' and deprecated != 'X'", inplace=True)
     # add columns if they are not present
     str_cols = ['trx', 'rcv', 'sys', 'selAst']
     for col in str_cols:
@@ -149,6 +151,8 @@ def create_optical_obs_df(body_id, optical_obs_file=None, t_min_tdb=None,
     for col in ades_column_types:
         if col not in obs_df:
             obs_df[col] = np.nan
+    # remove any columns that are not in ades_column_types
+    obs_df = obs_df[ades_column_types.keys()]
     if verbose:
         print(f"Read in {len(obs_df)} observations from the MPC.")
     _ades_mode_check(obs_df)
@@ -213,6 +217,7 @@ def add_psv_obs(psv_obs_file, obs_df, t_min_tdb=None, t_max_tdb=None, verbose=Fa
     psv_df['sigCorr'] = psv_df['rmsCorr']
     times = Time(psv_df['obsTime'].to_list(), format='isot', scale='utc')
     psv_df['obsTimeMJD'] = times.utc.mjd
+    psv_df['obsTimeMJDTDB'] = times.tdb.mjd
     add_counter = 0
     if verbose:
         print(f"Read in {len(psv_df)} observations from the ADES PSV file.")
@@ -324,6 +329,7 @@ def add_gaia_obs(obs_df, t_min_tdb=None, t_max_tdb=None, gaia_dr='gaiadr3', verb
         obs_df.loc[idx, 'provID'] = prov_id
         obs_df.loc[idx, 'obsTime'] = f'{obs_time.utc.iso}Z'
         obs_df.loc[idx, 'obsTimeMJD'] = data['epoch_utc'] + 55197.0
+        obs_df.loc[idx, 'obsTimeMJDTDB'] = obs_time.tdb.mjd
         obs_df.loc[idx, 'mode'] = 'CCD'
         obs_df.loc[idx, 'stn'] = '258'
         obs_df.loc[idx, 'ra'] = data['ra']
@@ -974,5 +980,4 @@ def get_optical_obs(body_id, optical_obs_file=None, t_min_tdb=None,
         obs_df = deweight_obs(obs_df, num_obs_per_night, verbose)
     elif eliminate:
         obs_df = eliminate_obs(obs_df, num_obs_per_night, verbose)
-    # obs_df.sort_values(by='obsTimeMJD', inplace=True, ignore_index=True)
     return obs_df
diff --git a/grss/fit/fit_radar.py b/grss/fit/fit_radar.py
@@ -93,6 +93,7 @@ def add_radar_obs(obs_df, t_min_tdb=None, t_max_tdb=None, verbose=False):
         obs_df.loc[idx,'provID'] = prov_id
         obs_df.loc[idx,'obsTime'] = f'{date.utc.isot}Z'
         obs_df.loc[idx,'obsTimeMJD'] = date.utc.mjd
+        obs_df.loc[idx,'obsTimeMJDTDB'] = date.tdb.mjd
         obs_df.loc[idx,'mode'] = 'RAD'
         obs_df.loc[idx,'trx'] = tx_code
         obs_df.loc[idx,'rcv'] = rx_code

diff --git a/grss/fit/fit_simulation.py b/grss/fit/fit_simulation.py
@@ -515,7 +515,6 @@ def __init__(self, x_init, obs_df, cov_init=None, n_iter_max=10,
         body_id = perm_id if isinstance(perm_id, str) else prov_id
         self.name = body_id
         self.t_sol = None
-        self.t_sol_utc = None
         self.x_init = None
         self.x_nom = None
         self.covariance_init = cov_init
@@ -610,7 +609,6 @@ def _check_initial_solution(self, x_init, cov_init):
             if key in x_init and x_init[key] != 0.0:
                 self.fit_nongrav = True
         self.t_sol = x_init['t']
-        self.t_sol_utc = Time(self.t_sol, format='mjd', scale='tdb').utc.mjd
         self.x_init = {key: x_init[key] for key in x_init if key != 't'}
         self.x_nom = self.x_init.copy()
         self.n_fit = len(self.x_nom)
@@ -663,6 +661,7 @@ def _add_simulated_obs(self):
             self.obs.loc[idx, 'provID'] = 'SIM_'+str(self.obs['provID'][0])
             self.obs.loc[idx, 'obsTime'] = f'{time.utc.isot}Z'
             self.obs.loc[idx, 'obsTimeMJD'] = time.utc.mjd
+            self.obs.loc[idx, 'obsTimeMJDTDB'] = time.tdb.mjd
             self.obs.loc[idx, 'mode'] = mode
             if mode in {'SIM_CCD', 'SIM_OCC'}:
                 self.obs.loc[idx, 'stn'] = 'SIM'
@@ -734,9 +733,9 @@ def _parse_observation_arrays(self, obs_df):
         self.radar_idx = self.delay_idx + self.doppler_idx
         self.n_obs = 2*len(self.optical_idx)+len(self.radar_idx)
 
-        self.past_obs_idx = self.obs.query('obsTimeMJD < @self.t_sol_utc').index
+        self.past_obs_idx = self.obs.query('obsTimeMJDTDB < @self.t_sol').index
         self.past_obs_exist = len(self.past_obs_idx) > 0
-        self.future_obs_idx = self.obs.query('obsTimeMJD >= @self.t_sol_utc').index
+        self.future_obs_idx = self.obs.query('obsTimeMJDTDB >= @self.t_sol').index
         self.future_obs_exist = len(self.future_obs_idx) > 0
         return None
 
@@ -1259,7 +1258,7 @@ def _get_analytic_partials(self, prop_sim_past, prop_sim_future):
             future_optical_partials = np.array(prop_sim_future.opticalPartials)
             future_radar_partials = np.array(prop_sim_future.radarPartials)
             future_light_time = np.array(prop_sim_future.lightTimeEval)
-        t_eval_tdb = Time(self.obs['obsTimeMJD'].values, format='mjd', scale='utc').tdb.mjd
+        t_eval_tdb = self.obs.obsTimeMJDTDB.values
         cos_dec = self.obs.cosDec.values
         for i, obs_info_len in enumerate(self.observer_info_lengths):
             if obs_info_len in {4, 7}:

diff --git a/grss/fit/fit_utils.py b/grss/fit/fit_utils.py
@@ -427,6 +427,17 @@ def get_sbdb_info(tdes):
             nongrav_params[nongrav_key_map[key]] = nongrav_default_dict[key]
         if key in cov_keys:
             elements[nongrav_key_map[key]] = nongrav_data_dict[key]
+    abs_mag = 0
+    albedo = 0.125
+    for par in raw_data['phys_par']:
+        if par['name'] == 'H':
+            abs_mag = float(par['value'])
+        if par['name'] == 'albedo':
+            albedo = float(par['value'])
+    body_radius = 1329*10**(-0.2*abs_mag)/(2*albedo**0.5)*1e3
+    if abs_mag == 0:
+        body_radius = 0
+    nongrav_params['radius'] = body_radius
     return [elements, cov_mat, nongrav_params]
 
 def get_similarity_stats(sol_1, cov_1, sol_2, cov_2):
@@ -464,7 +475,13 @@ def get_similarity_stats(sol_1, cov_1, sol_2, cov_2):
     det_2 = np.linalg.det(cov_2)
     det_big = np.linalg.det(big_cov)
     term_1 = 1/8 * (sol_1-sol_2) @ big_cov @ (sol_1-sol_2).T
-    term_2 = 1/2 * np.log(det_big/np.sqrt(det_1*det_2))
+    # simplify the natural log in term_2 because sometimes the determinant
+    # is too small and the product of the two determinants is beyond the lower
+    # limit of the float64 type
+    # term_2 = 1/2 * np.log(det_big/np.sqrt(det_1*det_2))
+    # term_2 = 1/2 * (np.log(det_big) - np.log(np.sqrt(det_1*det_2)))
+    # term_2 = 1/2 * (np.log(det_big) - 1/2 * np.log(det_1*det_2))
+    term_2 = 1/2 * np.log(det_big) - 1/4 * (np.log(det_1) + np.log(det_2))
     bhattacharya = term_1 + term_2
     bhatt_coeff = np.exp(-bhattacharya)
     return maha_dist_1, maha_dist_2, bhattacharya, bhatt_coeff
diff --git a/grss/version.txt b/grss/version.txt
@@ -1 +1 @@
-4.0.2
+4.0.3
diff --git a/include/interpolate.h b/include/interpolate.h
@@ -1,7 +1,7 @@
 #ifndef INTERPOLATE_H
 #define INTERPOLATE_H
 
-#include "force.h"
+#include "observe.h"
 
 /**
  * @brief Compute the sum of two real numbers with a compensated summation.
@@ -44,13 +44,6 @@ void approx_xInteg(const std::vector<real> &xInteg0,
  */
 void interpolate_on_the_fly(PropSimulation *propSim, const real &t, const real &dt);
 
-/**
- * @brief Interpolate the integrator state for one evaluation time.
- */
-void evaluate_one_interpolation(const PropSimulation *propSim, const real &t,
-                                const real &dt, const real &tInterp,
-                                std::vector<real> &xInterp);
-
 /**
  * @brief Determine whether the next interpolation index is within the window
  * of the time step that was just completed.
@@ -80,70 +73,4 @@ void get_lightTimeOneBody(PropSimulation *propSim, const size_t &i,
                           const bool bouncePointAtCenterOfMass, const real &t,
                           const real &dt, real &lightTimeOneBody);
 
-/**
- * @brief Compute the correction to the apparent state of the body due to the
- * gravitational light bending.
- */
-void get_glb_correction(PropSimulation *propSim, const real &tInterpGeom,
-                        std::vector<real> &xInterpApparentBary);
-
-/**
- * @brief Get the relevant measurement (optical/radar) for a given measurement time.
- */
-void get_measurement(PropSimulation *propSim, const size_t &interpIdx,
-                     const real &t, const real &dt, const real tInterpGeom,
-                     const std::vector<real> &xInterpGeom,
-                     const std::vector<real> &xInterpApparent);
-
-/**
- * @brief Get the optical measurement and partials.
- */
-void get_optical_measurement(PropSimulation *propSim,
-                             const std::vector<real> &xInterpApparent,
-                             std::vector<real> &opticalMeasurement,
-                             std::vector<real> &opticalPartials);
-
-/**
- * @brief Get the radar measurement and partials.
- */
-void get_radar_measurement(PropSimulation *propSim, const size_t &interpIdx,
-                           const real &t, const real &dt,
-                           const real tInterpGeom,
-                           const std::vector<real> &xInterpGeom,
-                           std::vector<real> &radarMeasurement,
-                           std::vector<real> &radarPartials);
-
-/**
- * @brief Get the radar delay measurement and partials.
- */
-void get_delay_measurement(PropSimulation *propSim, const size_t &interpIdx,
-                           const real &t, const real &dt, const size_t &i,
-                           const real tInterpGeom,
-                           const std::vector<real> &xInterpGeom,
-                           const real &receiveTimeTDB, real &transmitTimeTDB,
-                           std::vector<real> &xObsBaryRcv,
-                           std::vector<real> &xTrgtBaryBounce,
-                           std::vector<real> &xObsBaryTx, real &delayMeasurement,
-                           std::vector<real> &delayPartials);
-
-/**
- * @brief Get the relativistic delay measurement correction.
- */
-void get_delta_delay_relativistic(PropSimulation *propSim,
-                                  const real &tForSpice,
-                                  const std::vector<real> &targetState,
-                                  real &deltaDelayRelativistic);
-
-/**
- * @brief Get the Doppler measurement and partials.
- */
-void get_doppler_measurement(PropSimulation *propSim, const size_t &i,
-                             const real receiveTimeTDB,
-                             const real transmitTimeTDB,
-                             const std::vector<real> xObsBaryRcv,
-                             const std::vector<real> xTrgtBaryBounce,
-                             const std::vector<real> xObsBaryTx,
-                             const real transmitFreq, real &dopplerMeasurement,
-                             std::vector<real> &dopplerPartials);
-
 #endif
diff --git a/include/observe.h b/include/observe.h
@@ -0,0 +1,80 @@
+#ifndef OBSERVE_H
+#define OBSERVE_H
+
+#include "force.h"
+
+/**
+ * @brief Compute the correction to the apparent state of the body due to the
+ * gravitational light bending.
+ */
+void get_glb_correction(PropSimulation *propSim, const real &tInterpGeom,
+                        std::vector<real> &xInterpApparentBary);
+
+/**
+ * @brief Get the relevant measurement (optical/radar) for a given measurement time.
+ */
+void get_measurement(PropSimulation *propSim, const size_t &interpIdx,
+                     const real &t, const real &dt, const real tInterpGeom,
+                     const std::vector<real> &xInterpGeom,
+                     const std::vector<real> &xInterpApparent);
+
+/**
+ * @brief Get the optical measurement and partials.
+ */
+void get_optical_measurement(PropSimulation *propSim,
+                             const std::vector<real> &xInterpApparent,
+                             std::vector<real> &opticalMeasurement,
+                             std::vector<real> &opticalPartials);
+
+/**
+ * @brief Get the radar measurement and partials.
+ */
+void get_radar_measurement(PropSimulation *propSim, const size_t &interpIdx,
+                           const real &t, const real &dt,
+                           const real tInterpGeom,
+                           const std::vector<real> &xInterpGeom,
+                           std::vector<real> &radarMeasurement,
+                           std::vector<real> &radarPartials);
+
+/**
+ * @brief Get the radar delay measurement and partials.
+ */
+void get_delay_measurement(PropSimulation *propSim, const size_t &interpIdx,
+                           const real &t, const real &dt, const size_t &i,
+                           const real tInterpGeom,
+                           const std::vector<real> &xInterpGeom,
+                           const real &receiveTimeTDB, real &transmitTimeTDB,
+                           std::vector<real> &xObsBaryRcv,
+                           std::vector<real> &xTrgtBaryBounce,
+                           std::vector<real> &xObsBaryTx, real &delayMeasurement,
+                           std::vector<real> &delayPartials);
+
+/**
+ * @brief Get the relativistic delay measurement correction.
+ */
+void get_delta_delay_relativistic(PropSimulation *propSim,
+                                  const real &tForSpice,
+                                  const std::vector<real> &targetState,
+                                  real &deltaDelayRelativistic);
+
+/**
+ * @brief Get the Doppler measurement and partials.
+ */
+void get_doppler_measurement(PropSimulation *propSim, const size_t &i,
+                             const real receiveTimeTDB,
+                             const real transmitTimeTDB,
+                             const std::vector<real> xObsBaryRcv,
+                             const std::vector<real> xTrgtBaryBounce,
+                             const std::vector<real> xObsBaryTx,
+                             const real transmitFreq, real &dopplerMeasurement,
+                             std::vector<real> &dopplerPartials);
+
+/**
+ * @brief Interpolate the integrator state for one evaluation time.
+ */
+void evaluate_one_interpolation(
+    const PropSimulation *propSim, const real &t, const real &dt,
+    const real &tInterp,
+    std::vector<real> &xInterp);  // defined in interpolate.cpp
+
+#endif
diff --git a/pyproject.toml b/pyproject.toml
@@ -3,7 +3,7 @@ name = "grss"
 dynamic = ["version", "dependencies"]
 description = "GRSS: Gauss-Radau Small-body Simulator"
 readme = "README.md"
-requires-python = ">=3.8"
+requires-python = ">=3.10"
 license = {file = "LICENSE"}
 keywords = [
     "astronomy",
-Original file line number
+Diff line change
@@ Expand Up / @@ -5,6 +5,7 @@ on: @@
           - main
           - dev
         paths:
+          - '.github/workflows/joss.yml'
           - 'joss/**'
     jobs:
@@ Expand Down @@