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.

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[i] = x_earth[i] - x_sun[i]               # x-coordinates r-vector
    ry[i] = y_earth[i] - y_sun[i]               # y-coordinates r-vector
    r_v = (rx[i], ry[i])                        # r-vector
    radius[i] = linalg.norm(r_v)                # Length r-vector

    F_x = ((-G*m_s*m_e)/radius[i]**3) * rx[i]     # Forces x-direction
    F_y = ((-G*m_s*m_e)/radius[i]**3) * ry[i]     # 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[i] + ax_s*dt       # Velocity, sun
    vy_sun[i+1] = vy_sun[i] + ay_s*dt       # Velocity, sun

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

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

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

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

 

[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.