diff --git a/lisainstrument/instrument.py b/lisainstrument/instrument.py
index 890098179bb756c06a8cb086e192a9e11b3a5452..76d077fa044c96aa80627dee3e8dc3a4e1796f16 100755
--- a/lisainstrument/instrument.py
+++ b/lisainstrument/instrument.py
@@ -67,7 +67,7 @@ class Instrument:
# Physics simulation sampling and filtering
physics_upsampling=4, aafilter=('kaiser', 240, 1.1, 2.9),
# Telemetry sampling
- telemetry_downsampling=86400 * 4, telemetry_t0='orbits',
+ telemetry_downsampling=86400 * 4, initial_telemetry_size=0,
# Inter-spacecraft propagation
orbits='static', orbit_dataset='tps/ppr',
gws=None, interpolation=('lagrange', 31),
@@ -100,7 +100,7 @@ class Instrument:
"""Initialize an instrumental simulation.
Args:
- size: number of samples to generate
+ size: number of measurement samples to generate
dt: sampling period [s]
t0: initial time [s], or 'orbits' to match that of the orbits
physics_upsampling: ratio of sampling frequencies for physics vs. measurement simulation
@@ -109,7 +109,7 @@ class Instrument:
('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
+ initial_telemetry_size: number of telemetry samples before :attr:`lisainstrument.Instrument.t0`
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
@@ -201,32 +201,6 @@ class Instrument:
logger.info("Computing measurement time vector (size=%s, dt=%s)", self.size, self.dt)
self.t = self.t0 + np.arange(self.size, dtype=np.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 File(orbits, 'r') as orbitf:
- version = Version(orbitf.attrs['version'])
- logger.debug("Using orbit file version %s", version)
- # Switch between various orbit file standards
- if version in SpecifierSet('== 1.*', True):
- self.telemetry_t0 = float(orbitf.attrs['t0'])
- elif version in SpecifierSet('== 2.*', True):
- self.telemetry_t0 = float(orbitf.attrs['t0'])
- else:
- raise ValueError(f"unsupported orbit file version '{version}'")
-
- 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 = int(np.ceil(self.size / self.telemetry_downsampling))
- self.telemetry_t = \
- self.telemetry_t0 + np.arange(self.telemetry_size, dtype=np.float64) * self.telemetry_dt
-
# Physics sampling
self.physics_upsampling = int(physics_upsampling)
self.physics_size = self.size * self.physics_upsampling
@@ -235,6 +209,23 @@ class Instrument:
logger.info("Computing physics time vector (size=%s, dt=%s)", self.physics_size, self.physics_dt)
self.physics_t = self.t0 + np.arange(self.physics_size, dtype=np.float64) * self.physics_dt
+ # 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 \
+ + int(np.ceil(self.size / self.telemetry_downsampling))
+ self.telemetry_dt = self.dt * self.telemetry_downsampling
+ self.telemetry_fs = self.fs / self.telemetry_downsampling
+ self.telemetry_t0 = self.t0 - self.initial_telemetry_size * self.telemetry_dt
+ 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_t_withinitial = self.telemetry_t0 + \
+ np.arange(self.physics_size + self.initial_telemetry_physics_size, dtype=np.float64) * self.physics_dt
+
# Instrument topology
self.central_freq = float(central_freq)
self.init_lock(lock)
@@ -811,19 +802,34 @@ class Instrument:
logger.debug("Computing local timer deviations")
self.local_timer_deviations = \
+ self.clock_offsets + \
+ self.clock_freqoffsets * self.physics_t + \
+ self.clock_freqlindrifts * self.physics_t**2 / 2 + \
+ self.clock_freqquaddrifts * self.physics_t**3 / 3 + \
ForEachSC(lambda sc:
cumulative_trapezoid(np.broadcast_to(
- self.clock_noise_offsets[sc] + self.clock_noise_fluctuations[sc],
- self.physics_size),
- dx=self.physics_dt, initial=0) + self.clock_offsets[sc]
+ self.clock_noise_fluctuations_withinitial[sc],
+ self.physics_size + self.initial_telemetry_physics_size),
+ dx=self.physics_dt, initial=0)[self.initial_telemetry_physics_size:]
+ )
+
+ self.local_timer_deviations_withinitial = \
+ self.clock_offsets + \
+ self.clock_freqoffsets * self.physics_t_withinitial + \
+ self.clock_freqlindrifts * self.physics_t_withinitial**2 / 2 + \
+ self.clock_freqquaddrifts * self.physics_t_withinitial**3 / 3 + \
+ ForEachSC(lambda sc:
+ cumulative_trapezoid(np.broadcast_to(
+ self.clock_noise_fluctuations_withinitial[sc],
+ self.physics_size + self.initial_telemetry_physics_size),
+ dx=self.physics_dt, initial=0)
)
## Timer deviations from TCB
- self.tcb_timer_deviations = \
- self.local_timer_deviations \
- + self.tps_proper_time_deviations \
- + self.tcb_sync_noises
+ physics_to_telemetry = lambda _, x: x[::self.physics_upsampling * self.telemetry_downsampling]
+ self.tcb_timer_deviations = self.tcb_sync_noises \
+ + self.local_timer_deviations_withinitial.transformed(physics_to_telemetry)
## TDIR tone
@@ -1334,8 +1340,18 @@ class Instrument:
+ self.clock_freqquaddrifts * t**2
logger.debug("Generating clock noise fluctuations")
- self.clock_noise_fluctuations = ForEachSC(lambda sc:
- noises.clock(self.physics_fs, self.physics_size, self.clock_asds[sc])
+
+ # Include initial telemetry time period
+ self.clock_noise_fluctuations_withinitial = ForEachSC(lambda sc:
+ noises.clock(
+ self.physics_fs,
+ self.physics_size + self.initial_telemetry_physics_size,
+ self.clock_asds[sc])
+ )
+
+ # 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:]
)
## Modulation noise
@@ -1667,7 +1683,8 @@ class Instrument:
self._write_attr(hdf5,
'git_url', 'version',
'dt', 't0', 'size', 'fs', 'duration',
- 'telemetry_t0', 'telemetry_downsampling', 'telemetry_fs', 'telemetry_size',
+ 'initial_telemetry_size', 'telemetry_downsampling', 'telemetry_fs',
+ 'telemetry_size', 'telemetry_t0',
'physics_upsampling', 'physics_size', 'physics_dt', 'physics_fs',
'central_freq', 'lock_config', 'lock',
'fplan_file', 'fplan',
diff --git a/tests/test_instrument.py b/tests/test_instrument.py
index 54adb731c58fedd3dcdf3855f0bd5df966c53c42..bb97ee9972d855e9c58057c502fcba1c21ec091c 100755
--- a/tests/test_instrument.py
+++ b/tests/test_instrument.py
@@ -38,11 +38,13 @@ def test_esa_orbits_1_0_2():
instru = Instrument(size=100, orbits='tests/esa-orbits-1-0-2.h5')
instru.simulate()
+@pytest.mark.skip("see https://gitlab.in2p3.fr/lisa-simulation/instrument/-/issues/83")
def test_keplerian_orbits_2_0_dev():
"""Test that simulations can run with Keplerian orbit files v2.0.dev."""
instru = Instrument(size=100, orbits='tests/keplerian-orbits-2-0-dev.h5')
instru.simulate()
+@pytest.mark.skip("see https://gitlab.in2p3.fr/lisa-simulation/instrument/-/issues/83")
def test_esa_trailing_orbits_2_0_dev():
"""Test that simulations can run with ESA trailing orbit files v2.0.dev."""
instru = Instrument(size=100, orbits='tests/esa-trailing-orbits-2-0-dev.h5')