Planetary Motion

#Planetary Motion Program 1 
from pylab import*
gm=10.0
x=20.0
y=0.0
r=sqrt(x*x+y*y)
#velocity=sqrt(gm/r), we get circular path
# for less or more velocities we get elliptical or hyperbolic path
f=1 #choose different f values, for example f=0.8, 1.0, 1.2, 1.4 ,see difference
vx=0
vy=f*sqrt(gm/r)
h=0.05
xpos=[x]
ypos=[y]
for i in range(30000):
    r=(x*x+y*y)**0.5
    vx=vx-(gm*x*h/(r*r*r))
    vy=vy-(gm*y*h/(r*r*r))
    x=x+vx*h
    y=y+vy*h
    xpos.append(x)
    ypos.append(y)
plot(xpos,ypos,'.b')
xlabel('x')
ylabel('Y')
axis('equal')
show()

 Out put

 


#-------------------PLANETARY MOTION program 2--------------------------#
#Gravitation is a conservative fore: E = T + V
#The total energy of the system can be expressed as two coupled 1st order odes:
#dr/dt = v              Where v is the velocity
#dv/dt = F/m            Where F is the force and m is the mass
#F = (G*m1*m2)/(r**2)   m1 and m2 are the mass of the sun and mars

from pylab import *

#Set parameters:
r=1 # Astronomical unit
N = 25                   # No of days in Earth's year
dt = 1.00/N                 # Time Step: Fractions of a year
mu = 4 * pi**2           # mu=4pi^2 is the Gravitational Parameter: mu = GM #where G=6.67e-11 is the Universal Gravitational Constant and M is the mass of #the body

#Create an array, for all variables, of size N with all entries equal to zero:
x =zeros((N,))
y=zeros((N,))
vx=zeros((N,))
vy=zeros((N,))

# Initial Conditions:
x[0] = r                   # (x0 = r, y0 = 0) in AU
vy[0] = sqrt(mu/r)      # (vx0 = 0, v_y0 = sqrt(mu/r)) AU/yr

#Implement Verlet Algorithm:
for k in range(0,N-1):
    vx[k+1]=vx[k]-mu*x[k]*dt/r**3
    x[k+1]=x[k]+vx[k+1]*dt
    vy[k+1]=vy[k]-mu*y[k]*dt/r**3
    y[k+1]=y[k]+vy[k+1]*dt

#Plot:
plot(x,y,'go')
title ('Circular Orbit of Earth')
xlabel ('x')
ylabel ('y')
axis('equal')
show()