from pylab import *
from random import *

d=0.1
y1=d/2.
y2=-d/2.

L=4.
N=100 # Number of bins
lam=0.01 # Wavelength of light (in meters)
omega=3.e8*2.*3.14159/lam

ymax=2.
dy=2*ymax/N

Iave=[]
Atot=[]
ypos=[]

for i in range(N):
	Iave.append(0)
	Atot.append(0)	
	ypos.append(dy*(i-N/2.))

photonrate=50000
dt=1
dn=int(photonrate*dt)
print dn
tmax=15.
t=0.

ion()
myplot=plot(ypos,Iave)

while t < tmax:
	t=t+dt
	print t
	for j in range(N):
		Atot[j]=0

	for i in range(dn):
		theta=1.57*random()-0.78
		y=L*tan(theta)
		r=sqrt(L*L+y*y)
		t_arr=t+r/3.e8
		A=1.*cos(omega*t_arr)
		dum=random()
		if (dum < 0.5) :	
			j=int((y+y1)/dy+N/2)
		else :
			j=int((y+y2)/dy+N/2)
		if (j >= 0) & (j < N):
			Atot[j]=Atot[j]+A
	

	for j in range(N):
		Iave[j]=Iave[j]+Atot[j]*Atot[j]

	myplot[0].set_ydata(Iave)
	ylim(0,max(Iave))
	draw()