.gitattributes

@@ -5,4 +5,6 @@ tests/esa-trailing-orbits-2-0.h5 filter=lfs diff=lfs merge=lfs -text
 tests/keplerian-fplan-1-1.h5 filter=lfs diff=lfs merge=lfs -text
 tests/esa-trailing-fplan-1-1.h5 filter=lfs diff=lfs merge=lfs -text
 tests/gws-1-1.h5 filter=lfs diff=lfs merge=lfs -text
-tests/gws-2-0-dev.h5 filter=lfs diff=lfs merge=lfs -text
+tests/gws-2-0.h5 filter=lfs diff=lfs merge=lfs -text
+tests/glitch-1-1.h5 filter=lfs diff=lfs merge=lfs -text
+tests/glitch-1-0.h5 filter=lfs diff=lfs merge=lfs -text

lisainstrument/hexagon.py

@@ -45,7 +45,8 @@ class Hexagon():
                  primary_laser_asd=100,
                  locked_laser_asd=60,
                  offset_freqs='default',
-                 central_freq=2.816E14):
+                 central_freq=2.816E14,
+                 laser_shape='white+infrared'):
 
         self.simulated = False
 
@@ -60,6 +61,7 @@ class Hexagon():
 
         self.primary_laser_asd = float(primary_laser_asd)
         self.locked_laser_asd = float(locked_laser_asd)
+        self.laser_shape = laser_shape
 
         self.central_freq = float(central_freq)
         if offset_freqs == 'default':
@@ -90,7 +92,7 @@ class Hexagon():
         logger.debug("Generating laser noise")
         self.laser_noises = np.empty((self.size, 3)) # (size, laser) [Hz]
         # Laser 1 has its own stability
-        self.laser_noises[:, 0] = noises.laser(self.fs, self.size, self.primary_laser_asd)
+        self.laser_noises[:, 0] = noises.laser(self.fs, self.size, self.primary_laser_asd, self.laser_shape)
         # Laser 2 replicated laser 1 with added locking noise
         self.laser_noises[:, 1] = self.laser_noises[:, 0] \
             + noises.white(self.fs, self.size, self.locked_laser_asd)

lisainstrument/instrument.py

@@ -75,8 +75,8 @@ class Instrument:
                  glitches=None,
                  # Laser locking and frequency plan
                  lock='N1-12', fplan='static',
-                 laser_asds=28.2, central_freq=2.816E14,
-                 offset_freqs='default',
+                 laser_asds=30, laser_shape='white+infrared',
+                 central_freq=2.816E14, offset_freqs='default',
                  # Laser phase modulation
                  modulation_asds='default', modulation_freqs='default', tdir_tone=None,
                  # Clocks
@@ -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
@@ -127,6 +129,7 @@ class Instrument:
             fplan: path to frequency-plan file, dictionary of locking beatnote frequencies [Hz],
                 or 'static' for a default set of constant locking beatnote frequencies
             laser_asds: dictionary of amplitude spectral densities for laser noise [Hz/sqrt(Hz)]
+            laser_shape: laser noise spectral shape, either 'white' or 'white+infrared'
             central_freq: laser central frequency from which all offsets are computed [Hz]
             offset_freqs: dictionary of laser frequency offsets for unlocked lasers [Hz],
                 defined with respect to :attr:`lisainstrument.Instrument.central_freq`,
@@ -156,13 +159,13 @@ class Instrument:
             ttl_coeffs: tuple (local_phi, distant_phi, local_eta, distant_eta) of dictionaries of
                 tilt-to-length coefficients on each MOSA [m/rad], 'default' for a default set of
                 coefficients, or 'random' to randomly generate a set of coefficients in a uniform
-                distribution [-1.6, 1.6] µm/rad (for local) and [-1.65, 1.65] µm/rad (for distant)
+                distribution [-2.3, 2.3] mm/rad
             sc_jitter_asds: tuple of dictionaries of angular jitter amplitude spectral densities
                 for spacecraft, ordered as (yaw, pitch, roll) [rad/sqrt(Hz)]
             mosa_jitter_asds: tuple of dictionaries of angular jitter amplitude spectral densities
                 for MOSA, ordered as (yaw, pitch) [rad/sqrt(Hz)]
             mosa_angles: dictionary of oriented MOSA opening angles [deg], or 'default'
