program LangevinDynamics2d_active
use TPS2d_active_def_enhanced
implicit none
	integer oldtraiet
	integer i,ntpt,l,j,k,nstep1,nstep2,m,u,n,ns,tr,tcount,cc,xint,yint
	integer rix,riy
	real*8 x,y,phi,xnew,ynew,phinew,xstart,xfinal,tpt,DW,DW2,s,mm,phi2,dist,inphi,rdist1
	real*8 rdistc,dx,dy,rdistr,thetai,theta_o,theta_n,dummy,a,vn,w,b
	real*8 x1,x2,y1,y2,r1,r2,r3,r4,xi,yi
	integer :: istarts(0:29999)
	integer :: istops(0:29999)
	integer :: T_reach(0:29999)
	
	idum = -13
	iduv = -7
!~ 	idum = -7
!~ 	iduv = -13
	phi = 0.d0
	file_name4 = 'start_stop_tpd.dat'
	file_name5 = 'trajs_tpd.dat'
	file_name6 = 'tpd1.dat'
	file_name7 = 'tpd2.dat'
	file_name8 = 'tpd3.dat'
	l = 0
	tr = 0
	j = 0
	tcount = 0
	
	do i = 0,29999
		istarts(i) = 0
		istops(i) = 0
		T_reach(i) = 0
	enddo
	
	do i = 0,999999
		trajx(i) = 0
		trajy(i) = 0
	enddo
	reactiveprobdensity1(:,:)=0.d0
	reactiveprobdensity2(:,:)=0.d0
	reactiveprobdensity3(:,:)=0.d0
	
!~ 	open(unit=15,file=file_name5,status='unknown')
	
	i = 0
	do while (i .lt. 10000)
		x = xc1
		y = yc1
		l = 0
		a = ranv()
		theta_o = 2*pi*a
		b = ranv()
		call GaussianVariable (DW,DW2)
		do while (l .lt. 3)
!~ 			write(15,*) x,y
			call GaussianVariable (DW,DW2)
			xnew = x + v*dcos(theta_o)*dt + fwx2d(x,y)/zeta
			ynew = y + v*dsin(theta_o)*dt + fwy2d(x,y)/zeta
			call GaussianVariable (DW,DW2)
			theta_n = theta_o + sigmar * DW
			rdistc = dsqrt((xnew-xc1)**2+(ynew-yc1)**2)
			rdistr = dsqrt(xnew**2+ynew**2)
			dx = xnew-x
			dy = ynew-y
			if ((rdistc .le. rc) .and. (l .eq. 0)) then
				inphi = theta_o
				l = 1
			else if ((rdistc .le. rc) .and. (l .eq. 1)) then
				inphi = theta_o
			else if ((rdistc .gt. rc) .and. ((l .eq. 1) .or. (l .eq. 2))) then
				if (l .eq. 1) tcount = 0
				if (l .eq. 1) dummy = j
				trajx(tcount) = x
				trajy(tcount) = y
				l = 2
				if (rdistr .le. rr) then
					l = 3
				else if (xnew .le. xtloc) then
					l = 3
					istarts(tr) = dummy
					istops(tr) = j
					T_reach(tr) = 1
				endif
				tcount = tcount + 1
			else if ((rdistc .le. rc) .and. (l .eq. 2)) then
				trajx(tcount) = x
				trajy(tcount) = y
				tcount = tcount + 1
			endif
			x = xnew
			y = ynew
			theta_o = theta_n
			j = j + 1
		enddo
!~ 		if (T_reach(tr) .eq. 1) then
!~ 			do cc = 0,tcount-1
!~ 				xint = Nint(trajx(cc)/grid_dx)
!~ 				yint = Nint(trajy(cc)/grid_dy)
!~ 				xi = trajx(cc)
!~ 				yi = trajy(cc)
!~ 				if (xi.ge.grid_x_left .and. xi.le.grid_x_right .and. yi.ge.grid_y_left .and. yi.le.grid_y_right) then
!~ 					if ((xi.ge.(xint*grid_dx)).and.(xi.lt.((xint+1)*grid_dx))) then
!~ 						rix = xint
!~ 					else if ((xi.ge.((xint-1)*grid_dx)).and.(xi.lt.(xint*grid_dx))) then
!~ 						rix = xint-1
!~ 					endif
!~ 					if ((yi.ge.(yint*grid_dy)).and.(yi.lt.((yint+1)*grid_dy))) then
!~ 						riy = yint
!~ 					else if ((yi.ge.((yint-1)*grid_dy)).and.(yi.lt.(yint*grid_dy))) then
!~ 						riy = yint-1
!~ 					endif
!~ 					reactiveprobdensity1(rix,riy) = reactiveprobdensity1(rix,riy) + 1.d0
!~ 				endif
!~ 			enddo
!~ 		endif
		if (T_reach(tr) .eq. 1) then
			do cc = 0,tcount-1
				xint = Nint(trajx(cc)/grid_dx)
				yint = Nint(trajy(cc)/grid_dy)
				if (xint.ge.n_grid_x_left .and. xint.le.n_grid_x_right .and. yint.ge.n_grid_y_left .and. yint.le.n_grid_y_right) then
					reactiveprobdensity1(xint,yint) = reactiveprobdensity1(xint,yint) + 1.d0
				endif
			enddo
		endif
		do cc = 0,999999
			trajx(cc) = 0
			trajy(cc) = 0
		enddo
		tr = tr + 1
		if (mod(i,100) .eq. 0) print*, i, xnew, ynew
		i = i+1
	enddo
	
	l = 0
	tcount = 0
	
	i = 0
	do while (i .lt. 10000)
		x = xc2
		y = yc2
		l = 0
		a = ranv()
		theta_o = 2*pi*a
		b = ranv()
		call GaussianVariable (DW,DW2)
		do while (l .lt. 3)
!~ 			write(15,*) x,y
			call GaussianVariable (DW,DW2)
			xnew = x + v*dcos(theta_o)*dt + fwx2d(x,y)/zeta
			ynew = y + v*dsin(theta_o)*dt + fwy2d(x,y)/zeta
			call GaussianVariable (DW,DW2)
			theta_n = theta_o + sigmar * DW
			rdistc = dsqrt((xnew-xc2)**2+(ynew-yc2)**2)
			rdistr = dsqrt(xnew**2+ynew**2)
			dx = xnew-x
			dy = ynew-y
			if ((rdistc .le. rc) .and. (l .eq. 0)) then
				inphi = theta_o
				l = 1
			else if ((rdistc .le. rc) .and. (l .eq. 1)) then
				inphi = theta_o
			else if ((rdistc .gt. rc) .and. ((l .eq. 1) .or. (l .eq. 2))) then
				if (l .eq. 1) tcount = 0
				if (l .eq. 1) dummy = j
				trajx(tcount) = x
				trajy(tcount) = y
				l = 2
				if (rdistr .le. rr) then
					l = 3
				else if (xnew .le. xtloc) then
					l = 3
					istarts(tr) = dummy
					istops(tr) = j
					T_reach(tr) = 1
				endif
				tcount = tcount + 1
			else if ((rdistc .le. rc) .and. (l .eq. 2)) then
				trajx(tcount) = x
				trajy(tcount) = y
				tcount = tcount + 1
			endif
			x = xnew
			y = ynew
			theta_o = theta_n
			j = j + 1
		enddo
		if (T_reach(tr) .eq. 1) then
			do cc = 0,tcount-1
				xint = Nint(trajx(cc)/grid_dx)
				yint = Nint(trajy(cc)/grid_dy)
				if (xint.ge.n_grid_x_left .and. xint.le.n_grid_x_right .and. yint.ge.n_grid_y_left .and. yint.le.n_grid_y_right) then
					reactiveprobdensity2(xint,yint) = reactiveprobdensity2(xint,yint) + 1.d0
				endif
			enddo
		endif
!~ 		if (T_reach(tr) .eq. 1) then
!~ 			do cc = 0,tcount-1
!~ 				xint = Nint(trajx(cc)/grid_dx)
!~ 				yint = Nint(trajy(cc)/grid_dy)
!~ 				xi = trajx(cc)
!~ 				yi = trajy(cc)
!~ 				if (xi.ge.grid_x_left .and. xi.le.grid_x_right .and. yi.ge.grid_y_left .and. yi.le.grid_y_right) then
!~ 					if ((xi.ge.(xint*grid_dx)).and.(xi.lt.((xint+1)*grid_dx))) then
!~ 						rix = xint
!~ 					else if ((xi.ge.((xint-1)*grid_dx)).and.(xi.lt.(xint*grid_dx))) then
!~ 						rix = xint-1
!~ 					endif
!~ 					if ((yi.ge.(yint*grid_dy)).and.(yi.lt.((yint+1)*grid_dy))) then
!~ 						riy = yint
!~ 					else if ((yi.ge.((yint-1)*grid_dy)).and.(yi.lt.(yint*grid_dy))) then
!~ 						riy = yint-1
!~ 					endif
!~ 					reactiveprobdensity2(rix,riy) = reactiveprobdensity2(rix,riy) + 1.d0
!~ 				endif
!~ 			enddo
!~ 		endif
		do cc = 0,999999
			trajx(cc) = 0
			trajy(cc) = 0
		enddo
		tr = tr + 1
		if (mod(i,100) .eq. 0) print*, i, xnew, ynew
		i = i+1
	enddo

	l = 0
	tcount = 0
	
	i = 0
	do while (i .lt. 10000)
		x = xc3
		y = yc3
		l = 0
		a = ranv()
		theta_o = 2*pi*a
		b = ranv()
		call GaussianVariable (DW,DW2)
		do while (l .lt. 3)
