title: Elliptical orbit in vpython
body: <p>I have the following code. This code is simulation of orbiting objects around other objects, E.g. Solar system. As you run it, the objects orbit in circular trajectory. </p>
  
  <pre><code>import math
  from vpython import *
  lamp = local_light(pos=vector(0,0,0), color=color.yellow)
  # Data in units according to the International System of Units
  G = 6.67 * math.pow(10,-11)
  
  # Mass of the Earth
  ME = 5.973 * math.pow(10,24)
  # Mass of the Moon
  MM = 7.347 * math.pow(10,22)
  # Mass of the Mars
  MMa = 6.39 * math.pow(10,23)
  # Mass of the Sun
  MS = 1.989 * math.pow(10,30)
  # Radius Earth-Moon
  REM = 384400000
  # Radius Sun-Earth
  RSE = 149600000000
  RMS = 227900000000
  # Force Earth-Moon
  FEM = G*(ME*MM)/math.pow(REM,2)
  # Force Earth-Sun
  FES = G*(MS*ME)/math.pow(RSE,2)
  # Force Mars-Sun
  FEMa = G*(MMa*MS)/math.pow(RMS,2)
  
  # Angular velocity of the Moon with respect to the Earth (rad/s)
  wM = math.sqrt(FEM/(MM * REM))
  # Velocity v of the Moon (m/s)
  vM = wM * REM
  print("Angular velocity of the Moon with respect to the Earth: ",wM," rad/s")
  print("Velocity v of the Moon: ",vM/1000," km/s")
  
  # Angular velocity of the Earth with respect to the Sun(rad/s)
  wE = math.sqrt(FES/(ME * RSE))
  # Angular velocity of the Mars with respect to the Sun(rad/s)
  wMa = math.sqrt(FEMa/(MMa * RMS))
  
  # Velocity v of the Earth (m/s)
  vE = wE * RSE
  # Velocity v of the Earth (m/s)
  vMa = wMa * RMS
  print("Angular velocity of the Earth with respect to the Sun: ",wE," rad/s")
  print("Velocity v of the Earth: ",vE/1000," km/s")
  
  
  # Initial angular position
  theta0 = 0
  
  # Position at each time
  def positionMoon(t):                                     
      theta = theta0 + wM * t
      return theta
  
  def positionMars(t):                                     
      theta = theta0 + wMa * t
      return theta
  
  def positionEarth(t):
      theta = theta0 + wE * t
      return theta
  
  
  def fromDaysToS(d):
      s = d*24*60*60
      return s
  
  def fromStoDays(s):
      d = s/60/60/24
      return d
  
  def fromDaysToh(d):
      h = d * 24
      return h
  
  # Graphical parameters
  print("\nSimulation Earth-Moon-Sun motion\n")
  days = 365
  seconds = fromDaysToS(days)
  print("Days: ",days)
  print("Seconds: ",seconds)
  
  v = vector(384,0,0)
  E = sphere(pos = vector(1500,0,0), color = color.blue, radius = 60, make_trail=True)
  Ma = sphere(pos = vector(2300,0,0), color = color.orange, radius = 30, make_trail=True)
  M = sphere(pos = E.pos + v, color = color.white,radius = 10, make_trail=True)
  S = sphere(pos = vector(0,0,0), color = color.yellow, radius=700)
  
  t = 0
  thetaTerra1 = 0
  dt = 5000
  dthetaE = positionEarth(t+dt)- positionEarth(t)
  dthetaM = positionMoon(t+dt) - positionMoon(t)
  dthetaMa = positionMars(t+dt) - positionMars(t)
  print("delta t:",dt,"seconds. Days:",fromStoDays(dt),"hours:",fromDaysToh(fromStoDays(dt)),sep=" ")
  print("Variation angular position of the Earth:",dthetaE,"rad/s that's to say",degrees(dthetaE),"degrees",sep=" ")
  print("Variation angular position of the Moon:",dthetaM,"rad/s that's to say",degrees(dthetaM),"degrees",sep=" ")
  
  while t &lt; seconds:
      rate(500)
      thetaEarth = positionEarth(t+dt)- positionEarth(t)
      thetaMoon = positionMoon(t+dt) - positionMoon(t)
      thetaMars = positionMars(t+dt) - positionMars(t)
      # Rotation only around z axis (0,0,1)
      E.pos = rotate(E.pos,angle=thetaEarth,axis=vector(0,1,0))
      Ma.pos = rotate(Ma.pos,angle=thetaMars,axis=vector(0,1,0))
      v = rotate(v,angle=thetaMoon,axis=vector(0,1,0))
      M.pos = E.pos + v
  t += dt
  </code></pre>
  
  <p>I am wondering How to change the path of orbit to elliptical?
  I have tried several ways but I could not manage to find any solution.</p>
  
  <p>Thank you.
  Thank you</p>