-            dws_asds: dictionary of amplitude spectral densities for DWS measurement noise [rad/s/sqrt(Hz)]
+            dws_asds: dictionary of amplitude spectral densities for DWS measurement noise [rad/sqrt(Hz)]
             ranging_biases: dictionary of ranging noise bias [s]
             ranging_asds: dictionary of ranging noise amplitude spectral densities [s/sqrt(Hz)]
             prn_ambiguity: distance after which PRN code repeats itself [m] (reasonable value is 300 km),
@@ -248,6 +251,10 @@ class Instrument:
 
         # Laser and modulation noise
         self.laser_asds = ForEachMOSA(laser_asds)
+        if laser_shape not in ['white', 'white+infrared']:
+            raise ValueError(f"invalid laser noise spectral shape '{laser_shape}'")
+        self.laser_shape = laser_shape
+
         if modulation_asds == 'default':
             # Default based on the performance model, with 10x amplification for right-sided
             # MOSAs, to account for errors in the frequency distribution system
@@ -682,27 +689,76 @@ 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:
             logger.debug("No gravitational-wave responses")
             self.gw_file = None
+            self.gw_group = None
             self.gws = ForEachMOSA(0)
         else:
             logger.info("Using user-provided gravitational-wave responses")
             self.gw_file = None
+            self.gw_group = 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.
 
@@ -1508,7 +1564,11 @@ class Instrument:
             mosa: laser index
         """
         logger.info("Generating laser noise for laser %s", mosa)
-        self.laser_noises[mosa] = noises.laser(self.physics_fs, self.physics_size, self.laser_asds[mosa])
+        self.laser_noises[mosa] = noises.laser(
+            fs=self.physics_fs,
+            size=self.physics_size,
+            asd=self.laser_asds[mosa],
+            shape=self.laser_shape)
 
         logger.debug("Computing carrier offsets for primary local beam %s", mosa)
         self.local_carrier_offsets[mosa] = self.offset_freqs[mosa]
@@ -1723,7 +1783,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',

lisainstrument/meta.py

@@ -2,6 +2,6 @@
 # -*- coding: utf-8 -*-
 # pylint: disable=missing-module-docstring
 
-__version__ = '1.1.dev'
+__version__ = '1.2.dev'
 __author__ = 'Jean-Baptiste Bayle'
 __email__ = 'j2b.bayle@gmail.com'

lisainstrument/noises.py

@@ -134,7 +134,7 @@ def infrared(fs, size, asd):
     return np.cumsum(red_noise) * (2 * pi / fs)
 
 
-def laser(fs, size, asd=30, fknee=2E-3):
+def laser(fs, size, asd, shape):
     """Generate laser noise [Hz].
 
     This is a white noise with an infrared relaxation towards low frequencies,
@@ -143,16 +143,30 @@ def laser(fs, size, asd=30, fknee=2E-3):
         S_p(f) = asd^2 [ 1 + (fknee / f)^4 ]
                = asd^2 + asd^2 fknee^4 / f^4.
 
+    The low-frequency part (infrared relaxation) can be disabled, in which
+    case the noise shape becomes
+
+        S_p(f) = asd^2.
+
     Args:
         asd: amplitude spectral density [Hz/sqrt(Hz)]
         fknee: cutoff frequency [Hz]
