Resolve "Simulate the estimate for the SCET offset from TCB"

Closes #26 (closed)

Simulate the timer deviations from TCB to THE.

Changes:

• New simulation parameters telemetry_downsampling, telemetry_t0 and sync_asds
• New attributes tcb_timer_deviations and tcb_sync_noise
• New HDF5 output datasets tcb_timer_deviations and tcb_sync_noise

lisainstrument/instrument.py

@@ -60,6 +60,8 @@ class Instrument:
                  size=2592000, dt=1/4, t0='orbits',
                  # Physics simulation sampling and filtering
                  physics_upsampling=4, aafilter=('kaiser', 240, 1.1, 2.9),
+                 # Telemetry sampling
+                 telemetry_downsampling=86400 * 4, telemetry_t0='orbits',
                  # Inter-spacecraft propagation
                  orbits='static', neglect_tps=False,
                  gws=None, interpolation=('lagrange', 31),
@@ -79,6 +81,8 @@ class Instrument:
                  # Optical pathlength noises
                  backlink_asds=3E-12, backlink_fknees=2E-3, testmass_asds=2.4E-15, testmass_fknees=0.4E-3,
                  oms_asds=(6.35E-12, 1.25E-11, 1.42E-12, 3.38E-12, 3.32E-12, 7.90E-12), oms_fknees=2E-3,
+                 # TCB synchronization
+                 sync_asds=2.2E-3,
                  # Tilt-to-length (TTL)
                  ttl_coeffs='default',
                  sc_jitter_asds=(1E-8, 1E-8, 1E-8), mosa_jitter_asds=1E-8, mosa_angles='default',
@@ -96,6 +100,8 @@ class Instrument:
                 or None for no filter; to design a filter from a Kaiser window, use a tuple
                 ('kaiser', attenuation [dB], f1 [Hz], f2 [Hz]) with f1 < f2 the frequencies defining
                 the transition band
+            telemetry_downsampling: ratio of sampling frequencies for measurements vs. telemetry events
+            telemetry_t0: time of initial telemetry event [s], or 'orbits' to match that of the orbits
             orbits: path to orbit file, dictionary of constant PPRs for static arms, 'static'
                 for a set of static PPRs corresponding to a fit of Keplerian orbits around t = 0,
                 or dictionary of PPR time series
@@ -134,6 +140,7 @@ class Instrument:
             testmass_fknees: dictionary of cutoff frequencies for test-mass noise [Hz]
             oms_asds: tuple of dictionaries of amplitude spectral densities for OMS noise [m/sqrt(Hz)],
                 ordered as (isc_carrier, isc_usb, tm_carrier, tm_usb, ref_carrier, ref_usb)
+            sync_asds: dictionary of amplitude spectral densities for TCB synchronization noise [s/sqrt(Hz)]
             oms_fknees: dictionary of cutoff frequencies for OMS noise
             ttl_coeffs: tuple (local_phi, distant_phi, local_eta, distant_eta) of dictionaries of
                 tilt-to-length coefficients on each MOSA [m/rad], 'default' for a default set of
@@ -154,22 +161,40 @@ class Instrument:
         self.simulated = False
 
         # Measurement sampling
-        self.size = int(size)
-        self.dt = float(dt)
         if t0 == 'orbits':
             if isinstance(orbits, str) and orbits != 'static':
                 logger.debug("Reading initial time from orbit file '%s'", orbits)
                 with h5py.File(orbits, 'r') as orbitf:
-                    self.t0 = float(orbitf.attrs['tau0'])
+                    attr = 't0' if neglect_tps else 'tau0'
+                    self.t0 = float(orbitf.attrs[attr])
             else:
                 self.t0 = 0.0
         else:
             self.t0 = float(t0)
+        self.size = int(size)
+        self.dt = float(dt)
         self.fs = 1 / self.dt
         self.duration = self.dt * self.size
         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
 
