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

# # Timelike geodesic in Schwarzschild spacetime
# 
# This Jupyter/SageMath worksheet presents the numerical computation of a timelike geodesic in Schwarzschild spacetime, given an initial point and tangent vector. 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.
# 
# Click [here](https://raw.githubusercontent.com/sagemanifolds/SageManifolds/master/Worksheets/v1.1/SM_simple_goed_Schwarz.ipynb) to download the worksheet file (ipynb format). To run it, you must start SageMath within the Jupyter notebook, via the command `sage -n jupyter`

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


# Then, 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 pick an initial point and an initial tangent vector:

# In[6]:


p0 = M.point((0, 8*m, pi/2, 1e-12), name='p_0')
v0 = M.tangent_space(p0)((1.297513, 0, 0, 0.0640625/m), name='v_0')
v0.display()


# We declare a geodesic with such initial conditions, denoting by $s$ the affine parameter (proper time), with $(s_{\rm min}, s_{\rm max})=(0, 1500\,m)$:

# In[7]:


s = var('s')
geod = M.integrated_geodesic(g, (s, 0, 1500), v0); 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 `X3` of $\mathbb{R}^3$:

# In[8]:


sol = geod.solve(parameters_values={m: 1})  # numerical integration
interp = geod.interpolate()                 # interpolation of the solution for the plot
graph = geod.plot_integrated(chart=X3, mapping=to_R3, plot_points=500, 
                             thickness=2, label_axes=False)           # the geodesic
graph += p0.plot(chart=X3, mapping=to_R3, size=4, parameters={m: 1}) # the starting point
graph += sphere(size=2, color='grey')                                 # the event horizon
show(graph, viewer='threejs', online=True)


# A Tachyon view of the geodesic:

# In[9]:


show(graph, viewer='tachyon', aspect_ratio=1, figsize=10)


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

# In[10]:


geod.system(verbose=True)


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

# In[11]:


g.christoffel_symbols_display()


# In[12]:


g.riemann().display_comp()


# In[ ]: