From cf91651d61d25dc2f4331dc443f6f01bd8833ea6 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Bayle <j2b.bayle@gmail.com>
Date: Mon, 28 Mar 2022 09:55:59 +0200
Subject: [PATCH] Add support for orbit files v2.0
---
lisainstrument/instrument.py | 70 +++++++++++++++++++++++++++++++-----
1 file changed, 61 insertions(+), 9 deletions(-)
diff --git a/lisainstrument/instrument.py b/lisainstrument/instrument.py
index f5a63ad..54c7e71 100644
--- a/lisainstrument/instrument.py
+++ b/lisainstrument/instrument.py
@@ -169,18 +169,27 @@ class Instrument:
self.version = meta.__version__
self.simulated = False
+ # Check orbit dataset
+ if orbit_dataset not in ['tcb/ltt', 'tps/ppr']:
+ raise ValueError(f"invalid orbit dataset '{orbit_dataset}'")
+
# Measurement sampling
if t0 == 'orbits':
if isinstance(orbits, str) and orbits != 'static':
logger.debug("Reading initial time from orbit file '%s'", orbits)
with File(orbits, 'r') as orbitf:
- if orbit_dataset == 'tcb/ltt':
- attr = 't0'
- elif orbit_dataset == 'tps/ppr':
- attr = 'tau0'
+ version = Version(orbitf.attrs['version'])
+ logger.debug("Using orbit file version %s", version)
+ # Warn for orbit file development version
+ if version.is_devrelease:
+ logger.warning("You are using an orbit file in a development version")
+ # Switch between various orbit file standards
+ if version in SpecifierSet('== 1.*', True):
+ self.t0 = float(orbitf.attrs['t0' if orbit_dataset == 'tcb/ltt' else 'tau0'])
+ elif version in SpecifierSet('== 2.*', True):
+ self.t0 = float(orbitf.attrs['t0'])
else:
- raise ValueError(f"invalid orbit dataset '{orbit_dataset}'")
- self.t0 = float(orbitf.attrs[attr])
+ raise ValueError(f"unsupported orbit file version '{version}'")
else:
self.t0 = 0.0
else:
@@ -197,7 +206,19 @@ class Instrument:
if isinstance(orbits, str) and orbits != 'static':
logger.debug("Reading telemetry initial time from orbit file '%s'", orbits)
with File(orbits, 'r') as orbitf:
- self.telemetry_t0 = float(orbitf.attrs['t0'])
+ version = Version(orbitf.attrs['version'])
+ logger.debug("Using orbit file version %s", version)
+ # Warn for orbit file development version
+ if version.is_devrelease:
+ logger.warning("You are using an orbit file in a development 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:
@@ -501,10 +522,14 @@ class Instrument:
with File(self.orbit_file, 'r') as orbitf:
version = Version(orbitf.attrs['version'])
logger.debug("Using orbit file version %s", version)
+ # Warn for orbit file development version
if version.is_devrelease:
logger.warning("You are using an orbit file in a development version")
- if version in SpecifierSet('~= 1.0', True):
+ # Switch between various orbit file standards
+ if version in SpecifierSet('== 1.*', True):
self.init_orbits_file_1_0(orbitf)
+ elif version in SpecifierSet('== 2.*', True):
+ self.init_orbits_file_2_0(orbitf)
else:
raise ValueError(f"unsupported orbit file version '{version}'")
else:
@@ -518,7 +543,7 @@ class Instrument:
self.orbit_dataset = None
def init_orbits_file_1_0(self, orbitf):
- """Initialize orbits from an orbit file version ~= 1.0."""
+ """Initialize orbits from an orbit file version == 1.*."""
def pprs(mosa):
if self.orbit_dataset == 'tcb/ltt':
@@ -562,6 +587,33 @@ class Instrument:
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
+ def init_orbits_file_2_0(self, orbitf):
+ """Initialize orbits from an orbit file version == 2.*."""
+
+ # Prepare common interpolation method
+ times = orbitf.attrs['t0'] + np.arange(orbitf.attrs['size']) * orbitf.attrs['dt']
+ interpolate = lambda data, t: InterpolatedUnivariateSpline(times, data, ext='raise')(t)
+ index = {'12': 0, '23': 1, '31': 2, '13': 3, '32': 4, '21': 5}
+
+ # Interpolate necessary orbital quantities,
+ # show a helpful error message if orbit file is too short
+ try:
+ logger.debug("Interpolating proper pseudo-ranges")
+ dataset = orbitf['tcb/ltt'] if self.orbit_dataset == 'tcb/ltt' else orbitf['tps/ppr']
+ self.pprs = ForEachMOSA(lambda mosa: interpolate(dataset[:, index[mosa]], self.physics_t))
+ logger.debug("Interpolating proper pseudo-range derivatives")
+ dataset = orbitf['tcb/d_ltt'] if self.orbit_dataset == 'tcb/ltt' else orbitf['tps/d_ppr']
+ self.d_pprs = ForEachMOSA(lambda mosa: interpolate(dataset[:, index[mosa]], self.physics_t))
+ logger.debug("Interpolating proper time deviation from TCB")
+ if self.orbit_dataset == 'tcb/ltt':
+ self.tps_proper_time_deviations = ForEachSC(lambda sc: 0)
+ else:
+ dataset = orbitf['tcb/delta_tau']
+ self.tps_proper_time_deviations = ForEachSC(lambda sc: interpolate(dataset[:, sc - 1], self.physics_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
+
def init_gws(self, gws):
"""Initialize gravitational-wave responses."""
if isinstance(gws, str):
--
GitLab