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

# # Timelike orbits in Schwarzschild spacetime
# 
# This Jupyter/SageMath worksheet is relative to the lectures
# [Geometry and physics of black holes](https://relativite.obspm.fr/blackholes/).
# 
# Click [here](https://raw.githubusercontent.com/egourgoulhon/BHLectures/master/sage/ges_orbits.ipynb) to download the worksheet file (ipynb format). To run it, you must start SageMath with the Jupyter notebook, with the command `sage -n jupyter`
# 
# It uses the `integrated_geodesic` functionality introduced by Karim Van Aelst in **SageMath 8.1**, in the framework of the [SageManifolds](http://sagemanifolds.obspm.fr) project.

# A version of SageMath at least equal to 8.1 is required to run this worksheet:

# In[1]:


version()


# In[2]:


get_ipython().run_line_magic('display', 'latex  # LaTeX rendering turned on')


# We define first the spacetime manifold $M$ and the standard Schwarzschild-Droste coordinates on it:

# In[3]:


M = Manifold(4, 'M')
X.<t,r,th,ph> = M.chart(r't r:(0,+oo) th:(0,pi):\theta ph:\phi')
X


# For graphical purposes, we introduce $\mathbb{R}^3$ and some coordinate map $M\to \mathbb{R}^3$:

# In[4]:


R3 = Manifold(3, 'R^3', latex_name=r'\mathbb{R}^3')
X3.<x,y,z> = R3.chart()
to_R3 = M.diff_map(R3, {(X, X3): [r*sin(th)*cos(ph), 
                                  r*sin(th)*sin(ph), r*cos(th)]})
to_R3.display()


# Next, we define the Schwarzschild metric:

# In[5]:


g = M.lorentzian_metric('g')
m = var('m'); assume(m >= 0)
g[0,0], g[1,1] = -(1-2*m/r), 1/(1-2*m/r)
g[2,2], g[3,3] = r^2, (r*sin(th))^2
g.display()


# We set the specific conserved energy and angular momentum:

# In[6]:


#var('E L', domain='real')
E = 0.973
L = 4.2*m


# In[7]:


u = M.vector_field(name='u')
u[0] = E/(1-2*m/r)
u[1] = sqrt(E^2 - (1-2*m/r)*(1+L^2/r^2))
u[3] = L/(r^2*sin(th)^2)
u.display()


# In[8]:


g(u,u).display()


# Pericenter and apocenter:

# In[9]:


eq = (r^3*((1-2/r)*(1+L^2/m^2/r^2) - E^2)).simplify_full()
eq


# In[10]:


R.<w> = RR[]
print(R)


# In[11]:


eqp = R(eq.subs(r=w))
eqp


# In[12]:


roots = sorted(eqp.real_roots())
roots


# In[13]:


r_per = roots[1]*m
r_apo = roots[2]*m
r_per, r_apo


# We pick an initial point and an initial tangent vector:

# In[14]:


r0 = r_per
p0 = M.point((0, r0, pi/2, 1e-12), name='p_0')
Tp0 = M.tangent_space(p0) 
u0 = u.at(p0)
u0.set_name('u_0')
print(u0)


# In[15]:


u0.display()


# Let us check that the scalar square of $u_0$ is $-1$, by means of the metric tensor at $p_0$:

# In[16]:


g0 = g.at(p0)
print(g0)


# The scalar square is indeed equal to $-1$: 

# In[17]:


g0(u0,u0).simplify_full()


# ### Check
# 
# Let us compute the conserved energy and angular momentum along the geodesic by 
# taking scalar products of the 4-velocity $v_0$ with the Killing vectors $\xi=\frac{\partial}{\partial t}$ and $\eta=\frac{\partial}{\partial t}$:

# In[18]:


xi = X.frame()[0]
eta = X.frame()[3]
xi, eta


# In[19]:


xi0 = xi.at(p0)
eta0 = eta.at(p0)
print(xi0)
print(eta0)


# The specific conserved energy $\varepsilon$ is:

# In[20]:


- g0(xi0, u0)


# while the specific conserved angular momentum $\ell$ is

# In[21]:


g0(eta0, u0)


# We declare the geodesic through $p_0$ having $v_0$ as inital vector, denoting by $s$ the affine parameter (proper time):

# In[22]:


s = var('s')
geod = M.integrated_geodesic(g, (s, 0, 1500), u0); geod


# We ask for the numerical integration of the geodesic, providing some numerical value for the parameter $m$, and then plot it in terms of the Cartesian chart:

# In[23]:


sol = geod.solve(step=4, parameters_values={m: 1})  # numerical integration
interp = geod.interpolate()                 # interpolation of the solution for the plot


# In[24]:


graph = geod.plot_integrated(chart=X3, mapping=to_R3, ambient_coords=(x,y), plot_points=500, 
                             thickness=2, label_axes=False)   


# In[25]:


graph += circle((0,0), r_per/m, linestyle=':')
graph += circle((0,0), r_apo/m, linestyle=':')
graph += circle((0,0), 2, fill=True, edgecolor='grey', facecolor='grey')
show(graph, aspect_ratio=1, axes_labels=[r'$x/m$', r'$y/m$'])


# In[26]:


file_name = "ges_orbit_e{:.3f}_l{:.3f}.pdf".format(float(E), float(L/m))
graph.save(file_name, aspect_ratio=1, axes_labels=[r'$x/m$', r'$y/m$'])
file_name


# Some details about the system solved to get the geodesic:

# In[27]:


geod.system(verbose=True)


# We recognize in the above list the Christoffel symbols of the metric $g$:

# In[28]:


g.christoffel_symbols_display()


# In[ ]: