Analysing Earth using EinsteinPy!

Calculating the eccentricity and speed at apehelion of Earth’s orbit

Various parameters of Earth’s orbit considering Sun as schwarzschild body and solving geodesic equations are calculated

1. Defining the initial parameters

[1]:

from astropy import units as u
import numpy as np
from einsteinpy.metric import Schwarzschild

[2]:

# Define position and velocity vectors in spherical coordinates
# Earth is at perihelion
M = 1.989e30 * u.kg  # mass of sun
pos_vec = [147.09e6 * u.km, np.pi / 2 * u.rad, np.pi * u.rad]
speed_at_perihelion = 30.29 * u.km / u.s
omega = (u.rad * speed_at_perihelion) / pos_vec[0]
vel_vec = [0 * u.km / u.s, 0 * u.rad / u.s, omega]

  • Defining \(\lambda\) (or \(\tau\)) for which to calculate trajectory

    • \(\lambda\) is proper time and is approximately equal to coordinate time in non-relativistic limits

[3]:

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

2. Making a class instance to get the trajectory in cartesian coordinates

[4]:

starting_time = 0 * u.s
obj = Schwarzschild.from_spherical(pos_vec, vel_vec, starting_time, M)
ans = obj.calculate_trajectory(
    end_lambda=end_lambda, OdeMethodKwargs={"stepsize": stepsize}, return_cartesian=True
)

/home/shreyas/Softwares/anaconda3/lib/python3.6/site-packages/scipy/integrate/_ivp/common.py:32: UserWarning: The following arguments have no effect for a chosen solver: `first_step`.
  .format(", ".join("`{}`".format(x) for x in extraneous)))

  • Results are returned in SI units in cartesian coordinates

[5]:

obj.units_list

[5]:

[Unit("s"),
 Unit("m"),
 Unit("m"),
 Unit("m"),
 Unit(dimensionless),
 Unit("m / s"),
 Unit("m / s"),
 Unit("m / s")]

  • Return value is a tuple consisting of 2 numpy array

    • First one contains list of \(\lambda\)

    • Seconds is array of shape (n,8) where each component is:

    • t - coordinate time

    • x - position in m

    • y - position in m

    • z - position in m

    • dt/d\(\lambda\)

    • dx/d\(\lambda\)

    • dy/d\(\lambda\)

    • dz/d\(\lambda\)

[6]:

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

[6]:

((13164,), (13164, 8))

Calculating distance at apehelion

  • Should be 152.10 million km

[7]:

r = np.sqrt(np.square(ans[1][:, 1]) + np.square(ans[1][:, 2]))
i = np.argmax(r)
(r[i] * u.m).to(u.km)

[7]:
$1.5205967 \times 10^{8} \; \mathrm{km}$

Speed at perihelion should be 29.29 km/s and should be along y-axis

[8]:

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

[8]:
$29.300051 \; \mathrm{\frac{km}{s}}$

Calculating the eccentricity

  • Should be 0.0167

[9]:

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

[9]:

0.01648174751185976

Plotting the trajectory with time

[10]:

%matplotlib inline
# Changing Backend

[11]:

from einsteinpy.plotting import StaticGeodesicPlotter
sgp = StaticGeodesicPlotter(M)
sgp.plot(pos_vec, vel_vec, end_lambda, stepsize)

/home/shreyas/Softwares/anaconda3/lib/python3.6/site-packages/scipy/integrate/_ivp/common.py:32: UserWarning: The following arguments have no effect for a chosen solver: `first_step`.
  .format(", ".join("`{}`".format(x) for x in extraneous)))

[11]:

[<matplotlib.lines.Line2D at 0x7fdecba49ef0>,
 <matplotlib.lines.Line2D at 0x7fdecba5d470>]
../_images/examples_Analysing_Earth_using_EinsteinPy!_20_2.png

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