Lab 10: Modeling Gravitational Motion - Duplicate

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

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()

#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 = R_Earth((2*(M_Sun))/M_Earth)^1/3 = 5.57 * 10^8