# HW 1

## Imports


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


```python
import itertools
from astropy import units
import sympy
import pandas
from scipy import integrate
import numpy
import plotly.express as px

from sympy.diffgeom import Manifold, Patch
from pystein import coords, metric, curvature, geodesic
from pystein.utilities import tensor_pow as tpow, full_simplify, boundary_filter
```


```python
show_plots = True
```

## Exercises

### B1 - Polar Coords for $\mathbb{R}^2$

#### Setup Metric


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

##### Cartesian


```python
x, y = sympy.symbols('x, y', nonnegative=False)
cs = coords.CoordSystem('cartesian', P, [x, y])
dx, dy = cs.base_oneforms()
ds2 = tpow(dx, 2) + tpow(dy, 2)
g_cart = metric.Metric(twoform=ds2)
g_cart
```

##### Polar


```python
r, theta = sympy.symbols('r theta', nonnegative=True)
cs = coords.CoordSystem('polar', P, [r, theta])
dr, dtheta = cs.base_oneforms()
ds2 = tpow(dr, 2) + r ** 2 * tpow(dtheta, 2)
g1 = metric.Metric(twoform=ds2)
g1
```

#### Compute curvature components


```python
crs, rmns, rcs = curvature.compute_components(g1)
```


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


```python
curvature.display_components(rmns)
```


```python
curvature.display_components(rcs)
```

#### Parallel Transport


```python
# Unit circle parameterized such that x(0) == x(1)
param = sympy.symbols('lambda')
r_0 = sympy.symbols('r_0', real=True)
curve = [
    r_0,
    2 * sympy.pi * param,
]
```


```python
soln = geodesic.parallel_transport_soln(param, curve, g1)
```


```python
soln.soln[0]
```


```python
g1.angle(soln.vec(0), soln.vec(1))
```

#### Compute Geodesics


```python
full_simplify(geodesic.geodesic_equation(0, sympy.symbols('lambda'), g1))
```


```python
full_simplify(geodesic.geodesic_equation(1, sympy.symbols('lambda'), g1))
```

##### Visualize Cartesian


```python
ls = numpy.arange(0, 2, 0.01)
init = (0.0, 0.0, 0.1, 0.1) 
df1cart = geodesic.numerical_geodesic(g_cart, init, ls)
```


```python
ls = numpy.arange(0, 2, 0.001)
df1scart = geodesic.numerical_sampler(g_cart, ls, (0, 0), num_angles=6)
```


```python
if show_plots:
    fig = px.scatter(df1scart, x="x", y='y', 
                     color='theta_0',
                     title="B1 Cart Geodesic",
                     height=600, width=600)

    fig.show()
```

##### Visualize Polar


```python
ls = numpy.arange(0, 2, 0.001)
r_0 = 1.0 
dfs = [geodesic.numerical_geodesic(g1, (r_0, theta_0, 0.0, numpy.pi/2), ls).assign(theta_0=theta_0) 
       for theta_0 in numpy.arange(0.0, 2 * numpy.pi, 0.2)]
df = pandas.concat(dfs, axis=0)
```


```python
df['theta'] = numpy.mod(df['theta'], 2*numpy.pi)
# df['lam'] = ls
```


```python
if show_plots:
    fig = px.scatter(df, x="theta", y='r', 
                     color='theta_0',
                     title="Sample Geodesic")

    fig.show()
```

### B2 2-Sphere

#### Setup Metric


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

theta, phi, a = sympy.symbols('theta phi a', nonnegative=True)
cs = coords.CoordSystem('spherical', P, [theta, phi])
dtheta, dphi = cs.base_oneforms()
ds2 = a**2 * (tpow(dtheta, 2) + sympy.sin(theta)**2 * tpow(dphi, 2))
g2 = metric.Metric(twoform=ds2)
g2
```

#### Compute Components


```python
g2.matrix
```


```python
g2.inverse.matrix
```


```python
crs, rmn, rcs = curvature.compute_components(g2)
```


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


```python
curvature.display_components(rmn)
```


```python
curvature.display_components(rcs)
```


```python
full_simplify(curvature.ricci_scalar(metric=g2))
```

#### Parallel Transport


```python
param, theta_0 = sympy.symbols('lambda theta_0')

curve = [
    theta_0,
    param,
]
```


```python
geodesic.parallel_transport_equation(0, curve, param, g2)
```


```python
full_simplify(geodesic.parallel_transport_equation(1, curve, param, g2))
```

#### Geodesics

##### Geodesic Equations


```python
full_simplify(geodesic.geodesic_equation(0, sympy.symbols('lambda'), g2))
```


```python
full_simplify(geodesic.geodesic_equation(1, sympy.symbols('lambda'), g2))
```

##### Single Geodesic


```python
ls = numpy.arange(0, 8, 0.0001)
init = (0.01, 0.001, 3.14/4, 0.0)
st = 2
sp = 2
df2 = geodesic.numerical_geodesic(g2, (numpy.pi/2, 0.0, numpy.pi/st, numpy.pi/sp), ls)
#        for theta_0 in numpy.arange(0.01, 3.14, 0.5)]
```


```python
df2['theta'] = numpy.mod(df2['theta'], 2*numpy.pi)
df2['phi'] = numpy.mod(df2['phi'], numpy.pi)
df2 = df2.reset_index().rename(columns={'index': 'order'})
df2['lam'] = ls
```


```python
if show_plots:
    fig = px.scatter(df2, x="theta", y='phi', 
                     color='order',
                     title="B2 Geodesic")

    fig.show()
```

##### Multiple Initial Conditions


```python
ls = numpy.arange(0, .04, 0.00001)
st = 2
sp = .01
dfs2 = [geodesic.numerical_geodesic(g2, (theta_0, 0.0, numpy.pi/st, numpy.pi/sp), ls).assign(theta_0=theta_0) 
       for theta_0 in list(numpy.arange(numpy.pi/6, 5*numpy.pi/6, 0.1)) + 
       list(numpy.arange(numpy.pi/6 + numpy.pi, 5*numpy.pi/6 + numpy.pi, 0.1))]
#        for theta_0 in numpy.arange(0.01, 3.14, 0.5)]
df2b = pandas.concat(dfs2, axis=0)
```


```python
df2b['theta'] = numpy.mod(df2b['theta'], 2*numpy.pi)
df2b = df2b.reset_index().rename(columns={'index': 'order'})
```


```python
lim_df2b = boundary_filter(df2b, theta=(0, 2*3.14), phi=(0, 3.14))
```


```python
if show_plots:
    fig = px.scatter(lim_df2b, x="theta", y='phi', 
                     color='theta_0',
                     title="B2 Geodesics")

    fig.show()
```

### B3 2-Sphere Sch

#### Setup Metric


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

rho, phi, a = sympy.symbols('rho phi a', nonnegative=True)
cs = coords.CoordSystem('schw', P, [rho, phi])
drho, dphi = cs.base_oneforms()
ds2 = a**2 * ( (1 / (1 - rho**2)) * tpow(drho, 2) + rho ** 2 * tpow(dphi, 2))
g3 = metric.Metric(twoform=ds2)
g3
```

#### Geodesics


```python
full_simplify(geodesic.geodesic_equation(0, sympy.symbols('lambda'), g3))
```


```python
full_simplify(geodesic.geodesic_equation(1, sympy.symbols('lambda'), g3))
```


```python
# init = (numpy.sin(numpy.pi/4), 0.0, numpy.cos(numpy.pi/4), numpy.pi/4)
init = (numpy.sin(0), 0.0, numpy.cos(numpy.pi/4), numpy.pi/4)
lambdas = numpy.arange(0, 2.1, 0.001)
df3 = geodesic.numerical_geodesic(g3, init, lambdas)
```


```python
df3['phi'] = numpy.mod(df3['phi'], numpy.pi)
df3 = df3.reset_index().rename(columns={'index': 'order'})
df3['lam'] = lambdas
```


```python
if show_plots:
    fig = px.scatter(df3, x="rho", y='phi', 
    #                  color='init',
                         title="B3 Geodesic")

    fig.show()
```

### B4 - Embedded Surface - General Case

#### Setup Metric


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

x, y = sympy.symbols('x y', nonnegative=False)
_coords = [x, y]
cs = coords.CoordSystem('Cartesian', P, [x, y])
dx, dy = cs.base_oneforms()
_1forms = [dx, dy]

f = sympy.Function('f')(x, y)

ds2 = []
for i in range(2):
    for j in range(2):
        ds2.append(((1 if i == j else 0) + sympy.diff(f, _coords[i]) * sympy.diff(f, _coords[j])) 
                   * sympy.diffgeom.TensorProduct(_1forms[i],  _1forms[j]))