+        shape: spectral shape, either 'white' or 'white+infrared'
     """
-    logger.debug("Generating laser noise (fs=%s Hz, size=%s, asd=%s Hz/sqrt(Hz), fknee=%s Hz)",
-        fs, size, asd, fknee)
-    return white(fs, size, asd) + infrared(fs, size, asd * fknee**2)
+    fknee = 2E-3
+    logger.debug(
+        "Generating laser noise (fs=%s Hz, size=%s, asd=%s "
+        "Hz/sqrt(Hz), fknee=%s Hz, shape=%s)",
+        fs, size, asd, fknee, shape)
+
+    if shape == 'white':
+        return white(fs, size, asd)
+    if shape == 'white+infrared':
+        return white(fs, size, asd) + infrared(fs, size, asd * fknee**2)
+    raise ValueError(f"invalid laser noise spectral shape '{shape}'")
 
 
-def clock(fs, size, asd=6.32E-14):
+def clock(fs, size, asd):
     """Generate clock noise fluctuations [ffd].
 
     The power spectral density in fractional frequency deviations is a pink noise,
@@ -166,7 +180,7 @@ def clock(fs, size, asd=6.32E-14):
     return pink(fs, size, asd)
 
 
-def modulation(fs, size, asd=5.2E-14):
+def modulation(fs, size, asd):
     """Generate modulation noise [ffd].
 
     The power spectral density as fractional frequency deviations reads
@@ -182,7 +196,7 @@ def modulation(fs, size, asd=5.2E-14):
     return powerlaw(fs, size, asd, 1/3)
 
 
-def backlink(fs, size, asd=3E-12, fknee=2E-3):
+def backlink(fs, size, asd, fknee):
     """Generate backlink noise as fractional frequency deviation [ffd].
 
     The power spectral density in displacement is given by
@@ -207,7 +221,7 @@ def backlink(fs, size, asd=3E-12, fknee=2E-3):
         + red(fs, size, 2 * pi * asd * fknee**2 / c)
 
 
-def ranging(fs, size, asd=3E-9):
+def ranging(fs, size, asd):
     """Generate stochastic ranging noise [s].
 
     This is a white noise as a timing jitter,
@@ -221,7 +235,7 @@ def ranging(fs, size, asd=3E-9):
     return white(fs, size, asd)
 
 
-def testmass(fs, size, asd=2.4E-15, fknee=0.4E-3):
+def testmass(fs, size, asd, fknee):
     """Generate test-mass acceleration noise [m/s].
 
     Expressed in acceleration, the noise power spectrum reads
@@ -286,7 +300,7 @@ def jitter(fs, size, asd):
     logger.debug("Generating jitter (fs=%s Hz, size=%s, asd=%s rad/sqrt(Hz))", fs, size, asd)
     return violet(fs, size, 2 * pi * asd)
 
-def dws(fs, size, asd=7E-8/335):
+def dws(fs, size, asd):
     """Generate DWS measurement noise.
 
     The power spectral density in angle is given by

tests/__init__.py

+#! /usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""Test setup functions."""
+
+import os
+
+
+def teardown_module(_):
+    """Remove test measurement file after all tests are executed."""
+
+    try:
+        os.remove('measurements.h5')
+    except FileNotFoundError:
+        pass

tests/test_fplan.py

@@ -92,12 +92,15 @@ def test_static_fplan():
     # Check locking beatnotes
     instru = Instrument(size=100, lock='N1-12', fplan='static')
     instru.simulate()
+    instru.write(mode='w')
     assert _consistent_locking_beatnotes(instru, static)
     instru = Instrument(size=100, lock='N1-21', fplan='static')
     instru.simulate()
+    instru.write(mode='w')
     assert _consistent_locking_beatnotes(instru, static)
     instru = Instrument(size=100, lock='N4-12', fplan='static')
     instru.simulate()
+    instru.write(mode='w')
     assert _consistent_locking_beatnotes(instru, static)
 
 def test_static_fplan_valid_with_all_lock_configs():
@@ -112,6 +115,7 @@ def test_static_fplan_valid_with_all_lock_configs():
             instru = Instrument(size=100, lock=f'{topology}-{primary}', fplan='static')
             instru.disable_all_noises()
             instru.simulate()
+            instru.write(mode='w')
 
             for mosa in instru.MOSAS:
                 assert np.all(5E6 <= np.abs(instru.isi_carriers[mosa]))
@@ -146,12 +150,15 @@ def test_constant_unequal_fplan():
     # Check locking beatnotes
     instru = Instrument(size=100, lock='N1-12', fplan=fplan)
     instru.simulate()
+    instru.write(mode='w')
     assert _consistent_locking_beatnotes(instru, fplan)
     instru = Instrument(size=100, lock='N1-21', fplan=fplan)
     instru.simulate()
+    instru.write(mode='w')
     assert _consistent_locking_beatnotes(instru, fplan)
     instru = Instrument(size=100, lock='N4-12', fplan=fplan)
     instru.simulate()
+    instru.write(mode='w')
     assert _consistent_locking_beatnotes(instru, fplan)
 
 def test_varying_equal_fplan():
@@ -167,12 +174,15 @@ def test_varying_equal_fplan():
     fplan = {mosa: fplan for mosa in Instrument.MOSAS}
     instru = Instrument(size=100, lock='N1-12', fplan=fplan)
     instru.simulate()
+    instru.write(mode='w')
     assert _consistent_locking_beatnotes(instru, fplan)
     instru = Instrument(size=100, lock='N1-21', fplan=fplan)
     instru.simulate()
+    instru.write(mode='w')
     assert _consistent_locking_beatnotes(instru, fplan)
     instru = Instrument(size=100, lock='N4-12', fplan=fplan)
     instru.simulate()
+    instru.write(mode='w')
     assert _consistent_locking_beatnotes(instru, fplan)
 
 def test_varying_unequal_fplan():
@@ -190,12 +200,15 @@ def test_varying_unequal_fplan():
     # Check locking beatnotes
     instru = Instrument(size=100, lock='N1-12', fplan=fplan)
     instru.simulate()
+    instru.write(mode='w')
     assert _consistent_locking_beatnotes(instru, fplan)
     instru = Instrument(size=100, lock='N1-21', fplan=fplan)
     instru.simulate()
+    instru.write(mode='w')
     assert _consistent_locking_beatnotes(instru, fplan)
     instru = Instrument(size=100, lock='N4-12', fplan=fplan)
     instru.simulate()
+    instru.write(mode='w')
     assert _consistent_locking_beatnotes(instru, fplan)
 
 def test_keplerian_fplan_1_1():
@@ -216,12 +229,15 @@ def test_keplerian_fplan_1_1():
     # Check locking beatnotes
     instru = Instrument(size=100, lock='N1-12', fplan='tests/keplerian-fplan-1-1.h5')
     instru.simulate()
+    instru.write(mode='w')
     assert _consistent_locking_beatnotes(instru)
     instru = Instrument(size=100, lock='N1-21', fplan='tests/keplerian-fplan-1-1.h5')
     instru.simulate()
+    instru.write(mode='w')
     assert _consistent_locking_beatnotes(instru)
     instru = Instrument(size=100, lock='N4-12', fplan='tests/keplerian-fplan-1-1.h5')
     instru.simulate()
+    instru.write(mode='w')
     assert _consistent_locking_beatnotes(instru)
 
 def test_esa_trailing_fplan_1_1():
@@ -242,12 +258,15 @@ def test_esa_trailing_fplan_1_1():
     # Check locking beatnotes
     instru = Instrument(size=100, lock='N1-12', fplan='tests/esa-trailing-fplan-1-1.h5')
     instru.simulate()
+    instru.write(mode='w')
     assert _consistent_locking_beatnotes(instru)
     instru = Instrument(size=100, lock='N1-21', fplan='tests/esa-trailing-fplan-1-1.h5')
     instru.simulate()
+    instru.write(mode='w')
     assert _consistent_locking_beatnotes(instru)
     instru = Instrument(size=100, lock='N4-12', fplan='tests/esa-trailing-fplan-1-1.h5')
     instru.simulate()
+    instru.write(mode='w')
     assert _consistent_locking_beatnotes(instru)
 
 def test_offset_freqs():
