Goal

The goal of this lab is to model the motion of a planet orbiting a star.

Overview

In this lab we model the Earth's motion around the Sun, but we write the code in such a way that it can be easily adapted to model other 2-body systems. Our distances are in Astronomical units (Au) and our time is measured in years. We also assume that the Sun is fixed at a central position and Earth revolves around it.

Functions

The first function we make is the initialize function to create our position, time, and velocity arrays.

def initialize(x0,y0,v0x,v0y,tmax,dt):
  """This function takes in initial positions and velocities and times and returnes time, position and velocity arrays"""
  nsteps=int(tmax/dt)
  time=np.zeros(nsteps)
  x=np.zeros(nsteps)
  y=np.zeros(nsteps)
  vx=np.zeros(nsteps)
  vy=np.zeros(nsteps)
  x[0]=x0
  y[0]=y0
  vx[0]=v0x
  vy[0]=v0y
  return time,x,y,vx,vy

The next function was called distance. This calculated the distance between the Earth and Sun given the Earth's x and y positions. We assume throughout that the Sun is located at 𝑥=0,𝑦=0 , and that the Sun does not move.

def distance(x,y):
  """This function takes in x and y position and returns the distance between earth and the sun, r."""
  r=np.sqrt(x**2 + y**2)

  return r

We then make a calculate function to gives us time, position, and velocity arrays. We use an if statement to determine weather we use the Euler method or the Euler-Cromer method.

def calculate(time,x,y,vx,vy,dt,cromer=True):
  """This function calculates vx, vy, and t arrays. As well uses either the Eulor method or Euler-Cromer method to get x and y arrays"""

  for i in range(1,len(x)):
    ax=-(4*np.pi**2)/(distance(x[i-1],y[i-1]))**3*x[i-1]
    ay=-(4*np.pi**2)/(distance(x[i-1],y[i-1]))**3*y[i-1]
    vx[i]=vx[i-1]+ax*dt
    vy[i]=vy[i-1]+ay*dt
    time[i]=time[i-1]+dt
    if cromer==True:
      x[i]=x[i-1]+vx[i]*dt
      y[i]=y[i-1]+vy[i]*dt
    if cromer==False:
      x[i]=x[i-1]+vx[i-1]*dt
      y[i]=y[i-1]+vy[i-1]*dt

  return time,x,y,vx,vy

We then make a function to make our plots, we create them as subplots.

def make_plots(time,x,y,vx,vy):
  """This function will plot all of our arrays"""
  plt.figure(figsize=(15,4))
  plt.subplots_adjust(wspace=.35)

  plt.subplot(1,3,1)
  plt.scatter(x,y,c=time,s=10,label="X vs Y")
  plt.colorbar(label="Time (yr)")
  plt.xlabel('X-Position (Au)')
  plt.ylabel('Y-position (Au)')
  plt.axis('equal')

  plt.subplot(1,3,2)
  plt.plot(time,x,label="X vs Time")
  plt.plot(time,y,label="Y vs Time")
  plt.plot(time, np.sqrt(x**2+y**2),label="Radius vs Time")
  plt.ylabel('Position (Au)')
  plt.xlabel('Time (yrs)')
  plt.legend()


  plt.subplot(1,3,3)
  plt.plot(time,vx,label="X-Velocity vs Y-Velocity")
  plt.plot(time,vy,label="y-Velocity vs Y-Velocity")
  plt.plot(time, np.sqrt(vx**2+vy**2),label="Magnitude of Velocity vs Time")

  plt.xlabel('Time (yrs)')
  plt.ylabel('Velocity (Au/yr)')
  plt.legend()

We then make a main function which returns our final time, position, and velocity.

def main(x0,y0,v0x,v0y,tmax,dt,cromer=True):
  """This function calls our previus three functions"""
  t1,x1,y1,vx1,vy1=initialize(x0,y0,v0x,v0y,tmax,dt)
  t2,x2,y2,vx2,vy2=calculate(t1,x1,y1,vx1,vy1,dt,cromer=cromer)
  #t3,x3,y3,vx3,vy3=
  make_plots(t2,x2,y2,vx2,vy2)
  return t2,x2,y2,vx2,vy2

Graphs and Plots

WHen we call our main function we get the following graph. This shows how Earth's orbit stays the same and the magnitude of velocity stays the same as well as the radius.

When we increase the initial velocity by 1, we get the following plots.

When we decrease the initial velocity by one, we get the following plots.

The following models results show that with the total time of two years with a new initial velocity of 2$\pi$.

We then used the Euler method for the following plots.

Adapted model

The following code shows the model of motion for the planet Venus to show that it works for any planet.

def make_plots(time,x,y,vx,vy):
  """This function will plot all of our arrays"""
  plt.figure(figsize=(15,4))
  plt.subplots_adjust(wspace=.35)

  plt.subplot(1,3,1)
  plt.scatter(x,y,c=time,s=10,label="X vs Y")
  plt.colorbar(label="Time (yr)")
  plt.xlabel('X-Position (Au)')
  plt.ylabel('Y-position (Au)')
  plt.axis('equal')

  plt.subplot(1,3,2)
  r=np.zeros(len(np.sqrt(x**2+y**2)))+.72
  plt.plot(time,x,label="X vs Time")
  plt.plot(time,y,label="Y vs Time")
  plt.plot(time, r,label="Radius vs Time")
  plt.ylabel('Position (Au)')
  plt.xlabel('Time (yrs)')
  plt.legend()


  plt.subplot(1,3,3)
  plt.plot(time,vx,label="X-Velocity vs Y-Velocity")
  plt.plot(time,vy,label="y-Velocity vs Y-Velocity")
  plt.plot(time, np.sqrt(vx**2+vy**2),label="Magnitude of Velocity vs Time")

  plt.xlabel('Time (yrs)')
  plt.ylabel('Velocity (Au/yr)')
  plt.legend()
def main(x0,y0,v0x,v0y,tmax,dt,cromer=True):
  """This function calls our previus three functions"""
  t1,x1,y1,vx1,vy1=initialize(x0,y0,v0x,v0y,tmax,dt)
  t2,x2,y2,vx2,vy2=calculate(t1,x1,y1,vx1,vy1,dt,cromer=cromer)
  #t3,x3,y3,vx3,vy3=
  make_plots(t2,x2,y2,vx2,vy2)
  return t2,x2,y2,vx2,vy2

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 =.001 # time step in years


# call main
t2,x2,y2,vx2,vy2=main(x0,y0,v0x,v0y,tmax,dt,cromer=True)