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
@@ -78,6 +78,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 ASD:
+                 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',
@@ -165,6 +167,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(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
@@ -230,6 +237,9 @@ class Instrument:
         else:
             self.clock_freqquaddrifts = ForEachSC(clock_freqquaddrifts)
 
+        # synchronization
+        self.sync_asds = ForEachSC(sync_asds)
+
         # Clock-noise inversion
         self.clockinv_tolerance = float(clockinv_tolerance)
         self.clockinv_maxiter = int(clockinv_maxiter)
@@ -453,6 +463,10 @@ class Instrument:
                     self.d_pprs = ForEachMOSA(lambda mosa: scipy.interpolate.InterpolatedUnivariateSpline(
                         orbitf['tps']['tau'][:], orbitf[f'tps/l_{mosa}']['d_ppr'], k=5, ext='raise')(self.physics_t)
                     )
+                    logger.debug("Interpolating proper time deviation from TCB")
+                    self.tcb_proper_time_deviations = ForEachSC(lambda sc: scipy.interpolate.InterpolatedUnivariateSpline(
+                        orbitf['tcb']['t'][:], orbitf[f'tcb/sc_{sc}']['tau'], k=5, ext='raise')(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
@@ -600,8 +614,6 @@ class Instrument:
 
         self.simulate_noises()
 
-        ## TDIR tone
-
         logger.debug("Computing local timer deviations")
         self.local_timer_deviations = \
             ForEachSC(lambda sc:
@@ -611,6 +623,17 @@ 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 + \
+            1 / 2 * self.clock_freqlindrifts * self.telemetry_t**2 + \
+            1 / 3 * self.clock_freqquaddrifts * self.telemetry_t**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(
@@ -1075,6 +1098,8 @@ class Instrument:
         # Laser noise are only generated for lasers locked on cavities,
         # in `simulate_locking.lock_on_cavity()`
 
+        # simulate synchronization noise! ForEachSC
+
         self.laser_noises = ForEachMOSA(0)
 
         ## Clock noise
@@ -1174,6 +1199,11 @@ 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")
@@ -1699,6 +1729,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()