This is a TEST notebook for the functions h_plus_particle
and h_cross_particle
. It does not make use of the function plot_h_particle
from the kerrgeodesic_gw
package; instead the plots are performed by loops calling h_plus_particle
and h_cross_particle
at each step. This is much less efficient than using plot_h_particle
.
version()
'SageMath version 9.4, Release Date: 2021-08-22'
%display latex
# For benchmarking
import time
system_time0 = time.time()
CPU_time0 = time.perf_counter()
from kerrgeodesic_gw import (h_plus_particle, h_cross_particle)
Directory to store the figure files (created if it does not exist already):
import os
figdir = "figures_gw_single_particle/"
if not os.path.exists(figdir):
os.makedirs(figdir)
def ordital_period(a, r0):
return RDF(2*pi*(r0^(3/2) + a))
Plot for comparison with Fig. 1 of E. Poisson, PRD 47, 1497 (1993):
a = 0.
r0 = 6.2 # near ISCO value used by Poisson (1993)
theta, phi = 4*pi/9, 0
norm_Poisson93 = -0.065^(2/3) # the minus sign arises from a different metric signature
hp = lambda t: h_plus_particle(a, r0, t, theta, phi, l_max=5,
algorithm_Zinf='1.5PN') / norm_Poisson93
hc = lambda t: h_cross_particle(a, r0, t, theta, phi, l_max=5,
algorithm_Zinf='1.5PN') / norm_Poisson93
tmax = 2*ordital_period(a, r0)
graph = plot(hp, (0, tmax), thickness=1.5, legend_label=r'$h_+$',
axes_labels=[r'$t/M$', r'$r h / \mu/(M\Omega)^{2/3}$'],
gridlines=True, frame=True, axes=False,
title=r"$a=0,\quad r_0=6.2M,\quad \theta=\frac{4\pi}{9}$")
graph += plot(hc, (0, tmax), color='red', thickness=1.5, legend_label=r'$h_\times$')
graph.save(figdir + "compar_fig1_Poisson93.pdf")
graph
a = 0.
r0 = 6.
theta, phi = 4*pi/9, 0
hp = lambda t: h_plus_particle(a, r0, t, theta, phi)
hc = lambda t: h_cross_particle(a, r0, t, theta, phi)
tmax = 2*ordital_period(a, r0)
graph = plot(hp, (0, tmax), thickness=1.5, legend_label=r'$h_+$',
axes_labels=[r'$t/M$', r'$r h / \mu$'],
gridlines=True, frame=True, axes=False,
title=r"$a=0,\quad r_0=6M,\quad \theta=\frac{4\pi}{9}$")
graph += plot(hc, (0, tmax), color='red', thickness=1.5, legend_label=r'$h_\times$')
hp = lambda t: h_plus_particle(a, r0, t, theta, phi, l_max=5,
algorithm_Zinf='1.5PN')
hc = lambda t: h_cross_particle(a, r0, t, theta, phi, l_max=5,
algorithm_Zinf='1.5PN')
graph += plot(hp, (0, tmax), linestyle=':', thickness=1.5,
legend_label=r'$h_+$ (1.5PN)')
graph += plot(hc, (0, tmax), color='red', linestyle=':', thickness=1.5,
legend_label=r'$h_\times$ (1.5PN)')
graph.save(figdir + "waveform_a0_r6.pdf")
graph
show(graph, figsize=7)
a = 0.
r0 = 12.
theta, phi = 4*pi/9, 0
hp = lambda t: h_plus_particle(a, r0, t, theta, phi)
hc = lambda t: h_cross_particle(a, r0, t, theta, phi)
tmax = 2*ordital_period(a, r0)
graph = plot(hp, (0, tmax), thickness=1.5, legend_label=r'$h_+$',
axes_labels=[r'$t/M$', r'$r h / \mu$'],
gridlines=True, frame=True, axes=False,
title=r"$a=0,\quad r_0=12M,\quad \theta=\frac{4\pi}{9}$")
graph += plot(hc, (0, tmax), color='red', thickness=1.5, legend_label=r'$h_\times$')
hp = lambda t: h_plus_particle(a, r0, t, theta, phi, l_max=5,
algorithm_Zinf='1.5PN')
hc = lambda t: h_cross_particle(a, r0, t, theta, phi, l_max=5,
algorithm_Zinf='1.5PN')
graph += plot(hp, (0, tmax), linestyle=':', thickness=1.5,
legend_label=r'$h_+$ (1.5PN)')
graph += plot(hc, (0, tmax), linestyle=':', color='red', thickness=1.5,
legend_label=r'$h_\times$ (1.5PN)')
graph.save(figdir + "waveform_a0_r12.pdf")
graph
a = 0.
r0 = 50.
theta, phi = 4*pi/9, 0
hp = lambda t: h_plus_particle(a, r0, t, theta, phi, l_max=5)
hc = lambda t: h_cross_particle(a, r0, t, theta, phi, l_max=5)
tmax = 2*ordital_period(a, r0)
graph = plot(hp, (0, tmax), thickness=1.5, legend_label=r'$h_+$',
axes_labels=[r'$t/M$', r'$r h / \mu$'],
gridlines=True, frame=True, axes=False,
title=r"$a=0,\quad r_0=50M,\quad \theta=\frac{4\pi}{9}$")
graph += plot(hc, (0, tmax), color='red', thickness=1.5, legend_label=r'$h_\times$')
hp = lambda t: h_plus_particle(a, r0, t, theta, phi, l_max=5,
algorithm_Zinf='1.5PN')
hc = lambda t: h_cross_particle(a, r0, t, theta, phi, l_max=5,
algorithm_Zinf='1.5PN')
graph += plot(hp, (0, tmax), linestyle=':', thickness=1.5,
legend_label=r'$h_+$ (1.5PN)')
graph += plot(hc, (0, tmax), linestyle=':', color='red',thickness=1.5,
legend_label=r'$h_\times$ (1.5PN)')
graph.save(figdir + "waveform_a0_r50.pdf")
graph
a = 0.5
r0 = 4.234
theta, phi = 4*pi/9, 0
hp = lambda t: h_plus_particle(a, r0, t, theta, phi, l_max=7)
hc = lambda t: h_cross_particle(a, r0, t, theta, phi, l_max=7)
tmax = 2*ordital_period(a, r0)
graph = plot(hp, (0, tmax), thickness=1.5, legend_label=r'$h_+$',
axes_labels=[r'$t/M$', r'$r h / \mu$'],
gridlines=True, frame=True, axes=False,
title=r"$a=0.5M,\quad r_0=4.234M,\quad \theta=\frac{4\pi}{9}$")
graph += plot(hc, (0, tmax), color='red', thickness=1.5, legend_label=r'$h_\times$')
graph.save(figdir + "waveform_a0.5_r4.234.pdf")
graph
a = 0.9
r0 = 2.321
theta, phi = 4*pi/9, 0
hp = lambda t: h_plus_particle(a, r0, t, theta, phi)
hc = lambda t: h_cross_particle(a, r0, t, theta, phi)
tmax = 2*ordital_period(a, r0)
graph = plot(hp, (0, tmax), thickness=1.5, legend_label=r'$h_+$', plot_points=200,
axes_labels=[r'$t/M$', r'$r h / \mu$'],
gridlines=True, frame=True, axes=False,
title=r"$a=0.9M,\quad r_0=2.321M,\quad \theta=\frac{4\pi}{9}$")
graph += plot(hc, (0, tmax), color='red', thickness=1.5, legend_label=r'$h_\times$',
plot_points=200)
graph.save(figdir + "waveform_a0.9_r2.321.pdf")
graph
Comparison with Fig. 8 of S.L. Detweiler, ApJ 225, 687 (1978) (given that Detweiler is using
$-u = r_*-t$ on the $x$-axis, we use $-t$ as argument of h_plus_particle
)
a = 0.5
r0 = 4.45
theta, phi = pi/2, 0
hp = lambda t: h_plus_particle(a, r0, -t, theta, phi, l_max=7) # NB: t --> -t
hc = lambda t: h_cross_particle(a, r0, -t, theta, phi, l_max=7)
graph = plot(hp, (0, 125), thickness=1.5, legend_label=r'$h_+$',
axes_labels=[r'$(r_*-t)/M$', r'$r h / \mu$'], ymin=-1, ymax=1,
gridlines=True, frame=True, axes=False,
title=r"$a=0.5M,\quad r_0=4.45M,\quad \theta=\frac{\pi}{2}$")
graph += plot(hc, (0, 125), color='red', thickness=1.5, legend_label=r'$h_\times$')
graph.save(figdir + "compar_fig8_Detweiler78_lmax7.pdf")
graph
Comparison with Fig. 9 of S.L. Detweiler, ApJ 225, 687 (1978):
a = 0.9
r0 = 2.4
theta, phi = pi/2, 0
hp = lambda t: h_plus_particle(a, r0, -t, theta, phi, l_max=7) # NB: t --> -t
hc = lambda t: h_cross_particle(a, r0, -t, theta, phi, l_max=7)
graph = plot(hp, (0, 60), thickness=1.5, legend_label=r'$h_+$',
axes_labels=[r'$(r_*-t)/M$', r'$r h / \mu$'], ymin=-1, ymax=1,
gridlines=True, frame=True, axes=False,
title=r"$a=0.9M,\quad r_0=2.4M,\quad \theta=\frac{\pi}{2}$")
graph += plot(hc, (0, 60), color='red', thickness=1.5, legend_label=r'$h_\times$')
graph.save(figdir + "compar_fig8_Detweiler78_lmax_7.pdf")
graph
a = 0.9
r0 = 2.4
theta, phi = pi/2, 0
graph = Graphics()
for l_max in [5, 7, 10]:
hp = lambda t: h_plus_particle(a, r0, -t, theta, phi, l_max=l_max) # NB: t --> -t
graph += plot(hp, (0, 60), thickness=1.5, color=hue((l_max-5)/8),
legend_label=r"$\ell_{{\rm max}} = {}$".format(l_max),
axes_labels=[r'$(r_*-t)/M$', r'$r h_+ / \mu$'], ymin=-1., ymax=1.,
gridlines=True, frame=True, axes=False,
title=r"$a=0.9M,\quad r_0=2.4M,\quad \theta=\frac{\pi}{2}$")
graph.save(figdir + "compar_fig9_Detweiler78.pdf")
graph
print("System time elapsed since the start of this notebook: {} s".format(
time.time() - system_time0))
print("CPU time elapsed since the start of this notebook: {} s".format(
time.perf_counter() - CPU_time0))
System time elapsed since the start of this notebook: 53.653218507766724 s CPU time elapsed since the start of this notebook: 53.653356110999994 s