From dc3110206341da49c34c03c770f00be98e36f2a0 Mon Sep 17 00:00:00 2001 From: Jean-Baptiste Bayle <j2b.bayle@gmail.com> Date: Fri, 22 Jan 2021 12:07:45 -0800 Subject: [PATCH] Add instrument simulation --- lisainstrument/__init__.py | 13 +- lisainstrument/dsp.py | 179 +++++ lisainstrument/lisa.py | 1346 +++++++++++++++++++++++++++++++++++- requirements.txt | 15 +- 4 files changed, 1537 insertions(+), 16 deletions(-) create mode 100644 lisainstrument/dsp.py diff --git a/lisainstrument/__init__.py b/lisainstrument/__init__.py index d8004ef..25ceeb8 100644 --- a/lisainstrument/__init__.py +++ b/lisainstrument/__init__.py @@ -6,15 +6,4 @@ from .meta import __version__ from .meta import __author__ from .meta import __email__ -from .noises import white -from .noises import violet -from .noises import pink -from .noises import red -from .noises import infrared - -from .noises import laser -from .noises import clock -from .noises import modulation -from .noises import backlink -from .noises import ranging -from .noises import testmass +from .lisa import Instrument diff --git a/lisainstrument/dsp.py b/lisainstrument/dsp.py new file mode 100644 index 0000000..c3506ae --- /dev/null +++ b/lisainstrument/dsp.py @@ -0,0 +1,179 @@ +#! /usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Digital signal processing. + +Authors: + Jean-Baptiste Bayle <j2b.bayle@gmail.com> +""" + +import logging +import scipy.sparse +import numpy + + +def timeshift(data, shift, order=31): + """Shift `data` in time by `shift` samples. + + Args: + data: array to be shited + shift: scalar or array of time shifts to be applied [samples] + order: interpolation order + """ + logging.debug("Time shifting data '%s' by '%s' samples (order=%d)", data, shift, order) + + if numpy.isscalar(data): + logging.debug("Data is a constant scalar and cannot be time shifted") + return data + if numpy.all(shift == 0): + logging.debug("Time shift is vanishing and cannot be applied") + return data + + data = numpy.asarray(data) + mode = 'timeseries' if isinstance(shift, numpy.ndarray) else 'constant' + logging.debug("Using mode '%s'", mode) + + if mode == 'timeseries' and data.size != shift.size: + raise ValueError(f"`data` and `shift` must be of the same size (got {data.size}, {shift.size})") + + num_tabs = order + 1 + p = num_tabs // 2 # pylint: disable=invalid-name + size = data.size + + k = numpy.floor(shift).astype(int) + eps = shift - k + coeffs = lagrange_coeffs(eps, num_tabs, p) + logging.debug("Using Lagrange coefficiens '%s'", coeffs) + + if mode == 'timeseries': + logging.debug("Computing Lagrange matrix") + indices = numpy.arange(size) + i_min = numpy.min(k - (p - 1) + indices) + i_max = numpy.max(k + p + indices + 1) + csr_ind = numpy.tile(numpy.arange(num_tabs), size) + numpy.repeat(k + indices, num_tabs) - (p - 1) + csr_ptr = num_tabs * numpy.arange(size + 1) + mat = scipy.sparse.csr_matrix((numpy.ravel(coeffs), csr_ind - i_min, csr_ptr), shape=(size, i_max - i_min)) + logging.debug("Padding data (left=%d, right=%d)", max(0, -i_min), max(0, i_max - size)) + data_padded = numpy.pad(data[max(0, i_min):min(size, i_max)], (max(0, -i_min), max(0, i_max - size))) + logging.debug("Computing matrix-vector product") + shifted = mat.dot(data_padded) + elif mode == 'constant': + i_min = k - (p - 1) + i_max = k + p + size + logging.debug("Padding data (left=%d, right=%d)", max(0, -i_min), max(0, i_max - size)) + data_padded = numpy.pad(data[max(0, i_min):min(size, i_max)], (max(0, -i_min), max(0, i_max - size))) + logging.debug("Computing correlation product") + shifted = numpy.correlate(data_padded, coeffs[0], mode="valid") + + return shifted + + +def lagrange_coeffs(eps, num_tabs, p): # pylint: disable=invalid-name + """Calculate coefficients for lagrange interpolation""" + coeffs = numpy.zeros([num_tabs, eps.size]) + + if p > 1: + factor = numpy.ones(eps.size, dtype=numpy.float64) + factor *= eps * (1 - eps) + + for j in range(1, p): + factor *= (-1) * (1 - j / p) / (1 + j / p) + coeffs[p - 1 - j] = factor / (j + eps) + coeffs[p + j] = factor / (j + 1 - eps) + + coeffs[p - 1] = 1 - eps + coeffs[p] = eps + + for j in range(2, p): + coeffs *= 1 - (eps / j)**2 + + coeffs *= (1 + eps) * (1 - eps / p) + else: + coeffs[p - 1] = 1 - eps + coeffs[p] = eps + + return coeffs.T + + +def time_shift(data, shift, order=31): + """Shift data in time by tau + + Args: + data (numpy.ndarray): data to be shited + tau (numpy.ndarray): amount of time each data point is to be shifted + (must be of same dimension as data) + fs (double): sampling frequency in (Hz) + order (int): interpolation order + """ + logging.debug("Time shifting data '%s' by '%s' (order=%d)", data, shift, order) + + if numpy.isscalar(data): + logging.debug("Data is a constant scalar and cannot be time shifted") + return data + if numpy.all(shift == 0): + logging.debug("Time shift is vanishing and cannot be applied") + return data + + data = numpy.asarray(data) + mode = "timeseries" if isinstance(shift, numpy.ndarray) else "constant" + logging.debug("Using mode '%s'", mode) + + if mode == "timeseries" and data.size != shift.size: + raise ValueError(f"`data` and `tau` must be of the same size (got {data.size}, {shift.size})") + + num_tabs = order + 1 + p = num_tabs // 2 # pylint: disable=invalid-name + size = data.size + + def lagrange_coeffs(eps): + """Calculate coefficients for lagrange interpolation""" + coeffs = numpy.zeros([num_tabs, eps.size]) + + if p > 1: + factor = numpy.ones(eps.size, dtype=numpy.float64) + factor *= eps * (1 - eps) + + for j in range(1, p): + factor *= (-1) * (1 - j / p) / (1 + j / p) + coeffs[p - 1 - j] = factor / (j + eps) + coeffs[p + j] = factor / (j + 1 - eps) + + coeffs[p - 1] = 1 - eps + coeffs[p] = eps + + for j in range(2, p): + coeffs *= 1 - (eps / j)**2 + + coeffs *= (1 + eps) * (1 - eps / p) + else: + coeffs[p - 1] = 1 - eps + coeffs[p] = eps + + return coeffs.T + + k = numpy.floor(shift).astype(int) + eps = shift - k + coeffs = lagrange_coeffs(eps) + logging.debug("Using Lagrange coefficiens '%s'", coeffs) + + if mode == "timeseries": + logging.debug("Computing Lagrange matrix") + indices = numpy.arange(size) + i_min = numpy.min(k - (p - 1) + indices) + i_max = numpy.max(k + p + indices + 1) + csr_ind = numpy.tile(numpy.arange(num_tabs), size) + numpy.repeat(k + indices, num_tabs) - (p - 1) + csr_ptr = num_tabs * numpy.arange(size + 1) + mat = scipy.sparse.csr_matrix((numpy.ravel(coeffs), csr_ind - i_min, csr_ptr), shape=(size, i_max - i_min)) + logging.debug("Padding data (left=%d, right=%d)", max(0, -i_min), max(0, i_max - size)) + data_padded = numpy.pad(data[max(0, i_min):min(size, i_max)], (max(0, -i_min), max(0, i_max - size))) + logging.debug("Computing matrix-vector product") + shifted = mat.dot(data_padded) + elif mode == "constant": + i_min = k - (p - 1) + i_max = k + p + size + logging.debug("Padding data (left=%d, right=%d)", max(0, -i_min), max(0, i_max - size)) + data_padded = numpy.pad(data[max(0, i_min):min(size, i_max)], (max(0, -i_min), max(0, i_max - size))) + logging.debug("Computing correlation product") + shifted = numpy.correlate(data_padded, coeffs[0], mode="valid") + + return shifted diff --git a/lisainstrument/lisa.py b/lisainstrument/lisa.py index 0db7853..3537a9b 100644 --- a/lisainstrument/lisa.py +++ b/lisainstrument/lisa.py @@ -1,5 +1,6 @@ #! /usr/bin/env python3 # -*- coding: utf-8 -*- +# pylint: disable=too-many-lines """ LISA Instrument module. @@ -7,5 +8,1346 @@ Authors: Jean-Baptiste Bayle <j2b.bayle@gmail.com> """ -from .meta import __version__ -from .meta import __author__ +import logging +import h5py +import scipy.interpolate +import numpy + +from . import meta +from . import dsp +from . import noises + + +# TODO: +# * modulated beam, 4 args (carrier_offsets, carrier_fluctuations, usb_offsets, usb_fluctuations) +# * args can be scalars (repeated in dict), dict (check and convert to float), callable (call with mosa as arg) + + +class Signal: + """Represent a signal expressed as frequency offsets and fluctuations.""" + + def __init__(self, offsets=0, fluctuations=0): + """Initialize a signal from frequency offsets and fluctuations. + + Args: + offsets: frequency offsets [Hz] + fluctuations: frequency fluctuations [Hz] + """ + if callable(offsets): + self.offsets = {mosa: offsets(mosa) for mosa in Instrument.MOSAS} + else: + self.offsets = Instrument.mosa_dict(offsets) + + if callable(fluctuations): + self.fluctuations = {mosa: fluctuations(mosa) for mosa in Instrument.MOSAS} + else: + self.fluctuations = Instrument.mosa_dict(fluctuations) + + def transformed(self, offsets=lambda x, mosa: x, fluctuations=lambda x, mosa: x): + """Return a new two-variable signal transforming offsets and fluctuations. + + Args: + offsets: function of (offsets, mosa) returning new offsets [Hz] + fluctuations: function of (fluctuations, mosa) returning new fluctuations [Hz] + + Returns: + A new `Signal` isntance where offsets and fluctuations have been transformed. + """ + return self.__class__( + offsets={ + mosa: offsets(self.offsets[mosa], mosa) + for mosa in Instrument.MOSAS + }, + fluctuations={ + mosa: fluctuations(self.fluctuations[mosa], mosa) + for mosa in Instrument.MOSAS + }, + ) + + def reduce(self, function=lambda offsets, fluctuations: 0): + """Compute a new MOSA dictionary from offsets and fluctuations. + + Args: + function: function of (offsets, fluctuations) returning new value + """ + return { + mosa: function(self.offsets[mosa], self.fluctuations[mosa]) + for mosa in Instrument.MOSAS + } + + @property + def totalfreq(self): + """Return total frequencies, as the sum of offsets and fluctuations.""" + return self.reduce(lambda offsets, fluctuations: offsets + fluctuations) + + +class ModulatedBeam: + """Represent a modulated beam, with a carrier, an upper sideband, and optionally timer deviations.""" + + def __init__(self, carrier, usb, timer_deviations=None): + """Initialize a modulated beam from carrier and upper sideband signals. + + Args: + carrier: carrier signal + usb: upper sideband signal + timer_deviations: timer deviations [s] + """ + if not isinstance(carrier) or not isinstance(usb, Signal): + raise TypeError("carrier and upper sideband should be instances of `Signal`") + + self.carrier = carrier + self.usb = usb + self.timer_deviations = Instrument.mosa_dict(timer_deviations) + + def transformed(self, + carrier_offsets=lambda x, mosa: x, carrier_fluctuations=lambda x, mosa: x, + usb_offsets=lambda x, mosa: x, usb_fluctuations=lambda x, mosa: x, + timer_deviations=lambda x, mosa: x): + """Return a new modulated beam after applying transformations. + + Args: + carrier_offsets: function of (offsets, mosa) returning new carrier offsets [Hz] + carrier_fluctuations: function of (fluctuations, mosa) returning new carrier fluctuations [Hz] + usb_offsets: function of (offsets, mosa) returning new upper sideband offsets [Hz] + usb_fluctuations: function of (fluctuations, mosa) returning new upper sideband fluctuations [Hz] + timer_deviations: function of (deviations, mosa) return new timer deviations [s] + + Returns: + A new `ModulatedBeam` isntance where signals have been transformed. + """ + return self.__class__( + carrier=self.carrier.transformed(carrier_offsets, carrier_fluctuations), + usb=self.usb.transformed(usb_offsets, usb_fluctuations), + timer_deviations={ + mosa: timer_deviations(self.timer_deviations[mosa], mosa) + for mosa in Instrument.MOSAS + } + ) + + +class Instrument: + """Represents an instrumental simulation.""" + + # Indexing conventions + SCS = ['1', '2', '3'] + MOSAS = ['12', '23', '31', '13', '32', '21'] + + def __init__(self, size=2592000, dt=0.3, t0=0, + # Inter-spacecraft propagation + orbits=None, pprs=None, d_pprs=None, gws=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, + # Clocks + clock_asds=6.32E-14, clock_offsets=None, clock_freqoffsets=None, clock_freqlindrifts=None, + clock_freqquaddrifts=None, clockinv_tolerance=1E-10, clockinv_maxiter=5, + # Optical pathlength noises + backlink_asds=3E-12, backlink_fknees=2E-3, testmass_asds=2.4E-15, testmass_fknees=0.4E-3, + # Pseudo-ranging + ranging_biases=0, ranging_asds=3E-9, + # Physics simulation sampling and filtering + physics_upsampling=10, aafilter=('kaiser', 240, 0.3, 1.35)): + """Initialize an instrumental simulation. + + Args: + size: number of samples to generate + dt: sampling period [s] + t0: initial time [s] + orbits: path to orbit file (used for proper pseudo-ranges and derivatives therefore) + pprs: dictionary of proper pseudo-ranges, overrides values from orbit file [s] + d_pprs: dictionary of proper pseudo-range derivatives, overrides values from orbit file [s/s] + gws: dictionary of gravitational-wave responses + interpolation: interpolation function or interpolation method and parameters; when using 'lagrande', + use a tuple ('lagrange', order) with `order` the odd interpolation order; an arbitrary function should + take (x, shift [number of samples]) as parameter + laser_asds: dictionary of amplitude spectral densities for laser noise [Hz/sqrt(Hz)] + modulation_asds: dictionary of amplitude spectral densities for modulation noise [s/sqrt(Hz)] + modulation_fknees: dictionary of cutoff frequencies for modulation noise [Hz] + central_freq: laser central frequency from which all offsets are computed [Hz] + modulation_freqs: dictionary of modulation frequencies [Hz] + offsets_freqs: dictionary of laser frequency offsets [Hz] + three_lasers: whether the only simulate one laser per spacecraft + clock_asds: dictionary of clock noise amplitude spectral densities + clock_offsets: dictionary of clock offsets + clock_freqoffsets: dictionary of clock frequency offsets [s^-1] + clock_freqlindrifts: dictionary of clock frequency linear drifts [s^-2] + clock_freqquaddrifts: dictionary of clock frequency quadratic drifts [s^-2] + clockinv_tolerance: convergence tolerance for timer deviation inversion [s] + clockinv_maxiter: maximum number of iterations for timer deviation inversion + backlink_asds: dictionary of amplitude spectral densities for backlink noise [m/sqrt(Hz)] + backlink_fknees: dictionary of cutoff frequencied for backlink noise [Hz] + testmass_asds: dictionary of amplitude spectral densities for test-mass noise [ms^(-2)/sqrt(Hz)] + testmass_fknees: dictionary of cutoff frequencied for test-mass noise [Hz] + ranging_biases: dictionary of ranging noise bias [s] + ranging_asds: dictionary of ranging noise amplitude spectral densities [s/sqrt(Hz)] + physics_upsampling: ratio of sampling frequencies for physics vs. measurement simulation + aafilter: antialiasing filter function or filter design method; when using 'kaiser', use a tuple + ('kaiser', attenuation [dB], f1 [Hz], f2 [Hz]) with attenuation the attenuation in the stopband, + and f1 < f2 the frequencies defining the transition band; use `None` for no filter + """ + # pylint: disable=too-many-arguments,too-many-locals + self.git_url = 'https://gitlab.in2p3.fr/lisa-simulation/instrument' + self.version = meta.__version__ + logging.info("Initializing instrumental simulation (dt=%.2f, t0=%.2f, size=%d)", dt, t0, size) + + self.size = int(size) + self.dt = float(dt) + self.t0 = float(t0) + self.fs = 1 / self.dt + + self.laser_asds = self.mosa_dict(laser_asds) + self.modulation_asds = self.mosa_dict(modulation_asds) + self.modulation_fknees = self.mosa_dict(modulation_fknees) + self.central_freq = float(central_freq) + self.clockinv_tolerance = float(clockinv_tolerance) + self.clockinv_maxiter = int(clockinv_maxiter) + self.three_lasers = bool(three_lasers) + + self.clock_asds = self.sc_dict(clock_asds) + self.clock_offsets = self.sc_dict(clock_offsets, default=0) + self.clock_freqoffsets = self.sc_dict(clock_freqoffsets, default={ + # Default based on LISANode + '1': 5E-8, '2': 6.25E-7, '3': -3.75E-7 + }) + self.clock_freqlindrifts = self.sc_dict(clock_freqlindrifts, default={ + # Default based on LISANode + '1': 1.6E-15, '2': 2E-14, '3': -1.2E-14 + }) + self.clock_freqquaddrifts = self.sc_dict(clock_freqquaddrifts, default={ + # Default based on LISANode + '1': 9E-24, '2': 6.75E-23, '3': -1.125e-22 + }) + + self.ranging_biases = self.mosa_dict(ranging_biases) + self.ranging_asds = self.mosa_dict(ranging_asds) + + self.backlink_asds = self.mosa_dict(backlink_asds) + self.backlink_fknees = self.mosa_dict(backlink_fknees) + self.testmass_asds = self.mosa_dict(testmass_asds) + self.testmass_fknees = self.mosa_dict(testmass_fknees) + + self.physics_upsampling = int(physics_upsampling) + self.physics_size = self.size * self.physics_upsampling + self.physics_dt = self.dt / self.physics_upsampling + self.physics_fs = self.fs * self.physics_upsampling + + if interpolation is None: + logging.info("Disabling interpolation") + designed_interpolation = lambda x: x + elif callable(interpolation): + logging.info("Using user-provided interpolation function") + designed_interpolation = interpolation + else: + logging.info("Designing interpolation function") + designed_interpolation = Instrument.design_interpolation(interpolation) + self.interpolate = lambda x, shift: x if numpy.isscalar(x) else designed_interpolation(x, shift) + + if aafilter is None: + logging.info("Disabling antialiasing filter") + designed_aafilter = lambda x: x + elif callable(aafilter): + logging.info("Using user-provided antialiasing filter function") + designed_aafilter = aafilter + else: + logging.info("Design antialiasing filter (%s)", aafilter) + designed_aafilter = self.design_aafilter(aafilter) + self.aafilter = lambda x: x if numpy.isscalar(x) else designed_aafilter(x) + + self.downsampled = lambda x: x if numpy.isscalar(x) else x[::self.physics_upsampling] + + if orbits is not None: + self.orbits_path = str(orbits) + self.orbits = h5py.File(self.orbits_path, 'r') + self.pprs = self.mosa_dict(pprs, default=None) + self.d_pprs = self.mosa_dict(d_pprs, default=None) + else: + self.orbits_path = None + self.orbits = None + self.pprs = self.mosa_dict(pprs, default={ + # Default PPRs based on first samples of Keplerian orbits (v1.0) + '12': 8.3324, '23': 8.3028, '31': 8.3324, + '13': 8.3315, '32': 8.3044, '21': 8.3315, + }) + self.d_pprs = self.mosa_dict(d_pprs, default={ + # Default PPR derivatives based on first samples of Keplerian orbits (v1.0) + '12': 3.1977E-9, '23': -5.0241E-9, '31': -3.1977E-9, + '13': -3.1977E-9, '32': -5.0252E-9, '21': 3.1977E-9, + }) + + # No gravitational-wave signal by default + self.gws = self.mosa_dict(gws, default=0) + + # Default based on mission baseline of 2.4 MHz and 2.401 MHz for left and right MOSAs respect. + self.modulation_freqs = self.mosa_dict(modulation_freqs, default={ + '12': 2.4E9, '23': 2.4E9, '31': 2.4E9, + '13': 2.401E9, '32': 2.401E9, '21': 2.401E9, + }) + + # Default based on default for LISANode + self.offsets_freqs = self.mosa_dict(offsets_freqs, default={ + '12': 8.1E6, '23': 9.2E6, '31': 10.3E6, '13': 1.4E6, '32': -11.6E6, '21': -9.5E6, + }) + + @classmethod + def mosa_dict(cls, value, default=None): + """Make sure that value is turned into a dictionary with the MOSAs as keys. + + Args: + value: a scalar or a dictionary + default: default value, or None + """ + if value is None and default is not None: + value = default + elif value is None and default is None: + return None + if isinstance(value, dict): + return {mosa: float(value[mosa]) for mosa in cls.MOSAS} + return {mosa: value for mosa in cls.MOSAS} + + @classmethod + def sc_dict(cls, value, default=None): + """Make sure that value is turned into a dictionary with the spacecraft as keys. + + Args: + value: a scalar or a dictionary + default: default value, or None + """ + if value is None and default is not None: + value = default + elif value is None and default is None: + return None + if isinstance(value, dict): + return {sc: float(value[sc]) for sc in cls.SCS} + return {sc: float(value) for sc in cls.SCS} + + def disable_all_noises(self, but=None): + """Turn off all instrumental noises. + + Args: + but: optional category of noises to keep on ['laser', 'modulation', 'clock', 'pathlength', 'ranging'] + """ + if but != 'laser': + self.laser_asds = self.mosa_dict(0) + if but != 'modulation': + self.modulation_asds = self.mosa_dict(0) + if but != 'clock': + self.disable_clock_noises() + if but != 'pathlength': + self.disable_pathlength_noises() + if but != 'ranging': + self.disable_ranging_noises() + return self + + def disable_clock_noises(self): + """Turn off all imperfections on clocks. + + This includes offsets, and frequency offsets and deviations. + """ + self.clock_asds = self.sc_dict(0) + self.clock_offsets = self.sc_dict(0) + self.clock_freqoffsets = self.sc_dict(0) + self.clock_freqlindrifts = self.sc_dict(0) + self.clock_freqquaddrifts = self.sc_dict(0) + return self + + def disable_pathlength_noises(self): + """Turn off all optical pathlength noises.""" + self.backlink_asds = self.mosa_dict(0) + self.testmass_asds = self.mosa_dict(0) + return self + + def disable_ranging_noises(self): + """Turn off all pseudo-ranging noises.""" + self.ranging_biases = self.mosa_dict(0) + self.ranging_asds = self.mosa_dict(0) + return self + + def disable_dopplers(self): + """Set proper pseudo-range derivatives to zero to turn off Doppler effects.""" + self.d_pprs = self.mosa_dict(0) + return self + + @classmethod + def distant(cls, mosa): + """Return index of distant MOSA. + + In practive, we invert the indices to swap emitter and receiver. + """ + if mosa not in cls.MOSAS: + raise ValueError(f"invalid MOSA index '{mosa}'") + return f'{mosa[1]}{mosa[0]}' + + @classmethod + def adjacent(cls, mosa): + """Return index of adjacent MOSA. + + In practice, we replace the second index by the only unused spacecraft index. + """ + if mosa not in cls.MOSAS: + raise ValueError(f"invalid MOSA index '{mosa}'") + unused = [sc for sc in cls.SCS if sc not in mosa] + if len(unused) != 1: + raise RuntimeError(f"cannot find adjacent MOSA for '{mosa}'") + return f'{mosa[0]}{unused[0]}' + + @staticmethod + def design_interpolation(parameters): + """Design the interpolation function. + + We currently support Lagrande interpolation. + + Args: + parameters: see `interpolation` docstring in `__init__()` + + Returns: + A function which interpolates data. + """ + method = parameters[0] + if method == 'lagrange': + logging.debug("Using Lagrange interpolations") + order = parameters[1] + logging.debug("Lagrande interpolation order is %s", order) + return lambda x, shift: dsp.timeshift(x, shift, order=order) + + raise ValueError(f"invalid interpolation parameters '{parameters}'") + + def design_aafilter(self, parameters): + """Design the antialiasing filter. + + We currently support finite-impulse response filter designed from a Kaiser window. + The order and beta parameters of the Kaiser window are deduced from the desired attenuation + and transition bandwidth of the filter. + + Args: + parameters: see `aafilter` docstring in `__init__()` + + Returns: + A function that filters data. + """ + method = parameters[0] + nyquist = self.physics_fs / 2 + + if method == 'kaiser': + logging.debug("Designing finite-impulse response filter from Kaiser window") + attenuation, freq1, freq2 = parameters[1], parameters[2], parameters[3] + if attenuation == 0: + logging.debug("Vanishing filter attenuation, disabling filtering") + return lambda x: x + logging.debug("Filter attenuation is %s dB", attenuation) + logging.debug("Filter transition band is [%s Hz, %s Hz]", freq1, freq2) + numtaps, beta = scipy.signal.kaiserord(attenuation, (freq2 - freq1) / nyquist) + logging.debug("Kaiser window has %s taps and beta is %s", numtaps, beta) + taps = scipy.signal.firwin(numtaps, (freq1 + freq2) / (2 * nyquist), window=('kaiser', beta)) + logging.debug("Filter taps are %s", taps) + return lambda x: scipy.signal.lfilter(taps, 1, x) + + raise ValueError(f"invalid filter parameters '{parameters}'") + + def invert_timer_deviations(self, timer_deviations, sc): + """Invert timer deviations of a given spacecraft. + + We recursively solve the implicit equation dtau(tau) = dtau_hat(tau - dtau(tau)) until the + convergence criteria (tolerance) is met, or we exceed the maximum number of iterations. + + Args: + timer_deviations: array of timer deviations + sc: spacecraft index + """ + logging.info("Inverting timer deviations for spacecraft %s", sc) + logging.debug("Solving iteratively (tolerance=%s s, maxiter=%s)", + self.clockinv_tolerance, self.clockinv_maxiter) + + niter = 0 + error = 0 + next_inverse = timer_deviations + while not niter or error > self.clockinv_tolerance: + if niter >= self.clockinv_maxiter: + logging.warning("Maximum number of iterations '%s' reached (error=%.2E)", niter, error) + break + logging.debug("Starting iteration #%s", niter) + inverse = next_inverse + next_inverse = self.interpolate(timer_deviations, -inverse * self.physics_fs) + error = numpy.max(numpy.abs(inverse - next_inverse)) + logging.debug("End of iteration %s, with an error of %.2E s", niter, error) + niter += 1 + logging.debug("End of timer deviation inversion after %s iterations with an error of %.2E s", niter, error) + return inverse + + def write_metadata(self, hdf5): + """Set all properties as HDF5 attributes on `object`. + + Args: + hdf5: an HDF5 file, or a dataset + """ + for key, value in self.__dict__.items(): + try: + hdf5.attrs[key] = value + except (TypeError, RuntimeError): + try: + hdf5.attrs[key] = str(value) + except RuntimeError: + logging.warning("Cannot write metadata '%s' on '%s'", key, hdf5) + + @staticmethod + def write_dset(hdf5, dname, data): + """Write a dataset `dname` on `hdf5` for each MOSA. + + Args: + hdf5: HDF5 file in which to write + dname: dataset name or path + data: dictionary of data to write + """ + # No-op if no data + if not data: + return + # Retreive the maximum size of data + size = 1 + for datum in data.values(): + if numpy.isscalar(datum): + continue + if size != 1 and len(datum) != size: + raise ValueError(f"incompatible sizes in dictionary '{size}' and '{len(datum)}'") + size = max(size, len(datum)) + logging.debug("Writing dataset of size '%s' with '%s' columns", size, len(data)) + # Write dataset + dtype = numpy.dtype({'names': list(data.keys()), 'formats': len(data) * [numpy.float64]}) + hdf5.create_dataset(dname, (size,), dtype=dtype) + for key, datum in data.items(): + hdf5[dname][key] = datum + + def write(self, output='measurements.h5', mode='x', write_all=False): + """Run a simulation. + + Args: + output: path to measurement file + mode: measurement file opening mode + write_all: whether to write all quantities to file + """ + # pylint: disable=too-many-locals,too-many-statements,too-many-branches + + logging.info("Starting simulation") + hdf5 = h5py.File(output, mode) + + ## Physics time vector + + logging.info("Computing physics time vector (size=%s, dt=%s)", self.physics_size, self.physics_dt) + t = self.t0 + numpy.arange(self.physics_size, dtype=numpy.float64) * self.physics_dt + + logging.info("Writing metadata and physics time dataset") + self.write_metadata(hdf5) + hdf5['physics_t'] = t + + ## Laser noise + + logging.info("Generating laser noise time series") + laser_noises = { + mosa: noises.laser(self.physics_fs, self.physics_size, self.laser_asds[mosa]) + for mosa in (self.MOSAS[:3] if self.three_lasers else self.MOSAS) + } + + if self.three_lasers: + for mosa in self.MOSAS[:3]: + laser_noises[self.adjacent(mosa)] = laser_noises[mosa] + + if write_all: + logging.info("Writing laser noise to '%s'", output) + self.write_dset(hdf5, 'laser_noises', laser_noises) + + ## Clock noise + + logging.info("Generating clock noise time series") + clock_noise_offsets = { + sc: self.clock_freqoffsets[sc] + self.clock_freqlindrifts[sc] * t + self.clock_freqquaddrifts[sc] * t**2 + for sc in self.SCS + } + + for sc in self.SCS: + if self.physics_size and numpy.all(clock_noise_offsets[sc] == clock_noise_offsets[sc][0]): + clock_noise_offsets[sc] = clock_noise_offsets[sc][0] + + logging.info("Generating clock noise fluctuations time series") + clock_noise_fluctuations = { + sc: noises.clock(self.physics_fs, self.physics_size, self.clock_asds[sc]) + for sc in self.SCS + } + + if write_all: + logging.info("Writing clock noise to '%s'", output) + self.write_dset(hdf5, 'clock_noise_offsets', clock_noise_offsets) + self.write_dset(hdf5, 'clock_noise_fluctuations', clock_noise_fluctuations) + + ## Modulation noise + + logging.info("Generating modulation noise time series") + modulation_noises = { + mosa: noises.modulation(self.physics_fs, self.physics_size, + self.modulation_asds[mosa], self.modulation_fknees[mosa]) + for mosa in self.MOSAS + } + + if write_all: + logging.info("Writing modulation noise to '%s'", output) + self.write_dset(hdf5, 'modulation_noises', modulation_noises) + + ## Local beams + + local_carrier_offsets = self.offsets_freqs + + logging.info("Computing carrier frequency fluctuations for local beams") + local_carrier_fluctuations = laser_noises + + logging.info("Computing upper sideband frequency offsets for local beams") + local_usb_offsets = { + mosa: self.offsets_freqs[mosa] + self.modulation_freqs[mosa] * (1 + clock_noise_offsets[mosa[0]]) # 2.67 + for mosa in self.MOSAS + } + + logging.info("Computing upper sideband frequency fluctuations for local beams") + local_usb_fluctuations = { + mosa: laser_noises[mosa] \ + + self.modulation_freqs[mosa] * (clock_noise_fluctuations[mosa[0]] + modulation_noises[mosa]) # 2.68 + for mosa in self.MOSAS + } + + logging.info("Computing local timer deviations") + local_timer_deviations = { + sc: self.clock_offsets[sc] + numpy.cumsum( + numpy.broadcast_to(clock_noise_offsets[sc] + clock_noise_fluctuations[sc], \ + (self.physics_size,)) * self.dt) + for sc in self.SCS + } + + if write_all: + logging.info("Writing local beams to '%s'", output) + self.write_dset(hdf5, 'local_carrier_offsets', local_carrier_offsets) + self.write_dset(hdf5, 'local_carrier_fluctuations', local_carrier_fluctuations) + self.write_dset(hdf5, 'local_usb_offsets', local_usb_offsets) + self.write_dset(hdf5, 'local_usb_fluctuations', local_usb_fluctuations) + self.write_dset(hdf5, 'local_timer_deviations', local_timer_deviations) + + ## Proper pseudo-ranges + if self.pprs is None: + + logging.info("Interpolating proper pseudo-ranges from orbit file '%s'", self.orbits_path) + self.pprs = { + mosa: scipy.interpolate.InterpolatedUnivariateSpline( + self.orbits['tcb']['t'][:], self.orbits['tcb'][f'l_{mosa}']['ppr'], k=1, ext='raise')(t) + for mosa in self.MOSAS + } + + if write_all: + logging.info("Writing proper pseudo-ranges to '%s'", output) + self.write_dset(hdf5, 'pprs', self.pprs) + + ## Proper pseudo-range derivatives + if self.d_pprs is None: + + logging.info("Interpolating proper pseudo-range derivatives from orbit file '%s'", self.orbits_path) + self.d_pprs = { + mosa: scipy.interpolate.InterpolatedUnivariateSpline( + self.orbits['tcb']['t'][:], self.orbits['tcb'][f'l_{mosa}']['d_ppr'], k=1, ext='raise')(t) + for mosa in self.MOSAS + } + + if write_all: + logging.info("Writing proper pseudo-range derivatives to '%s'", output) + self.write_dset(hdf5, 'd_pprs', self.d_pprs) + + ## Propagation to distant MOSA + + logging.info("Propagating carrier frequency offsets to distant MOSA") + distant_carrier_offsets = { + mosa: (1 - self.d_pprs[mosa]) * \ + self.interpolate(local_carrier_offsets[self.distant(mosa)], -self.pprs[mosa] * self.physics_fs) \ + - self.d_pprs[mosa] * self.central_freq + for mosa in self.MOSAS + } + + logging.info("Propagating carrier frequency fluctuations to distant MOSA") + distant_carrier_fluctuations = { + mosa: self.gws[mosa] + (1 - self.d_pprs[mosa]) * \ + self.interpolate(local_carrier_fluctuations[self.distant(mosa)], -self.pprs[mosa] * self.physics_fs) + for mosa in self.MOSAS + } + + logging.info("Propagating upper sideband frequency offsets to distant MOSA") + distant_usb_offsets = { + mosa: (1 - self.d_pprs[mosa]) * \ + self.interpolate(local_usb_offsets[self.distant(mosa)], -self.pprs[mosa] * self.physics_fs) \ + - self.d_pprs[mosa] * self.central_freq + for mosa in self.MOSAS + } + + logging.info("Propagating upper sideband frequency fluctuations to distant MOSA") + distant_usb_fluctuations = { + mosa: self.gws[mosa] + (1 - self.d_pprs[mosa]) * \ + self.interpolate(local_usb_fluctuations[self.distant(mosa)], -self.pprs[mosa] * self.physics_fs) + for mosa in self.MOSAS + } + + logging.info("Propagating timer deviations to distant MOSA") + distant_timer_deviations = { + mosa: self.interpolate(local_timer_deviations[mosa[1]], -self.pprs[mosa] * self.physics_fs) \ + - self.pprs[mosa] + for mosa in self.MOSAS + } + + if write_all: + logging.info("Writing propagated distant beams to '%s'", output) + self.write_dset(hdf5, 'distant_carrier_offsets', distant_carrier_offsets) + self.write_dset(hdf5, 'distant_carrier_fluctuations', distant_carrier_fluctuations) + self.write_dset(hdf5, 'distant_usb_offsets', distant_usb_offsets) + self.write_dset(hdf5, 'distant_usb_fluctuations', distant_usb_fluctuations) + self.write_dset(hdf5, 'distant_timer_deviations', distant_timer_deviations) + + ## Backlink noise + + logging.info("Generating backlink noise time series") + backlink_noises = { + mosa: noises.backlink(self.physics_fs, self.physics_size, + self.backlink_asds[mosa], self.backlink_fknees[mosa]) + for mosa in self.MOSAS + } + + if write_all: + logging.info("Writing backlink noise to '%s'", output) + self.write_dset(hdf5, 'backlink_noises', backlink_noises) + + ## Propagation to adjacent MOSA + + logging.info("Propagating carrier frequency offsets to adjacent MOSA") + adjacent_carrier_offsets = { + mosa: local_carrier_offsets[self.adjacent(mosa)] + for mosa in self.MOSAS + } + + logging.info("Propagating carrier frequency fluctuations to adjacent MOSA") + adjacent_carrier_fluctuations = { + mosa: local_carrier_fluctuations[self.adjacent(mosa)] + self.central_freq * backlink_noises[mosa] + for mosa in self.MOSAS + } + + logging.info("Propagating upper sideband frequency offsets to adjacent MOSA") + adjacent_usb_offsets = { + mosa: local_usb_offsets[self.adjacent(mosa)] + for mosa in self.MOSAS + } + + logging.info("Propagating upper sideband frequency fluctuations to adjacent MOSA") + adjacent_usb_fluctuations = { + mosa: local_usb_fluctuations[self.adjacent(mosa)] + self.central_freq * backlink_noises[mosa] + for mosa in self.MOSAS + } + + if write_all: + logging.info("Writing propagated adjacent beams to '%s'", output) + self.write_dset(hdf5, 'adjacent_carrier_offsets', adjacent_carrier_offsets) + self.write_dset(hdf5, 'adjacent_carrier_fluctuations', adjacent_carrier_fluctuations) + self.write_dset(hdf5, 'adjacent_usb_offsets', adjacent_usb_offsets) + self.write_dset(hdf5, 'adjacent_usb_fluctuations', adjacent_usb_fluctuations) + + ## Inter-spacecraft interferometer local beams + + logging.info("Propagating local carrier frequency offsets to inter-spacecraft interferometer") + local_isc_carrier_offsets = local_carrier_offsets + + logging.info("Propagating local carrier frequency fluctuations to inter-spacecraft interferometer") + local_isc_carrier_fluctuations = local_carrier_fluctuations + + logging.info("Propagating local upper sideband frequency offsets to inter-spacecraft interferometer") + local_isc_usb_offsets = local_usb_offsets + + logging.info("Propagating local upper sideband frequency fluctuations to inter-spacecraft interferometer") + local_isc_usb_fluctuations = local_usb_fluctuations + + if write_all: + logging.info("Writing local beams at inter-spacecraft interferometer to '%s'", output) + self.write_dset(hdf5, 'local_isc_carrier_offsets', local_isc_carrier_offsets) + self.write_dset(hdf5, 'local_isc_carrier_fluctuations', local_isc_carrier_fluctuations) + self.write_dset(hdf5, 'local_isc_usb_offsets', local_isc_usb_offsets) + self.write_dset(hdf5, 'local_isc_usb_fluctuations', local_isc_usb_fluctuations) + + ## Inter-spacecraft interferometer distant beams + + logging.info("Propagating distant carrier frequency offsets to inter-spacecraft interferometer") + distant_isc_carrier_offsets = distant_carrier_offsets + + logging.info("Propagating distant carrier frequency fluctuations to inter-spacecraft interferometer") + distant_isc_carrier_fluctuations = distant_carrier_fluctuations + + logging.info("Propagating distant upper sideband frequency offsets to inter-spacecraft interferometer") + distant_isc_usb_offsets = distant_usb_offsets + + logging.info("Propagating distant upper sideband frequency fluctuations to inter-spacecraft interferometer") + distant_isc_usb_fluctuations = distant_usb_fluctuations + + if write_all: + logging.info("Writing distant beams at inter-spacecraft interferometer to '%s'", output) + self.write_dset(hdf5, 'distant_isc_carrier_offsets', distant_isc_carrier_offsets) + self.write_dset(hdf5, 'distant_isc_carrier_fluctuations', distant_isc_carrier_fluctuations) + self.write_dset(hdf5, 'distant_isc_usb_offsets', distant_isc_usb_offsets) + self.write_dset(hdf5, 'distant_isc_usb_fluctuations', distant_isc_usb_fluctuations) + + ## Inter-spacecraft interferometer beatnotes on TPS (high-frequency) + + logging.info("Computing inter-spacecraft carrier beatnote frequency offsets on TPS") + tps_isc_carrier_offsets = { + mosa: distant_isc_carrier_offsets[mosa] - local_isc_carrier_offsets[mosa] + for mosa in self.MOSAS + } + + logging.info("Computing inter-spacecraft carrier beatnote frequency fluctuations on TPS") + tps_isc_carrier_fluctuations = { + mosa: distant_isc_carrier_fluctuations[mosa] - local_isc_carrier_fluctuations[mosa] + for mosa in self.MOSAS + } + + logging.info("Computing inter-spacecraft upper sideband beatnote frequency offsets on TPS") + tps_isc_usb_offsets = { + mosa: distant_isc_usb_offsets[mosa] - local_isc_usb_offsets[mosa] + for mosa in self.MOSAS + } + + logging.info("Computing inter-spacecraft upper sideband beatnote frequency fluctuations on TPS") + tps_isc_usb_fluctuations = { + mosa: distant_isc_usb_fluctuations[mosa] - local_isc_usb_fluctuations[mosa] + for mosa in self.MOSAS + } + + if write_all: + logging.info("Writing inter-spacecraft beatnotes on TPS to '%s'", output) + self.write_dset(hdf5, 'tps_isc_carrier_offsets', tps_isc_carrier_offsets) + self.write_dset(hdf5, 'tps_isc_carrier_fluctuations', tps_isc_carrier_fluctuations) + self.write_dset(hdf5, 'tps_isc_usb_offsets', tps_isc_usb_offsets) + self.write_dset(hdf5, 'tps_isc_usb_fluctuations', tps_isc_usb_fluctuations) + + ## Ranging noise + + logging.info("Generating ranging noise time series") + ranging_noises = { + mosa: self.ranging_biases[mosa] \ + + noises.ranging(self.physics_fs, self.physics_size, self.ranging_asds[mosa]) + for mosa in self.MOSAS + } + + if write_all: + logging.info("Writing ranging noise to '%s'", output) + self.write_dset(hdf5, 'ranging_noises', ranging_noises) + + ## Measured pseudo-ranging on TPS grid (high-frequency) + + logging.info("Computing measured pseudo-ranges on TPS") + tps_mprs = { + mosa: local_timer_deviations[mosa[0]] - distant_timer_deviations[mosa] + ranging_noises[mosa] + for mosa in self.MOSAS + } + + if write_all: + logging.info("Writing measured pseudo-ranges on TPS to '%s'", output) + self.write_dset(hdf5, 'tps_mprs', tps_mprs) + + ## Test-mass acceleration noise + + logging.info("Generating test-mass acceleration noise time series") + testmass_noises = { + mosa: noises.testmass(self.physics_fs, self.physics_size, + self.testmass_asds[mosa], self.testmass_fknees[mosa]) + for mosa in self.MOSAS + } + + if write_all: + logging.info("Writing test-mass acceleration noise to '%s'", output) + self.write_dset(hdf5, 'tm_noises', testmass_noises) + + ## Test-mass interferometer local beams + + logging.info("Propagating local carrier frequency offsets to test-mass interferometer") + local_tm_carrier_offsets = local_carrier_offsets + + logging.info("Propagating local carrier frequency fluctuations to test-mass interferometer") + local_tm_carrier_fluctuations = { + mosa: local_carrier_fluctuations[mosa] + self.central_freq * testmass_noises[mosa] + for mosa in self.MOSAS + } + + logging.info("Propagating local upper sideband frequency offsets to test-mass interferometer") + local_tm_usb_offsets = local_usb_offsets + + logging.info("Propagating local upper sideband frequency fluctuations to test-mass interferometer") + local_tm_usb_fluctuations = { + mosa: local_usb_fluctuations[mosa] + self.central_freq * testmass_noises[mosa] + for mosa in self.MOSAS + } + + if write_all: + logging.info("Writing local beams at test-mass interferometer to '%s'", output) + self.write_dset(hdf5, 'local_tm_carrier_offsets', local_tm_carrier_offsets) + self.write_dset(hdf5, 'local_tm_carrier_fluctuations', local_tm_carrier_fluctuations) + self.write_dset(hdf5, 'local_tm_usb_offsets', local_tm_usb_offsets) + self.write_dset(hdf5, 'local_tm_usb_fluctuations', local_tm_usb_fluctuations) + + ## Test-mass interferometer adjacent beams + + logging.info("Propagating adjacent carrier frequency offsets to test-mass interferometer") + adjacent_tm_carrier_offsets = adjacent_carrier_offsets + + logging.info("Propagating adjacent carrier frequency fluctuations to test-mass interferometer") + adjacent_tm_carrier_fluctuations = adjacent_carrier_fluctuations + + logging.info("Propagating adjacent upper sideband frequency offsets to test-mass interferometer") + adjacent_tm_usb_offsets = adjacent_usb_offsets + + logging.info("Propagating adjacent upper sideband frequency fluctuations to test-mass interferometer") + adjacent_tm_usb_fluctuations = adjacent_usb_fluctuations + + if write_all: + logging.info("Writing adjacent beams at test-mass interferometer to '%s'", output) + self.write_dset(hdf5, 'adjacent_tm_carrier_offsets', adjacent_tm_carrier_offsets) + self.write_dset(hdf5, 'adjacent_tm_carrier_fluctuations', adjacent_tm_carrier_fluctuations) + self.write_dset(hdf5, 'adjacent_tm_usb_offsets', adjacent_tm_usb_offsets) + self.write_dset(hdf5, 'adjacent_tm_usb_fluctuations', adjacent_tm_usb_fluctuations) + + ## Test-mass interferometer beatnotes on TPS (high-frequency) + + logging.info("Computing test-mass carrier beatnote frequency offsets on TPS") + tps_tm_carrier_offsets = { + mosa: adjacent_tm_carrier_offsets[mosa] - local_tm_carrier_offsets[mosa] + for mosa in self.MOSAS + } + + logging.info("Computing test-mass carrier beatnote frequency fluctuations on TPS") + tps_tm_carrier_fluctuations = { + mosa: adjacent_tm_carrier_fluctuations[mosa] - local_tm_carrier_fluctuations[mosa] + for mosa in self.MOSAS + } + + logging.info("Computing test-mass upper sideband beatnote frequency offsets on TPS") + tps_tm_usb_offsets = { + mosa: adjacent_tm_usb_offsets[mosa] - local_tm_usb_offsets[mosa] + for mosa in self.MOSAS + } + + logging.info("Computing test-mass upper sideband beatnote frequency fluctuations on TPS") + tps_tm_usb_fluctuations = { + mosa: adjacent_tm_usb_fluctuations[mosa] - local_tm_usb_fluctuations[mosa] + for mosa in self.MOSAS + } + + if write_all: + logging.info("Writing test-mass beatnotes on TPS to '%s'", output) + self.write_dset(hdf5, 'tps_tm_carrier_offsets', tps_tm_carrier_offsets) + self.write_dset(hdf5, 'tps_tm_carrier_fluctuations', tps_tm_carrier_fluctuations) + self.write_dset(hdf5, 'tps_tm_usb_offsets', tps_tm_usb_offsets) + self.write_dset(hdf5, 'tps_tm_usb_fluctuations', tps_tm_usb_fluctuations) + + ## Reference interferometer local beams + + logging.info("Propagating local carrier frequency offsets to reference interferometer") + local_ref_carrier_offsets = local_carrier_offsets + + logging.info("Propagating local carrier frequency fluctuations to reference interferometer") + local_ref_carrier_fluctuations = local_carrier_fluctuations + + logging.info("Propagating local upper sideband frequency offsets to reference interferometer") + local_ref_usb_offsets = local_usb_offsets + + logging.info("Propagating local upper sideband frequency fluctuations to reference interferometer") + local_ref_usb_fluctuations = local_usb_fluctuations + + if write_all: + logging.info("Writing local beams at reference interferometer to '%s'", output) + self.write_dset(hdf5, 'local_ref_carrier_offsets', local_ref_carrier_offsets) + self.write_dset(hdf5, 'local_ref_carrier_fluctuations', local_ref_carrier_fluctuations) + self.write_dset(hdf5, 'local_ref_usb_offsets', local_ref_usb_offsets) + self.write_dset(hdf5, 'local_ref_usb_fluctuations', local_ref_usb_fluctuations) + + ## Reference interferometer adjacent beams + + logging.info("Propagating adjacent carrier frequency offsets to reference interferometer") + adjacent_ref_carrier_offsets = adjacent_carrier_offsets + + logging.info("Propagating adjacent carrier frequency fluctuations to reference interferometer") + adjacent_ref_carrier_fluctuations = adjacent_carrier_fluctuations + + logging.info("Propagating adjacent upper sideband frequency offsets to reference interferometer") + adjacent_ref_usb_offsets = adjacent_usb_offsets + + logging.info("Propagating adjacent upper sideband frequency fluctuations to reference interferometer") + adjacent_ref_usb_fluctuations = adjacent_usb_fluctuations + + if write_all: + logging.info("Writing adjacent beams at reference interferometer to '%s'", output) + self.write_dset(hdf5, 'adjacent_ref_carrier_offsets', adjacent_ref_carrier_offsets) + self.write_dset(hdf5, 'adjacent_ref_carrier_fluctuations', adjacent_ref_carrier_fluctuations) + self.write_dset(hdf5, 'adjacent_ref_usb_offsets', adjacent_ref_usb_offsets) + self.write_dset(hdf5, 'adjacent_ref_usb_fluctuations', adjacent_ref_usb_fluctuations) + + ## Reference interferometer beatnotes on TPS (high-frequency) + + logging.info("Computing reference carrier beatnote frequency offsets on TPS") + tps_ref_carrier_offsets = { + mosa: adjacent_ref_carrier_offsets[mosa] - local_ref_carrier_offsets[mosa] + for mosa in self.MOSAS + } + + logging.info("Computing reference carrier beatnote frequency fluctuations on TPS") + tps_ref_carrier_fluctuations = { + mosa: adjacent_ref_carrier_fluctuations[mosa] - local_ref_carrier_fluctuations[mosa] + for mosa in self.MOSAS + } + + logging.info("Computing reference upper sideband beatnote frequency offsets on TPS") + tps_ref_usb_offsets = { + mosa: adjacent_ref_usb_offsets[mosa] - local_ref_usb_offsets[mosa] + for mosa in self.MOSAS + } + + logging.info("Computing reference upper sideband beatnote frequency fluctuations on TPS") + tps_ref_usb_fluctuations = { + mosa: adjacent_ref_usb_fluctuations[mosa] - local_ref_usb_fluctuations[mosa] + for mosa in self.MOSAS + } + + if write_all: + logging.info("Writing reference beatnotes on TPS to '%s'", output) + self.write_dset(hdf5, 'tps_ref_carrier_offsets', tps_ref_carrier_offsets) + self.write_dset(hdf5, 'tps_ref_carrier_fluctuations', tps_ref_carrier_fluctuations) + self.write_dset(hdf5, 'tps_ref_usb_offsets', tps_ref_usb_offsets) + self.write_dset(hdf5, 'tps_ref_usb_fluctuations', tps_ref_usb_fluctuations) + + ## Sampling beatnotes and measured pseudo-ranges to THE grid + + inverse_timer_deviations = { + sc: self.invert_timer_deviations(local_timer_deviations[sc], sc) + for sc in self.SCS + } + + timestamp = lambda x, sc: self.interpolate(x, -inverse_timer_deviations[sc] * self.physics_fs) + + logging.info("Sampling inter-spacecraft carrier beatnote frequency offsets to THE grid") + the_isc_carrier_offsets = { + mosa: timestamp(tps_isc_carrier_offsets[mosa] / (1 + clock_noise_offsets[mosa[0]]), mosa[0]) + for mosa in self.MOSAS + } + logging.info("Sampling inter-spacecraft carrier beatnote frequency fluctuations to THE grid") + the_isc_carrier_fluctuations = { + mosa: timestamp(tps_isc_carrier_fluctuations[mosa] / (1 + clock_noise_offsets[mosa[0]]) + - tps_isc_carrier_offsets[mosa] * clock_noise_fluctuations[mosa[0]] + / (1 + clock_noise_offsets[mosa[0]])**2, mosa[0]) + for mosa in self.MOSAS + } + + logging.info("Sampling inter-spacecraft upper sideband beatnote frequency offsets to THE grid") + the_isc_usb_offsets = { + mosa: timestamp(tps_isc_usb_offsets[mosa] / (1 + clock_noise_offsets[mosa[0]]), mosa[0]) + for mosa in self.MOSAS + } + + logging.info("Sampling inter-spacecraft upper sideband beatnote frequency fluctuations to THE grid") + the_isc_usb_fluctuations = { + mosa: timestamp(tps_isc_usb_fluctuations[mosa] / (1 + clock_noise_offsets[mosa[0]]) + - tps_isc_usb_offsets[mosa] * clock_noise_fluctuations[mosa[0]] + / (1 + clock_noise_offsets[mosa[0]])**2, mosa[0]) + for mosa in self.MOSAS + } + + if write_all: + logging.info("Writing inter-spacecraft beatnotes sampled to THE grid to '%s'", output) + self.write_dset(hdf5, 'the_isc_carrier_offsets', the_isc_carrier_offsets) + self.write_dset(hdf5, 'the_isc_carrier_fluctuations', the_isc_carrier_fluctuations) + self.write_dset(hdf5, 'the_isc_usb_offsets', the_isc_usb_offsets) + self.write_dset(hdf5, 'the_isc_usb_fluctuations', the_isc_usb_fluctuations) + + logging.info("Sampling measured pseudo-ranges to THE grid") + the_mprs = { + mosa: timestamp(tps_mprs[mosa], mosa[0]) + for mosa in self.MOSAS + } + + if write_all: + logging.info("Writing measured pseudo-ranges sampled to THE grid to '%s'", output) + self.write_dset(hdf5, 'the_mprs', the_mprs) + + logging.info("Sampling test-mass beatnotes to THE grid") + the_tm_carrier_offsets = { + mosa: timestamp(tps_tm_carrier_offsets[mosa] / (1 + clock_noise_offsets[mosa[0]]), mosa[0]) + for mosa in self.MOSAS + } + the_tm_carrier_fluctuations = { + mosa: timestamp(tps_tm_carrier_fluctuations[mosa] / (1 + clock_noise_offsets[mosa[0]]) + - tps_tm_carrier_offsets[mosa] * clock_noise_fluctuations[mosa[0]] + / (1 + clock_noise_offsets[mosa[0]])**2, mosa[0]) + for mosa in self.MOSAS + } + the_tm_usb_offsets = { + mosa: timestamp(tps_tm_usb_offsets[mosa] / (1 + clock_noise_offsets[mosa[0]]), mosa[0]) + for mosa in self.MOSAS + } + the_tm_usb_fluctuations = { + mosa: timestamp(tps_tm_usb_fluctuations[mosa] / (1 + clock_noise_offsets[mosa[0]]) + - tps_tm_usb_offsets[mosa] * clock_noise_fluctuations[mosa[0]] + / (1 + clock_noise_offsets[mosa[0]])**2, mosa[0]) + for mosa in self.MOSAS + } + + if write_all: + logging.info("Writing test-mass beatnotes sampled to THE grid to '%s'", output) + self.write_dset(hdf5, 'the_tm_carrier_offsets', the_tm_carrier_offsets) + self.write_dset(hdf5, 'the_tm_carrier_fluctuations', the_tm_carrier_fluctuations) + self.write_dset(hdf5, 'the_tm_usb_offsets', the_tm_usb_offsets) + self.write_dset(hdf5, 'the_tm_usb_fluctuations', the_tm_usb_fluctuations) + + logging.info("Sampling reference beatnotes to THE grid") + the_ref_carrier_offsets = { + mosa: timestamp(tps_ref_carrier_offsets[mosa] / (1 + clock_noise_offsets[mosa[0]]), mosa[0]) + for mosa in self.MOSAS + } + the_ref_carrier_fluctuations = { + mosa: timestamp(tps_ref_carrier_fluctuations[mosa] / (1 + clock_noise_offsets[mosa[0]]) + - tps_ref_carrier_offsets[mosa] * clock_noise_fluctuations[mosa[0]] + / (1 + clock_noise_offsets[mosa[0]])**2, mosa[0]) + for mosa in self.MOSAS + } + the_ref_usb_offsets = { + mosa: timestamp(tps_ref_usb_offsets[mosa] / (1 + clock_noise_offsets[mosa[0]]), mosa[0]) + for mosa in self.MOSAS + } + the_ref_usb_fluctuations = { + mosa: timestamp(tps_ref_usb_fluctuations[mosa] / (1 + clock_noise_offsets[mosa[0]]) + - tps_ref_usb_offsets[mosa] * clock_noise_fluctuations[mosa[0]] + / (1 + clock_noise_offsets[mosa[0]])**2, mosa[0]) + for mosa in self.MOSAS + } + + if write_all: + logging.info("Writing reference beatnotes sampled to THE grid to '%s'", output) + self.write_dset(hdf5, 'the_ref_carrier_offsets', the_ref_carrier_offsets) + self.write_dset(hdf5, 'the_ref_carrier_fluctuations', the_ref_carrier_fluctuations) + self.write_dset(hdf5, 'the_ref_usb_offsets', the_ref_usb_offsets) + self.write_dset(hdf5, 'the_ref_usb_fluctuations', the_ref_usb_fluctuations) + + ## Measurement time vector + + logging.info("Computing measurement time vector (size=%s, dt=%s)", self.size, self.dt) + logging.info("Writing measurement time dataset to '%s'", output) + hdf5['t'] = t[::self.physics_upsampling] + + ## Antialiasing filtering + + logging.info("Filtering inter-spacecraft beatnotes") + filtered_isc_carrier_offsets = { + mosa: self.aafilter(the_isc_carrier_offsets[mosa]) + for mosa in self.MOSAS + } + filtered_isc_carrier_fluctuations = { + mosa: self.aafilter(the_isc_carrier_fluctuations[mosa]) + for mosa in self.MOSAS + } + filtered_isc_usb_offsets = { + mosa: self.aafilter(the_isc_usb_offsets[mosa]) + for mosa in self.MOSAS + } + filtered_isc_usb_fluctuations = { + mosa: self.aafilter(the_isc_usb_fluctuations[mosa]) + for mosa in self.MOSAS + } + + if write_all: + logging.info("Writing filtered inter-spacecraft beatnotes to '%s'", output) + self.write_dset(hdf5, 'filtered_isc_carrier_offsets', filtered_isc_carrier_offsets) + self.write_dset(hdf5, 'filtered_isc_carrier_fluctuations', filtered_isc_carrier_fluctuations) + self.write_dset(hdf5, 'filtered_isc_usb_offsets', filtered_isc_usb_offsets) + self.write_dset(hdf5, 'filtered_isc_usb_fluctuations', filtered_isc_usb_fluctuations) + + logging.info("Filtering measured pseudo-ranges") + filtered_mprs = { + mosa: self.aafilter(the_mprs[mosa]) + for mosa in self.MOSAS + } + + if write_all: + logging.info("Writing filtered measured pseudo-ranges to '%s'", output) + self.write_dset(hdf5, 'filtered_mprs', filtered_mprs) + + logging.info("Filtering test-mass beatnotes") + filtered_tm_carrier_offsets = { + mosa: self.aafilter(the_tm_carrier_offsets[mosa]) + for mosa in self.MOSAS + } + filtered_tm_carrier_fluctuations = { + mosa: self.aafilter(the_tm_carrier_fluctuations[mosa]) + for mosa in self.MOSAS + } + filtered_tm_usb_offsets = { + mosa: self.aafilter(the_tm_usb_offsets[mosa]) + for mosa in self.MOSAS + } + filtered_tm_usb_fluctuations = { + mosa: self.aafilter(the_tm_usb_fluctuations[mosa]) + for mosa in self.MOSAS + } + + if write_all: + logging.info("Writing filtered test-mass beatnotes to '%s'", output) + self.write_dset(hdf5, 'filtered_tm_carrier_offsets', filtered_tm_carrier_offsets) + self.write_dset(hdf5, 'filtered_tm_carrier_fluctuations', filtered_tm_carrier_fluctuations) + self.write_dset(hdf5, 'filtered_tm_usb_offsets', filtered_tm_usb_offsets) + self.write_dset(hdf5, 'filtered_tm_usb_fluctuations', filtered_tm_usb_fluctuations) + + logging.info("Filtering reference beatnotes") + filtered_ref_carrier_offsets = { + mosa: self.aafilter(the_ref_carrier_offsets[mosa]) + for mosa in self.MOSAS + } + filtered_ref_carrier_fluctuations = { + mosa: self.aafilter(the_ref_carrier_fluctuations[mosa]) + for mosa in self.MOSAS + } + filtered_ref_usb_offsets = { + mosa: self.aafilter(the_ref_usb_offsets[mosa]) + for mosa in self.MOSAS + } + filtered_ref_usb_fluctuations = { + mosa: self.aafilter(the_ref_usb_fluctuations[mosa]) + for mosa in self.MOSAS + } + + if write_all: + logging.info("Writing filtered reference beatnotes to '%s'", output) + self.write_dset(hdf5, 'filtered_ref_carrier_offsets', filtered_ref_carrier_offsets) + self.write_dset(hdf5, 'filtered_ref_carrier_fluctuations', filtered_ref_carrier_fluctuations) + self.write_dset(hdf5, 'filtered_ref_usb_offsets', filtered_ref_usb_offsets) + self.write_dset(hdf5, 'filtered_ref_usb_fluctuations', filtered_ref_usb_fluctuations) + + ## Downsampling filtering + + logging.info("Downsampling inter-spacecraft beatnotes") + isc_carrier_offsets = { + mosa: self.downsampled(filtered_isc_carrier_offsets[mosa]) + for mosa in self.MOSAS + } + isc_carrier_fluctuations = { + mosa: self.downsampled(filtered_isc_carrier_fluctuations[mosa]) + for mosa in self.MOSAS + } + isc_usb_offsets = { + mosa: self.downsampled(filtered_isc_usb_offsets[mosa]) + for mosa in self.MOSAS + } + isc_usb_fluctuations = { + mosa: self.downsampled(filtered_isc_usb_fluctuations[mosa]) + for mosa in self.MOSAS + } + + logging.info("Writing downsampled inter-spacecraft beatnotes to '%s'", output) + self.write_dset(hdf5, 'isc_carrier_offsets', isc_carrier_offsets) + self.write_dset(hdf5, 'isc_carrier_fluctuations', isc_carrier_fluctuations) + self.write_dset(hdf5, 'isc_usb_offsets', isc_usb_offsets) + self.write_dset(hdf5, 'isc_usb_fluctuations', isc_usb_fluctuations) + + logging.info("Filtering and downsampling measured pseudo-ranges") + mprs = { + mosa: self.downsampled(filtered_mprs[mosa]) + for mosa in self.MOSAS + } + + logging.info("Writing downsampled measured pseudo-ranges to '%s'", output) + self.write_dset(hdf5, 'mprs', mprs) + + logging.info("Filtering and downsampling test-mass beatnotes") + tm_carrier_offsets = { + mosa: self.downsampled(filtered_tm_carrier_offsets[mosa]) + for mosa in self.MOSAS + } + tm_carrier_fluctuations = { + mosa: self.downsampled(filtered_tm_carrier_fluctuations[mosa]) + for mosa in self.MOSAS + } + tm_usb_offsets = { + mosa: self.downsampled(filtered_tm_usb_offsets[mosa]) + for mosa in self.MOSAS + } + tm_usb_fluctuations = { + mosa: self.downsampled(filtered_tm_usb_fluctuations[mosa]) + for mosa in self.MOSAS + } + + logging.info("Writing downsampled test-mass beatnotes to '%s'", output) + self.write_dset(hdf5, 'tm_carrier_offsets', tm_carrier_offsets) + self.write_dset(hdf5, 'tm_carrier_fluctuations', tm_carrier_fluctuations) + self.write_dset(hdf5, 'tm_usb_offsets', tm_usb_offsets) + self.write_dset(hdf5, 'tm_usb_fluctuations', tm_usb_fluctuations) + + logging.info("Filtering and downsampling reference beatnotes") + ref_carrier_offsets = { + mosa: self.downsampled(filtered_ref_carrier_offsets[mosa]) + for mosa in self.MOSAS + } + ref_carrier_fluctuations = { + mosa: self.downsampled(filtered_ref_carrier_fluctuations[mosa]) + for mosa in self.MOSAS + } + ref_usb_offsets = { + mosa: self.downsampled(filtered_ref_usb_offsets[mosa]) + for mosa in self.MOSAS + } + ref_usb_fluctuations = { + mosa: self.downsampled(filtered_ref_usb_fluctuations[mosa]) + for mosa in self.MOSAS + } + + logging.info("Writing downsampled reference beatnotes to '%s'", output) + self.write_dset(hdf5, 'ref_carrier_offsets', ref_carrier_offsets) + self.write_dset(hdf5, 'ref_carrier_fluctuations', ref_carrier_fluctuations) + self.write_dset(hdf5, 'ref_usb_offsets', ref_usb_offsets) + self.write_dset(hdf5, 'ref_usb_fluctuations', ref_usb_fluctuations) + + ## Total frequencies + + logging.info("Computing inter-spacecraft beatnote total frequencies") + isc_carriers = { + mosa: isc_carrier_offsets[mosa] + isc_carrier_fluctuations[mosa] + for mosa in self.MOSAS + } + isc_usbs = { + mosa: isc_usb_offsets[mosa] + isc_usb_fluctuations[mosa] + for mosa in self.MOSAS + } + + logging.info("Writing inter-spacecraft beatnote total frequencies to '%s'", output) + self.write_dset(hdf5, 'isc_carriers', isc_carriers) + self.write_dset(hdf5, 'isc_usbs', isc_usbs) + + logging.info("Computing test-mass beatnote total frequencies") + tm_carriers = { + mosa: tm_carrier_offsets[mosa] + tm_carrier_fluctuations[mosa] + for mosa in self.MOSAS + } + tm_usbs = { + mosa: tm_usb_offsets[mosa] + tm_usb_fluctuations[mosa] + for mosa in self.MOSAS + } + + logging.info("Writing test-mass beatnote total frequencies to '%s'", output) + self.write_dset(hdf5, 'tm_carriers', tm_carriers) + self.write_dset(hdf5, 'tm_usbs', tm_usbs) + + logging.info("Computing reference beatnote total frequencies") + ref_carriers = { + mosa: ref_carrier_offsets[mosa] + ref_carrier_fluctuations[mosa] + for mosa in self.MOSAS + } + ref_usbs = { + mosa: ref_usb_offsets[mosa] + ref_usb_fluctuations[mosa] + for mosa in self.MOSAS + } + + logging.info("Writing reference beatnote total frequencies to '%s'", output) + self.write_dset(hdf5, 'ref_carriers', ref_carriers) + self.write_dset(hdf5, 'ref_usbs', ref_usbs) + + ## Closing simulation + + logging.info("Closing output file '%s'", output) + hdf5.close() + logging.info("Simulation complete") diff --git a/requirements.txt b/requirements.txt index 1fc1774..a3df11e 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,14 +1,25 @@ astroid==2.4.2 +attrs==20.3.0 +cycler==0.10.0 h5py==3.1.0 +iniconfig==1.1.1 isort==5.7.0 +kiwisolver==1.3.1 lazy-object-proxy==1.4.3 +lisaconstants==0.0.3 +matplotlib==3.3.3 mccabe==0.6.1 numpy==1.19.5 +packaging==20.8 +Pillow==8.1.0 +pluggy==0.13.1 +py==1.10.0 pylint==2.6.0 +pyparsing==2.4.7 pyplnoise==1.3 +pytest==6.2.1 +python-dateutil==2.8.1 scipy==1.6.0 six==1.15.0 toml==0.10.2 wrapt==1.12.1 - -git+https://lisainstrument-ci:kSRBLBj7s3R7zfgnQiZx@gitlab.in2p3.fr/lisa-simulation/python-constants.git@v0.0.3 -- GitLab