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

# # Numerical values
# 
# 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/numerical_values.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`
# 

# In[1]:


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


# ## Physical constants

# Values of the fundamental constants in SI units  (taken from the [Particle Data Group](http://pdg.lbl.gov/) (2017)):

# In[2]:


c = 2.99792458e8
G = 6.67408e-11
hbar = 1.054571800e-34


# The solar mass $M_\odot$ in SI units:

# In[3]:


M_sol = 1.98848e30


# The Planck mass in kg:

# In[4]:


sqrt(hbar*c/G)


# The solar mass $M_\odot$ in meters:

# In[5]:


M_sol_m = G*M_sol/c^2
M_sol_m


# The solar mass $M_\odot$ in seconds:

# In[6]:


M_sol_s = G*M_sol/c^3
M_sol_s


# The astronomical unit (au) in SI units:

# In[7]:


au = 1.49597870700e11


# ## A sample of black holes masses
# 
# Masses in $M_\odot$:

# In[8]:


masses = [15., 4.3e6, 6e9]
masses


# Schwarzschild radii:

# In[9]:


for m in masses:
    Rs = RDF(2*m * M_sol_m)
    Rs_km = RDF(Rs/1000)
    Rs_au = RDF(Rs/au)
    print("m = {} M_sol:  2m = {} km = {} au".format(m, Rs_km, Rs_au))    


# ## Free fall in a Schwarzschild black hole

# Proper time to reach the singularity starting from rest at the ISCO:

# In[10]:


var('r0')
tau_f(r0) = pi*sqrt(r0^3/8)
tau_f(r0)


# Proper time to reach the horizon starting from rest at the ISCO:

# In[11]:


tau_h(r0) = sqrt(r0^3/2)*(atan(sqrt(r0/2-1)) + sqrt(2/r0*(1-2/r0)))
tau_h(r0)


# In[12]:


for m in masses:
    t_s = RDF(tau_f(6) * m * M_sol_s)
    t_h = RDF(t_s/3600)
    t_d = RDF(t_h/24)
    print("m = {} M_sol:".format(m))    
    print("  tau_f = {} s = {} h = {} days".format(t_s, t_h, t_d))


# In[13]:


for m in masses:
    t_s = RDF(tau_h(6) * m * M_sol_s)
    t_h = RDF(t_s/3600)
    t_d = RDF(t_h/24)
    print("m = {} M_sol:".format(m))    
    print("  tau_h = {} s = {} h = {} days".format(t_s, t_h, t_d))


# Time spent into the black hole:

# In[14]:


for m in masses:
    t_s = RDF((tau_f(6) - tau_h(6)) * m * M_sol_s)
    t_h = RDF(t_s/3600)
    t_d = RDF(t_h/24)
    print("m = {} M_sol:".format(m))    
    print("  tau_inside = {} s = {} h = {} days".format(t_s, t_h, t_d))


# ## ISCO frequency

# ### Schwarzschild

# In[15]:


for m in masses:
    RI = RDF(6*m * M_sol_m)
    RI_km = RDF(RI/1000)
    RI_au = RDF(RI/au)
    f_Hz = RDF(1/(6*sqrt(6)* m * M_sol_s)/(2*pi))
    T_s = 1/f_Hz
    T_min = T_s/60
    print("m = {} M_sol:  r_ISCO = {} km = {} au".format(m, RI_km, RI_au)) 
    print("f_ISCO = {} Hz, T_ISCO = {} s = {} min".format(f_Hz, T_s, T_min))    


# ### Extreme Kerr

# In[16]:


for m in masses:
    RI = RDF(m * M_sol_m)
    RI_km = RDF(RI/1000)
    RI_au = RDF(RI/au)
    f_Hz = RDF(1/(2* m * M_sol_s)/(2*pi))
    T_s = 1/f_Hz
    T_min = T_s/60
    print("m = {} M_sol:  r_ISCO = {} km = {} au".format(m, RI_km, RI_au)) 
    print("f_ISCO = {} Hz, T_ISCO = {} s = {} min".format(f_Hz, T_s, T_min))    


# In[ ]: