split up observation and measurement routines; added obsTimeMJDTDB co…

…lumn to df

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
 

grss/fit/fit_optical.py

@@ -136,11 +136,10 @@ 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)
@@ -218,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.")
@@ -329,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']
@@ -979,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

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

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}:

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

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

src/CMakeLists.txt

@@ -3,6 +3,7 @@ approach.cpp
 elements.cpp
 force.cpp
 gr15.cpp
+observe.cpp
 interpolate.cpp
 parallel.cpp
 pck.cpp