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

# ## Fast Kerr geodesics

# In[1]:


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


# Using the manifold catalog (https://trac.sagemath.org/ticket/25869):

# In[2]:


K.<t, r, th, ph> = manifolds.Kerr(2, -1, coordinates="BL")


# Using the metric to normalize the initial 4-velocity.

# In[3]:


g = K.metric()
g[:]


# $\tau$ is the proper time.

# In[4]:


tau = var('tau')


# We choose a starting point $p$ and a normalized initial 4-speed $v$ in the tangent space $T_p$.

# In[5]:


p = K((0, 14.98, pi/2, 0))
Tp = K.tangent_space(p)
v = Tp((2, 0, 0.005, 0.05))
v = v / sqrt(-g.at(p)(v, v))


# Declaration of the geodesic:

# In[6]:


c = K.integrated_geodesic(g, (tau, 0, 3000), v)


# Fast integration:

# In[7]:


sol = c.solve(step = 1, method="ode_int")


# Interpolation:

# In[8]:


interp = c.interpolate()


# An embedding in $R^3$ is needed to plot the result:

# In[9]:


E.<x,y,z> = manifolds.Euclidean(3)
phi = K.diff_map(E, [r*sin(th)*cos(ph), r*sin(th)*sin(ph), r*cos(th)])


# Plotting the result:

# In[10]:


P = c.plot_integrated(mapping=phi, color="red", thickness=2, plot_points=3000)
P = P + sage.plot.plot3d.shapes.Sphere(4, color='grey')
P.show(aspect_ratio=[1, 1, 1], viewer='threejs', online=True)


# ### A few speed tests:

# In[11]:


get_ipython().run_line_magic('time', 'sol = c.solve(step=1, method="ode_int")')


# In[12]:


get_ipython().run_line_magic('time', 'sol = c.solve(step=1, method="rk4_maxima")')


# In[13]:


get_ipython().run_line_magic('timeit', 'K.default_chart().valid_coordinates(1, 1, 1, 1)')


# In[14]:


get_ipython().run_line_magic('timeit', 'K.default_chart().valid_coordinates_numerical(1, 1, 1, 1)')


# In[ ]: