# Homework 4

## Imports


```python
import nbtools
nbtools.setup_nb()
```


```python
from astropy import units, constants
import itertools
import functools

import numpy
from scipy import stats
import pandas
import plotly.express as px
import plotly.figure_factory as ff
import plotly.graph_objects as go

import sympy
from sympy.diffgeom import Manifold, Patch, TensorProduct as tprod

from pystein import coords, metric, curvature, geodesic
from pystein.utilities import tensor_pow as tpow, full_simplify
```

## Utilities


```python
import plotly.colors

def get_continuous_color(colorscale, intermed):
    """
    Plotly continuous colorscales assign colors to the range [0, 1]. This function computes the intermediate
    color for any value in that range.

    Plotly doesn't make the colorscales directly accessible in a common format.
    Some are ready to use:
    
        colorscale = plotly.colors.PLOTLY_SCALES["Greens"]

    Others are just swatches that need to be constructed into a colorscale:

        viridis_colors, scale = plotly.colors.convert_colors_to_same_type(plotly.colors.sequential.Viridis)
        colorscale = plotly.colors.make_colorscale(viridis_colors, scale=scale)

    :param colorscale: A plotly continuous colorscale defined with RGB string colors.
    :param intermed: value in the range [0, 1]
    :return: color in rgb string format
    :rtype: str
    """
    if len(colorscale) < 1:
        raise ValueError("colorscale must have at least one color")

    if intermed <= 0 or len(colorscale) == 1:
        return colorscale[0][1]
    if intermed >= 1:
        return colorscale[-1][1]

    for cutoff, color in colorscale:
        if intermed > cutoff:
            low_cutoff, low_color = cutoff, color
        else:
            high_cutoff, high_color = cutoff, color
            break

    # noinspection PyUnboundLocalVariable
    return plotly.colors.find_intermediate_color(
        lowcolor=low_color, highcolor=high_color,
        intermed=((intermed - low_cutoff) / (high_cutoff - low_cutoff)),
        colortype="rgb")
```

## Theory

### A1

Part $(ii)$


```python
M = Manifold('M', dim=4)
P = Patch('origin', M)

t, r, theta, phi = sympy.symbols('t r theta phi')
A = sympy.Function('A')(t, r)
B = sympy.Function('B')(t, r)

cs = coords.CoordSystem('Schwarzschild', P, [t, r, theta, phi])
dt, dr, dtheta, dphi = cs.base_oneforms()

ds2 = (- sympy.exp(2 * A) * tpow(dt, 2) 
       + sympy.exp(2 * B) * tpow(dr, 2) 
       + r ** 2 * (tpow(dtheta, 2) + sympy.sin(theta) ** 2 * tpow(dphi, 2)))

g = metric.Metric(twoform=ds2)
g
```


```python
G_mu_nu = [(mu, nu, curvature.einstein_tensor_component(mu, nu, g)) 
           for mu, nu in itertools.product(range(2), range(2))]
```


```python
full_simplify(G_mu_nu[0][2].doit())
```


```python
full_simplify(G_mu_nu[1][2].doit())
```


```python
full_simplify(G_mu_nu[3][2].doit())
```


```python
full_simplify(curvature.christoffel_symbol_component(1, 0, 0, g))
```

### A2


```python
M = Manifold('M', dim=4)
P = Patch('origin', M)

v, r, theta, phi = sympy.symbols('v r theta phi')
alpha = sympy.Function('alpha')(v, r)
beta = sympy.Function('beta')(v, r)

cs = coords.CoordSystem('Eddington-Finkelstein', P, [v, r, theta, phi])
dv, dr, dtheta, dphi = cs.base_oneforms()

ds2 = (- alpha * tpow(dv, 2) 
       + beta * tprod(dv, dr) 
       + beta * tprod(dr, dv) 
       + r ** 2 * (tpow(dtheta, 2) + sympy.sin(theta) ** 2 * tpow(dphi, 2)))

g2 = metric.Metric(twoform=ds2)
g2
```


