Read orbits from orbit file by default

lisainstrument/instrument.py

@@ -34,7 +34,7 @@ class Instrument:
     SCS = ForEachSC.indices()
     MOSAS = ForEachMOSA.indices()
 
-    def __init__(self, size=2592000, dt=1/4, t0=0,
+    def __init__(self, size=2592000, dt=1/4, t0='orbits',
                  # Inter-spacecraft propagation
                  orbits=None, gws=None, glitches=None, interpolation=('lagrange', 31),
                  # Lasers
@@ -54,7 +54,7 @@ class Instrument:
         Args:
             size: number of samples to generate
             dt: sampling period [s]
-            t0: initial time [s]
+            t0: initial time [s], or 'orbits' to match that of the orbits
             orbits: path to orbit file, or dictionary of scalars or time series for PPRs,
                 or None for default set of PPRs corresponding to static arms fit from Keplerian orbits
                 (from LISA Orbits v1.0) around t = 0
@@ -89,7 +89,7 @@ class Instrument:
                 ('kaiser', attenuation [dB], f1 [Hz], f2 [Hz]) with f1 < f2 the frequencies defining
                 the transition band
         """
-        # pylint: disable=too-many-arguments,too-many-statements
+        # pylint: disable=too-many-arguments,too-many-statements,too-many-locals,too-many-branches
         logger.info("Initializing instrumental simulation")
         self.git_url = 'https://gitlab.in2p3.fr/lisa-simulation/instrument'
         self.version = meta.__version__
@@ -98,7 +98,15 @@ class Instrument:
         # Measurement sampling
         self.size = int(size)
         self.dt = float(dt)
-        self.t0 = float(t0)
+        if t0 == 'orbits':
+            if isinstance(orbits, str):
+                logger.debug("Reading initial time from orbit file '%s'", orbits)
+                with h5py.File(orbits, 'r') as orbitf:
+                    self.t0 = float(orbitf.attrs['tau0'])
+            else:
+                self.t0 = 0.0
+        else:
+            self.t0 = float(t0)
         self.fs = 1 / self.dt
         logger.info("Computing measurement time vector (size=%s, dt=%s)", self.size, self.dt)
         self.t = self.t0 + numpy.arange(self.size, dtype=numpy.float64) * self.dt