diff --git a/.pylintrc b/.pylintrc
index 8a1b76f74082030df0add1d4ab99205360b68ec0..4b7439a7ecbf61c934686b421f45814890899b26 100644
--- a/.pylintrc
+++ b/.pylintrc
@@ -19,3 +19,6 @@ max-locals=30
# Maximum number of arguments for function / method
max-args=10
+
+# Maximum number of public methods
+max-public-methods=30
diff --git a/lisainstrument/instrument.py b/lisainstrument/instrument.py
index 47ec2c45ad3d0c2ac1b01d18388bdabd85a41dc1..8413f7bc01ee64eef2d83ed42c5af326d2ce8649 100644
--- a/lisainstrument/instrument.py
+++ b/lisainstrument/instrument.py
@@ -12,6 +12,7 @@ import logging
import h5py
import scipy.interpolate
import numpy
+import matplotlib.pyplot
from .containers import ForEachSC
from .containers import ForEachMOSA
@@ -31,7 +32,7 @@ class Instrument:
def __init__(self, size=2592000, dt=0.3, t0=0,
# Inter-spacecraft propagation
- orbits=None, gws=None, interpolation=('lagrange', 31),
+ orbits=None, gws=None, glitches=None, interpolation=('lagrange', 31),
# Lasers
laser_asds=28.2, modulation_asds=1E-14, modulation_fknees=1.5E-2,
central_freq=2.816E14, modulation_freqs=None, offsets_freqs=None, three_lasers=False,
@@ -52,6 +53,7 @@ class Instrument:
t0: initial time [s]
orbits: path to orbit file or None for default set of proper pseudo-ranges and derivatives thereof
gws: dictionary of gravitational-wave responses or path to gravitational-wave file
+ orbits: dictionary of glitch signals for injection points or path to glitch file
interpolation: interpolation function or interpolation method and parameters; when using 'lagrange',
use a tuple ('lagrange', order) with `order` the odd interpolation order; an arbitrary function should
take (x, shift [number of samples]) as parameter
@@ -152,12 +154,10 @@ class Instrument:
'12': 8.1E6, '23': 9.2E6, '31': 10.3E6, '13': 1.4E6, '32': -11.6E6, '21': -9.5E6,
})
- # Orbits and gravitational waves
- self.gws = None
- self.pprs = None
- self.d_pprs = None
+ # Orbits, gravitational waves, glitches
self.init_orbits(orbits)
self.init_gws(gws)
+ self.init_glitches(glitches)
# Interpolation and antialiasing filter
self.init_interpolation(interpolation)
@@ -282,8 +282,8 @@ class Instrument:
self.gw_file = gws
logging.info("Interpolating gravitational-wave responses from GW file '%s'", self.gw_file)
gwf = h5py.File(self.gw_file, 'r')
- self.gws = ForEachMOSA(lambda mosa: scipy.interpolate.InterpolatedUnivariateSpline(
- gwf['t'][:], gwf[f'l_{mosa}'][:], k=1, ext='raise')(self.physics_t)
+ self.glitch_tms = ForEachMOSA(lambda mosa: scipy.interpolate.InterpolatedUnivariateSpline(
+ gwf['t'][:], gwf[f'l_{mosa}'][:], k=5, ext='raise')(self.physics_t)
)
gwf.close()
elif gws is None:
@@ -295,6 +295,26 @@ class Instrument:
self.gw_file = None
self.gws = ForEachMOSA(gws)
+ def init_glitches(self, glitches):
+ """Initialize glitches."""
+ self.glitch_tms = ForEachMOSA(0)
+ if isinstance(glitches, str):
+ self.glitch_file = glitches
+ logging.info("Interpolating glitch signals from glitch file '%s'", self.glitch_file)
+ glitchf = h5py.File(self.glitch_file, 'r')
+ for mosa in self.MOSAS:
+ inj_point = f'tm_{mosa}'
+ if inj_point in glitchf:
+ logging.debug("Interpolating injection point '%s'", inj_point)
+ self.glitch_tms[mosa] = scipy.interpolate.InterpolatedUnivariateSpline(
+ glitchf['t'][:], glitchf[inj_point][:], k=5, ext='const')(self.physics_t)
+ glitchf.close()
+ elif glitches is None:
+ logging.debug("No glitch")
+ self.glitch_file = None
+ else:
+ raise ValueError(f"invalid value '{glitches}' for glitches")
+
def disable_all_noises(self, but=None):
"""Turn off all instrumental noises.
@@ -540,6 +560,7 @@ class Instrument:
logging.info("Computing test-mass carrier beatnote frequency fluctuations on TPS")
self.tps_tm_carrier_fluctuations = ForEachMOSA(lambda mosa:
self.adjacent_tm_carrier_fluctuations[mosa] - self.local_tm_carrier_fluctuations[mosa]
+ + self.glitch_tms[mosa]
)
logging.info("Computing test-mass upper sideband beatnote frequency offsets on TPS")
@@ -1058,3 +1079,137 @@ class Instrument:
logging.info("Closing output file '%s'", output)
hdf5.close()
+
+ def plot_fluctuations(self, output=None, skip=0):
+ """Plot measurements generated by the simulation.
+
+ Args:
+ output: output file, None to show the plots
+ skip: number of initial samples to skip [samples]
+ """
+ # Run simulation if needed
+ if not self.simulated:
+ self.simulate()
+ # Plot signals
+ logging.info("Plotting beatnote frequency fluctuations")
+ _, axes = matplotlib.pyplot.subplots(3, 1, figsize=(12, 12))
+ plot = lambda axis, x, label: axis.plot(self.t[skip:], numpy.broadcast_to(x, self.size)[skip:], label=label)
+ for mosa in self.MOSAS:
+ plot(axes[0], self.isc_carrier_fluctuations[mosa], mosa)
+ plot(axes[1], self.tm_carrier_fluctuations[mosa], mosa)
+ plot(axes[2], self.ref_carrier_fluctuations[mosa], mosa)
+ # Format plot
+ axes[0].set_title("Beatnote frequency fluctuations")
+ axes[2].set_xlabel("Time [s]")
+ axes[0].set_ylabel("Inter-spacecraft frequency [Hz]")
+ axes[1].set_ylabel("Test-mass frequency [Hz]")
+ axes[2].set_ylabel("Reference frequency [Hz]")
+ for axis in axes:
+ axis.grid()
+ axis.legend()
+ # Save or show glitch
+ if output is not None:
+ logging.info("Saving plot to %s", output)
+ matplotlib.pyplot.savefig(output, bbox_inches='tight')
+ else:
+ matplotlib.pyplot.show()
+
+ def plot_offsets(self, output=None, skip=0):
+ """Plot measurements generated by the simulation.
+
+ Args:
+ output: output file, None to show the plots
+ skip: number of initial samples to skip [samples]
+ """
+ # Run simulation if needed
+ if not self.simulated:
+ self.simulate()
+ # Plot signals
+ logging.info("Plotting beatnote frequency offsets")
+ _, axes = matplotlib.pyplot.subplots(3, 1, figsize=(12, 12))
+ plot = lambda axis, x, label: axis.plot(self.t[skip:], numpy.broadcast_to(x, self.size)[skip:], label=label)
+ for mosa in self.MOSAS:
+ plot(axes[0], self.isc_carrier_offsets[mosa], mosa)
+ plot(axes[1], self.tm_carrier_offsets[mosa], mosa)
+ plot(axes[2], self.ref_carrier_offsets[mosa], mosa)
+ # Format plot
+ axes[0].set_title("Beatnote frequency offsets")
+ axes[2].set_xlabel("Time [s]")
+ axes[0].set_ylabel("Inter-spacecraft frequency [Hz]")
+ axes[1].set_ylabel("Test-mass frequency [Hz]")
+ axes[2].set_ylabel("Reference frequency [Hz]")
+ for axis in axes:
+ axis.grid()
+ axis.legend()
+ # Save or show glitch
+ if output is not None:
+ logging.info("Saving plot to %s", output)
+ matplotlib.pyplot.savefig(output, bbox_inches='tight')
+ else:
+ matplotlib.pyplot.show()
+
+ def plot_totals(self, output=None, skip=0):
+ """Plot measurements generated by the simulation.
+
+ Args:
+ output: output file, None to show the plots
+ skip: number of initial samples to skip [samples]
+ """
+ # Run simulation if needed
+ if not self.simulated:
+ self.simulate()
+ # Plot signals
+ logging.info("Plotting beatnote total frequencies")
+ _, axes = matplotlib.pyplot.subplots(3, 1, figsize=(12, 12))
+ plot = lambda axis, x, label: axis.plot(self.t[skip:], numpy.broadcast_to(x, self.size)[skip:], label=label)
+ for mosa in self.MOSAS:
+ plot(axes[0], self.isc_carriers[mosa], mosa)
+ plot(axes[1], self.tm_carriers[mosa], mosa)
+ plot(axes[2], self.ref_carriers[mosa], mosa)
+ # Format plot
+ axes[0].set_title("Beatnote total frequencies")
+ axes[2].set_xlabel("Time [s]")
+ axes[0].set_ylabel("Inter-spacecraft frequency [Hz]")
+ axes[1].set_ylabel("Test-mass frequency [Hz]")
+ axes[2].set_ylabel("Reference frequency [Hz]")
+ for axis in axes:
+ axis.grid()
+ axis.legend()
+ # Save or show glitch
+ if output is not None:
+ logging.info("Saving plot to %s", output)
+ matplotlib.pyplot.savefig(output, bbox_inches='tight')
+ else:
+ matplotlib.pyplot.show()
+
+ def plot_mprs(self, output=None, skip=0):
+ """Plot measurements generated by the simulation.
+
+ Args:
+ output: output file, None to show the plots
+ skip: number of initial samples to skip [samples]
+ """
+ # Run simulation if needed
+ if not self.simulated:
+ self.simulate()
+ # Plot signals
+ logging.info("Plotting measured pseudo-ranges")
+ _, axes = matplotlib.pyplot.subplots(2, 1, figsize=(12, 8))
+ plot = lambda axis, x, label: axis.plot(self.t[skip:], numpy.broadcast_to(x, self.size)[skip:], label=label)
+ for mosa in self.MOSAS:
+ plot(axes[0], self.mprs[mosa], mosa)
+ plot(axes[1], numpy.gradient(self.mprs[mosa], self.dt), mosa)
+ # Format plot
+ axes[0].set_title("Measured pseudo-ranges")
+ axes[1].set_xlabel("Time [s]")
+ axes[0].set_ylabel("Pseudo-range [s]")
+ axes[1].set_ylabel("Pseudo-range derivative [s/s]")
+ for axis in axes:
+ axis.grid()
+ axis.legend()
+ # Save or show glitch
+ if output is not None:
+ logging.info("Saving plot to %s", output)
+ matplotlib.pyplot.savefig(output, bbox_inches='tight')
+ else:
+ matplotlib.pyplot.show()