@@ -257,6 +276,7 @@ def test_offset_freqs():
     instru = Instrument(size=100, lock='six', offset_freqs='default')
     instru.disable_all_noises()
     instru.simulate()
+    instru.write(mode='w')
     for mosa in instru.MOSAS:
         assert np.all(5E6 <= np.abs(instru.isi_carriers[mosa]))
         assert np.all(5E6 <= np.abs(instru.isi_usbs[mosa]))
@@ -274,15 +294,18 @@ def test_offset_freqs():
     }
     instru = Instrument(size=100, lock='six', offset_freqs=offset_freqs)
     instru.simulate()
+    instru.write(mode='w')
     for mosa in instru.MOSAS:
         assert instru.local_carrier_offsets[mosa] == offset_freqs[mosa]
 
     # Check that frequency offsets are only set on unlocked lasers
     instru = Instrument(size=100, lock='N1-12', offset_freqs=offset_freqs)
-    instru.simulate()
     instru.disable_all_noises()
+    instru.simulate()
+    instru.write(mode='w')
     assert instru.local_carrier_offsets['12'] == offset_freqs['12']
     instru = Instrument(size=100, lock='N1-23', offset_freqs=offset_freqs)
-    instru.simulate()
     instru.disable_all_noises()
+    instru.simulate()
+    instru.write(mode='w')
     assert instru.local_carrier_offsets['23'] == offset_freqs['23']

tests/test_glitch.py

+#! /usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""Test reading of glitch files."""
+
+import numpy as np
+from pytest import approx
+from h5py import File
+from lisaconstants import c
+from lisainstrument import Instrument
+
+
+def test_no_glitches():
+    """Test that we can simulate no glitches."""
+
+    instru = Instrument(size=100, glitches=None)
+    instru.disable_all_noises()
+
+    assert instru.glitch_file is None
+    for mosa in instru.MOSAS:
+        assert instru.glitch_tms[mosa] == 0
+        assert instru.glitch_lasers[mosa] == 0
+
+    instru.simulate()
+    instru.write(mode='w')
+
+    for mosa in instru.MOSAS:
+        assert instru.local_carrier_fluctuations[mosa] == approx(0)
+        assert instru.local_tmi_carrier_fluctuations[mosa] == approx(0)
+        assert instru.local_tmi_usb_fluctuations[mosa] == approx(0)
+
+def test_glitch_file_1_0():
+    """Test that we can read test-mass and laser glitches from files v1.0.
+
+    Test glitch file can be generated using LISA Glitch and the following script.
+
+    from lisaglitch import ShapeletGlitch
+    path = 'tests/glitch-1-0.h5'
+    ShapeletGlitch(beta=2, quantum_n=1, inj_point='tm_12', t_inj=10.0, size=200).write(path, mode='w')
+    ShapeletGlitch(beta=1, quantum_n=2, inj_point='tm_23', t_inj=0.0).write(path)
+    ShapeletGlitch(beta=3, quantum_n=2, inj_point='tm_31', t_inj=5.0).write(path)
+    ShapeletGlitch(beta=4, quantum_n=1, inj_point='laser_12', t_inj=0.0).write(path)
+    ShapeletGlitch(beta=4, quantum_n=3, inj_point='laser_31', t_inj=0.0).write(path)
+    """
+
+    glitch_file = 'tests/glitch-1-0.h5'
+    glitchf = File(glitch_file, 'r')
+
+    instru = Instrument(
+        size=glitchf.attrs['size'],
+        t0=glitchf.attrs['t0'],
+        dt=glitchf.attrs['dt'],
+        physics_upsampling=1,
+        aafilter=None,
+        lock='six',
+        glitches=glitch_file,
+        orbits='tests/keplerian-orbits-1-0-2.h5')
+    instru.disable_all_noises()
+
+    assert instru.glitch_file is glitch_file
+    for mosa in instru.MOSAS:
+
+        # Missing datasets should be zero
+        tm_dataset = f'tm_{mosa}'
+        if tm_dataset in glitchf:
+            assert instru.glitch_tms[mosa] == approx(glitchf[tm_dataset][:])
+        else:
+            assert instru.glitch_tms[mosa] == 0
+
+        laser_dataset = f'laser_{mosa}'
+        if laser_dataset in glitchf:
+            assert instru.glitch_lasers[mosa] == approx(glitchf[laser_dataset][:])
+        else:
+            assert instru.glitch_lasers[mosa] == 0
+
+    instru.simulate()
+    instru.write(mode='w')
+
+    for mosa in instru.MOSAS:
+
+        laser_dataset = f'laser_{mosa}'
+        if laser_dataset in glitchf:
+            assert approx(instru.local_carrier_fluctuations[mosa]) == \
+                glitchf[laser_dataset][:]
+        else:
+            assert approx(instru.local_carrier_fluctuations[mosa]) == 0
+
+        tm_dataset = f'tm_{mosa}'
+
+        if tm_dataset in glitchf:
+            assert approx(instru.local_tmi_carrier_fluctuations[mosa], abs=1E-6) == \
+                instru.local_carrier_fluctuations[mosa] \
+                + 2 * glitchf[tm_dataset][:] / c \
+                * (instru.central_freq + instru.local_tmi_carrier_offsets[mosa])
+        else:
+            assert approx(instru.local_tmi_carrier_fluctuations[mosa]) == \
+                instru.local_carrier_fluctuations[mosa]
+
+        if tm_dataset in glitchf:
+            assert approx(instru.local_tmi_usb_fluctuations[mosa], abs=1E-6) == \
+                instru.local_usb_fluctuations[mosa] \
+                + 2 * glitchf[tm_dataset][:] / c \
+                * (instru.central_freq + instru.local_tmi_usb_offsets[mosa])
+        else:
+            assert approx(instru.local_tmi_usb_fluctuations[mosa]) == \
+                instru.local_usb_fluctuations[mosa]
+
+def test_glitch_file_1_1():
+    """Test that we can read test-mass and laser glitches from files v1.1.
+
+    Test glitch file can be generated using LISA Glitch and the following script.
+
+    from lisaglitch import ShapeletGlitch
+    path = 'tests/glitch-1-1.h5'
+    ShapeletGlitch(beta=2, quantum_n=1, inj_point='tm_12', t_inj=10.0, size=200).write(path, mode='w')
+    ShapeletGlitch(beta=1, quantum_n=2, inj_point='tm_23', t_inj=0.0).write(path)
+    ShapeletGlitch(beta=3, quantum_n=2, inj_point='tm_31', t_inj=5.0).write(path)
+    ShapeletGlitch(beta=4, quantum_n=1, inj_point='laser_12', t_inj=0.0).write(path)
+    ShapeletGlitch(beta=4, quantum_n=3, inj_point='laser_31', t_inj=0.0).write(path)
+    """
+
+    glitch_file = 'tests/glitch-1-1.h5'
+    glitchf = File(glitch_file, 'r')
+
+    instru = Instrument(
+        size=glitchf.attrs['size'],
+        t0=glitchf.attrs['t0'],
+        dt=glitchf.attrs['dt'],
+        physics_upsampling=1,
+        aafilter=None,
+        lock='six',
+        glitches=glitch_file,
+        orbits='tests/keplerian-orbits-1-0-2.h5')
+    instru.disable_all_noises()
+
+    assert instru.glitch_file is glitch_file
+    for mosa in instru.MOSAS:
+
+        # Missing datasets should be zero
+        tm_dataset = f'tm_{mosa}'
+        if tm_dataset in glitchf:
+            assert instru.glitch_tms[mosa] == approx(glitchf[tm_dataset][:])
+        else:
+            assert instru.glitch_tms[mosa] == 0
+
+        laser_dataset = f'laser_{mosa}'
+        if laser_dataset in glitchf:
+            assert instru.glitch_lasers[mosa] == approx(glitchf[laser_dataset][:])
+        else:
+            assert instru.glitch_lasers[mosa] == 0
+
+    instru.simulate()
+    instru.write(mode='w')
+
+    for mosa in instru.MOSAS:
+
+        laser_dataset = f'laser_{mosa}'
+        if laser_dataset in glitchf:
+            assert approx(instru.local_carrier_fluctuations[mosa]) == \
+                glitchf[laser_dataset][:]
+        else:
+            assert approx(instru.local_carrier_fluctuations[mosa]) == 0
+
+        tm_dataset = f'tm_{mosa}'
+
+        if tm_dataset in glitchf:
+            assert approx(instru.local_tmi_carrier_fluctuations[mosa], abs=1E-6) == \
+                instru.local_carrier_fluctuations[mosa] \
+                + 2 * glitchf[tm_dataset][:] / c \
+                * (instru.central_freq + instru.local_tmi_carrier_offsets[mosa])
+        else:
+            assert approx(instru.local_tmi_carrier_fluctuations[mosa]) == \
+                instru.local_carrier_fluctuations[mosa]
+
+        if tm_dataset in glitchf:
+            assert approx(instru.local_tmi_usb_fluctuations[mosa], abs=1E-6) == \
+                instru.local_usb_fluctuations[mosa] \
+                + 2 * glitchf[tm_dataset][:] / c \
+                * (instru.central_freq + instru.local_tmi_usb_offsets[mosa])
+        else:
+            assert approx(instru.local_tmi_usb_fluctuations[mosa]) == \
+                instru.local_usb_fluctuations[mosa]