!~ 			write(15,*) x,y
			call GaussianVariable (DW,DW2)
			xnew = x + v*dcos(theta_o)*dt + fwx2d(x,y)/zeta
			ynew = y + v*dsin(theta_o)*dt + fwy2d(x,y)/zeta
			call GaussianVariable (DW,DW2)
			theta_n = theta_o + sigmar * DW
			rdistc = dsqrt((xnew-xc3)**2+(ynew-yc3)**2)
			rdistr = dsqrt(xnew**2+ynew**2)
			dx = xnew-x
			dy = ynew-y
			if ((rdistc .le. rc) .and. (l .eq. 0)) then
				inphi = theta_o
				l = 1
			else if ((rdistc .le. rc) .and. (l .eq. 1)) then
				inphi = theta_o
			else if ((rdistc .gt. rc) .and. ((l .eq. 1) .or. (l .eq. 2))) then
				if (l .eq. 1) tcount = 0
				if (l .eq. 1) dummy = j
				trajx(tcount) = x
				trajy(tcount) = y
				l = 2
				if (rdistr .le. rr) then
					l = 3
				else if (xnew .le. xtloc) then
					l = 3
					istarts(tr) = dummy
					istops(tr) = j
					T_reach(tr) = 1
				endif
				tcount = tcount + 1
			else if ((rdistc .le. rc) .and. (l .eq. 2)) then
				trajx(tcount) = x
				trajy(tcount) = y
				tcount = tcount + 1
			endif
			x = xnew
			y = ynew
			theta_o = theta_n
			j = j + 1
		enddo
		if (T_reach(tr) .eq. 1) then
			do cc = 0,tcount-1
				xint = Nint(trajx(cc)/grid_dx)
				yint = Nint(trajy(cc)/grid_dy)
				if (xint.ge.n_grid_x_left .and. xint.le.n_grid_x_right .and. yint.ge.n_grid_y_left .and. yint.le.n_grid_y_right) then
					reactiveprobdensity3(xint,yint) = reactiveprobdensity3(xint,yint) + 1.d0
				endif
			enddo
		endif
!~ 		if (T_reach(tr) .eq. 1) then
!~ 			do cc = 0,tcount-1
!~ 				xint = Nint(trajx(cc)/grid_dx)
!~ 				yint = Nint(trajy(cc)/grid_dy)
!~ 				xi = trajx(cc)
!~ 				yi = trajy(cc)
!~ 				if (xi.ge.grid_x_left .and. xi.le.grid_x_right .and. yi.ge.grid_y_left .and. yi.le.grid_y_right) then
!~ 					if ((xi.ge.(xint*grid_dx)).and.(xi.lt.((xint+1)*grid_dx))) then
!~ 						rix = xint
!~ 					else if ((xi.ge.((xint-1)*grid_dx)).and.(xi.lt.(xint*grid_dx))) then
!~ 						rix = xint-1
!~ 					endif
!~ 					if ((yi.ge.(yint*grid_dy)).and.(yi.lt.((yint+1)*grid_dy))) then
!~ 						riy = yint
!~ 					else if ((yi.ge.((yint-1)*grid_dy)).and.(yi.lt.(yint*grid_dy))) then
!~ 						riy = yint-1
!~ 					endif
!~ 					reactiveprobdensity3(rix,riy) = reactiveprobdensity3(rix,riy) + 1.d0
!~ 				endif
!~ 			enddo
!~ 		endif
		do cc = 0,999999
			trajx(cc) = 0
			trajy(cc) = 0
		enddo
		tr = tr + 1
		if (mod(i,100) .eq. 0) print*, i, xnew, ynew
		i = i+1
	enddo
	
