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

# # Periastron/apoastron of null geodesics in Schwarzschild spacetime and related quantities

# This Jupyter/SageMath notebook 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_null_periastron.ipynb) to download the notebook file (ipynb format). To run it, you must start SageMath with `sage -n jupyter`.

# In[2]:


version()


# In[3]:


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


# ## The cubic polynomial

# In[4]:


r, b, m = var('r b m', domain='real')


# In[5]:


P(r, b, m) = r^3 - b^2*r + 2*m*b^2
P(r, b, m)


# ## Roots via Viète's trigonometric method

# In[6]:


r0 = 2/sqrt(3)*b*cos(pi/3 - arccos(3*sqrt(3)*m/b)/3)
r0


# In[7]:


P(r0, b, m).simplify_full().trig_reduce()


# In[8]:


r1 = 2/sqrt(3)*b*cos(pi - arccos(3*sqrt(3)*m/b)/3)
r1


# In[9]:


P(r1, b, m).simplify_full().trig_reduce()


# In[10]:


r2 = 2/sqrt(3)*b*cos(5*pi/3 - arccos(3*sqrt(3)*m/b)/3)
r2


# In[11]:


P(r2, b, m).simplify_full().trig_reduce()


# In[12]:


rm(b) = r1.subs({m: 1})
rp(b) = r0.subs({m: 1})
ra(b) = r2.subs({m: 1})


# In[13]:


g = plot(rp(b), (b, 3*sqrt(3), 8), color='red', legend_label=r'$r_{\rm per}$',
         axes_labels=[r'$b/m$', r'$r$']) \
    + plot(ra(b), (b, 3*sqrt(3), 8), color='green', legend_label=r'$r_{\rm apo}$') \
    + plot(rm(b), (b, 3*sqrt(3), 8), color='blue', legend_label=r'$r_{\rm neg}$')
g


# In[14]:


g = plot(rp(b), (b, 3*sqrt(3), 10), color='red', thickness=1.5,
         legend_label=r'$r_{\rm per}$', axes_labels=[r'$b/m$', r'$r/m$'], 
         frame=True, gridlines=True, aspect_ratio=1) \
    + plot(ra(b), (b, 3*sqrt(3), 10), color='green', thickness=1.5,
           legend_label=r'$r_{\rm apo}$')
g


# In[15]:


show(g, xmin=5, ymin=2)


# In[16]:


g.save('ges_null_per_apo.pdf', xmin=5, ymin=2)


# ## Expansions close to $b_{\rm crit}$

# ### Expansion of the form $b = b_{\rm crit} + x$

# In[17]:


bcrit = 3*sqrt(3)


# In[18]:


rpx = rp(bcrit + x).taylor(x, 0, 2)
rpx


# In[19]:


rax = ra(bcrit + x).taylor(x, 0, 2)
rax


# In[20]:


upx = (1/rpx).taylor(x, 0, 2)
upx


# In[21]:


uax = (1/rax).taylor(x, 0, 2)
uax


# ### Expansion of the form $u_{\rm per} = \frac{1}{3} (1 - x)$

# In[22]:


up = (1 - x)/3


# In[23]:


bx = 1/(up*sqrt(1 - 2*up)).simplify_full()
bx


# In[24]:


upx = (1/rp(bx)).simplify_full().taylor(x, 0, 3)
upx


# In[25]:


P(1/upx, bx, 1).simplify_full()


# In[26]:


uax = (1/ra(bx)).simplify_full().taylor(x, 0, 3)
uax


# In[27]:


P(1/uax, bx, 1).simplify_full().taylor(x, 0, 3)


# In[28]:


umx = (1/rm(bx)).simplify_full().taylor(x, 0, 3)
umx


# In[29]:


P(1/umx, bx, 1).simplify_full().taylor(x, 0, 3)


# Another check:

# In[30]:


umx + upx + uax


# ## The cubic polynomial $P_b(u)$

# In[31]:


P(u, b) = 2*u^3 - u^2 + 1/b^2
P


# In[32]:


b_crit = 3*sqrt(3)
b_sel = [3, 4, b_crit, 7, 20]
b_sel


# In[33]:


graph = Graphics()
for b in b_sel:
    if b == b_crit:
        legend_label = r'$b = 3\sqrt{3}\, m$'
    else:
        legend_label=r'$b = {:.0f} \, m$'.format(float(b))
    graph += plot(P(u, b), (u, -0.5, 0.8), color = hue((b - b_crit)/6), 
                  legend_label=legend_label, thickness=1.5,
                  axes_labels=[r"$u$", r"$P_b(u)$"])
show(graph, ymin=-0.2, ymax=0.2)


# In[34]:


graph.save('ges_polynomial_b_u.pdf', ymin=-0.2, ymax=0.2)


# ### Roots of $P_b(u)$
# 
# The depressed polynomial:

# In[35]:


b = var('b')
pd = (P(x + 1/6, b)/2).expand()
pd


# In[36]:


p = -1/12
q = 1/(2*b^2) - 1/108
p, q


# Roots via Viète formula:

# In[37]:


3*q/(2*p)*sqrt(-3/p)


# In[38]:


psi = arcsin(b_crit/b)
un(b) = 1/3*cos(2*psi/3 + 2*pi/3) + 1/6
un(b)