ds2 = sum(ds2)
g4 = metric.Metric(twoform=ds2)
g4
```

#### Compute Components


```python
g4.matrix
```


```python
sympy.det(g4.matrix)
```


```python
g4.inverse.matrix
```


```python
crs, rms, rcs = curvature.compute_components(g4)
```


```python
# Simplify notation
subscript_notation = [
    (sympy.diff(sympy.diff(f, x), x), sympy.symbols('f_xx')),
    (sympy.diff(sympy.diff(f, y), y), sympy.symbols('f_yy')),
    (sympy.diff(sympy.diff(f, x), y), sympy.symbols('f_xy')),
    (sympy.diff(f, x), sympy.symbols('f_x')),
    (sympy.diff(f, y), sympy.symbols('f_y')),
]
```


```python
crs = [(c[0], c[1].subs(subscript_notation)) for c in crs]
rms = [(c[0], c[1].subs(subscript_notation)) for c in rms]
rcs = [(c[0], c[1].subs(subscript_notation)) for c in rcs]
```


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


```python
curvature.display_components(rms)
```


```python
curvature.display_components(rcs)
```

Clean Expr


```python
simplified = full_simplify(curvature.ricci_scalar(g4)).subs(subscript_notation)
simplified
```

#### Fun Examples

Paraboloid


```python
a, b = sympy.symbols('a b')
subs = {f: x**2/a + y**2/b}
```


```python
full_simplify(full_simplify(curvature.ricci_scalar(g4)).subs(subs))
```

### B5 - Embedded Surface - Sphere

#### Substitution


```python
a_0 = sympy.symbols('a_0')
subs = {f: sympy.sqrt(a_0**2 - x**2 - y**2)}
g5 = metric.Metric(twoform=g4.twoform.subs(subs))
g5
```

#### Curvature


```python
full_simplify(curvature.ricci_scalar(g5, simplify_intermediate=True))
```

#### Parallel Transport


```python
param = sympy.symbols('\\lambda')

curve = [
    sympy.cos(2 * sympy.pi * param),
    sympy.sin(2 * sympy.pi * param),
]
```


```python
full_simplify(geodesic.parallel_transport_equation(0, curve, param, g5))
```


```python
full_simplify(geodesic.parallel_transport_equation(1, curve, param, g5))
```

#### Geodesic

##### Geodesic Equations


```python
full_simplify(geodesic.geodesic_equation(0, sympy.symbols('lambda'), g5))
```


```python
full_simplify(geodesic.geodesic_equation(1, sympy.symbols('lambda'), g5))
```

##### Visual


```python
a_0_val = 1.0
g5_num = metric.Metric(twoform=g5.twoform.subs({a_0: a_0_val}))
```


```python
ls = numpy.arange(0.0, 8.0, 0.001)
df5 = geodesic.numerical_sampler(g5_num, ls, (0.5, 0.5), tangent_scale=0.1)
```


```python
lim_df5 = boundary_filter(df5, x=(-1, 1), y=(-1, 1))
```


```python
if show_plots:
    fig = px.scatter(lim_df5, x="x", y='y', 
                     color='theta_0',
                     title="B5 Geodesic", 
                     height=600, width=600)
    fig.show()
```

### B6 - Embedded Surface - Hyperboloid

#### Substitution


```python
a_0 = sympy.symbols('a_0')
subs = {f: sympy.sqrt(a_0**2 - x**2 + y**2)}
g6 = metric.Metric(twoform=g4.twoform.subs(subs))
g6
```


```python
g6.matrix.doit()
```

#### Curvature


```python
full_simplify(curvature.ricci_scalar(g6, simplify_intermediate=True))
```

#### Parallel Transport


```python
full_simplify(geodesic.parallel_transport_equation(0, curve, param, g6))
```


```python
full_simplify(geodesic.parallel_transport_equation(0, curve, param, g6))
```

#### Geodesic


```python
a_0_val = 1.0
g6_num = metric.Metric(twoform=g6.twoform.subs({a_0: a_0_val}))
```


```python
ls = numpy.arange(0.0, 8.0, 0.001)
df6 = geodesic.numerical_sampler(g6_num, ls, (0.5, 0.5), tangent_scale=0.1)
```


```python
lim_df6 = boundary_filter(df6, x=(-1, 1), y=(-1, 1))
```


```python
if show_plots:
    fig = px.scatter(lim_df6, x="x", y='y', 
                     color='theta_0',
                     title="B6 Geodesic",
                     height=600, width=600)

    fig.show()
```

### B7 - Embedded Surface - Cylinder

#### Substitution


```python
a_0 = sympy.symbols('a_0')
subs = {f: sympy.sqrt(a_0**2 - y**2)}
g7 = metric.Metric(twoform=g4.twoform.subs(subs))
g7
```


```python
g7.matrix.doit()
```

#### Curvature


```python
crs, rmn, rcs = curvature.compute_components(g7)
```


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


```python
full_simplify(curvature.ricci_scalar(g7))
```

#### Parallel Transport


```python
full_simplify(geodesic.parallel_transport_equation(0, curve, param, g7))
```


```python
full_simplify(geodesic.parallel_transport_equation(1, curve, param, g7))
```

#### Geodesic


```python
a_0_val = 1.0
g7_num = metric.Metric(twoform=g7.twoform.subs({a_0: a_0_val}))
```


```python
ls = numpy.arange(0.0, 8.0, 0.001)
df7 = geodesic.numerical_sampler(g7_num, ls, (0.5, 0.5), tangent_scale=0.1)
```


```python
lim_df7 = boundary_filter(df7, x=(-1, 1), y=(-1, 1))
```


```python
if show_plots:
    fig = px.scatter(lim_df7, x="x", y='y', 
                     color='theta_0',
                     title="B7 Geodesic",
                     height=600, width=600)
    fig.show()
```

### B8 - Embedded Surface - Cone

#### Substitution


```python
a_0 = sympy.symbols('a_0')
subs = {f: sympy.sqrt(a_0 * (x**2 + y**2))}
g8 = metric.Metric(twoform=g4.twoform.subs(subs).doit())
g8
```


```python
full_simplify(g8.matrix)
```

#### Curvature


```python
crs, rmn, rcs = curvature.compute_components(g8)
```


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


```python

```


```python
full_simplify(curvature.ricci_scalar(g8))
```

#### Parallel Transport


```python
full_simplify(geodesic.parallel_transport_equation(0, curve, param, g8))
```


```python
full_simplify(geodesic.parallel_transport_equation(1, curve, param, g8))
```

#### Geodesic


```python
a_0_val = 1.0
g8_num = metric.Metric(twoform=g8.twoform.subs({a_0: a_0_val}))
```


```python
ls = numpy.arange(0.0, 8.0, 0.001)
df8 = geodesic.numerical_sampler(g8_num, ls, (0.5, 0.5), tangent_scale=0.1)
```


```python
lim_df8 = boundary_filter(df8, x=(-1, 1), y=(-1, 1))
```


```python
if show_plots:
    fig = px.scatter(lim_df8, x="x", y='y', 
                     color='theta_0',
                     title="B8 Geodesic", 
                     height=600, width=600)

    fig.show()
```

## Applications

### C1 - Lower bound on size of universe


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

r, theta, phi = sympy.symbols('r theta phi', nonnegative=True)
a_0 = sympy.symbols('a_0', real=True, nonnegative=True)

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

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

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


```python
gc1.matrix
```


```python
R_0 = 1.7e-54 * (1 / units.meter) ** 2
```


```python
a0_val = numpy.sqrt(6 / R_0)
```


```python
V_bound = 2 * numpy.pi ** 2 * (a0_val) ** (3)
```


```python
V_bound
```

### C2 - Radial Geodesic

$$ L = \int_{\gamma} \sqrt{g_{\mu\nu}\frac{d x^{\mu}}{d\lambda}\frac{d x^{\nu}}{d\lambda}} d \lambda$$

Assume that radial geodesic implies that $$\frac{d}{d\lambda}\theta = \frac{d}{d\lambda}\phi = 0$$

$$ L_0 = \int_{\gamma} \sqrt{g_{rr}\left(\frac{d r}{d\lambda}\right)^2} d \lambda = \int_{0}^{r_0} \sqrt{g_{rr}}d r $$


```python
r_0 = sympy.symbols('r_0', real=True, nonnegative=True)
```


```python
soln = sympy.integrate(sympy.sqrt(gc1.matrix[0,0]), (r, 0, r_0))
soln
```


```python
soln.args[0][0]
```

$$ r_0 = \sin\left(\frac{L_0}{a_0}\right)$$


```python
L_0 = (46e9 * units.lightyear).to(units.meter)
L_0
```


```python
a0_val
```


```python
L_0 / a0_val
```


```python
r0_val = numpy.sin((L_0/a0_val).value)
r0_val
```


```python
v_obs = 4 / 3 * numpy.pi * r0_val ** 3
v_obs_scaled = v_obs * (a0_val) ** 3
v_obs_scaled
```


```python
v_obs_scaled / V_bound
```

88ecb2812829eb8eb730d88499f025e0e7f83db9
ipynb
Jupyter Notebook
notebooks/homework-1.ipynb
JWKennington/CourseGR1
ba81a1bf0427210c14d4e97c8f2cf8a345297812
MIT
notebooks/homework-1.ipynb
JWKennington/CourseGR1
ba81a1bf0427210c14d4e97c8f2cf8a345297812
MIT
notebooks/homework-1.ipynb
JWKennington/CourseGR1
ba81a1bf0427210c14d4e97c8f2cf8a345297812
MIT
Qwen/Qwen-72B
1. YES
2. YES