```python
g2.matrix
```


```python
G_mu_nu = [(mu, nu, curvature.einstein_tensor_component(mu, nu, g2)) 
           for mu, nu in itertools.product(range(2), range(2))]
```


```python
full_simplify(G_mu_nu[0][2].doit())
```


```python
full_simplify(G_mu_nu[1][2].doit())
```


```python
full_simplify(G_mu_nu[3][2].doit())
```

The above the says that $\beta = \beta(v)$, which can simplify $G_{vv}$ and $G_{vr}$


```python
beta_simp = sympy.Function('beta')(v)
ds2_simp = (- alpha * tpow(dv, 2) 
       + beta_simp * tprod(dv, dr) 
       + beta_simp * tprod(dr, dv) 
       + r ** 2 * (tpow(dtheta, 2) + sympy.sin(theta) ** 2 * tpow(dphi, 2)))

g2_simp = metric.Metric(twoform=ds2_simp)
```


```python
G_mu_nu_simp = [(mu, nu, curvature.einstein_tensor_component(mu, nu, g2_simp)) 
                 for mu, nu in itertools.product(range(2), range(2))]
```


```python
full_simplify(G_mu_nu_simp[0][2].doit())
```


```python
full_simplify(G_mu_nu_simp[1][2].doit())
```

The above gives equation for $\alpha$ as $G_{vr} = 0$
$$ \frac{1}{r^2 \beta} \left[ r\partial_r \alpha + \alpha - \beta^2 \right] = 0 \implies  r\partial_r \alpha + \alpha - \beta^2 = 0 \implies \alpha = \left[1 + \frac{f(v)}{r}\right] \beta^2(v)$$


```python
f = sympy.Function('f')(v)

# Second simplification to metric
alpha_simp = (1 + f / r) * beta_simp ** 2
ds2_simp2 = (- alpha_simp * tpow(dv, 2) 
             + beta_simp * tprod(dv, dr) 
             + beta_simp * tprod(dr, dv) 
             + r ** 2 * (tpow(dtheta, 2) + sympy.sin(theta) ** 2 * tpow(dphi, 2)))

g2_simp2 = metric.Metric(twoform=ds2_simp2)
```


```python
G_mu_nu_simp2 = [(mu, nu, curvature.einstein_tensor_component(mu, nu, g2_simp2)) 
                  for mu, nu in itertools.product(range(2), range(2))]
```


```python
full_simplify(G_mu_nu_simp2[0][2].doit())
```


```python
Phi = sympy.Function('Phi')(r)
ds2_temp = (- (1 + 2 * Phi) * tpow(dv, 2) 
             + tprod(dv, dr) 
             + tprod(dr, dv) 
             + r ** 2 * (tpow(dtheta, 2) + sympy.sin(theta) ** 2 * tpow(dphi, 2)))

g2_temp = metric.Metric(twoform=ds2_temp)
g2_temp
```


```python
g2_temp.matrix
```


```python
g2_temp.inverse.matrix
```


```python
tmp_00 = curvature.einstein_tensor_component(0, 0, g2_temp)
```


```python
full_simplify(tmp_00.doit())
```

### A3


```python
M = Manifold('M', dim=4)
P = Patch('origin', M)

t, r, theta, phi = sympy.symbols('t r theta phi')
A = sympy.Function('A')(r)
B = sympy.Function('B')(r)

cs = coords.CoordSystem('Schwarzschild', P, [t, r, theta, phi])
dt, dr, dtheta, dphi = cs.base_oneforms()

ds2 = (- sympy.exp(2 * A) * tpow(dt, 2) 
       + sympy.exp(2 * B) * tpow(dr, 2) 
       + r ** 2 * (tpow(dtheta, 2) + sympy.sin(theta) ** 2 * tpow(dphi, 2)))

g = metric.Metric(twoform=ds2)
g
```


```python
crs = curvature.compute_components(g, riemann=False, ricci=False)
```


```python
curvature.display_components(crs)
```


```python
full_simplify(sympy.exp(-2 * A) * curvature.einstein_tensor_component(0, 0, g).doit())
```


```python
full_simplify(sympy.exp(-2 * B) * curvature.einstein_tensor_component(1, 1, g).doit())
```


```python
full_simplify(sympy.exp(-2 * A) * curvature.einstein_tensor_component(0, 0, g).doit() +
              sympy.exp(-2 * B) * curvature.einstein_tensor_component(1, 1, g).doit())
```

## Exercises

### B2


```python
M = Manifold('M', dim=4)
P = Patch('origin', M)

t, r, theta, phi, r_s = sympy.symbols('t r theta phi r_s')

cs = coords.CoordSystem('Schwarzschild', P, [t, r, theta, phi])
dt, dr, dtheta, dphi = cs.base_oneforms()

schwarz = (1 - r_s / r)

ds2 = (- schwarz * tpow(dt, 2) 
       + (schwarz ** -1)  * tpow(dr, 2) 
       + r ** 2 * (tpow(dtheta, 2) + sympy.sin(theta) ** 2 * tpow(dphi, 2)))

g = metric.Metric(twoform=ds2)
g
```


```python
full_simplify(curvature.ricci_scalar(g, simplify_intermediate=False).doit())
```


```python
crs = curvature.compute_components(g, riemann=False, ricci=False)
```


```python
1e8 / (365.25 * 24 * 60 * 60)
```


```python
curvature.display_components(crs)
```


```python
def a_norm(r, r_s):
    return 0.5 * (r_s / (r ** 2) ) * (1 / numpy.sqrt(1 - (r_s / r)) )
```


```python
rs = numpy.arange(1 + 0.0001, 2 + 0.001, 0.0001)
df = pandas.DataFrame({'radius': rs, 'a_norm': a_norm(rs, 1)})
fig = px.line(df, x='radius', y='a_norm')
fig.update_layout(yaxis_range=[0,5],
                  xaxis_range=[0.8,2],
                  width=700,
                  height=500,
                  showlegend=False,
                  title_text=r'$\text{Stationary Observer Acceleration in Schwarzschild }(t, r)$', 
                  title_x=0.5,
                  xaxis_title=r'$r / r_s$',
                  yaxis_title=r'$r_s a \equiv r_s \sqrt{g_{\mu\nu}a^\mu a^\nu}$')
fig.show()
```


```python
ep = sympy.symbols('varepsilon')
a_ep = sympy.Rational(1, 2) * r_s / ( (r_s + ep) ** 2 ) * sympy.sqrt((r_s + ep) / ep)
a_ep
```


```python
sympy.series(a_ep, x=ep, n=0)
```

### B3


```python
r_s, r = sympy.symbols('r_s r', real=True, nonnegative=True)
dt_dr = 1 / (1 - r_s / r)
dt_dr
```

Ingoing


```python
sympy.integrate(dt_dr, r)
```


```python
sympy.integrate(-dt_dr, r)
```


```python
def ring(r, r_s, c):
    return -r - r_s * numpy.log(numpy.abs(r - r_s)) + c

def rong(r, r_s, c):
    return r + r_s * numpy.log(numpy.abs(r - r_s)) + c
```


```python
_r_s = 1.0
_ep = 1e-3
_r = numpy.arange(0*_r_s + _ep, 5 * _r_s + _ep, _ep)

_c = numpy.arange(-1.0, 2.0, 1.0)

dfs = []
for c in _c:
    df_i = pandas.DataFrame({'r': _r, 't': ring(_r, _r_s, c)}).assign(c=c, kind='In')
    df_o = pandas.DataFrame({'r': _r, 't': rong(_r, _r_s, c)}).assign(c=c, kind='Out')
    dfs.extend([df_i, df_o])
df = pandas.concat(dfs, axis=0)
```


```python
fig = px.line(df, x='r', y='t', color='c', line_dash='kind')
fig.update_layout(yaxis_range=[-5,5],
                  xaxis_range=[0,3],
                  width=800,
                  height=600,
                  showlegend=True,
                  title_text=r'$\text{Null Radial Geodesics in Schwarzschild }(t, r)$', 
                  title_x=0.5,
                  xaxis_title=r'$r / r_s$',
                  yaxis_title=r'$t(r)$')
fig.show()
```


```python
# fig.write_image('/Users/jim/repos/tex/homework/courses/PHYS510 - GR I/hw4/figures/fig-b3.pdf')
```

### B4


```python
def z(r_e, r_r, r_s = 1.0):
    return numpy.sqrt(1 - (r_s / r_r)) / numpy.sqrt(1 - (r_s / r_e)) - 1
```


```python
z(1.01, 10)
```


```python
r_e , r_r = numpy.meshgrid(numpy.arange(1.0, 2.0, 0.01),
                           numpy.arange(1.0, 2.0, 0.01))
zs = z(r_e, r_r)
```


```python
df = pandas.DataFrame({'r_e': r_e.ravel(), 'r_r': r_r.ravel(), 'z':zs.ravel()})
df = df[df['r_r'] >= df['r_e']]
df['logz'] = numpy.log(df['z'])
fig = px.scatter(df, x='r_e', y='r_r', color='logz')
fig_ep = 0.25
fig.update_layout(
#     yaxis_range=[1.0 - fig_ep, 5.0 + fig_ep],
#                   xaxis_range=[1.0 - fig_ep, 5.0 + fig_ep],
                  width=600,
                  height=600,
#                   showlegend=False,
                  title_text=r'$\text{Gravitational Redshift in Schwarzschild }(t, r)$', 
                  title_x=0.5,
                  xaxis_title=r'$r_E/r_s$',
                  yaxis_title=r'$r_R/r_s$')
fig.show()
```

### B5


```python
r, r_s, r_o, R = sympy.symbols('r r_s r_o R')
r_dot = sympy.sqrt(r_s / r - r_s / r_o)
r_dot
```


```python
sympy.integrate(-1/r_dot, (r, r_o, R))
```

Solved for t(r) by integrating... ugly : $a=r_s$, $b=r_o$, $x=r$

$$t(r) = \frac{1}{\sqrt{x} \sqrt{b-a}} \sqrt{\frac{x}{b-x}}\left(2 a^{3 / 2} \sqrt{b-x} \tanh ^{-1}\left(\frac{\sqrt{a} \sqrt{b-x}}{\sqrt{x} \sqrt{b-a}}\right)-\right. 
\left.\sqrt{b-a}\left(\sqrt{b}(2 a+b) \sqrt{1-\frac{x}{b}} \sin ^{-1}\left(\frac{\sqrt{x}}{\sqrt{b}}\right)+\sqrt{x}(x-b)\right)\right)$$


```python
# Previous Attempt
# 
# def t(r, r_s, r_o):
#     coeff = numpy.sqrt(numpy.abs(r_o - r_s) / r_s)
    
#     f1 = 1.0 / numpy.sqrt(r * numpy.abs(r_s - r_o))
#     f2 = numpy.sqrt(r / numpy.abs(r - r_o))
#     f3 = numpy.sqrt(numpy.abs(r_s - r_o))
    
#     g1 = numpy.sqrt(r_o) * (2 * r_s + r_o) * numpy.sqrt(1.0 - r / r_o) * numpy.arcsin(numpy.sqrt(r / r_o))
#     g2 = numpy.sqrt(r) * (r - r_o)
    
#     h1 = 2 * (r_s) ** (3.0/2) * numpy.sqrt(numpy.abs(r - r_o))
#     h2 = numpy.arctanh(numpy.sqrt(r_s * numpy.abs(r - r_o)) / numpy.sqrt(r * numpy.abs(r_s - r_o)))
    
#     return coeff * f1 * f2 * (h1 * h2 - f3 * (g1 + g2))
```


