Source code for einsteinpy.hypersurface.schwarzschildembedding
import warnings
import numpy as np
from astropy import units as u
from matplotlib import pyplot as plt
[docs]class SchwarzschildEmbedding:
"""
Class for Utility functions for Schwarzschild Embedding surface to
implement gravitational lensing
Attributes
----------
input_units : list
list of input units of M
units_list : list
customized units to handle values of M and render plots
within grid range
r_init : ~astropy.units.m
"""
def __init__(self, M):
"""
Constructor
Initialize mass and embedding initial radial coordinate in appropiate units
in order to render the plots of the surface in finite grid. The initial r
is taken to be just greater than schwarzschild radius but it is important
to note that the embedding breaks at r < 9m/4.
Parameters
----------
M : ~astropy.units.kg
Mass of the body
"""
self.input_units = [M.unit]
self.units_list = [u.kg * 10e22, u.m / M.to(u.kg * 10e22).value]
M = M.to(self.units_list[0])
self.M = M
self.r_init = (((3 * self.M.value + 0.0001) / self.M.value) * u.m).to(
self.units_list[1]
)
[docs] def gradient(self, r):
"""
Calculate gradient of Z coordinate w.r.t r to update the value of r and
thereby get value of spherical radial coordinate R.
Parameters
----------
r : float
schwarzschild coordinate at which gradient is supposed to be obtained
Returns
-------
float
gradient of Z w.r.t r at the point r (passed as argument)
"""
R = r / np.sqrt(1 - (2 * self.M.value / r))
num_one = 1 - (3 * self.M.value / r)
num_two = np.sqrt(
((4 * self.M.value * r - 9 * self.M.value * self.M.value) * R)
/ (r - 3 * self.M.value) ** 2
)
deno = np.sqrt(1 - (2 * self.M.value / r)) ** 3
return num_one * num_two / deno
[docs] def radial_coord(self, r):
"""
Returns spherical radial coordinate (of the embedding) from given schwarzschild
coordinate.
Parameters
----------
r : float
Returns
-------
float
spherical radial coordinate of the 3d embedding
"""
return r / np.sqrt(1 - (2 * self.M.value / r))
[docs] def get_values(self, alpha):
"""
Obtain the Z coordinate values and corrosponding R values for range of
r as 9m/4 < r < 9m.
Parameters
----------
alpha : float
scaling factor to obtain the step size for incrementing r
Returns
-------
tuple
(list, list) : values of R (x_axis) and Z (y_axis)
"""
x_axis = []
y_axis = []
r_initial = self.r_init.value
r_step = self.M.value / alpha
z = 0
r = r_initial
while r < 9 * self.M.value:
x_axis.append(self.radial_coord(r))
y_axis.append(z)
z = z + self.gradient(r) * r_step
r = r + r_step
z = 0
r = r_initial
while r > (9 * self.M.value / 4):
x_axis.append(self.radial_coord(r))
y_axis.append(z)
z = z + self.gradient(r) * r_step
r = r - r_step
return x_axis, y_axis
[docs] def get_values_surface(self, alpha):
"""
Obtain the same values as of the get_values function but reshapes them to obtain
values for all points on the solid of revolution about Z axis (as the
embedding is symmetric in angular coordinates).
Parameters
----------
alpha : float
scaling factor to obtain the step size for incrementing r
Returns
-------
tuple
(~numpy.array of X, ~numpy.array of Y, ~numpy.array of Z) values in cartesian coordinates
obtained after applying solid of revolution
"""
r_initial = self.r_init.value
r_step = self.M.value / alpha
phi_values = np.linspace(0, 2 * np.pi, 60)
R_values = []
z_values = []
z = 0
r = r_initial
while r < 20 * self.M.value:
R_values.append(self.radial_coord(r))
z_values.append(z)
z = z + self.gradient(r) * r_step
r = r + r_step
R_values = np.array(R_values)
R_values, phi_values = np.meshgrid(R_values, phi_values)
X = R_values * np.cos(phi_values)
Y = R_values * np.sin(phi_values)
x_len = X.shape[0]
Z = np.array(z_values)
z_values = np.array(z_values)
for i in range(0, x_len - 1):
Z = np.concatenate((Z, z_values), axis=0)
Z.reshape((X.shape[0], X.shape[1]))
return X, Y, Z