from astropy import units as u
import numpy as np
from einsteinpy.metric import Schwarzschild
from einsteinpy.coordinates import SphericalDifferential
# 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 = (u.rad * speed_at_perihelion) / distance
sph_obj = SphericalDifferential(distance, np.pi / 2 * u.rad, np.pi * u.rad,
0 * u.km / u.s, 0 * u.rad / u.s, omega)
# 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
obj = Schwarzschild.from_coords(sph_obj, M)
ans = obj.calculate_trajectory(
end_lambda=end_lambda, OdeMethodKwargs={"stepsize": stepsize}, return_cartesian=True
)
ans[0].shape, ans[1].shape
((13150,), (13150, 8))
r = np.sqrt(np.square(ans[1][:, 1]) + np.square(ans[1][:, 2]))
i = np.argmax(r)
(r[i] * u.m).to(u.km)
((ans[1][i][6]) * u.m / u.s).to(u.km / u.s)
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
0.016618341586578822
from einsteinpy.bodies import Body
from einsteinpy.geodesic import Geodesic
Sun = Body(name="Sun", mass=M, parent=None)
Object = Body(name="Earth", differential=sph_obj, parent=Sun)
geodesic = Geodesic(body=Object, time=0 * u.s, end_lambda=end_lambda, step_size=stepsize)
from einsteinpy.plotting import GeodesicPlotter
sgp = GeodesicPlotter()
sgp.plot(geodesic)
sgp.show()
FigureWidget({ 'data': [{'marker': {'color': '#ffcc00', 'line': {'width': 2}, 'size': 10}, '…
All data regarding earth's orbit is taken from https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html