!~ 	close(15)
	
!~ 	open(unit=10,file=file_name4,status='unknown')
!~ 	do i = 0,29999
!~ 		write(10,*) istarts(i),istops(i),T_reach(i)
!~ 	enddo
!~ 	close(10)
	
!~ 	do i=n_grid_x_left,n_grid_x_right
!~ 		do j=n_grid_y_left,n_grid_y_right
!~ 			x1 = i*grid_dx
!~ 			y1 = j*grid_dy
!~ 			x2 = i*grid_dx+grid_dx
!~ 			y2 = j*grid_dy+grid_dy
!~ 			r1 = dsqrt(x1**2+y1**2)
!~ 			r2 = dsqrt(x1**2+y2**2)
!~ 			r3 = dsqrt(x2**2+y1**2)
!~ 			r4 = dsqrt(x2**2+y2**2)
!~ 			if ((r1.gt.1) .or. (r2.gt.1) .or. (r3.gt.1) .or. (r4.gt.1)) then
!~ 				reactiveprobdensity1(i,j) = 0
!~ 				reactiveprobdensity2(i,j) = 0
!~ 				reactiveprobdensity3(i,j) = 0
!~ 			endif
!~ 			if ((r1.lt.0.1d0) .or. (r2.lt.0.1d0) .or. (r3.lt.0.1d0) .or. (r4.lt.0.1d0)) then
!~ 				reactiveprobdensity1(i,j) = 0
!~ 				reactiveprobdensity2(i,j) = 0
!~ 				reactiveprobdensity3(i,j) = 0
!~ 			endif
!~ 		enddo
!~ 	enddo

	do i=n_grid_x_left,n_grid_x_right
		do j=n_grid_y_left,n_grid_y_right
			x1 = i*grid_dx-(grid_dx/2.d0)
			y1 = j*grid_dy-(grid_dy/2.d0)
			x2 = i*grid_dx+(grid_dx/2.d0)
			y2 = j*grid_dy+(grid_dy/2.d0)
			r1 = dsqrt(x1**2+y1**2)
			r2 = dsqrt(x1**2+y2**2)
			r3 = dsqrt(x2**2+y1**2)
			r4 = dsqrt(x2**2+y2**2)
			if ((r1.gt.1) .or. (r2.gt.1) .or. (r3.gt.1) .or. (r4.gt.1)) then
				reactiveprobdensity1(i,j) = 0
				reactiveprobdensity2(i,j) = 0
				reactiveprobdensity3(i,j) = 0
			endif
			if ((r1.lt.0.1d0) .or. (r2.lt.0.1d0) .or. (r3.lt.0.1d0) .or. (r4.lt.0.1d0)) then
				reactiveprobdensity1(i,j) = 0
				reactiveprobdensity2(i,j) = 0
				reactiveprobdensity3(i,j) = 0
			endif
		enddo
	enddo
	
	open(unit=10,file=file_name6,status='unknown')
	write(10,*) '# ',grid_x_left,grid_x_right,grid_dx
	write(10,*) '# ',grid_y_left,grid_y_right,grid_dy
	do i=n_grid_x_left,n_grid_x_right
		do j=n_grid_y_left,n_grid_y_right
			write(10,*) i*grid_dx,j*grid_dy,reactiveprobdensity1(i,j)
		enddo
	enddo
	close(10)
	
	open(unit=10,file=file_name7,status='unknown')
	write(10,*) '# ',grid_x_left,grid_x_right,grid_dx
	write(10,*) '# ',grid_y_left,grid_y_right,grid_dy
	do i=n_grid_x_left,n_grid_x_right
		do j=n_grid_y_left,n_grid_y_right
			write(10,*) i*grid_dx,j*grid_dy,reactiveprobdensity2(i,j)
		enddo
	enddo
	close(10)
	
	open(unit=10,file=file_name8,status='unknown')
	write(10,*) '# ',grid_x_left,grid_x_right,grid_dx
	write(10,*) '# ',grid_y_left,grid_y_right,grid_dy
	do i=n_grid_x_left,n_grid_x_right
		do j=n_grid_y_left,n_grid_y_right
			write(10,*) i*grid_dx,j*grid_dy,reactiveprobdensity3(i,j)
		enddo
	enddo
	close(10)

end program LangevinDynamics2d_active