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
	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
	integer :: istarts(0:29999)
	integer :: istops(0:29999)
	integer :: angular_slice(0:29999)
	integer :: T_reach(0:29999)
	integer :: point(0:29999)
	
	idum = -13
	iduv = -7
	phi = 0.d0
	file_name1 = 'comm_left_sim_circle.dat'
	file_name2 = 'comm_up_sim_circle.dat'
	file_name3 = 'comm_right_sim_circle.dat'
	file_name4 = 'start_stop_circle.dat'
	file_name5 = 'trajs_circle.dat'
	l = 0
	tr = 0
	j = 0
	
!~ 	do i = 0,19
!~ 		a = ranv()
!~ 		print*, a
!~ 	enddo
	
	do i = 0,7
		q1(i) = 0.d0
		q2(i) = 0.d0
		q3(i) = 0.d0
		trajang1(i) = 0.d0
		trajang2(i) = 0.d0
		trajang3(i) = 0.d0
	enddo
	
	do i = 0,3
		q1f(i) = 0.d0
		q2f(i) = 0.d0
		q3f(i) = 0.d0
		trajangf1(i) = 0.d0
		trajangf2(i) = 0.d0
		trajangf3(i) = 0.d0
		nr1(i) = 0.d0
		nr2(i) = 0.d0
		nr3(i) = 0.d0
		nt1(i) = 0.d0
		nt2(i) = 0.d0
		nt3(i) = 0.d0
	enddo
	
	do i = 0,29999
		istarts(i) = 0
		istops(i) = 0
		angular_slice(i) = 0
		T_reach(i) = 0
		point(i) = 0
	enddo
	
	open(unit=15,file=file_name5,status='unknown')
	
	i = 0
	do while (i .lt. 10000)
		x = xc1
		y = yc1
		l = 0
!~ 		j = 0
		a = ranv()
		theta_o = 2*pi*a
		b = ranv()
		vn = v + (b*0.1d0)
		call GaussianVariable (DW,DW2)
		w = DW*0.4d0
		do while (l .lt. 3)
			write(15,*) x,y
			call GaussianVariable (DW,DW2)
			xnew = x + vn*dcos(theta_o)*dt + fwx2d(x,y)/zeta
			ynew = y + vn*dsin(theta_o)*dt + fwy2d(x,y)/zeta
			call GaussianVariable (DW,DW2)
			theta_n = theta_o + w*dt + sigmar * DW
			rdistc = dsqrt((xnew-xc1)**2+(ynew-yc1)**2)
			rdistr = dsqrt(xnew**2+ynew**2)
			dx = xnew-x
			dy = ynew-y
!~ 			if ((dx .gt. 0.d0) .and. (dy .ge. 0.d0)) then
!~ 				thetai = atan(dy/dx)
!~ 			else if ((dx .lt. 0.d0) .and. (dy .ge. 0.d0)) then
!~ 				thetai = pi + atan(dy/dx)
!~ 			else if ((dx .lt. 0.d0) .and. (dy .lt. 0.d0)) then
!~ 				thetai = pi + atan(dy/dx)
!~ 			else if ((dx .gt. 0.d0) .and. (dy .lt. 0.d0)) then
!~ 				thetai = 2*pi + atan(dy/dx)
!~ 			else if (dx .eq. 0.d0) then
!~ 				if (dy .gt. 0.d0) then
!~ 					thetai = pi/2.d0
!~ 				else if (dy .lt. 0.d0) then
!~ 					thetai = 3*(pi/2.d0)
!~ 				else
!~ 					print*, 'An error occurred!', dx, dy
!~ 				endif
!~ 			else
!~ 				print*, 'An error occurred!'
!~ 			endif
			if ((rdistc .le. rc) .and. (l .eq. 0)) then
!~ 				inphi = thetai
				inphi = theta_o
				l = 1
			else if ((rdistc .le. rc) .and. (l .eq. 1)) then
