#!/usr/bin/env python
# coding: utf-8

# # Gravitational wave emission by a single orbiting particle
# 
# 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`. 

# In[1]:


version()


# In[2]:


get_ipython().run_line_magic('display', 'latex')


# In[3]:


# For benchmarking
import time
system_time0 = time.time()
CPU_time0 = time.perf_counter()


# In[4]:


from kerrgeodesic_gw import (h_plus_particle, h_cross_particle)


# Directory to store the figure files (created if it does not exist already):

# In[5]:


import os
figdir = "figures_gw_single_particle/"
if not os.path.exists(figdir):
    os.makedirs(figdir)


# In[6]:


def ordital_period(a, r0):
    return RDF(2*pi*(r0^(3/2) + a))


# ### Waveform from near the ISCO for $a=0$

# Plot for comparison with Fig. 1 of E. Poisson, PRD **47**, 1497 (1993):

# In[7]:


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


# ### Waveform from the ISCO for $a=0$

# In[8]:


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


# In[9]:


show(graph, figsize=7)


# ### Waveform from $r_0=12M$ for $a=0$

# In[10]:


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


# ### Waveform from $r_0=50M$ for $a=0$

# In[11]:


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


# ### Waveform from the ISCO for $a=0.5 M$

# In[12]:


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


# ### Waveform from the ISCO for $a=0.9 M$

# In[13]:


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 Detweiler (1978)
# 
# 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`)

# In[14]:


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):

# In[15]:


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


# In[16]:


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


# In[17]:


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))