# In[39]:


P(un(b), b).simplify_full().trig_reduce().trig_expand()


# In[40]:


up(b) = 1/3*cos(2*psi/3 + 4*pi/3) + 1/6
up(b)


# In[41]:


P(up(b), b).simplify_full().trig_reduce().trig_expand()


# In[42]:


ua(b) = 1/3*cos(2*psi/3) + 1/6
ua(b)


# In[43]:


P(ua(b), b).simplify_full().trig_reduce().trig_expand()


# In[44]:


g = plot(ua, (b_crit, 16), color='green', thickness=1.5,
         legend_label=r'$u_{\rm a}$', axes_labels=[r'$b/m$', r'$u$'], 
         frame=True, gridlines=True) \
    + plot(up, (b_crit, 16), color='red', thickness=1.5,
           legend_label=r'$u_{\rm p}$') \
    + plot(un, (b_crit, 16), color='maroon', thickness=1.5,
           legend_label=r'$u_{\rm n}$')
g


# ### Adding $u_{\rm n}$ for $b< b_{\rm c}$:

# In[45]:


xi = b_crit/b - sqrt((b_crit/b)^2 - 1)
u_neg(b) = (1 - xi^(2/3) - xi^(-2/3))/6
u_neg(b)


# Check:

# In[46]:


P(u_neg(b), b).simplify_full().factor()


# In[47]:


g += plot(u_neg, (0.1, b_crit), color='maroon', thickness=1.5) \
     + line([(b_crit, -3.5), (b_crit, 0.6)], linestyle='dashed',
            color='black')


# In[48]:


show(g, xmax=15, ymin=-0.5, ymax=0.5)


# In[49]:


g.save('gis_u_per_apo_neg.pdf', xmax=15, ymin=-0.5, ymax=0.5)


# ## Elliptic modulus $k$

# In[50]:


k(b) = sqrt((up(b) - un(b))/(ua(b) - un(b))).simplify_full()
k(b)


# In[51]:


g = plot(k(b), (b, 3*sqrt(3), 20), color='red', thickness=1.5,
         axes_labels=[r'$b/m$', r'$k$'], 
         frame=True, gridlines=True)
g


# In[52]:


g.save("gis_elliptic_mod.pdf")


# In[53]:


k1(b) = sqrt(2) / sqrt(sqrt(3)*cot(2/3*arcsin(3*sqrt(3)/b)) + 1)
k1(b)


# In[54]:


(k(b)^2 - k1(b)^2).simplify_full()


# ## Variable $t$

# In[59]:


tb(u, b) = sqrt((up(b) - u)/(ua(b) - u))/k(b)
tb(u, b)


# In[60]:


bc = b_crit
g = Graphics()
for b1 in [bc, 5.3, 5.5, 6, 8, 10, 100]:
    g += plot(tb(u, b1), (u, 0, up(b1)), 
              legend_label="b={}".format(float(b1)), 
              color=hue(b1/bc), axes_labels=[r'$u$', r'$t$'], 
              frame=True, gridlines=True)
g


# ## Function $\phi(u, b)$

# In[61]:


g = Graphics()
phib(u, b) = arcsin(tb(u, b))
for b1 in [bc, 5.3, 5.5, 6, 8, 10, 100]:
    g += plot(phib(u, b1), (u, 0, up(b1)), 
              legend_label="b={}".format(float(b1)), 
              color=hue(b1/bc), axes_labels=[r'$u$', r'$\phi(u, b)$'], 
              frame=True, gridlines=True)
g


# ### Series expansion for $b$ close to $b_{\rm crit}$

# In[54]:


upx = 1/3 - x


# In[55]:


bx = 1/(upx*sqrt(1 - 2*upx)).simplify_full()
bx


# In[56]:


up(bx).taylor(x, 0, 2)


# In[57]:


bx.taylor(x, 0, 2)


# In[58]:


uax = ua(bx).taylor(x, 0, 2)
uax


# In[59]:


unx = un(bx).taylor(x, 0, 2)
unx


# In[60]:


kx = k(bx).taylor(x, 0, 2)
kx


# In[61]:


sin_ax = (sqrt(upx/uax)/kx).taylor(x, 0, 2)
sin_ax


# In[62]:


ax = (arcsin(sqrt(upx/uax)/kx)).taylor(x, 0, 2)
ax


# In[63]:


tan_ax = tan(arcsin(sqrt(upx/uax)/kx)).taylor(x, 0, 2)
tan_ax


# In[64]:


(tan_ax * sqrt(1 - k(bx)^2)).taylor(x, 0, 2)


# In[65]:


(1/_).taylor(x, 0, 2)


# In[66]:


sin(arctan(1/sqrt(2))).simplify_full()


# In[67]:


arctan(1/sqrt(2)).simplify_full()


# In[68]:


sqrt(1 - k(bx)^2).taylor(x, 0, 2)


# In[69]:


aax = arcsin(sqrt(upx/uax)/kx)
(tan(aax) + sec(aax)).taylor(x, 0, 2)


# In[70]:


(3*sqrt(3)/bx).taylor(x, 0, 2)


# In[71]:


(2/3)*arcsin((3*sqrt(3)/bx)).taylor(x, 0, 2)


# In[ ]: