Animate the signal model for GW150914 including highpassed and filtered LIGO data

gw150914-animate.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# The MIT License (MIT)
# Copyright (c) 2017 Duncan Macleod
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
# IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
# DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
# OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE
# OR OTHER DEALINGS IN THE SOFTWARE.

"""Animate a gravitational-wave signal, and with detectior noise, and
with filtered detector noise
"""

from __future__ import division

import numpy

from matplotlib import animation
from matplotlib import pyplot

from gwpy.timeseries import TimeSeries
from gwpy.signal import filter_design
from gwpy.plotter import (rcParams, TimeSeriesPlot)

from pycbc.waveform import get_td_waveform
from pycbc.detector import Detector

pyplot.style.use('dark_background')
rcParams.update({
    'figure.dpi': 600,
    'figure.subplot.left': 0.,
    'figure.subplot.right': 1.,
    'figure.subplot.bottom': 0.,
    'figure.subplot.top': 1.,
    'text.usetex': False,
    'font.family': 'sans-serif',
    'font.sans-serif': ['Helvetica Neue', 'Helvetica'],
})

# -- animation helpers --------------------------------------------------------

def create_figure(xlim=None, ylim=None):
    plot = TimeSeriesPlot(figsize=[12, 4])
    ax = plot.gca()
    ax.set_axis_off()
    ax.arrow(0.02, 0.05, 0.02, 0, color='white',
             head_width=0.02, head_length=0.01, transform=ax.transAxes)
    ax.text(0.06, 0.05, 'Time',
            transform=ax.transAxes, color='white', va='center')
    return plot, ax


def animate(target, fig, animate_func, init_func, fps=30., duration=10.):
    nf = fps * duration
    delay = int(1 / fps * 1000)
    anim = animation.FuncAnimation(fig, animate_func, init_func=init_func,
                                   frames=int(nf), interval=delay, blit=True)
    anim.save(target, fps=fps, extra_args=['-vcodec', 'libxvid'])
    print('{} written'.format(target))
    return anim


def legend(ax, artists, labels):
    leg = ax.legend(artists, labels, loc='upper left', frameon=False,
                    fontsize=16)
    for text in leg.get_texts():  # change font colour
        text.set_color('white')
    return leg


# -- get data and simulate waveform -------------------------------------------

# create waveform
gw150914 = {
    'mass1': 40,
    'mass2': 32,
    'distance': 410,
    'coa_phase': 3.12,
}
hp, hx = get_td_waveform(approximant='IMRPhenomD', delta_t=1/4096.,
                         f_lower=30, **gw150914)

# project waveform onto LHO at the right time (and sky location)
ra = 1.78
dec = -1.22
psi = 1.56935358208
gps = 1126259462.425
Fp, Fx = Detector('H1').antenna_pattern(ra, dec, psi, gps)
sig = TimeSeries.from_pycbc(Fp * hp + Fx * hx).pad((0, 512))
print('Generated waveform')

# get real detector data
lho = TimeSeries.fetch_open_data('H1', 1126259446, 1126259478)
llo = TimeSeries.fetch_open_data('L1', 1126259446, 1126259478)
print('Downloaded LIGO data')

# shift l1 data to account for time-delay and orientation
llo.t0 += 0.0069 * llo.t0.unit
llo *= -1 * llo.unit

# shift detector data to merge at t=0
t0 = gps * lho.t0.unit
lho.t0 -= t0
llo.t0 -= t0

# rescale to sensible amplitude (makes plotting easier)
lho *= 1e21
llo *= 1e21
sig *= 1e21

# -- animate waveform ---------------------------------------------------------

print('Animating waveform...')

# crop signal to interval we want to animate
sigc = sig.crop(-.4, .2)

# set animation parameters
fps = 30.
duration = 10.
times = numpy.arange(sigc.size)

# create figure and animation methods
plot, ax = create_figure()
line, = ax.plot([], [], color='white')
ax.set_xlim(0, times.size)
ax.set_ylim(-2, 2)
npf = sig.size / (fps * duration)

def init():
    line.set_data([], [])
    return line,

def update(n):
    n = int(n * npf)
    line.set_data(times[:n], sigc[:n])
    return line,

# write animation
animate('gw150914-signal.mp4', plot, update, init, fps=fps, duration=duration)
plot.close()

# -- animate filtered data ----------------------------------------------------

filter1 = filter_design.highpass(10, lho.sample_rate)
bp = filter_design.bandpass(50, 250, lho.sample_rate)
notches = [filter_design.notch(f, lho.sample_rate) for f in (60, 120, 180)]
filter2 = filter_design.concatenate_zpks(bp, *notches)

for zpk, ylim, tag in [
       (filter1, (-300, 300), 'highpass'),
       (filter2, (-2, 2), 'filtered'),
]:
    print('Animating {} data...'.format(tag))

    lhof = lho.filter(zpk, filtfilt=True).crop(-.4, .2)
    llof = llo.filter(zpk, filtfilt=True).crop(-.4, .2)
    sigf = sig.filter(zpk, filtfilt=True).crop(-.4, .2)[:lhof.size]
    times = numpy.arange(sigf.size)

    # create figure and animation methods
    plot, ax = create_figure()
    lholine, = ax.plot([], [], color='gwpy:ligo-hanford', linewidth=1)
    lloline, = ax.plot([], [], color='gwpy:ligo-livingston', linewidth=1)
    sigline, = ax.plot([], [], color='white')
    ax.set_xlim(0, times.size)
    ax.set_ylim(*ylim)

    # add a legend
    legend(ax, [sigline, lholine, lloline],
           ['GW150914 signal', 'LIGO-Hanford', 'LIGO-Livingston'])

    def init():
        sigline.set_data([], [])
        lholine.set_data([], [])
        lloline.set_data([], [])
        return sigline, lholine, lloline

    def update(n):
        n = int(n * npf)
        sigline.set_data(times[:n], sigf.value[:n])
        lholine.set_data(times[:n], lhof.value[:n])
        lloline.set_data(times[:n], llof.value[:n])
        return sigline, lholine, lloline

    # write animation
    animate('gw150914-{}.mp4'.format(tag), plot, update, init,
            fps=fps, duration=duration)
    plot.close()

Hi again!
So i've managed to reach my goal and make this works on all detections + all the available detectors.. polarization looks like is not influencing the final result as you said. Thanks for the huge work, gonna push the videos on that host i gave you among with other plots and data <3
Kudos!!

https://gist.github.com/iGio90/b109b02914fb460eeb6e7633256922c7