!~ 				inphi = thetai
				inphi = theta_o
			else if ((rdistc .gt. rc) .and. ((l .eq. 1) .or. (l .eq. 2))) then
				if (l .eq. 1) dummy = j
				l = 2
				if (rdistr .le. rr) then
					l = 3
					istarts(tr) = dummy
					istops(tr) = j
					mm=mod(inphi,2*pi)
					if (mm.lt.0) mm=mm+2*pi
					k = int(mm/grid_dtheta)
					angular_slice(tr) = k
					point(tr) = 1
					trajang1(k) = trajang1(k) + 1.d0
				else if (xnew .le. xtloc) then
					l = 3
					istarts(tr) = dummy
					istops(tr) = j
					mm=mod(inphi,2*pi)
					if (mm.lt.0) mm=mm+2*pi
					k = int(mm/grid_dtheta)
					angular_slice(tr) = k
					T_reach(tr) = 1
					point(tr) = 1
					trajang1(k) = trajang1(k) + 1.d0
					q1(k) = q1(k) + 1.d0
				endif
			endif
			x = xnew
			y = ynew
			theta_o = theta_n
			j = j + 1
		enddo
		tr = tr + 1
		if (mod(i,100) .eq. 0) print*, i, xnew, ynew
		i = i+1
	enddo
	
	trajangf1(0) = trajang1(0)+trajang1(7)
	trajangf1(1) = trajang1(1)+trajang1(2)
	trajangf1(2) = trajang1(3)+trajang1(4)
	trajangf1(3) = trajang1(5)+trajang1(6)
	q1f(0) = q1(0)+q1(7)
	q1f(1) = q1(1)+q1(2)
	q1f(2) = q1(3)+q1(4)
	q1f(3) = q1(5)+q1(6)
	
	do i = 0,3
		nt1(i) = q1f(i)
		nr1(i) = trajangf1(i)-q1f(i)
	enddo
	
	i = 0
	do while (i .lt. 10000)
		x = xc2
		y = yc2
		l = 0
!~ 		j = 0
		a = ranv()
		theta_o = 2*pi*a
		b = ranv()
		vn = v + (b*0.1d0)
		call GaussianVariable (DW,DW2)
		w = DW*0.4d0
		do while (l .lt. 3)
			write(15,*) x,y
			call GaussianVariable (DW,DW2)
			xnew = x + vn*dcos(theta_o)*dt + fwx2d(x,y)/zeta
			ynew = y + vn*dsin(theta_o)*dt + fwy2d(x,y)/zeta
			call GaussianVariable (DW,DW2)
			theta_n = theta_o + w*dt + sigmar * DW
			rdistc = dsqrt((xnew-xc2)**2+(ynew-yc2)**2)
			rdistr = dsqrt(xnew**2+ynew**2)
			dx = xnew-x
			dy = ynew-y
!~ 			if ((dx .gt. 0.d0) .and. (dy .ge. 0.d0)) then
!~ 				thetai = atan(dy/dx)
!~ 			else if ((dx .lt. 0.d0) .and. (dy .ge. 0.d0)) then
!~ 				thetai = pi + atan(dy/dx)
!~ 			else if ((dx .lt. 0.d0) .and. (dy .lt. 0.d0)) then
!~ 				thetai = pi + atan(dy/dx)
!~ 			else if ((dx .gt. 0.d0) .and. (dy .lt. 0.d0)) then
!~ 				thetai = 2*pi + atan(dy/dx)
!~ 			else if (dx .eq. 0.d0) then
!~ 				if (dy .gt. 0.d0) then
!~ 					thetai = pi/2.d0
!~ 				else if (dy .lt. 0.d0) then
!~ 					thetai = 3*(pi/2.d0)
!~ 				else
!~ 					print*, 'An error occurred!'
!~ 				endif
!~ 			else
!~ 				print*, 'An error occurred!'
!~ 			endif
			if ((rdistc .le. rc) .and. (l .eq. 0)) then
!~ 				inphi = thetai
				inphi = theta_o
				l = 1
			else if ((rdistc .le. rc) .and. (l .eq. 1)) then
