Modeling Planetary Motion Lab

• Goal The goal of this lab is to model planetary motion orbiting a star. In the code below we will start by using the model of the Earth orbiting the Sun, the code can be adjusted for any mass orbiting another.

Part 1 Creating functions

*Initialize function The goal of this function will do the following to set all our initial condition arrays to zero by using np.zero

def initialize(x0, y0, v0x, v0y, tmax, dt):
  
  nsteps = int(tmax/dt)
  tarr = np.zeros(nsteps)
  xarr = np.zeros(nsteps)
  yarr = np.zeros(nsteps)
  vxarr = np.zeros(nsteps)
  vyarr = np.zeros(nsteps)
  xarr[0] = x0
  yarr[0] = y0
  vxarr[0] = v0x
  vyarr[0] = v0y
  return tarr, xarr, yarr, vxarr, vyarr

• Distance function The purpose of this function is to constantly know how far the Earth is from the sun from any given point of x,y cords. Assuming the sun is set at 0,0

def distance(x, y):
  r = np.sqrt(x**2 + y**2)
  return r

• Calculate function The calculate function should complete the following compute vx,vy, and time. Depending on if Euler-cromer true or Euler-cromer false will cause the Euler method to be computed using standard Euler compared to Euler-cromer

GM = 4*np.pi**2 # G*Msun in AU^3/yr^2
def calculate(tarr, xarr, yarr, vxarr, vyarr, dt, cromer=True):
  for i in range(len(tarr)-1):
    tarr[i+1] = tarr[i] + dt
    if cromer:
      vxarr[i+1] = vxarr[i] - GM*xarr[i]/(distance(xarr[i],yarr[i])**3)*dt
      vyarr[i+1] = vyarr[i] - GM*yarr[i]/(distance(xarr[i],yarr[i])**3)*dt
      xarr[i+1] = xarr[i] + vxarr[i+1]*dt
      yarr[i+1] = yarr[i] + vyarr[i+1]*dt
    else:
      vxarr[i+1] = vxarr[i] - GM*xarr[i]/(distance(xarr[i],yarr[i])**3)*dt
      vyarr[i+1] = vyarr[i] - GM*yarr[i]/(distance(xarr[i],yarr[i])**3)*dt
      xarr[i+1] = xarr[i] + vxarr[i]*dt
      yarr[i+1] = yarr[i] + vyarr[i]*dt
  return tarr, xarr, yarr, vxarr, vyarr

• Plot function This is a first in the Wiki so I will explain the purpose of this code it is to allow us to save more time than copying and pasting a series of single plot functions for every test we will run in this Lab. This function will plot the same plot for every test preformed

def make_plots(tarr, xarr, yarr, vxarr, vyarr):
  plt.figure(figsize=(15,4))
  plt.subplots_adjust(wspace=.35)
  plt.subplot(1,3,1)
  plt.scatter(xarr,yarr,c=tarr,s=10)
  plt.colorbar(label="Time (yr)")
  plt.axis('equal')
  plt.plot(0,0,'ro',label='Sun')
  plt.legend()
  plt.xlabel('x (AU)')
  plt.ylabel('y (AU)')
  plt.subplot(1,3,2)
  plt.plot(tarr,xarr,label='x')
  plt.plot(tarr,yarr,label='y')
  plt.plot(tarr,distance(xarr,yarr),label='r')
  plt.xlabel('Time (yr)')
  plt.ylabel('Position (AU)')
  plt.legend()
  plt.subplot(1,3,3)
  plt.plot(tarr,vxarr,label='vx')
  plt.plot(tarr,vyarr,label='vy')
  plt.plot(tarr,np.sqrt(vxarr**2+vyarr**2),label='v')
  plt.xlabel('Time (yr)')
  plt.ylabel('Velocity (AU/yr)')
  plt.legend()
  plt.show()

• Main program The purpose of the main program is that it will use all the other functions within one function saving time when computing and creating graphs for the series of test conditions we will be preforming later in the lab.

