C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C
C   the stochastic predator prey model 
C   using the Gillespie algorithm
C
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	integer		  seed,y1,y2,N
	double precision  x,c1,c2,c3,dt,t
	double precision  a1,a2,a3,a0
	
	x = 100000
	c1 = .0002
	c2 = .01
	c3 = 10
	N=100000
	
	seed=12423
	y1 = 1000
	y2 = 1000
	t=0.0
	
	call srand(seed)
	open(unit=10, file='lotka.txt')
	
	Do i = 1, N
		a1 = x*c1*y1
		a2 = y1*y2*c2
		a3 = y2*c3
		a0=a1+a2+a3
		
		dt = -log(rand())/a0
		r = rand()*a0
		
		if (r .lt. a1) then
			y1 = y1+1
		else
		if (r .lt. (a1+a2)) then
			y1=y1-1
			y2=y2+1
		else
			y2=y2-1
		end if
		end if
		
		t=t+dt
		
		write(10,*) t, y1, y2
	END DO
	close(10)
	End