```python

```


```python
def t(r, r_s, r_o):
    c = numpy.sqrt(r * (r_o - r_s) / r_s)
    
    f1 = 

    
    return 
```


```python
_rs = numpy.arange(0.01, 3, 0.01)
_ts = t(_rs, 1.0, 30.0)
```


```python
df = pandas.DataFrame({'Time': _ts, 'Radius': _rs})
fig = px.line(df, x='Radius', y='Time')
fig_ep = 0.25
fig.update_layout(
#                   yaxis_range=[1.0 - fig_ep, 5.0 + fig_ep],
#                   xaxis_range=[1.0 - fig_ep, 5.0 + fig_ep],
                  width=800,
                  height=400,
#                   showlegend=False,
#                   title_text=r'$\text{Gravitational Redshift in Schwarzschild }(t, r)$', 
#                   title_x=0.5,
#                   xaxis_title=r'$r_E$',
#                   yaxis_title=r'$r_R$',
)
fig.show()
```

### B6


```python
def v_out(r, r_o, r_s = 1):
    return 2 * (r - r_o) + r_s * numpy.log(numpy.abs((r - r_s) / (r_o - r_s)))
```


```python
# Hack method
# _ros = numpy.arange(0.0, 3.0, 0.1)
# narrow_ep = 0.01
# _ros = numpy.concatenate((_ros, numpy.arange(1.0 - narrow_ep, 1.0 + narrow_ep, narrow_ep / 2)))

# Wrinkly method
# _quants = numpy.arange(0.34, .977, 0.03)
_quants = numpy.arange(0.25, .97, 0.03)
_ros = stats.norm.ppf(_quants, 1, .5)


df = []
for _ro in _ros:
    _rs = numpy.arange(0.01, 2.0, 0.0001)
    _vo = v_out(_rs, r_o=_ro)
    
    if _ro > 1:
        _vo = numpy.where(_rs < _ro, numpy.full(_rs.shape, numpy.nan), _vo)
    if _ro < 1:
        _vo = numpy.where(_rs > _ro, numpy.full(_rs.shape, numpy.nan), _vo)
    
    df.append(pandas.DataFrame({'v': _vo, 'r': _rs}).assign(r_o=_ro))
df = pandas.concat(df, axis=0)
```


```python
fig = px.scatter(df, x='r', y='v', color='r_o', color_continuous_scale=px.colors.sequential.Agsunset)

# fig.update_layout(
#               yaxis_range=[0.0, 2.0],
#               xaxis_range=[0.0, 2.0],
#                   width=700,
#                   height=700,
# #                   showlegend=False,
#                   title_text=r'$\text{Outgoing Null Geodesics in Schwarzschild }(v, r)$', 
#                   title_x=0.5,
#                   xaxis_title=r'$r / r_s$',
#                   yaxis_title=r'$v$')
# fig.show()
```


```python
# fig.write_image('/Users/jim/repos/tex/homework/courses/PHYS510 - GR I/hw4/figures/fig-b6.pdf')
```

### B7


```python

```


```python

```

### B8


```python

```


```python

```


```python

```

## Problem


```python

```

3f7b579162db704334152e5c9a3753b0369b04f1
ipynb
Jupyter Notebook
notebooks/homework-4.ipynb
JWKennington/CourseGR1
ba81a1bf0427210c14d4e97c8f2cf8a345297812
MIT
notebooks/homework-4.ipynb
JWKennington/CourseGR1
ba81a1bf0427210c14d4e97c8f2cf8a345297812
MIT
notebooks/homework-4.ipynb
JWKennington/CourseGR1
ba81a1bf0427210c14d4e97c8f2cf8a345297812
MIT
Qwen/Qwen-72B
1. YES
2. YES