diff --git a/CHANGELOG.md b/CHANGELOG.md index 5bf71af4c39ce868f2627bcc57b343e833d43294..97ff0e8ebfbd7cc5c3716823cf58303cad899462 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,6 +22,7 @@ Changes currently on master, but not under a tag. - Add functionality to sample in redshift and reconstruction of source frame masses. - Add functionality to combine result objects. - Enable initial values for emcee to be specified. +- Added Interferometer.plot_time_domain_strain ### Changed - Changed to using `setuptools` for installation. diff --git a/examples/injection_examples/plot_time_domain_data.py b/examples/injection_examples/plot_time_domain_data.py new file mode 100644 index 0000000000000000000000000000000000000000..a0b1a61a889cafb1569030680fec478f12d4da5b --- /dev/null +++ b/examples/injection_examples/plot_time_domain_data.py @@ -0,0 +1,38 @@ +#!/bin/python +""" +""" +from __future__ import division, print_function +import numpy as np +import tupak + +np.random.seed(1) + +duration = 4 +sampling_frequency = 2048. + +outdir = 'outdir' +label = 'example' + +injection_parameters = dict( + mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0, + phi_12=1.7, phi_jl=0.3, luminosity_distance=1000., iota=0.4, psi=2.659, + phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108) + +waveform_arguments = dict(waveform_approximant='IMRPhenomPv2', + reference_frequency=50.) + +waveform_generator = tupak.gw.WaveformGenerator( + duration=duration, sampling_frequency=sampling_frequency, + frequency_domain_source_model=tupak.gw.source.lal_binary_black_hole, + parameters=injection_parameters, waveform_arguments=waveform_arguments) +hf_signal = waveform_generator.frequency_domain_strain() + +H1 = tupak.gw.detector.get_interferometer_with_fake_noise_and_injection( + 'H1', injection_polarizations=hf_signal, + injection_parameters=injection_parameters, duration=duration, + sampling_frequency=sampling_frequency, outdir=outdir) + +t0 = injection_parameters['geocent_time'] +H1.plot_time_domain_data(outdir=outdir, label=label, notches=[50], + bandpass_frequencies=(50, 200), start_end=(-0.5, 0.5), + t0=t0) diff --git a/tupak/gw/detector.py b/tupak/gw/detector.py index 83f74aa020f0f16abb1fb8c8d8fdc5c13d96912e..b571d7cc55e8912e3948db7d8dba62ad338e9798 100644 --- a/tupak/gw/detector.py +++ b/tupak/gw/detector.py @@ -14,6 +14,7 @@ from .calibration import Recalibrate try: import gwpy + import gwpy.signal except ImportError: logger.warning("You do not have gwpy installed currently. You will " " not be able to use some of the prebuilt functions.") @@ -1356,6 +1357,15 @@ class Interferometer(object): """ The frequency domain strain in units of strain / Hz """ return self.strain_data.frequency_domain_strain + @property + def time_domain_strain(self): + """ The time domain strain in units of s """ + return self.strain_data.time_domain_strain + + @property + def time_array(self): + return self.strain_data.time_array + def time_delay_from_geocenter(self, ra, dec, time): """ Calculate the time delay from the geocenter for the interferometer. @@ -1491,6 +1501,73 @@ class Interferometer(object): '{}/{}_{}_frequency_domain_data.png'.format( outdir, self.name, label)) + def plot_time_domain_data( + self, outdir='.', label=None, bandpass_frequencies=(50, 250), + notches=None, start_end=None, t0=None): + """ Plots the strain data in the time domain + + Parameters + ---------- + outdir, label: str + Used in setting the saved filename. + bandpass: tuple, optional + A tuple of the (low, high) frequencies to use when bandpassing the + data, if None no bandpass is applied. + notches: list, optional + A list of frequencies specifying any lines to notch + start_end: tuple + A tuple of the (start, end) range of GPS times to plot + t0: float + If given, the reference time to subtract from the time series before + plotting. + + """ + + # We use the gwpy timeseries to perform bandpass and notching + if notches is None: + notches = list() + timeseries = gwpy.timeseries.TimeSeries( + data=self.time_domain_strain, times=self.time_array) + zpks = [] + if bandpass_frequencies is not None: + zpks.append(gwpy.signal.filter_design.bandpass( + bandpass_frequencies[0], bandpass_frequencies[1], + self.strain_data.sampling_frequency)) + if notches is not None: + for line in notches: + zpks.append(gwpy.signal.filter_design.notch( + line, self.strain_data.sampling_frequency)) + if len(zpks) > 0: + zpk = gwpy.signal.filter_design.concatenate_zpks(*zpks) + strain = timeseries.filter(zpk, filtfilt=True) + else: + strain = timeseries + + fig, ax = plt.subplots() + + if t0: + x = self.time_array - t0 + xlabel = 'GPS time [s] - {}'.format(t0) + else: + x = self.time_array + xlabel = 'GPS time [s]' + + ax.plot(x, strain) + ax.set_xlabel(xlabel) + ax.set_ylabel('Strain') + + if start_end is not None: + ax.set_xlim(*start_end) + + fig.tight_layout() + + if label is None: + fig.savefig( + '{}/{}_time_domain_data.png'.format(outdir, self.name)) + else: + fig.savefig( + '{}/{}_{}_time_domain_data.png'.format(outdir, self.name, label)) + class TriangularInterferometer(InterferometerList):