Resolve "Error when using orbit files v2.0"

Closes #82 (closed)

• Fix the issue with orbit files v2.0
• Avoid printing multiple warnings for orbit files in a development version
• Add unit tests for various orbit files

lisainstrument/instrument.py

@@ -180,9 +180,6 @@ class Instrument:
                 with File(orbits, 'r') as orbitf:
                     version = Version(orbitf.attrs['version'])
                     logger.debug("Using orbit file version %s", version)
-                    # Warn for orbit file development version
-                    if version.is_devrelease:
-                        logger.warning("You are using an orbit file in a development version")
                     # Switch between various orbit file standards
                     if version in SpecifierSet('== 1.*', True):
                         self.t0 = float(orbitf.attrs['t0' if orbit_dataset == 'tcb/ltt' else 'tau0'])
@@ -208,9 +205,6 @@ class Instrument:
                 with File(orbits, 'r') as orbitf:
                     version = Version(orbitf.attrs['version'])
                     logger.debug("Using orbit file version %s", version)
-                    # Warn for orbit file development version
-                    if version.is_devrelease:
-                        logger.warning("You are using an orbit file in a development version")
                     # Switch between various orbit file standards
                     if version in SpecifierSet('== 1.*', True):
                         self.telemetry_t0 = float(orbitf.attrs['t0'])
@@ -593,23 +587,24 @@ class Instrument:
         # Prepare common interpolation method
         times = orbitf.attrs['t0'] + np.arange(orbitf.attrs['size']) * orbitf.attrs['dt']
         interpolate = lambda data, t: InterpolatedUnivariateSpline(times, data, ext='raise')(t)
-        index = {'12': 0, '23': 1, '31': 2, '13': 3, '32': 4, '21': 5}
+        link_index = {'12': 0, '23': 1, '31': 2, '13': 3, '32': 4, '21': 5}
+        sc_index = {'1': 0, '2': 1, '3': 2}
 
         # Interpolate necessary orbital quantities,
         # show a helpful error message if orbit file is too short
         try:
             logger.debug("Interpolating proper pseudo-ranges")
             dataset = orbitf['tcb/ltt'] if self.orbit_dataset == 'tcb/ltt' else orbitf['tps/ppr']
-            self.pprs = ForEachMOSA(lambda mosa: interpolate(dataset[:, index[mosa]], self.physics_t))
+            self.pprs = ForEachMOSA(lambda mosa: interpolate(dataset[:, link_index[mosa]], self.physics_t))
             logger.debug("Interpolating proper pseudo-range derivatives")
             dataset = orbitf['tcb/d_ltt'] if self.orbit_dataset == 'tcb/ltt' else orbitf['tps/d_ppr']
-            self.d_pprs = ForEachMOSA(lambda mosa: interpolate(dataset[:, index[mosa]], self.physics_t))
+            self.d_pprs = ForEachMOSA(lambda mosa: interpolate(dataset[:, link_index[mosa]], self.physics_t))
             logger.debug("Interpolating proper time deviation from TCB")
             if self.orbit_dataset == 'tcb/ltt':
                 self.tps_proper_time_deviations = ForEachSC(lambda sc: 0)
             else:
                 dataset = orbitf['tcb/delta_tau']
-                self.tps_proper_time_deviations = ForEachSC(lambda sc: interpolate(dataset[:, sc - 1], self.physics_t))
+                self.tps_proper_time_deviations = ForEachSC(lambda sc: interpolate(dataset[:, sc_index[sc]], self.physics_t))
         except ValueError as error:
             logger.error("Missing orbit information at \n%s", self.physics_t)
             raise ValueError("missing orbit information, use longer orbit file or adjust sampling") from error
@@ -1552,9 +1547,11 @@ class Instrument:
         logger.debug("Solving iteratively (tolerance=%s s, maxiter=%s)",
             self.clockinv_tolerance, self.clockinv_maxiter)
 
-        niter = 0
+        # Drop samples at the edges to compute error
+        edge = min(100, len(timer_deviations) // 2 - 1)
         error = 0
-        edge = 100 # drop samples at the edges to check convergence
+
+        niter = 0
         next_inverse = timer_deviations
         while not niter or error > self.clockinv_tolerance:
             if niter >= self.clockinv_maxiter: