Central field motion -motion of satellite around earth

from pylab import*

gm=10.0
x=20.0
y=0.0
#Satellite
cir=Circle((x,y),radius=2,  fc='g')
gca().add_patch(cir)
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=0.8 #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]
#Earth
cir=Circle((0,0),radius=5,  fc='y')
gca().add_patch(cir)
 
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')

axis('scaled')
xlabel('x')
ylabel('Y')
show()

Simulation of standing waves

from pylab import*
a=2.0
lam=1.0
v=1.0
tmax=10
xmax=3
dt=0.05
x=linspace(0,xmax,1000)
t=0
while t<=tmax:
    y=a*(sin(2*pi*(x-v*t)/lam)+sin(2*pi*(x+v*t)/lam))
    plot(x,y,'b')
    t=t+dt

show()

Driven LCR circuit

from pylab import*
from math import*
c=1
l=1
r=1.5
q0=0
i0=0
v0=1
w=1.5
t=0.0
h=0.01
def f1(t,i,q):
    return((-1/l)*((q/c)+r*i-v0*sin(w*t)))
def f2(t,i,q):
    return(i)
i=[i0]
q=[q0]
time=[t]

while t<200:
    k1=h*f1(t,i0,q0)
    l1=h*f2(t,i0,q0)
    k2=h*f1(t + h/2.0,i0+k1/2.0,q0+l1/2.0 )
    l2=h*f2(t + h/2.0,i0+k1/2.0,q0+l1/2.0 )
    k3=h*f1(t + h/2.0,i0+k2/2.0,q0+l2/2.0 )
    l3=h*f2(t + h/2.0,i0+k2/2.0,q0+l2/2.0 )
    k4=h*f1(t + h, i0+k3, q0+l3 )
    l4=h*f2(t + h, i0+k3, q0+l3 )
      
    i1=i0+(k1 + 2*k2+ 2*k3 + k4)/6.0
    q1=q0+(l1 + 2*l2+ 2*l3 + l4)/6.0
    i.append(i1)
    q.append(q1)
    t=t+h
    time.append(t)
    i0=i1
    q0=q1

figure(1)
title("Current in driven LCR circuit")
xlabel(" Time")
ylabel("Current")
plot(time,i)
grid(True)

figure(2)
title(" Charge in driven LCR Circuit")
xlabel("Time")
ylabel("Charge")
plot(time,q)
grid(True)

#figure(3)
#title("Harmonic oscillator motion")
#xlabel("x")
#ylabel("Velocity")
#plot(i,q)
#grid(True)

show()

Radio active decay - Monte carlo method

from pylab import*
N0=1000
la=2
dt=0.05
tmax=50
N=N0
t=0
while t<=tmax:
    plot(t,N,'.b')
    for i in range(1,N):
        x=random()
       
        if x<la*dt:
            N=N-1
   
    t=t+dt
td=linspace(0,50,1000)
nd=N0*exp(-la*td)
plot(td, nd, 'r-', linewidth=2)
show()

Value of Pi hit and miss method Monte carlo integration

from pylab import*
darts=1000
hits=0
for i in range(1,darts):
    x=random()
    y=random()
    plot(x,y,'.b')
   
    dist=(x**2+y**2)
    if dist<=1.0:
        hits=hits+1.0
pi=4*(hits/darts)
print pi
cx = linspace(0,1,100)
cy = sqrt(1 - cx**2)
plot(cx, cy, 'r-', linewidth=2)
show()