Analysing Earth using EinsteinPy!

Importing required modules

import numpy as np
from astropy import units as u

from einsteinpy.coordinates import SphericalDifferential
from einsteinpy.geodesic import Timelike
from einsteinpy.metric import Schwarzschild
from einsteinpy.plotting import GeodesicPlotter

Defining the Metric and Initial Parameters

# Define position and velocity vectors in spherical coordinates
# Earth is at perihelion
M = 1.989e30 * u.kg  # Mass of Sun
distance = 147.09e6 * u.km
speed_at_perihelion = 30.29 * u.km / u.s
omega = (speed_at_perihelion.value) / distance.value * u.rad / u.s # Angular Speed

# Defining the initial coordinates of the test particle
# in SI
sph = SphericalDifferential(
    t=0. * u.s,
    r=distance,
    theta=np.pi / 2 * u.rad,
    phi=np.pi * u.rad,
    v_r=0. * u.m / u.s,
    v_th=0. * u.rad / u.s,
    v_p=omega,
)

# Schwarzschild Metric Object
ms = Schwarzschild(coords=sph, M=M)

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.

# Set lambda to complete an year
# Lambda is always specified in seconds
end_lambda = ((1 * u.year).to(u.s)).value
# Choosing stepsize for ODE solver to be 5 minutes
step_size = ((5 * u.min).to(u.s)).value

Calculating Geodesic

geod = Timelike(metric=ms, coords=sph, end_lambda=end_lambda, step_size=step_size)
geod

Geodesic Object:
            Type = (Time-like)
            Metric = ((
            Name: (Schwarzschild Metric),
            Coordinates: (Spherical Polar Coordinates:
            t = (0.0 s), r = (147090000.0 km), theta = (1.5707963267948966 rad), phi = (3.141592653589793 rad)
            v_t: 1.0000000151460953 m / s, v_r: 0.0 m / s, v_th: 0.0 rad / s, v_p: 2.0592834319124346e-07 rad / s),
            Mass: (1.989e+30 kg),
            Spin parameter: (0.0),
            Charge: (0.0 C),
            Schwarzschild Radius: (2954.126555055405)
)),
            Initial State = ([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.15506405e+07 -1.47089888e+11 -1.82615574e+08 ...  3.69925133e+01
  -3.02899770e+04 -1.16558621e-11]
 [ 3.15530405e+07 -1.47089782e+11 -2.55311495e+08 ...  5.17185190e+01
  -3.02899551e+04 -1.16567636e-11]
 [ 3.15554405e+07 -1.47089640e+11 -3.28007353e+08 ...  6.64445117e+01
  -3.02899259e+04 -1.16576651e-11]])

ans = geod.trajectory

ans[0].shape, ans[1].shape

((8,), (8,))

Calculating distance and speed at aphelion

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

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

$1.5204946 \times 10^{8} \; \mathrm{km}$

• Speed should be ~\(29.29\) km/s and should be along Y-axis

((ans[i][6]) * u.m / u.s).to(u.km / u.s)

$29.302018 \; \mathrm{\frac{km}{s}}$

Calculating the eccentricity of the orbit

• Eccentricity should be ~\(0.0167\)

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

0.016709193287269494

Plotting the trajectory with time

sgp = GeodesicPlotter()
sgp.plot(geod)
sgp.show()

Data type cannot be displayed: application/vnd.plotly.v1+json

All data regarding earth’s orbit is taken from https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html