Analysing Earth using EinsteinPy!

Importing required modules

[1]:
import numpy as np
from astropy import units as u

from einsteinpy.coordinates.utils import four_position, stacked_vec
from einsteinpy.geodesic import Geodesic
from einsteinpy.metric import Schwarzschild
from einsteinpy.plotting import GeodesicPlotter

Defining the Metric and Initial Parameters

[2]:
# Define position and velocity vectors in spherical coordinates
# Earth is at perihelion
t = 0.
M = 1.989e30  # Mass of Sun
distance = 147.09e9 # m
speed_at_perihelion = 30290 # m/s
omega = (speed_at_perihelion) / distance # Angular Speed

x_vec = np.array([distance, np.pi / 2, np.pi])  # 3-Pos
v_vec = np.array([0.0, 0.0, omega])  # 3-Vel

# Schwarzschild Metric Object
ms_cov = Schwarzschild(M=M)
# Getting Position 4-Vector
x_4vec = four_position(t, x_vec)
# Calculating Schwarzschild Metric at x_4vec
ms_cov_mat = ms_cov.metric_covariant(x_4vec)
# Getting stacked (Length-8) initial vector, containing 4-Pos and 4-Vel
init_vec = stacked_vec(ms_cov_mat, t, x_vec, v_vec, time_like=True)

Defining constraints on Affine Parameter \(\lambda\), that parameterizes the geodesic

For timelike geodesics (like here), \(\lambda\) can be taken to be proper time (\(\tau\)) and it becomes approximately equal to Coordinate Time (\(t\)) in non-relativistic limits.

[3]:
# Set lambda to complete an year.
# Lambda is always specified in secs
end_lambda = 3.154e7
# Choosing stepsize for ODE solver to be 5 minutes
step_size = 300

Calculating Geodesic

[4]:
geod = Geodesic(metric=ms_cov, init_vec=init_vec, end_lambda=end_lambda, step_size=step_size)
geod
[4]:
Geodesic Object:

Metric = ((
Name: (Schwarzschild Metric),
Coordinates: (S),
Mass: (1.989e+30),
Spin parameter: (0),
Charge: (0),
Schwarzschild Radius: (2954.126555055405)
)),

Initial Vector = ([0.00000000e+00 1.47090000e+11 1.57079633e+00 3.14159265e+00
 1.00000002e+00 0.00000000e+00 0.00000000e+00 2.05928343e-07]),

Trajectory = ([[ 0.00000000e+00 -1.47090000e+11  1.80133298e-05 ... -3.70945515e-12
  -3.02900000e+04  0.00000000e+00]
 [ 2.40000004e+02 -1.47090000e+11 -7.26960000e+06 ...  1.47260227e+00
  -3.02900000e+04 -9.01708829e-17]
 [ 2.64000004e+03 -1.47089979e+11 -7.99655961e+07 ...  1.61986242e+01
  -3.02899956e+04 -9.91879708e-16]
 ...
 [ 3.15338405e+07 -1.47089644e+11  3.26256116e+08 ... -6.60895799e+01
  -3.02899267e+04 -1.16495498e-11]
 [ 3.15362405e+07 -1.47089785e+11  2.53560256e+08 ... -5.13635868e+01
  -3.02899557e+04 -1.16504517e-11]
 [ 3.15386405e+07 -1.47089891e+11  1.80864335e+08 ... -3.66375809e+01
  -3.02899775e+04 -1.16513535e-11]])
[5]:
ans = geod.trajectory

ans[0].shape, ans[1].shape
[5]:
((8,), (8,))

Calculating distance and speed at aphelion

  • Distance should be ~\(152.10\) million km

[6]:
r = np.sqrt(np.square(ans[:, 1]) + np.square(ans[:, 2]))
i = np.argmax(r)
(r[i] * u.m).to(u.km) # Converting m to km
[6]:
$1.5204946 \times 10^{8} \; \mathrm{km}$
  • Speed should be ~\(29.29\) km/s and should be along Y-axis

[7]:
((ans[i, 6]) * u.m / u.s).to(u.km / u.s)
[7]:
$29.302018 \; \mathrm{\frac{km}{s}}$

Calculating the eccentricity of the orbit

  • Eccentricity should be ~\(0.0167\)

[8]:
xlist, ylist = ans[:, 1], ans[:, 2]
i = np.argmax(ylist)
x, y = xlist[i], ylist[i]
eccentricity = x / (np.sqrt(x ** 2 + y ** 2))
eccentricity
[8]:
0.016709193287269494

Plotting the trajectory with time

[9]:
sgp = GeodesicPlotter()
sgp.plot(geod)
sgp.show()