Analytic solution of the Earth's orbit around the Sun

[Nasa123123]

Homework Statement

Find the analytical solution of earths orbit around the sun, the sun is centered in origo. Include the pull from earth on the sun, and sun on the earth.

We know that
mass of earth = 5.972*10^24 kg
mass of sun = 1.989*10^30 kg
G = 6.67 * 10^-11 Nm^2/kg^2
distance at aphelion = 1.521*10^11 m
velocity (only in y-direction at aphelion) = 29290 m/s

Relevant Equations

F = - G*m1*m2/r^2

The exercise is to compare numerical and analytical solution. I have worked out the code from earlier exercise (see code under this text), but I don't understand how the analytical solution works. I have tried to use the equation r(theta) = a(1-e^2)/(1+e*cos(theta)), which is OK but I don't think this is how we are supposed to do it. I think we are supposed to do some calculations with the center of mass but I just don't understand it.
Ps. I know my code probably can be more efficient with two dimensional arrays etc, but it works so I am just asking for help with the problem.

[CODE lang="python" title="Numerical solution"]from numpy import *
import matplotlib.pyplot as plt

G = 6.67E-11 # Gravity constant
m_s = 1.989E30 # Mass of sun
m_e = 5.972E24 # Mass of earth

t = 60*60*24*365*2
dt = 60*60*12

n = int(t/dt)
t = zeros(n)

x_sun = zeros(n, float)
y_sun = zeros(n, float)
x_earth = zeros(n, float)
y_earth = zeros(n, float)

vx_sun = zeros(n, float)
vy_sun = zeros(n, float)
vx_earth = zeros(n, float)
vy_earth = zeros(n, float)

radius = zeros(n, float) # Distance sun and earth
rx = zeros(n, float) # r-vector x
ry = zeros(n, float) # r-vector y

x_earth[0] = 1.521E11 # Start position x-coordinate
y_earth[0] = 0

x_sun[0] = 0
y_sun[0] = 0

vx_earth[0] = 0
vy_earth[0] = 29290 # Start velocity, y-coordinate Earth at aphelion

vx_sun[0] = 0 # Sun centered at origo
vy_sun[0] = 0for i in range(n-1):
rx = x_earth - x_sun # x-coordinates r-vector
ry = y_earth - y_sun # y-coordinates r-vector
r_v = (rx, ry) # r-vector
radius = linalg.norm(r_v) # Length r-vector

F_x = ((-G*m_s*m_e)/radius**3) * rx # Forces x-direction
F_y = ((-G*m_s*m_e)/radius**3) * ry # Forces y-direction

ax_e = F_x/m_e # Acceleration x-direction, Earth on sun
ay_e = F_y/m_e # Acceleration y-direction, Earth on sun

ax_s = -F_x/m_s # Acceleration x-direction, sun on earth
ay_s = -F_y/m_s # Acceleration y-direction, sun on earth

vx_sun[i+1] = vx_sun + ax_s*dt # Velocity, sun
vy_sun[i+1] = vy_sun + ay_s*dt # Velocity, sun

vx_earth[i+1] = vx_earth + ax_e*dt # Velocity, earth
vy_earth[i+1] = vy_earth + ay_e*dt # Velocity, earth

x_sun[i+1] = x_sun + vx_sun[i+1]*dt # x-coordinate, sun
y_sun[i+1] = y_sun + vy_sun[i+1]*dt # y-coordinate, sun

x_earth[i+1] = x_earth + vx_earth[i+1]*dt # x-coordinate, earth
y_earth[i+1] = y_earth + vy_earth[i+1]*dt # y-coordinate, earth

t[i+1] = t + dt # Update timeplt.plot(x_earth, y_earth)
plt.plot(x_sun, y_sun)
plt.show()
[/CODE]

 

[mfb]

The Sun can't stay at the origin in an inertial reference frame. With your starting parameters the whole system will translate and move away from the origin over time. Putting the origin at the center of mass is much, much more useful.