tests/test_gws.py

@@ -19,6 +19,7 @@ def test_no_gws():
         assert instru.gws[mosa] == 0
 
     instru.simulate()
+    instru.write(mode='w')
 
     for mosa in instru.MOSAS:
         assert instru.distant_carrier_fluctuations[mosa] == approx(0)
@@ -44,6 +45,7 @@ def test_dict_of_gws():
         assert np.all(instru.gws[mosa] == responses[mosa])
 
     instru.simulate()
+    instru.write(mode='w')
 
     for mosa in instru.MOSAS:
         assert approx(instru.distant_carrier_fluctuations[mosa]) == \
@@ -90,7 +92,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
@@ -98,6 +101,7 @@ def test_gw_file_1_1():
         assert instru.gws[mosa] == approx(gwf[f'tcb/l_{mosa}'][:])
 
     instru.simulate()
+    instru.write(mode='w')
 
     for mosa in instru.MOSAS:
         assert approx(instru.distant_carrier_fluctuations[mosa]) == \
@@ -109,8 +113,37 @@ 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()
+    instru.write(mode='w')
+
+    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.
+    """Test that we can read GW files v2.0 with orbit files v2.0.
 
     Test GW file can be generated using LISA GW Response and the following script.
 
@@ -130,10 +163,10 @@ def test_gw_file_2_0_dev():
             gw_beta=0, gw_lambda=0,
         )
 
-        galbin.write('tests/gws-2-0-dev.h5')
+        galbin.write('tests/gws-2-0.h5')
     """
 
-    gw_file = 'tests/gws-2-0-dev.h5'
+    gw_file = 'tests/gws-2-0.h5'
     gwf = File(gw_file, 'r')
 
     instru = Instrument(
@@ -144,7 +177,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
@@ -152,6 +186,7 @@ def test_gw_file_2_0_dev():
         assert instru.gws[mosa] == approx(gwf['tcb/y'][:, imosa])
 
     instru.simulate()
+    instru.write(mode='w')
 
     for imosa, mosa in enumerate(instru.MOSAS):
         assert approx(instru.distant_carrier_fluctuations[mosa]) == \
@@ -162,3 +197,32 @@ 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()
+    instru.write(mode='w')
+
+    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]

tests/test_hexagon.py

@@ -15,6 +15,6 @@ def test_write():
     """Test that simulations can run and be written."""
     hexagon = Hexagon()
     hexagon.simulate()
-    hexagon.write('test.h5', mode='w')
+    hexagon.write('measurements.h5', mode='w')
 
-    assert os.path.isfile('test.h5')
+    assert os.path.isfile('measurements.h5')

