Animation of impact of beam scattering noise on spectral sensitivity of LIGO-Hanford gravitational wave detector

ligo-scattering-animation.py

#!/usr/bin/env python

"""Load, filter, and plot LIGO-Hanford data including optical scattering noise
"""

import glob
import os
import subprocess

from matplotlib import (pyplot, rcParams)
from matplotlib.animation import (FuncAnimation, writers)
from matplotlib.gridspec import GridSpec

from gwpy.timeseries import TimeSeries
from gwpy.plotter import TimeSeriesPlot

rcParams.update({  # separate subplots more than the default
    'figure.subplot.hspace': .7,
    'figure.subplot.top': .95,
})

# set GPS times for data access, and for plotting
gps = 1134294143.940
start = int(gps) - 32
end = int(gps) + 32
pstart = int(gps) - 10
pend = int(gps) + 10

# get data
data = TimeSeries.fetch_open_data('H1', start, end, cache=True, verbose=True)

# clean data using bandpass and notch filtering
clean = data.bandpass(10, 200).notch(60).crop(pstart, pend)

# calculate ASD of quiet time
quiet = (int(gps) - 7, int(gps) - 1)
qasd = data.crop(*quiet).asd(fftlength=.5)

# initialise figure
plot = pyplot.figure(figsize=(12, 7))
grid = GridSpec(4, 1)
tax = plot.add_subplot(grid[:1], projection='timeseries')
fax = plot.add_subplot(grid[1:], projection='frequencyseries')

# plot timeseries
tax.plot(clean)
tax.set_ylim(-3e-19, 3e-19)
tax.set_ylabel('Strain amplitude', y=.38)
marker = tax.add_patch(pyplot.Polygon([[0, 0], [0, 0], [0, 0], [0, 0]],
                                      fc='lightgray', ec='r'))
tax.text(1., 1., '10--200 Hz bandpass, 60 Hz notch', ha='right', va='bottom',
         fontsize=12, transform=tax.transAxes)

# plot reference spectrum and initialise line for moving spectrum
line, = fax.plot([], [],)
fax.plot(qasd, zorder=0)
fax.set_xscale('log')
fax.set_xlim(12, 1000)
fax.set_yscale('log')
fax.set_ylim(3e-24, 5e-21)
fax.set_xlabel('Frequency [Hz]')
fax.set_ylabel(r'Strain noise amplitude [Hz$^{-1/2}$]')
fax.set_xticks([100, 1000])
fax.set_xticklabels(['100', '1000'])

start = float(int(gps) - 10)


def update(n):
    """Draw ASD and timestamp marker at each timestep
    """
    y1, y2 = tax.get_ylim()
    t = start + n * .5
    ts = data.crop(t, t+1.5)
    asd = ts.asd(fftlength=.5)
    marker.set_xy([(t, y1), (t+1.5, y1), (t+1.5, y2), (t, y2)])
    line.set_data(asd.xindex.value, asd.value)
    return marker, line


# make animation
nframes = 40
fps = 5.
delay = int(1 / fps * 1000)
anim = FuncAnimation(plot, update, init_func=init, frames=nframes,
                     interval=delay, blit=True)
writer = writers['ffmpeg'](fps=fps)
anim.save('scattering.mp4', writer, dpi=600)

Dear Duncan,

Thank you for this interesting code.
However, I needed to manipulate it a little bit.

1. glob, os, subprocess are unused in this script, so I deleted them.
2. gwpy has no module gwpy.plotter hence the error: ModuleNotFoundError: No module named 'gwpy.plotter'
   Anyway, it is also unused, so I deleted its import, too.
3. Then I had this error: ValueError: Unknown projection 'timeseries', so I commented out the following two lines:
   tax = plot.add_subplot(grid[:1], projection='timeseries')
   fax = plot.add_subplot(grid[1:], projection='frequencyseries')
4. Then I got: NameError: name 'init' is not defined, so I commented out the following line:
   anim = FuncAnimation(plot, update, init_func=init, frames=nframes,interval=delay, blit=True)
5. Now I have: RuntimeError: Requested MovieWriter (ffmpeg) not available

But it works and I have the beautiful and very instructive animation!
I use miniconda3 (Python 3.7) and gwpy 0.15.0

Best wishes

Kamran Derakhshani