Lab 06 - Curlss/Computational_physics GitHub Wiki

Modeling planetary motion—more especially, Earth's orbit around the Sun—was our aim in this experiment. We intended to do this by using the Euler-Cromer technique and the Euler method for our arrangement. We also wanted to be able to modify this configuration for any two-body system, such as additional planets. Additionally, we wanted to explain how Earth moves, how altering its velocity impacts it, and how Kepler's first and second laws connect to this.

The foundation of this lab is the ability to comprehend and characterize planetary motion. We modeled circular or elliptical orbits using Python in order to better comprehend planetary motion, and we developed a model to illustrate the differences between the Euler-Cromer approach and the conventional Euler method for updating the system's state over time. Kepler's first law, which states that orbits are ellipses with circles being an exception, and Kepler's second law, which states that a planet moves more slowly when it is farther from the Sun and more quickly when it is closer, were both supported by our planetary motion models for Earth's motion around the Sun. It was observed that the orbital form is influenced by the velocity value; an orbit with the proper velocity is round and has a constant radius and speed; orbits with variations from this value are elliptical. In addition to demonstrating how beginning circumstances impact orbital motion, this code helps us comprehend the mathematical modeling and physical principles influencing planetary motion by implementing simple functions and interpreting the results via visual charts.

The lab

First, I had to create the starting circumstances for position and velocity in part 1 and allocate arrays for time, location (x and y), and velocity (vx and vy) in order to mimic planetary motion.

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

To create a function that calculates the radial distance between the Earth and the Sun, the following section was required.

# write your distance function here
def distance(x, y):
  r = np.sqrt(x**2 + y**2)
  return r

In this instance, we had to update the position and velocity arrays using the Euler-Cromer technique when cromer=True and the Euler-Cromer method when cromer=False in order to calculate the Euler-Cromer method.

def calculate(tarr, xarr, yarr, vxarr, vyarr, dt, cromer=True):
  
  for i in range(1,len(tarr)):
    r = distance(xarr[i-1], yarr[i-1])
    ax = -(4*np.pi**2*xarr[i-1])/r**3
    ay = -(4*np.pi**2*yarr[i-1])/r**3
    if cromer:
      vxarr[i] = vxarr[i-1] + ax*dt
      vyarr[i] = vyarr[i-1] + ay*dt
      xarr[i] = xarr[i-1] + vxarr[i]*dt
      yarr[i] = yarr[i-1] + vyarr[i]*dt
      tarr[i] = tarr[i-1] + dt
    else:
      xarr[i] = xarr[i-1] + vxarr[i-1]*dt
      yarr[i] = yarr[i-1] + vyarr[i-1]*dt
      vxarr[i] = vxarr[i-1] + ax*dt
      vyarr[i] = vyarr[i-1] + ay*dt
      tarr[i] = tarr[i-1] + dt


  return tarr, xarr, yarr, vxarr, vyarr

A scatter plot of the orbit (y vs. x) with time-coded points, a plot of the x and y locations vs time, coupled with the computed radial distance, and a plot displaying the x and y components of velocity along with the overall speed are all shown by this portion, which specifies a function.

# write your make_plots function here
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.subplot(1,3,2)
  plt.plot(tarr, xarr, label='x')
  plt.plot(tarr, yarr, label='y')
  plt.plot(tarr, np.sqrt(xarr**2+yarr**2), 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()

In order to create a main function that calls our initialize, calculate, and make_plots methods to illustrate the results, we must lastly have this section. The calculated arrays are returned by this primary function, which permits further examination of the simulation results. By changing the beginning circumstances, we may alter our velocity, allowing us to describe planetary motion and see the stability of the orbital model.

# write your main function here
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=True)
  make_plots(tarr, xarr, yarr, vxarr, vyarr)
  return tarr, xarr, yarr, vxarr, vyarr

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 = 1 # time in years
dt = .01# time step in years


# call main
main(x0, y0, v0x, v0y, tmax, dt,cromer=True)
plt.show()

image