diff --git a/lisainstrument/instrument.py b/lisainstrument/instrument.py
index 672d04bd531c336633661b3868384f343f3b5444..143d3c2087b2dcea3176a6ec380c368cb713bc6d 100755
--- a/lisainstrument/instrument.py
+++ b/lisainstrument/instrument.py
@@ -231,22 +231,28 @@ class Instrument:
# Telemetry sampling
self.telemetry_downsampling = int(telemetry_downsampling)
- self.initial_telemetry_size = int(initial_telemetry_size)
- self.initial_telemetry_measurement_size = \
- self.initial_telemetry_size * self.telemetry_downsampling
- self.initial_telemetry_physics_size = \
- self.initial_telemetry_measurement_size * self.physics_upsampling
- self.telemetry_size = self.initial_telemetry_size \
- + math.ceil(self.size / self.telemetry_downsampling) + 1
self.telemetry_dt = self.dt * self.telemetry_downsampling
self.telemetry_fs = self.fs / self.telemetry_downsampling
+ # Extra telemetry samples before t0
+ self.initial_telemetry_size = int(initial_telemetry_size)
self.telemetry_t0 = self.t0 - self.initial_telemetry_size * self.telemetry_dt
+ # Total telemetry size, includes initial telemetry samples
+ # plus telemetry samples covering the entire measurement time vector,
+ # hence the use of ``math.ceil`` -- the +1 is for the sample at t0
+ self.telemetry_size = self.initial_telemetry_size + 1 \
+ + math.ceil(self.size / self.telemetry_downsampling)
logger.info("Computing telemetry time vector (size=%s, dt=%s)", self.telemetry_size, self.telemetry_dt)
- self.telemetry_t = self.telemetry_t0 + np.arange(self.telemetry_size, dtype=np.float64) * self.telemetry_dt
-
- self.physics_et_withinitial = self.telemetry_t0 - self.t0 + \
- np.arange(self.physics_size + self.initial_telemetry_physics_size, dtype=np.float64) * self.physics_dt
- self.physics_t_withinitial = self.t0 + self.physics_et_withinitial
+ self.telemetry_t = self.telemetry_t0 \
+ + np.arange(self.telemetry_size, dtype=np.float64) * self.telemetry_dt
+ # Physics time vector covering telemetry samples
+ self.telemetry_to_physics_dt = self.telemetry_downsampling * self.physics_upsampling
+ self.physics_size_covering_telemetry = self.telemetry_size * self.telemetry_to_physics_dt
+ self.physics_et_covering_telemetry = self.telemetry_t0 - self.t0 + \
+ np.arange(self.physics_size_covering_telemetry, dtype=np.float64) * self.physics_dt
+ self.physics_t_covering_telemetry = self.t0 + self.physics_et_covering_telemetry
+ # How to cut physics data covering telemetry to regular size
+ first_sample = self.initial_telemetry_size * self.telemetry_to_physics_dt
+ self.telemetry_physics_slice = slice(first_sample, first_sample + self.physics_size)
# Orbits, gravitational waves, glitches
self.init_orbits(orbits, orbit_dataset)
@@ -654,7 +660,7 @@ class Instrument:
logger.debug("Interpolating proper pseudo-range derivatives")
self.d_pprs = ForEachMOSA(lambda mosa: pprs(mosa).derivative()(self.physics_t))
logger.debug("Interpolating TPSs with respect to TCB")
- self.tps_wrt_tcb = ForEachSC(lambda sc: tps_wrt_tcb(sc)(self.physics_t_withinitial))
+ self.tps_wrt_tcb = ForEachSC(lambda sc: tps_wrt_tcb(sc)(self.physics_t_covering_telemetry))
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
@@ -687,9 +693,10 @@ class Instrument:
self.tps_wrt_tcb = ForEachSC(lambda sc: 0)
else:
dataset = orbitf['tcb/delta_tau']
- self.tps_wrt_tcb = ForEachSC(lambda sc: interpolate(dataset[:, sc_index[sc]], self.physics_t_withinitial))
+ self.tps_wrt_tcb = ForEachSC(lambda sc:
+ interpolate(dataset[:, sc_index[sc]], self.physics_t_covering_telemetry))
except ValueError as error:
- logger.error("Missing orbit information at \n%s", self.physics_t_withinitial)
+ logger.error("Missing orbit information at \n%s", self.physics_t_covering_telemetry)
raise ValueError("missing orbit information, use longer orbit file or adjust sampling") from error
def init_gws(self, gws):
@@ -962,18 +969,18 @@ class Instrument:
+ self.integrated_clock_noise_fluctuations
logger.debug("Computing THE with respect to TCB")
- t = self.physics_et_withinitial
+ t = self.physics_et_covering_telemetry
self.the_wrt_tcb_withinitial = \
self.tps_wrt_tcb \
+ self.clock_offsets \
+ self.clock_freqoffsets * (t + self.tps_wrt_tcb) \
+ self.clock_freqlindrifts * (t + self.tps_wrt_tcb)**2 / 2 \
+ self.clock_freqquaddrifts * (t + self.tps_wrt_tcb)**3 / 3 \
- + self.tps_wrt_tcb * self.clock_noise_fluctuations_withinitial \
- + self.integrated_clock_noise_fluctuations_withinitial
+ + self.tps_wrt_tcb * self.clock_noise_fluctuations_covering_telemetry \
+ + self.integrated_clock_noise_fluctuations_covering_telemetry
logger.debug("Computing MOC time correlations")
- physics_to_telemetry = lambda _, x: x[::self.physics_upsampling * self.telemetry_downsampling]
+ physics_to_telemetry = lambda _, x: x[::self.telemetry_to_physics_dt]
self.moc_time_correlations = self.moc_time_correlation_noises \
+ self.the_wrt_tcb_withinitial.transformed(physics_to_telemetry)
@@ -1521,35 +1528,36 @@ class Instrument:
logger.debug("Generating clock noise fluctuations")
# Include initial telemetry time period
- self.clock_noise_fluctuations_withinitial = ForEachSC(lambda sc:
+ self.clock_noise_fluctuations_covering_telemetry = ForEachSC(lambda sc:
noises.clock(
self.physics_fs,
- self.physics_size + self.initial_telemetry_physics_size,
+ self.physics_size_covering_telemetry,
self.clock_asds[sc]),
concurrent=self.concurrent
)
# Slice to only select physics period
- self.clock_noise_fluctuations = self.clock_noise_fluctuations_withinitial.transformed(
- lambda _, x: x if np.isscalar(x) else x[self.initial_telemetry_physics_size:]
- )
+ self.clock_noise_fluctuations = \
+ self.clock_noise_fluctuations_covering_telemetry.transformed(
+ lambda _, x: x if np.isscalar(x) else x[self.telemetry_physics_slice]
+ )
logger.debug("Integrating clock noise fluctuations")
# Include initial telemetry time period
- self.integrated_clock_noise_fluctuations_withinitial = \
+ self.integrated_clock_noise_fluctuations_covering_telemetry = \
ForEachSC(lambda sc:
cumulative_trapezoid(np.broadcast_to(
- self.clock_noise_fluctuations_withinitial[sc],
- self.physics_size + self.initial_telemetry_physics_size),
+ self.clock_noise_fluctuations_covering_telemetry[sc],
+ self.physics_size_covering_telemetry),
dx=self.physics_dt, initial=0),
concurrent=self.concurrent
)
# Slice to only select physics period
self.integrated_clock_noise_fluctuations = \
- self.integrated_clock_noise_fluctuations_withinitial.transformed(
- lambda _, x: x if np.isscalar(x) else x[self.initial_telemetry_physics_size:]
+ self.integrated_clock_noise_fluctuations_covering_telemetry.transformed(
+ lambda _, x: x if np.isscalar(x) else x[self.telemetry_physics_slice]
)
## Modulation noise
@@ -1986,9 +1994,9 @@ class Instrument:
logger.debug("Writing clock noise to '%s'", output)
self.clock_noise_offsets.write(hdf5, 'clock_noise_offsets')
self.clock_noise_fluctuations.write(hdf5, 'clock_noise_fluctuations')
- self.clock_noise_fluctuations_withinitial.write(
+ self.clock_noise_fluctuations_covering_telemetry.write(
hdf5, 'clock_noise_fluctuations_withinitial')
- self.integrated_clock_noise_fluctuations_withinitial.write(
+ self.integrated_clock_noise_fluctuations_covering_telemetry.write(
hdf5, 'integrated_clock_noise_fluctuations_withinitial')
self.integrated_clock_noise_fluctuations.write(
hdf5, 'integrated_clock_noise_fluctuations')