Elliptical orbit in vpython

Asked by At

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.

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

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.

Thank you. Thank you

0 Answers