class ball:
def __init__(self, m=1, r=np.array([0,0,0]), v=np.array([0,0,0]), a=np.array([0,0,0])):
self.m = m
self.r = r
self.v = v
self.a = a
self.r_array = np.array([])
self.v_array = np.array([])
self.a_array = np.array([])
# Gravitational constant
G = 6.674e-11 # N m^2 / kg^2
# Force function: inputs are ball_1 and ball_2
def force(ball_1, ball_2):
'''
This is the force that ball_1 feels due to ball_2. This is dependent
only on the position and mass of the balls, and the gravitational constant (G).
The order in which we put the balls matters for calcualting the correct direction of the force.
'''
output_force = np.zeros(3) # Initialze the returned force array with zeros
r = ball_1.r - ball_2.r # Position vector from ball_2 to ball_1
# Quick check if r^2 is zero to avoid division by zero
if np.sqrt(r.dot(r)) == 0:
return output_force
r_hat = r / np.sqrt(r.dot(r)) # Unit vector pointing from ball_2 to ball_1
##### Begin Editing #####
# YOUR CODE HERE
#note if you wish to use the magnitude of r^2, use r.dot(r)
output_force = -((G*ball_1.m*ball_2.m)/r.dot(r))*r_hat
##### End Editing #####
return output_force
#Insert the initial velocity for the Earth
##### Begin Editing #####
v_initial_earth = 22000 # Your answer here
##### End Editing #####
# Physical constants: mass of the Sun, mass of the Earth, and average distance between the two
M_Sun = 1.988e30 # kg
M_Earth = 5.972e24 # kg
R_Earth = 1.496e11 # m
# Create instances of the ball class for the earth and sun
earth = ball(m=M_Earth, r=np.array([R_Earth,0,0]),
v=np.array([0, v_initial_earth, 0]))
sun = ball(m=M_Sun)
##### Begin Editing #####
T = 30000000# Your answer here
dt = 3600# Your answer here
##### End Editing #####
times = np.arange(0, T+dt, dt)
N = times.size
earth.a_array = np.empty((N,3))
earth.a_array[0] = earth.a
earth.v_array = np.empty((N,3))
earth.v_array[0] = earth.v
earth.r_array = np.empty((N,3))
earth.r_array[0] = earth.r
sun.a_array = np.empty((N,3))
sun.a_array[0] = sun.a
sun.v_array = np.empty((N,3))
sun.v_array[0] = sun.v
sun.r_array = np.empty((N,3))
sun.r_array[0] = sun.r
i = 1
for t in times[1:]:
# first update the accelerations
earth.a = force(earth, sun)/earth.m
sun.a = np.zeros(3) # Here we assume no force on the sun (stays static)
# then append the new value to the list
earth.a_array[i] = earth.a
sun.a_array[i] = sun.a
# second update the velocities
earth.v = earth.v + earth.a * dt
sun.v = sun.v + sun.a * dt
# then append the new value to the list
earth.v_array[i] = earth.v
sun.v_array[i] = sun.v
# third update the positions
earth.r = earth.r + earth.v * dt
sun.r = sun.r + sun.v * dt
# then append the new value to the list
earth.r_array[i] = earth.r
sun.r_array[i] = sun.r
# update iteration count
i = i + 1
plt.plot(earth.r_array[:,0], earth.r_array[:,1], '.')
plt.plot(sun.r_array[:,0], sun.r_array[:,1], '.')
plt.gca().set_aspect('equal')
plt.show()
#Use this code cell to calculate the eccentricty of your orbit.
#Here we use earth.r_array to determine the aphelion and perihelion distances of the orbit
rMag = np.sqrt(earth.r_array[:,0]**2 + earth.r_array[:,1]**2) #turns earth.r_array to list of distances (magnitudes)
A = max(rMag) #aphelion is maximum distance
P = min(rMag) #perihelion is minimum distance
##### Begin Editing #####
#Code in the formula for eccentricity here in terms of A and P, then run the code cell to print out result
ecc = (A-P)/(A+P) #Look at paragraph above for formula
##### End Editing #####
print('The eccentricity of our simulated orbit is {:.4f}.'.format(ecc))
def init():
line.set_data([], [])
return(line)
def animate_ball_list(ball_list):
def animation_function(i):
x = []
y = []
for j in range(0,len(ball_list)):
x.append(ball_list[j].r_array[i, 0])
y.append(ball_list[j].r_array[i, 1])
line.set_data(x, y)
return (line)
return animation_function
fig3, axs3 = plt.subplots()
axs3.set_xlim((np.min(earth.r_array[:,0]) - 0.1*np.abs(np.min(earth.r_array[:,0])),
np.max(earth.r_array[:,0]) + 0.1*np.abs(np.max(earth.r_array[:,0]))))
axs3.set_ylim((np.min(earth.r_array[:,1]) - 0.1*np.abs(np.min(earth.r_array[:,1])),
np.max(earth.r_array[:,1]) + 0.1*np.abs(np.max(earth.r_array[:,1])))
)
axs3.set_aspect('equal')
axs3.plot([-1,1], [0,0], color='C1', lw=2)
line, = axs3.plot([], [], '.', markersize=24, lw=2)
plt.close()
anim_2 = animation.FuncAnimation(fig3, animate_ball_list([sun, earth]), init_func=init, frames=N, interval=40)
anim_2
rMag = np.sqrt(earth.r_array[:,0]**2 + earth.r_array[:,1]**2)
rMaxIndex = np.argmax(rMag)
rMinIndex = np.argmin(rMag)
vMag = np.sqrt(earth.v_array[:,0]**2 + earth.v_array[:,1]**2)
print('At its maximum distance, the speed of earth is {:.3f} m/s'.format(vMag[rMaxIndex]))
print('At its minimum distance, the speed of earth is {:.3f} m/s'.format(vMag[rMinIndex]))
def initialize_planet(planet, N):
'''
This takes in a planet and initializes it with
empty arrays for the position, velocity, and
acceleration histories. N denotes the total number
of time steps these histories will contain.
'''
planet.a_array = np.empty((N,3))
planet.a_array[0] = planet.a
planet.v_array = np.empty((N,3))
planet.v_array[0] = planet.v
planet.r_array = np.empty((N,3))
planet.r_array[0] = planet.r
def update_planet(planet, force_felt, dt):
'''
This takes in a planet and the force it felt and
updates the planet by one time step.
'''
# first update the acceleration
planet.a = force_felt/planet.m
# then append the new value to the list
planet.a_array[i] = planet.a
# second update the velocities
planet.v = planet.v + planet.a * dt
# then append the new value to the list
planet.v_array[i] = planet.v
# third update the positions
planet.r = planet.r + planet.v * dt
# then append the new value to the list
planet.r_array[i] = planet.r
#Insert the initial speed, T, and dt for the simulation
##### Begin Editing #####
v_initial_earth = 29750 # Your answer here
T = 40000000# Your answer here
dt = 3600# Your answer here
##### End Editing #####
times = np.arange(0, T+dt, dt)
N = times.size
earth = ball(m=M_Earth, r=np.array([R_Earth,0,0]),
v=np.array([0, v_initial_earth, 0]))
sun = ball(m=M_Sun)
initialize_planet(earth, N)
initialize_planet(sun, N)
i = 1
for t in times[1:]:
# First find the force
earth_force = force(earth, sun)
sun_force = np.zeros(3) # Here we assume no force on the sun still (stays static)
# Update the planet(s)
update_planet(earth, earth_force, dt)
update_planet(sun, sun_force, dt)
i = i + 1
# Plotting
plt.plot(earth.r_array[:,0], earth.r_array[:,1], '.')
plt.plot(sun.r_array[:,0], sun.r_array[:,1], '.')
plt.gca().set_aspect('equal')
plt.show()
def length_of_year(planet, dt, tolerance=R_Earth*0.01):
'''
This takes in a planet after it has orbitted and calculates
the length of the year based on how long it takes for it to
return back to its starting point. Because our simulation
is not perfect, some tolerance is allowed in case it misses
its starting point by a little
'''
def distance_from_start(index):
'''
This calculates the distance at a given time the planet
is from its start location
'''
start_loc = planet.r_array[0]
curr_loc = planet.r_array[index]
distance = curr_loc - start_loc
return np.sqrt(distance.dot(distance))
# This loop seems a little complicated, but really it is just first
# checking that the planet got away from start and then came back,
# so it doesn't just think it has finished a loop immediately
near_start = True
for i in range(len(planet.r_array)):
if near_start and distance_from_start(i) > tolerance:
near_start = False
if not near_start and distance_from_start(i) < tolerance:
return i*dt
return 0 # If something went wrong, it will say the year was instant
v_initial_earth = 29800 # m/s
earth = ball(m=M_Earth, r=np.array([R_Earth,0,0]),
v=np.array([0, v_initial_earth, 0]))
sun = ball(m=M_Sun)
planets = [earth, sun]
for planet in planets:
initialize_planet(planet, N)
i = 1
for t in times[1:]:
# first find the force
for planet in planets:
planet.force = 0
for second_planet in planets:
planet.force += force(planet, second_planet)
# Update the planet(s)
for planet in planets:
update_planet(planet, planet.force, dt)
i = i + 1
print('The length of one year is ~'+str(length_of_year(earth, dt))+' s.')
plt.plot(earth.r_array[:,0], earth.r_array[:,1], '.')
plt.plot(sun.r_array[:,0], sun.r_array[:,1], '.')
plt.gca().set_aspect('equal')
plt.show()
#Code cell provided for calculations of the force. Show your work to recieve credit.
$F_G = $ (G*earth_mass*sun_mass)/r^2
$F_G = $ ((6.67*10^-11)(5.97*10^24)(1.989*10^30))/(1.496*10^11)^2 = 3.54*10^22
#Code cell provided for calculations of the accelerations. Show your work to recieve credit.
$a_{Earth} = $ (G*sun_mass)/r^2 = ((6.67*10^-11)(1.989*10^30))/(1.496*10^11)^2
$a_{Earth} = $ 0.00593 m/s^2
$a_{Sun} = $ (G*earth_mass)/r^2 = ((6.67*10^-11)(5.97*10^24))/(1.496*10^11)^2
$a_{Sun} = $ 1.779*10^-8 m/s^2
#Use this code cell to determine the ratio of Sun's distance traveled during the simulation to
#the radius of the Sun
#If you are redoing this section after tackling the CAT, re-run all code cells from Part 3 to ensure
#your values are updated properly for this code.
R_Sun = 6.957e8 #[m]
FinalPositionSun = sun.r_array[-1,:]
#This indexes the position of the sun at the last time step (i.e. -1 first index)
##### Begin Editing #####
#Determine the magnitude of the position to be the norm of the final position of the Sun
#Hint: Use np.linalg.norm()
DisplacementSun = np.linalg.norm(FinalPositionSun)
# Calculate Ratio here
RatioDistoRadius = (DisplacementSun/R_Sun)
##### End Editing #####
print('Ratio of the displacement of the Sun to the radius of the Sun is '+str(RatioDistoRadius))
T = 4e7
dt = 1e4
#Insert the above conditions into the arrays below
##### Begin Editing #####
suns = [ball(m=M_Sun, r=np.array([-R_Earth, 0, 0]), v=np.array([0, -9000, 0])),
ball(m=M_Sun, r=np.array([R_Earth, 0, 0]), v=np.array([0, 9000, 0])),
ball(m=M_Earth, r=np.array([-(0.9*R_Earth), 0, 0]), v=np.array([0, 100000, 0]))]
##### End Editing #####
times = np.arange(0, T+dt, dt)
N = times.size
for sun in suns:
initialize_planet(sun, N)
i = 1
for t in times[1:]:
# first find the force
for planet in suns:
planet.force = 0
for second_planet in suns:
planet.force += force(planet, second_planet)
# Update the planet(s)
for planet in suns:
update_planet(planet, planet.force, dt)
i = i + 1
print('The length of one year is ~'+str(length_of_year(suns[2], dt, tolerance=R_Earth*0.1))+' s.')
for sun in suns:
plt.plot(sun.r_array[:,0], sun.r_array[:,1], '.')
plt.xlim(-1.5*R_Earth, 1.5*R_Earth)
plt.ylim(-R_Earth, R_Earth)
plt.gca().set_aspect('equal')
plt.show()
fig, axs = plt.subplots()
xMin = 0
xMax = 0
yMin = 0
yMax = 0
for sun in suns:
xMin = min(xMin, np.min(sun.r_array[:,0]))
yMin = min(yMin, np.min(sun.r_array[:,1]))
xMax = max(xMax, np.max(sun.r_array[:,0]))
yMax = max(yMax, np.max(sun.r_array[:,1]))
axs.set_xlim(xMin, xMax)
axs.set_ylim(yMin, yMax)
plt.gca().set_aspect('equal')
axs.plot([-1,1], [0,0], color='C1', lw=2)
line, = axs.plot([], [], '.', markersize=24, lw=2)
plt.close()
anim_many = animation.FuncAnimation(fig, animate_ball_list(suns), init_func=init, frames=N, interval=20)
anim_many
# Code cell provided for calculations with values recopied from above to use.
R_Earth = 6.38e6 # m Earth's radius
M_Sun = 1.988e30 # kg
M_Earth = 5.972e24 # kg
D_avg_Earth_Sun = 1.496e11 # m The average distance between Earth and Sun
$F_{G,closest}~=~$ (6.67*10^-11*M_Earth*M_Sun)/(D_avg_Earth_Sun-R_Earth) = 3.5386*10^22 N
# Code cell provided for calculations with values recopied from above to use.
R_Earth = 6.38e6 # m Earth's radius
M_Sun = 1.988e30 # kg
M_Earth = 5.972e24 # kg
D_avg_Earth_Sun = 1.496e11 # m The average distance between Earth and Sun
$F_{G,furthest}~=~$ (6.67*10^-11*M_Earth*M_Sun)/(D_avg_Earth_Sun+R_Earth) = 3.5386*10^22 N
Difference = $F_{G,closest}~=~$ - $F_{G,furthest}~=~$ = 6.036 * 10^18 N
# Code cell provided for calculations with values recopied from above to use.
R_Earth = 6.38e6 # m Earth's radius
M_Sun = 1.988e30 # kg
M_Earth = 5.972e24 # kg
M_Moon = 7.346e22 #kg
D_avg_Earth_Moon = 3.84467e8
F_sunTidal = (6.67*10^-11*M_Sun*M_Earth)*((2*R_Earth)/(D_avg_Earth_Sun)^3) = 3.017 * 10^18 N
F_moonTidal = (6.67*10^-11*M_Moon*M_Earth)*((2*R_Earth)/(D_avg_Earth_Moon)^3) = 6.568 * 10^18 N
Ratio_Tides = F_sunTidal/F_moonTidal = 0.459
# Code cell provided for calculations with values recopied from above to use.
R_Earth = 6.38e6 # m Earth's radius
M_Sun = 1.988e30 # kg
M_Earth = 5.972e24 # kg
d = R_Earth((2*(M_Sun))/M_Earth)^1/3 = 5.57 * 10^8