+        # Sampling of telemetry events given in TCB
+        if telemetry_t0 == 'orbits':
+            if isinstance(orbits, str) and orbits != 'static':
+                logger.debug("Reading telemetry initial time from orbit file '%s'", orbits)
+                with h5py.File(orbits, 'r') as orbitf:
+                    self.telemetry_t0 = float(orbitf.attrs['t0'])
+            else:
+                self.telemetry_t0 = 0.0
+        else:
+            self.telemetry_t0 = float(telemetry_t0)
+        self.telemetry_downsampling = int(telemetry_downsampling)
+        self.telemetry_dt = self.dt * self.telemetry_downsampling
+        self.telemetry_fs = 1 / self.telemetry_dt
+        self.telemetry_size = self.size // self.telemetry_downsampling
+        self.telemetry_t = \
+            self.telemetry_t0 + numpy.arange(self.telemetry_size, dtype=numpy.float64) * self.telemetry_dt
+
         # Physics sampling
         self.physics_upsampling = int(physics_upsampling)
         self.physics_size = self.size * self.physics_upsampling
@@ -235,6 +260,9 @@ class Instrument:
         else:
             self.clock_freqquaddrifts = ForEachSC(clock_freqquaddrifts)
 
+        # TCB synchronization
+        self.sync_asds = ForEachSC(sync_asds)
+
         # Clock-noise inversion
         self.clockinv_tolerance = float(clockinv_tolerance)
         self.clockinv_maxiter = int(clockinv_maxiter)
@@ -433,14 +461,17 @@ class Instrument:
         """Initialize orbits."""
         if orbits == 'static':
             logger.info("Using default set of static proper pseudo-ranges")
+            if not neglect_tps:
+                logger.warning("Neglecting TPS despite configuration because static PPRs are used")
             self.orbit_file = None
-            self.neglect_tps = None
             self.pprs = ForEachMOSA({
                 # Default PPRs based on first samples of Keplerian orbits (v1.0)
                 '12': 8.3324, '23': 8.3028, '31': 8.3324,
                 '13': 8.3315, '32': 8.3044, '21': 8.3315,
             })
             self.d_pprs = ForEachMOSA(0)
+            self.neglect_tps = True
+            self.tps_proper_time_deviations = ForEachSC(0)
         elif isinstance(orbits, str):
             logger.info("Using orbit file '%s'", orbits)
             self.orbit_file = orbits
@@ -456,15 +487,19 @@ class Instrument:
                     raise ValueError(f"unsupported orbit file version '{orbits_version}'")
         else:
             logger.info("Using user-provided proper pseudo-ranges and derivatives thereof")
+            if not self.neglect_tps:
+                logger.warning("Neglecting TPS despite configuration because user-provided PPRs are used")
             self.orbit_file = None
-            self.neglect_tps = None
+            self.neglect_tps = True
             self.pprs = ForEachMOSA(orbits)
             self.d_pprs = self.pprs.transformed(lambda sc, x:
                 0 if numpy.isscalar(x) else numpy.gradient(x, self.physics_dt)
             )
+            self.tps_proper_time_deviations = ForEachSC(0)
 
     def init_orbits_file_1_0(self, orbitf):
         """Initialize orbits from an orbit file version 1.0, 1.0.1 or 1.0.2."""
+
         def pprs(mosa):
             if self.neglect_tps:
                 times = orbitf['tcb']['t'][:]
@@ -473,6 +508,7 @@ class Instrument:
                 times = orbitf['tps']['tau'][:]
                 values = orbitf[f'tps/l_{mosa}']['ppr']
             return scipy.interpolate.InterpolatedUnivariateSpline(times, values, k=5, ext='raise')
+
         def d_pprs(mosa):
             if self.neglect_tps:
                 times = orbitf['tcb']['t'][:]