tests/test_instrument.py

@@ -10,16 +10,19 @@ def test_run():
     """Test that simulations can run."""
     instru = Instrument(size=100)
     instru.simulate()
+    instru.write(mode='w')
 
 def test_run_no_aafilter():
     """Test that simulations can run with no filter."""
     instru = Instrument(size=100, aafilter=None)
     instru.simulate()
+    instru.write(mode='w')
 
 def test_run_no_upsampling():
     """Test that simulations can run with no filter."""
     instru = Instrument(size=100, physics_upsampling=1, aafilter=None)
     instru.simulate()
+    instru.write(mode='w')
 
 def test_no_orbit_file():
     """Test that simulations fail with an invalid orbit file."""
@@ -32,30 +35,38 @@ def test_keplerian_orbits_1_0_2():
     """Test that simulations can run with Keplerian orbit files v1.0.2."""
     instru = Instrument(size=100, orbits='tests/keplerian-orbits-1-0-2.h5')
     instru.simulate()
+    instru.write(mode='w')
 
 def test_esa_orbits_1_0_2():
     """Test that simulations can run with ESA orbit files v1.0.2."""
     instru = Instrument(size=100, orbits='tests/esa-orbits-1-0-2.h5')
     instru.simulate()
+    instru.write(mode='w')
 
 def test_keplerian_orbits_2_0():
     """Test that simulations can run with Keplerian orbit files v2.0."""
     instru = Instrument(size=100, orbits='tests/keplerian-orbits-2-0.h5')
     instru.simulate()
+    instru.write(mode='w')
 
 def test_esa_trailing_orbits_2_0():
     """Test that simulations can run with ESA trailing orbit files v2.0."""
     instru = Instrument(size=100, orbits='tests/esa-trailing-orbits-2-0.h5')
     instru.simulate()
+    instru.write(mode='w')
 
 def test_locking():
     """Test that simulations can run with various lock configurations."""
     # Test six free-running lasers
-    Instrument(size=100, lock='six').simulate()
+    instru = Instrument(size=100, lock='six')
+    instru.simulate()
+    instru.write(mode='w')
     # Test non-swap configurations
     for i in range(1, 6):
         for primary in ['12', '23', '31', '13', '32', '21']:
-            Instrument(size=100, lock=f'N{i}-{primary}').simulate()
+            instru = Instrument(size=100, lock=f'N{i}-{primary}')
+            instru.simulate()
+            instru.write(mode='w')
     # Test that any other raises an error
     with pytest.raises(ValueError):
         Instrument(size=100, lock='whatever')

tests/test_ranging.py

@@ -21,6 +21,7 @@ def test_no_prn_ambiguity():
     instru.disable_clock_noises()
     instru.disable_ranging_noises()
     instru.simulate()
+    instru.write(mode='w')
     for mosa in instru.MOSAS:
         assert instru.mprs[mosa][50:] == approx(pprs[mosa])
 
@@ -41,6 +42,7 @@ def test_prn_ambiguity_with_static_orbits():
     instru.disable_clock_noises()
     instru.disable_ranging_noises()
     instru.simulate()
+    instru.write(mode='w')
     for mosa in instru.MOSAS:
         assert instru.mprs[mosa][50:] == approx(pprs[mosa] % (300E3 / c))
 
@@ -49,6 +51,7 @@ def test_prn_ambiguity_with_static_orbits():
     instru.disable_clock_noises()
     instru.disable_ranging_noises()
     instru.simulate()
+    instru.write(mode='w')
     for mosa in instru.MOSAS:
         assert instru.mprs[mosa][50:] == approx(pprs[mosa] % (100E3 / c))
 
@@ -64,6 +67,7 @@ def test_prn_ambiguity_with_esa_orbits():
     instru.disable_clock_noises()
     instru.disable_ranging_noises()
     instru.simulate()
+    instru.write(mode='w')
 
     for mosa in instru.MOSAS:
         assert np.all(instru.mprs[mosa][50:] <= 300E3 / c)