Add plotting methods to Instrument

.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

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()