!~ 				inphi = thetai
				inphi = theta_o
			else if ((rdistc .gt. rc) .and. ((l .eq. 1) .or. (l .eq. 2))) then
				if (l .eq. 1) dummy = j
				l = 2
				if (rdistr .le. rr) then
					l = 3
					istarts(tr) = dummy
					istops(tr) = j
					mm=mod(inphi,2*pi)
					if (mm.lt.0) mm=mm+2*pi
					k = int(mm/grid_dtheta)
					angular_slice(tr) = k
					point(tr) = 2
					trajang2(k) = trajang2(k) + 1.d0
				else if (xnew .le. xtloc) then
					l = 3
					istarts(tr) = dummy
					istops(tr) = j
					mm=mod(inphi,2*pi)
					if (mm.lt.0) mm=mm+2*pi
					k = int(mm/grid_dtheta)
					angular_slice(tr) = k
					T_reach(tr) = 1
					point(tr) = 2
					trajang2(k) = trajang2(k) + 1.d0
					q2(k) = q2(k) + 1.d0
				endif
			endif
			x = xnew
			y = ynew
			theta_o = theta_n
			j = j + 1
		enddo
		tr = tr + 1
		if (mod(i,100) .eq. 0) print*, i, xnew, ynew
		i = i+1
	enddo
	
	trajangf2(0) = trajang2(0)+trajang2(7)
	trajangf2(1) = trajang2(1)+trajang2(2)
	trajangf2(2) = trajang2(3)+trajang2(4)
	trajangf2(3) = trajang2(5)+trajang2(6)
	q2f(0) = q2(0)+q2(7)
	q2f(1) = q2(1)+q2(2)
	q2f(2) = q2(3)+q2(4)
	q2f(3) = q2(5)+q2(6)
	
	do i = 0,3
		nt2(i) = q2f(i)
		nr2(i) = trajangf2(i)-q2f(i)
	enddo
	
	i = 0
	do while (i .lt. 10000)
		x = xc3
		y = yc3
		l = 0
!~ 		j = 0
		a = ranv()
		theta_o = 2*pi*a
		b = ranv()
		vn = v + (b*0.1d0)
		call GaussianVariable (DW,DW2)
		w = DW*0.4d0
		do while (l .lt. 3)
			write(15,*) x,y
			call GaussianVariable (DW,DW2)
			xnew = x + vn*dcos(theta_o)*dt + fwx2d(x,y)/zeta
			ynew = y + vn*dsin(theta_o)*dt + fwy2d(x,y)/zeta
			call GaussianVariable (DW,DW2)
			theta_n = theta_o + w*dt + sigmar * DW
			rdistc = dsqrt((xnew-xc3)**2+(ynew-yc3)**2)
			rdistr = dsqrt(xnew**2+ynew**2)
			dx = xnew-x
			dy = ynew-y
!~ 			if ((dx .gt. 0.d0) .and. (dy .ge. 0.d0)) then
!~ 				thetai = atan(dy/dx)
!~ 			else if ((dx .lt. 0.d0) .and. (dy .ge. 0.d0)) then
!~ 				thetai = pi + atan(dy/dx)
!~ 			else if ((dx .lt. 0.d0) .and. (dy .lt. 0.d0)) then
!~ 				thetai = pi + atan(dy/dx)
!~ 			else if ((dx .gt. 0.d0) .and. (dy .lt. 0.d0)) then
!~ 				thetai = 2*pi + atan(dy/dx)
!~ 			else if (dx .eq. 0.d0) then
!~ 				if (dy .gt. 0.d0) then
!~ 					thetai = pi/2.d0
!~ 				else if (dy .lt. 0.d0) then
!~ 					thetai = 3*(pi/2.d0)
!~ 				else
!~ 					print*, 'An error occurred!'
!~ 				endif
!~ 			else
!~ 				print*, 'An error occurred!'
!~ 			endif
			if ((rdistc .le. rc) .and. (l .eq. 0)) then
!~ 				inphi = thetai
				inphi = theta_o
				l = 1
			else if ((rdistc .le. rc) .and. (l .eq. 1)) then
