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

@@ -57,7 +57,7 @@ class Instrument:
 
     def __init__(self,
                  # Sampling parameters
-                 size=2592000, dt=1/4, t0='orbits',
+                 size=2592000, dt=1/4, t0='orbits', telemetry_dt=86400,
                  # Physics simulation sampling and filtering
                  physics_upsampling=4, aafilter=('kaiser', 240, 1.1, 2.9),
                  # Inter-spacecraft propagation
@@ -79,6 +79,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',
@@ -134,6 +136,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
@@ -170,6 +173,11 @@ class Instrument:
         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
 
+        self.telemetry_dt = telemetry_dt
+        self.telemetry_fs = 1 / telemetry_dt
+        self.telemetry_size = int(numpy.ceil(size / telemetry_dt * self.dt))
+        self.telemetry_t = self.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 +243,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)
@@ -465,6 +476,7 @@ class Instrument:
 
     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 +485,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 +495,18 @@ class Instrument:
                 values = orbitf[f'tps/l_{mosa}']['d_ppr']
             return scipy.interpolate.InterpolatedUnivariateSpline(times, values, k=5, ext='raise')
 
+        def tcb_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.tcb_proper_time_deviations = ForEachSC(lambda sc: tcb_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 +644,6 @@ class Instrument:
 
         self.simulate_noises()
 
-        ## TDIR tone
-
         logger.debug("Computing local timer deviations")
         self.local_timer_deviations = \
             ForEachSC(lambda sc:
@@ -635,6 +653,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.tcb_proper_time_deviations + \
+            self.tcb_sync_noise
+
+        ## 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 +1228,12 @@ class Instrument:
             noises.dws(self.physics_fs, self.physics_size, self.dws_asds[mosa])
         )
 
+        ## TCB synchronization noise
+
+        self.tcb_sync_noise = ForEachSC(lambda sc:
+            noises.white(self.telemetry_fs, self.telemetry_size, self.sync_asds[sc])
+        )
+
         ## Angular jitters
 
         logger.info("Generating spacecraft angular jitters")
@@ -1723,6 +1759,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()