Nasa123123
- 4
- 0
- 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[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()
[/CODE]
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[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()
[/CODE]