!~ 				inphi = thetai
				inphi = theta_o
			else if ((rdistc .gt. rc) .and. ((l .eq. 1) .or. (l .eq. 2))) then
				if (l .eq. 1) dummy = j
				l = 2
				if (rdistr .le. rr) then
					l = 3
					istarts(tr) = dummy
					istops(tr) = j
					mm=mod(inphi,2*pi)
					if (mm.lt.0) mm=mm+2*pi
					k = int(mm/grid_dtheta)
					angular_slice(tr) = k
					point(tr) = 3
					trajang3(k) = trajang3(k) + 1.d0
				else if (xnew .le. xtloc) then
					l = 3
					istarts(tr) = dummy
					istops(tr) = j
					mm=mod(inphi,2*pi)
					if (mm.lt.0) mm=mm+2*pi
					k = int(mm/grid_dtheta)
					angular_slice(tr) = k
					T_reach(tr) = 1
					point(tr) = 3
					trajang3(k) = trajang3(k) + 1.d0
					q3(k) = q3(k) + 1.d0
				endif
			endif
			x = xnew
			y = ynew
			theta_o = theta_n
			j = j + 1
		enddo
		tr = tr + 1
		if (mod(i,100) .eq. 0) print*, i, xnew, ynew
		i = i+1
	enddo
	
	close(15)
	
	trajangf3(0) = trajang3(0)+trajang3(7)
	trajangf3(1) = trajang3(1)+trajang3(2)
	trajangf3(2) = trajang3(3)+trajang3(4)
	trajangf3(3) = trajang3(5)+trajang3(6)
	q3f(0) = q3(0)+q3(7)
	q3f(1) = q3(1)+q3(2)
	q3f(2) = q3(3)+q3(4)
	q3f(3) = q3(5)+q3(6)
	
	do i = 0,3
		nt3(i) = q3f(i)
		nr3(i) = trajangf3(i)-q3f(i)
	enddo
	
	do i = 0,3
		succ1(i) = q1f(i)
		if (trajangf1(i) .eq. 0d0) then
			q1f(i) = 2.d0
		else
			q1f(i) = q1f(i)/trajangf1(i)
		endif
		std1(i) = dsqrt( ( nt1(i)*(1-q1f(i))**2 + nr1(i)*(0-q1f(i))**2 )/(nt1(i)+nr1(i)) )			!standard deviation
		e1(i) = std1(i)/dsqrt(nt1(i)+nr1(i))														!error on the mean
		e_r1(i) = e1(i)/q1f(i)																		!relative error on the mean, must be smaller than 0.05
	enddo
	
	do i = 0,3
		succ2(i) = q2f(i)
		if (trajangf2(i) .eq. 0d0) then
			q2f(i) = 2.d0
		else
			q2f(i) = q2f(i)/trajangf2(i)
		endif
		std2(i) = dsqrt( ( nt2(i)*(1-q2f(i))**2 + nr2(i)*(0-q2f(i))**2 )/(nt2(i)+nr2(i)) )			!standard deviation
		e2(i) = std2(i)/dsqrt(nt2(i)+nr2(i))														!error on the mean
		e_r2(i) = e2(i)/q2f(i)																		!relative error on the mean, must be smaller than 0.05
	enddo
	
	do i = 0,3
		succ3(i) = q3f(i)
		if (trajangf3(i) .eq. 0d0) then
			q3f(i) = 2.d0
		else
			q3f(i) = q3f(i)/trajangf3(i)
		endif
		std3(i) = dsqrt( ( nt3(i)*(1-q3f(i))**2 + nr3(i)*(0-q3f(i))**2 )/(nt3(i)+nr3(i)) )			!standard deviation
		e3(i) = std3(i)/dsqrt(nt3(i)+nr3(i))														!error on the mean
		e_r3(i) = e3(i)/q3f(i)																		!relative error on the mean, must be smaller than 0.05
	enddo
	
	open(unit=10,file=file_name1,status='unknown')
	do i = 0,3
		write(10,*) q1f(i), std1(i), e1(i), e_r1(i), succ1(i), trajangf1(i)
	enddo
	close(10)
	
	open(unit=10,file=file_name2,status='unknown')
	do i = 0,3
		write(10,*) q2f(i), std2(i), e2(i), e_r2(i), succ2(i), trajangf2(i)
	enddo
	close(10)
	
	open(unit=10,file=file_name3,status='unknown')
	do i = 0,3
		write(10,*) q3f(i), std3(i), e3(i), e_r3(i), succ3(i), trajangf3(i)
	enddo
	close(10)
	
	open(unit=10,file=file_name4,status='unknown')
	do i = 0,29999
		write(10,*) istarts(i),istops(i),point(i),angular_slice(i),T_reach(i)
	enddo
	close(10)

end program LangevinDynamics2d_active