From d5068c67ef062e2042495435114c376e00950396 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Bayle <j2b.bayle@gmail.com>
Date: Fri, 1 Jul 2022 12:56:33 +0200
Subject: [PATCH] Read TPS GW responses
---
lisainstrument/instrument.py | 71 ++++++++++++++++++++++++++++++------
tests/test_gws.py | 62 ++++++++++++++++++++++++++++++-
2 files changed, 120 insertions(+), 13 deletions(-)
diff --git a/lisainstrument/instrument.py b/lisainstrument/instrument.py
index 144789e..595b17b 100755
--- a/lisainstrument/instrument.py
+++ b/lisainstrument/instrument.py
@@ -113,11 +113,13 @@ class Instrument:
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
- orbit_dataset: datasets to read from the orbit file, must be 'tps/ppr' or 'tcb/ltt';
+ orbit_dataset: datasets to read from orbit file, must be 'tps/ppr' or 'tcb/ltt';
if set to 'tps/ppr', read proper pseudo-ranges (PPRs) in TPSs (proper times),
if set to 'tcb/ltt', read light travel times (LTTs) in TCB (coordinate time);
ignored if no orbit files are used
- gws: path to gravitational-wave file, or dictionary of gravitational-wave responses
+ gws: path to gravitational-wave file, or dictionary of gravitational-wave responses;
+ if ``orbit_dataset`` is ``'tps/ppr'``, we try to read link responses as functions
+ of the TPS instead of link responses in the TCB (fallback behavior)
interpolation: interpolation function or interpolation method and parameters;
use a tuple ('lagrange', order) with `order` the odd Lagrange interpolation order;
an arbitrary function should take (x, shift [number of samples]) as parameter
@@ -683,16 +685,11 @@ class Instrument:
if version.is_devrelease:
logger.warning("You are using a GW file in a development version")
if version in SpecifierSet('< 1.1', True):
- interpolate = lambda data, t: InterpolatedUnivariateSpline(gwf['t'][:], data, k=5, ext='zeros')(t)
- self.gws = ForEachMOSA(lambda mosa: interpolate(gwf[f'l_{mosa}'][:], self.physics_t))
+ self._init_gw_file_before_1_1(gwf)
elif version in SpecifierSet('== 1.1', True):
- interpolate = lambda data, t: InterpolatedUnivariateSpline(gwf['t'][:], data, k=5, ext='zeros')(t)
- self.gws = ForEachMOSA(lambda mosa: interpolate(gwf[f'tcb/l_{mosa}'][:], self.physics_t))
+ self._init_gw_file_1_1(gwf)
elif version in SpecifierSet('== 2.*', True):
- times = gwf.attrs['t0'] + np.arange(gwf.attrs['size']) * gwf.attrs['dt']
- interpolate = lambda data, t: InterpolatedUnivariateSpline(times, data, k=5, ext='zeros')(t)
- link_index = {'12': 0, '23': 1, '31': 2, '13': 3, '32': 4, '21': 5}
- self.gws = ForEachMOSA(lambda mosa: interpolate(gwf['tcb/y'][:, link_index[mosa]], self.physics_t))
+ self._init_gw_file_2(gwf)
else:
raise ValueError(f"unsupported GW file version '{version}'")
elif gws is None:
@@ -704,6 +701,58 @@ class Instrument:
self.gw_file = None
self.gws = ForEachMOSA(gws)
+ def _init_gw_file_before_1_1(self, gwf):
+ """Initialize GW responses from GW file version < 1.1.
+
+ Args:
+ gwf (:obj:`h5py.File`): GW file object
+ """
+ self.gw_group = None
+ interpolate = lambda data, t: InterpolatedUnivariateSpline(gwf['t'][:], data, k=5, ext='zeros')(t)
+ self.gws = ForEachMOSA(lambda mosa: interpolate(gwf[f'l_{mosa}'][:], self.physics_t))
+
+ def _init_gw_file_1_1(self, gwf):
+ """Initialize GW responses from GW file version == 1.1.
+
+ Args:
+ gwf (:obj:`h5py.File`): GW file object
+ """
+ if self.orbit_dataset == 'tps/ppr' and 'tps' in gwf:
+ logger.debug("Using link responses in TPS (following orbit dataset)")
+ self.gw_group = 'tps'
+ elif self.orbit_dataset == 'tps/ppr' and 'tcb' in gwf:
+ logger.warning("TPS link responses not found on '%s', fall back to TCB responses", self.gw_file)
+ logger.debug("Using link responses in TCB (inconsistent with orbit dataset)")
+ self.gw_group = 'tcb'
+ else:
+ logger.debug("Using link responses in TCB (following orbit dataset)")
+ self.gw_group = 'tcb'
+
+ interpolate = lambda data, t: InterpolatedUnivariateSpline(gwf['t'][:], data, k=5, ext='zeros')(t)
+ self.gws = ForEachMOSA(lambda mosa: interpolate(gwf[f'{self.gw_group}/l_{mosa}'][:], self.physics_t))
+
+ def _init_gw_file_2(self, gwf):
+ """Initialize GW responses from GW file version == 2.*.
+
+ Args:
+ gwf (:obj:`h5py.File`): GW file object
+ """
+ if self.orbit_dataset == 'tps/ppr' and 'tps' in gwf:
+ logger.debug("Using link responses in TPS (following orbit dataset)")
+ self.gw_group = 'tps'
+ elif self.orbit_dataset == 'tps/ppr' and 'tcb' in gwf:
+ logger.warning("TPS link responses not found on '%s', fall back to TCB responses", self.gw_file)
+ logger.debug("Using link responses in TCB (inconsistent with orbit dataset)")
+ self.gw_group = 'tcb'
+ else:
+ logger.debug("Using link responses in TCB (following orbit dataset)")
+ self.gw_group = 'tcb'
+
+ link_index = {'12': 0, '23': 1, '31': 2, '13': 3, '32': 4, '21': 5}
+ times = gwf.attrs['t0'] + np.arange(gwf.attrs['size']) * gwf.attrs['dt']
+ interpolate = lambda data, t: InterpolatedUnivariateSpline(times, data, k=5, ext='zeros')(t)
+ self.gws = ForEachMOSA(lambda mosa: interpolate(gwf[f'{self.gw_group}/y'][:, link_index[mosa]], self.physics_t))
+
def init_glitches(self, glitches):
"""Initialize glitches.
@@ -1724,7 +1773,7 @@ class Instrument:
'sc_jitter_theta_asds', 'mosa_jitter_phi_asds',
'dws_asds', 'mosa_angles',
'orbit_file', 'orbit_dataset',
- 'gw_file',
+ 'gw_file', 'gw_group',
'glitch_file',
'interpolation_order',
'electro_delays_isis', 'electro_delays_tmis', 'electro_delays_rfis',
diff --git a/tests/test_gws.py b/tests/test_gws.py
index 5951f04..9691aea 100644
--- a/tests/test_gws.py
+++ b/tests/test_gws.py
@@ -90,7 +90,8 @@ def test_gw_file_1_1():
aafilter=None,
lock='six',
gws=gw_file,
- orbits='tests/keplerian-orbits-1-0-2.h5')
+ orbits='tests/keplerian-orbits-1-0-2.h5',
+ orbit_dataset='tcb/ltt')
instru.disable_all_noises()
assert instru.gw_file is gw_file
@@ -109,6 +110,34 @@ def test_gw_file_1_1():
assert approx(instru.tps_isi_usb_fluctuations[mosa]) == \
-(instru.central_freq + instru.local_usb_offsets[mosa]) * gwf[f'tcb/l_{mosa}'][:]
+ instru = Instrument(
+ size=100,
+ t0=gwf.attrs['t0'],
+ dt=gwf.attrs['dt'],
+ physics_upsampling=1,
+ aafilter=None,
+ lock='six',
+ gws=gw_file,
+ orbits='tests/keplerian-orbits-1-0-2.h5',
+ orbit_dataset='tps/ppr')
+ instru.disable_all_noises()
+
+ assert instru.gw_file is gw_file
+ for mosa in instru.MOSAS:
+ assert instru.gws[mosa] == approx(gwf[f'tcb/l_{mosa}'][:])
+
+ instru.simulate()
+
+ for mosa in instru.MOSAS:
+ assert approx(instru.distant_carrier_fluctuations[mosa]) == \
+ -(instru.central_freq + instru.local_carrier_offsets[mosa]) * gwf[f'tps/l_{mosa}'][:]
+ assert approx(instru.distant_usb_fluctuations[mosa]) == \
+ -(instru.central_freq + instru.local_usb_offsets[mosa]) * gwf[f'tps/l_{mosa}'][:]
+ assert approx(instru.tps_isi_carrier_fluctuations[mosa]) == \
+ -(instru.central_freq + instru.local_carrier_offsets[mosa]) * gwf[f'tps/l_{mosa}'][:]
+ assert approx(instru.tps_isi_usb_fluctuations[mosa]) == \
+ -(instru.central_freq + instru.local_usb_offsets[mosa]) * gwf[f'tps/l_{mosa}'][:]
+
def test_gw_file_2_0_dev():
"""Test that we can read GW files v2.0.dev with orbit files v2.0.
@@ -144,7 +173,8 @@ def test_gw_file_2_0_dev():
aafilter=None,
lock='six',
gws=gw_file,
- orbits='tests/keplerian-orbits-2-0.h5')
+ orbits='tests/keplerian-orbits-2-0.h5',
+ orbit_dataset='tcb/ltt')
instru.disable_all_noises()
assert instru.gw_file is gw_file
@@ -162,3 +192,31 @@ def test_gw_file_2_0_dev():
-(instru.central_freq + instru.local_carrier_offsets[mosa]) * gwf['tcb/y'][:, imosa]
assert approx(instru.tps_isi_usb_fluctuations[mosa]) == \
-(instru.central_freq + instru.local_usb_offsets[mosa]) * gwf['tcb/y'][:, imosa]
+
+ instru = Instrument(
+ size=100,
+ t0=gwf.attrs['t0'],
+ dt=gwf.attrs['dt'],
+ physics_upsampling=1,
+ aafilter=None,
+ lock='six',
+ gws=gw_file,
+ orbits='tests/keplerian-orbits-2-0.h5',
+ orbit_dataset='tps/ppr')
+ instru.disable_all_noises()
+
+ assert instru.gw_file is gw_file
+ for imosa, mosa in enumerate(instru.MOSAS):
+ assert instru.gws[mosa] == approx(gwf['tcb/y'][:, imosa])
+
+ instru.simulate()
+
+ for imosa, mosa in enumerate(instru.MOSAS):
+ assert approx(instru.distant_carrier_fluctuations[mosa]) == \
+ -(instru.central_freq + instru.local_carrier_offsets[mosa]) * gwf['tps/y'][:, imosa]
+ assert approx(instru.distant_usb_fluctuations[mosa]) == \
+ -(instru.central_freq + instru.local_usb_offsets[mosa]) * gwf['tps/y'][:, imosa]
+ assert approx(instru.tps_isi_carrier_fluctuations[mosa]) == \
+ -(instru.central_freq + instru.local_carrier_offsets[mosa]) * gwf['tps/y'][:, imosa]
+ assert approx(instru.tps_isi_usb_fluctuations[mosa]) == \
+ -(instru.central_freq + instru.local_usb_offsets[mosa]) * gwf['tps/y'][:, imosa]
--
GitLab