# Copyright (C) 2023 Leo Tsukada
#
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 2 of the License, or (at your
# option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
## @file
## @package plotdtdphi
#
# =============================================================================
#
# Preamble
#
# =============================================================================
#
import matplotlib
from matplotlib import figure
from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
matplotlib.rcParams.update({
"font.size": 10.0,
"axes.titlesize": 10.0,
"axes.labelsize": 10.0,
"xtick.labelsize": 8.0,
"ytick.labelsize": 8.0,
"legend.fontsize": 8.0,
"figure.dpi": 300,
"savefig.dpi": 300,
"text.usetex": True
})
from matplotlib import pyplot
import numpy
import lal
from lalburst.snglcoinc import light_travel_time
from gstlal.stats import inspiral_extrinsics
TPDPDF = inspiral_extrinsics.InspiralExtrinsics()
R = {"H1":lal.CachedDetectors[lal.LHO_4K_DETECTOR].response,"L1":lal.CachedDetectors[lal.LLO_4K_DETECTOR].response,"V1":lal.CachedDetectors[lal.VIRGO_DETECTOR].response,"K1":lal.CachedDetectors[lal.KAGRA_DETECTOR].response}
X = {"H1":lal.CachedDetectors[lal.LHO_4K_DETECTOR].location,"L1":lal.CachedDetectors[lal.LLO_4K_DETECTOR].location,"V1":lal.CachedDetectors[lal.VIRGO_DETECTOR].location,"K1":lal.CachedDetectors[lal.KAGRA_DETECTOR].location}
[docs]def dtdphi2prob(dt, dphi, ifo1, ifo2, refsnr, refhorizon):
"""function that takes dt,dphi samples and return the probability (density) values derived from the new pdf
:dt: the time delay between detector pair
:dphi: the difference of the coalescence phase between detector pair
:returns: dtdpipdf(dt, dphi) * snr^4 (assuming snr = {"H1": 5., "V1": 2.25} and horizon = {"H1": 110., "V1": 45.})
"""
t = {ifo1: 0, ifo2: dt}
phi = {ifo1: 0, ifo2: dphi}
snr = {ifo1: refsnr[ifo1], ifo2: refsnr[ifo2]}
horizon = {ifo1: refhorizon[ifo1], ifo2: refhorizon[ifo2]}
# signature is (time, phase, snr, horizon)
p = TPDPDF.time_phase_snr(t, phi, snr, horizon) * (sum(s**2 for s in snr.values())**.5)**4
return float(p)
[docs]def init_plot(figsize):
fig = figure.Figure()
FigureCanvas(fig)
fig.set_size_inches(figsize)
axes = fig.gca()
return fig, axes
[docs]def plots_dtdphi(ifo1, ifo2, snrs, horizons, sngl=None):
len_arr = 100
dt_max = light_travel_time(ifo1, ifo2)
dtvec = numpy.linspace(-dt_max, dt_max, len_arr)
dphivec = numpy.linspace(-2 * numpy.pi, 2 * numpy.pi, len_arr)
DTProb = numpy.zeros(len_arr)
DPProb = numpy.zeros(len_arr)
Prob = numpy.zeros((len_arr,len_arr))
for j, dt in enumerate(dtvec):
for i, dphi in enumerate(dphivec):
p = dtdphi2prob(dt, dphi, ifo1, ifo2, snrs, horizons)
Prob[j,i] = p
DPProb[i] += p
DTProb[j] += p
fig = pyplot.figure(figsize=(7.5,7.5))
ax1 = fig.add_subplot(221)
ax1.plot(dphivec, DPProb / numpy.sum(DPProb) / (dphivec[1] - dphivec[0]), label ="Direct Calculation")
ax1.set_ylabel(r"""$P(\Delta\phi | s, \{D_{%s}, D_{%s}\}, \{\rho_{%s}, \rho_{%s}\})$""" % (ifo1, ifo2, ifo1, ifo2))
ax1.legend(bbox_to_anchor=(1, 1), loc='upper left', ncol=1)
ax2 = fig.add_subplot(223)
ax2.pcolor(dphivec, dtvec, Prob)
ax2.set_xlabel(r"""$\phi_{%s} - \phi_{%s}$""" % (ifo2, ifo1))
ax2.set_ylabel(r"""$t_{%s} - t_{%s}$""" % (ifo2, ifo1))
ax3 = fig.add_subplot(224)
ax3.plot(DTProb / numpy.sum(DTProb) / (dtvec[1] - dtvec[0]), dtvec)
ax3.set_xlabel(r"""$P(\Delta t | s, \{D_{%s}, D_{%s}\}, \{\rho_{%s}, \rho_{%s}\})$""" % (ifo1, ifo2, ifo1, ifo2))
if sngl is not None:
dt_ref = sngl["dt"]
dphi_ref = sngl["dphi"]
ax1.axvline(dphi_ref, ls="--", color = "r", lw=4)
ax2.plot(dphi_ref, dt_ref, "ko", mfc = "None", mec = "r", ms = 14, mew=4)
ax3.axhline(dt_ref, ls="--", color = "r", lw=4)
fig.tight_layout(pad = .8)
return fig