@@ -482,11 +518,18 @@ class Instrument:
                 values = orbitf[f'tps/l_{mosa}']['d_ppr']
             return scipy.interpolate.InterpolatedUnivariateSpline(times, values, k=5, ext='raise')
 
+        def tps_proper_time_deviations(sc):
+            times = orbitf['tcb']['t'][:]
+            values = orbitf[f'tcb/sc_{sc}']['tau']
+            return scipy.interpolate.InterpolatedUnivariateSpline(times, values, k=5, ext='raise')
+
         try:
             logger.debug("Interpolating proper pseudo-ranges")
             self.pprs = ForEachMOSA(lambda mosa: pprs(mosa)(self.physics_t))
             logger.debug("Interpolating proper pseudo-range derivatives")
             self.d_pprs = ForEachMOSA(lambda mosa: d_pprs(mosa)(self.physics_t))
+            logger.debug("Interpolating proper time deviation from TCB")
+            self.tps_proper_time_deviations = ForEachSC(lambda sc: tps_proper_time_deviations(sc)(self.telemetry_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
@@ -624,8 +667,6 @@ class Instrument:
 
         self.simulate_noises()
 
-        ## TDIR tone
-
         logger.debug("Computing local timer deviations")
         self.local_timer_deviations = \
             ForEachSC(lambda sc:
@@ -635,6 +676,18 @@ class Instrument:
                 dx=self.physics_dt, initial=self.clock_offsets[sc])
             )
 
+        ## Timer deviations from TCB
+
+        self.tcb_timer_deviations = \
+            self.clock_offsets + \
+            self.clock_freqoffsets * self.telemetry_t + \
+            self.clock_freqlindrifts * self.telemetry_t**2 / 2 + \
+            self.clock_freqquaddrifts * self.telemetry_t**3 / 3 + \
+            self.tps_proper_time_deviations + \
+            self.tcb_sync_noises
+
+        ## TDIR tone
+
         self.tdir_tones = ForEachMOSA(lambda mosa:
             0 if self.tdir_tone_amplitudes[mosa] == 0 \
             else self.tdir_tone_amplitudes[mosa] * numpy.sin(
@@ -1198,6 +1251,12 @@ class Instrument:
             noises.dws(self.physics_fs, self.physics_size, self.dws_asds[mosa])
         )
 
+        ## TCB synchronization noise
+
+        self.tcb_sync_noises = ForEachSC(lambda sc:
+            noises.white(self.telemetry_fs, self.telemetry_size, self.sync_asds[sc])
+        )
+
         ## Angular jitters
 
         logger.info("Generating spacecraft angular jitters")
@@ -1480,6 +1539,7 @@ class Instrument:
 
         self.write_metadata(hdf5)
         hdf5['physics_t'] = self.physics_t
+        hdf5['telemetry_t'] = self.telemetry_t
 
         if keep_all:
 
@@ -1524,6 +1584,9 @@ class Instrument:
             self.oms_ref_carrier_noises.write(hdf5, 'oms_ref_carrier_noises')
             self.oms_ref_usb_noises.write(hdf5, 'oms_ref_usb_noises')
 
+            logger.debug("Writing TCB synchronization noise to '%s'", output)
+            self.tcb_sync_noises.write(hdf5, 'tcb_sync_noises')
+
             logger.debug("Writing spacecraft angular jitter to '%s'", output)
             self.sc_jitter_phis.write(hdf5, 'sc_jitter_phis')
             self.sc_jitter_etas.write(hdf5, 'sc_jitter_etas')
@@ -1723,6 +1786,9 @@ class Instrument:
         self.ref_usb_fluctuations.write(hdf5, 'ref_usb_fluctuations')
         self.ref_usbs.write(hdf5, 'ref_usbs')
 
+        logger.debug("Write timer deviations from TCB to '%s'", output)
+        self.tcb_timer_deviations.write(hdf5, 'tcb_timer_deviations')
+
         logger.info("Closing measurement file '%s'", output)
         hdf5.close()