def main(x0, y0, v0x, v0y, tmax, dt, cromer=True):
  tarr, xarr, yarr, vxarr, vyarr = initialize(x0, y0, v0x, v0y, tmax, dt)
  tarr, xarr, yarr, vxarr, vyarr = calculate(tarr, xarr, yarr, vxarr, vyarr, dt, cromer)
  make_plots(tarr, xarr, yarr, vxarr, vyarr)
  return tarr, xarr, yarr, vxarr, vyarr

• Running the program This is the example of the function being ran with the following conditions set below

x0 = 1  # initial x position in AU
y0 = 0 # initial y position in AU
v0x = 1   # initial x velocity in AU/yr
v0y = 2* np.pi # initial y velocity in AU/yr
tmax = 1 # time in years
dt =0.01 # time step in years


# call main
main_1 = main(x0, y0, v0x, v0y, tmax, dt)

• Outcome

What if scenarios

In this section we will changing various variables to see the difference in the outcome we receive -Increasing initial velocity by 1 code below

# write your code here to increase the initial velocity by one

x0 =1  # initial x position in AU
y0 =0  # initial y position in AU
v0x =2  # initial x velocity in AU/yr
v0y =2*np.pi  # initial y velocity in AU/yr
tmax =3 # time in years
dt = .01 # time step in years


# call main
main_2 = main(x0, y0, v0x, v0y, tmax, dt)

-Outcome

• decreasing the initial velocity by 1

# write code to decrease initial velocity by one
x0 =1  # initial x position in AU
y0 =0  # initial y position in AU
v0x =0  # initial x velocity in AU/yr
v0y =2*np.pi  # initial y velocity in AU/yr
tmax = 3# time in years
dt =.01 # time step in years


# call main
main_3 = main(x0, y0, v0x, v0y, tmax, dt)

-Creating more stability in our models In this model we will be changing the time interval, as well as setting v0y = 2 * np.pi

x0 =1 # initial x position in AU
y0 = 0 # initial y position in AU
v0x = 1 # initial x velocity in AU/yr
v0y =2*np.pi  # initial y velocity in AU/yr
tmax =3 # time in years
dt = 0.01# time step in years


# call main
main_4 = main(x0, y0, v0x, v0y, 2, dt)
main_5 = main(x0, y0, v0x, v0y, 3, dt)
main_6 = main(x0, y0, v0x, v0y, 4, dt)

-Outcome

Euler compared to the Euler-Cromer

In this test we will set Euler-cromer to false allowing the Standard method of Euler to be preformed

x0 =1  # initial x position in AU
y0 =0  # initial y position in AU
v0x =0  # initial x velocity in AU/yr
v0y =2*np.pi  # initial y velocity in AU/yr
tmax = 4# time in years
dt = 0.01# time step in years


# call main
main_7 = main(x0, y0, v0x, v0y, tmax, dt, cromer=False)

The results we got were that the Euler method was not as effective as the Euler-Cromer method when handling this type of modeling

• Below is code that can be adapted to be used on any body of mass/ planet in our solar system or beyond. This code will be for planets in our solar system sorry pluto

planets =np.array(['Mercury', 'Venus', 'Earth', 'Mars', 'Jupiter', 'Saturn', 'Uranus', 'Neptune', 'Pluto'])
orbital_radius = np.array([0.39, 0.72, 1.00, 1.52, 5.20, 9.54, 19.19, 30.06, 39.53])
orbital_eccentricity = np.array([0.206, 0.007, 0.017, 0.093, 0.048, 0.056, 0.046, 0.010, 0.248])
mass = np.array([2.4E23, 4.9E24, 6.0E24, 6.6E23, 1.9E27, 5.7E26, 8.8E25, 1.03E26, 1.3E22])
r_max = orbital_radius*(1+orbital_eccentricity)
M_sun = 1.989e33
e = orbital_eccentricity
M_P = mass
v_min = (2*np.pi)*np.sqrt((1-e)/r_max*(1+M_P/M_sun))
for r,v in zip(r_max,v_min):
  main(r,0,0,v,250,.001)