program LangevinDynamics2d_active
use TPS2d_active_def_enhanced
implicit none

	integer i,counts_v,counts_a,l,tr,k
	real*8 v_avg,rdist1,xnew,ynew,d_r,dx,dy,v_t,thetap,xp,yp,thetai,xn,yn,dtheta_sq,theta_avg,Dr_ext,cutoff,thetapp,thetaip
	real*8 rdistc,inphi,rdistr,mm,dumm,dummy, fpt
!~ 	real*8 :: x_arr (0:136349)							!map: the first block of the repeated stuff in the code refers to the parameters to analyze the first
!~ 	real*8 :: y_arr (0:136349)							!trajectory, the second to the second trajectory and so on.
!~ 	real*8 :: x_arr (0:757041)
!~ 	real*8 :: y_arr (0:757041)
!~ 	real*8 :: x_arr (0:487966)
!~ 	real*8 :: y_arr (0:487966)
!~ 	real*8 :: x_arr (0:278771)
!~ 	real*8 :: y_arr (0:278771)
	real*8 :: x_arr (0:271501)
	real*8 :: y_arr (0:271501)
	character (LEN=75) ::firstline
!~ 	real*8 :: theta_arr (0:20000000)
!~ 	integer :: istarts(0:282)
!~ 	integer :: istops(0:282)

	file_name1 = '5_traj.txt'
!~ 	file_name2 = 'start_stop_trajs_2old.dat'
	file_name3 = 'comm_left_exp_5.dat'
	file_name4 = 'comm_up_exp_5.dat'
	file_name5 = 'comm_right_exp_5.dat'
	file_name6 = 'fpt_exp_5.dat'

	open(unit=10,file=file_name1,status='unknown')
	read(10,*) firstline
!~ 	do i = 0,136349
!~ 	do i = 0,757041
!~ 	do i = 0,487966
!~ 	do i = 0,278771
	do i = 0,271501
!~ 		read(10,*) x_arr(i), y_arr(i), theta_arr(i)
		read(10,*) dumm, x_arr(i), y_arr(i)
	enddo
	close(10)
	
!~ 	do i = 0,136349
!~ 	do i = 0,757041
!~ 	do i = 0,487966
!~ 	do i = 0,278771
	do i = 0,271501
!~ 		x_arr(i) = x_arr(i) - 600.d0
!~ 		y_arr(i) = y_arr(i) - 518.d0
!~ 		x_arr(i) = x_arr(i) - 616.d0
!~ 		y_arr(i) = y_arr(i) - 476.d0
!~ 		x_arr(i) = x_arr(i) - 512.d0
!~ 		y_arr(i) = y_arr(i) - 530.d0
!~ 		x_arr(i) = x_arr(i) - 555.d0
!~ 		y_arr(i) = y_arr(i) - 548.d0
		x_arr(i) = x_arr(i) - 555.d0
		y_arr(i) = y_arr(i) - 569.d0
	enddo
	
!~ 	print*, x_arr(0), y_arr(0)

	open(unit=9,file=file_name6,status='unknown')
	fpt = 0.d0
	
	
	!UNCOMMENT THE NEXT THREE BLOCKS TO ESTIMATE V AND DROT
	
!~ 	v_avg = 0.d0
!~ 	theta_avg = 0.d0
!~ 	counts_v = 0
!~ 	counts_a = 0
!~ 	cutoff = 49.d0																			!49 is about 5 times the hard radius of the particle, to be sure that these points are not influenced by the potential
	
!~ 	do i = 0,136348
!~ 	do i = 0,757041
!~ 	do i = 1,487965
!~ 		xnew = x_arr(i)
!~ 		ynew = y_arr(i)
!~ 		xp = x_arr(i-1)
!~ 		yp = y_arr(i-1)
!~ 		xn = x_arr(i+1)
!~ 		yn = y_arr(i+1)
!~ 		rdist1 = radius-dsqrt(xnew**2+ynew**2)
!~ 		if ((rdist1 .gt. cutoff) .and. (((xnew .ne. xp) .or. (ynew .ne. yp)) .and. ((xnew .ne. xn) .or. (ynew .ne. yn)))) then
!~ 			dx = xnew-xp
!~ 			dy = ynew-yp
!~ 			d_r = dsqrt(dx**2+dy**2)
!~ 			v_t = d_r/dt
!~ 			v_avg = v_avg + v_t
!~ 			if (i .lt. 136349) then
!~ 			if (i .lt. 757041) then
!~ 			if (i .lt. 487966) then
!~ 				if ((dx .gt. 0.d0) .and. (dy .ge. 0.d0)) then
!~ 					thetap = atan(dy/dx)
!~ 				else if ((dx .lt. 0.d0) .and. (dy .ge. 0.d0)) then
!~ 					thetap = pi + atan(dy/dx)
!~ 				else if ((dx .lt. 0.d0) .and. (dy .lt. 0.d0)) then
!~ 					thetap = pi + atan(dy/dx)
!~ 				else if ((dx .gt. 0.d0) .and. (dy .lt. 0.d0)) then
!~ 					thetap = 2*pi + atan(dy/dx)
!~ 				else if (dx .eq. 0.d0) then
!~ 					if (dy .gt. 0.d0) then
!~ 						thetap = pi/2.d0
!~ 					else if (dy .lt. 0.d0) then
!~ 						thetap = 3*(pi/2.d0)
!~ 					else
!~ 						print*, 'An error occurred! 1'
!~ 					endif
!~ 				else
!~ 					print*, 'An error occurred! 2'
!~ 				endif
!~ 				if (((xn-xnew) .gt. 0.d0) .and. ((yn-ynew) .ge. 0.d0)) then
!~ 					thetai = atan((yn-ynew)/(xn-xnew))
!~ 				else if (((xn-xnew) .lt. 0.d0) .and. ((yn-ynew) .ge. 0.d0)) then
!~ 					thetai = pi + atan((yn-ynew)/(xn-xnew))
!~ 				else if (((xn-xnew) .lt. 0.d0) .and. ((yn-ynew) .lt. 0.d0)) then
!~ 					thetai = pi + atan((yn-ynew)/(xn-xnew))
!~ 				else if (((xn-xnew) .gt. 0.d0) .and. ((yn-ynew) .lt. 0.d0)) then
!~ 					thetai = 2*pi + atan((yn-ynew)/(xn-xnew))
!~ 				else if ((xn-xnew) .eq. 0.d0) then
!~ 					if ((yn-ynew) .gt. 0.d0) then
!~ 						thetai = pi/2.d0
!~ 					else if ((yn-ynew) .lt. 0.d0) then
!~ 						thetai = 3*(pi/2.d0)
!~ 					else
!~ 						print*, 'An error occurred! 3'
!~ 					endif
!~ 				else
!~ 					print*, 'An error occurred! 4'
!~ 				endif
!~ 				if (((thetap .ge. 0.d0) .and. (thetap .lt. (pi/2.d0))) .and. ((thetai .gt. (3*(pi/2.d0))) .and. (thetai .lt. (2*pi)))) then
!~ 					thetai = thetai - (2*pi)
!~ 				else if (((thetap .gt. (3*(pi/2.d0))) .and. (thetap .lt. (2*pi))) .and. ((thetai .ge. 0.d0) .and. (thetai .lt. (pi/2.d0)))) then
!~ 					thetap = thetap - (2*pi)
!~ 				endif
!~ 				dtheta_sq = (thetai-thetap)**2
!~ 				theta_avg = theta_avg + dtheta_sq
!~ 				counts_a = counts_a + 1
!~ 			endif
!~ 			counts_v = counts_v + 1
!~ 		else if ((rdist1 .gt. cutoff)) then
!~ 			v_t = 0.d0
!~ 			v_avg = v_avg + v_t
!~ 			counts_v = counts_v + 1
!~ 			print*, 'The particle did not move!'
!~ 		endif
!~ 	enddo
	
!~ 	v_avg = v_avg/counts_v
!~ 	theta_avg = theta_avg/counts_a
!~ 	Dr_ext = theta_avg/(2*dt)
!~ 	print*, v_avg, counts_v, Dr_ext, counts_a
	
	!the next three repeated parts are used to compute v and Dr at the beginning, middle and end of the trajectory, beginning and end are different
	!from middle because I used 3 blocks as in middle but also beg and end to estimate parameters roughly, while I used the block at beg and end only
	!to extract them more accurately but not at the middle
	
!~ 	v_avg = 0.d0
!~ 	theta_avg = 0.d0
!~ 	counts_v = 0
!~ 	counts_a = 0
!~ 	cutoff = 49.d0
!~ 	i = 0
	
!~ 	!do i = 0,100
!~ 	do while (counts_v < 1000)
!~ 		xnew = x_arr(i)
!~ 		ynew = y_arr(i)
!~ 		xp = x_arr(i-1)
!~ 		yp = y_arr(i-1)
!~ 		xn = x_arr(i+1)
!~ 		yn = y_arr(i+1)
!~ 		rdist1 = radius-dsqrt(xnew**2+ynew**2)
!~ 		if ((rdist1 .gt. cutoff) .and. (((xnew .ne. xp) .or. (ynew .ne. yp)) .and. ((xnew .ne. xn) .or. (ynew .ne. yn)))) then
!~ 			dx = xnew-xp
!~ 			dy = ynew-yp
!~ 			d_r = dsqrt(dx**2+dy**2)
!~ 			v_t = d_r/dt
!~ 			v_avg = v_avg + v_t
!~ 			if (i .lt. 136349) then
!~ 			if (i .lt. 757041) then
!~ 			if (i .lt. 487966) then
!~ 				if ((dx .gt. 0.d0) .and. (dy .ge. 0.d0)) then
!~ 					thetap = atan(dy/dx)
!~ 				else if ((dx .lt. 0.d0) .and. (dy .ge. 0.d0)) then
!~ 					thetap = pi + atan(dy/dx)
!~ 				else if ((dx .lt. 0.d0) .and. (dy .lt. 0.d0)) then
!~ 					thetap = pi + atan(dy/dx)
!~ 				else if ((dx .gt. 0.d0) .and. (dy .lt. 0.d0)) then
!~ 					thetap = 2*pi + atan(dy/dx)
!~ 				else if (dx .eq. 0.d0) then
!~ 					if (dy .gt. 0.d0) then
!~ 						thetap = pi/2.d0
!~ 					else if (dy .lt. 0.d0) then
!~ 						thetap = 3*(pi/2.d0)
!~ 					else
!~ 						!print*, 'An error occurred! 1'
!~ 					endif
!~ 				else
!~ 					!print*, 'An error occurred! 2'
!~ 				endif
!~ 				if (((xn-xnew) .gt. 0.d0) .and. ((yn-ynew) .ge. 0.d0)) then
!~ 					thetai = atan((yn-ynew)/(xn-xnew))
!~ 				else if (((xn-xnew) .lt. 0.d0) .and. ((yn-ynew) .ge. 0.d0)) then
!~ 					thetai = pi + atan((yn-ynew)/(xn-xnew))
!~ 				else if (((xn-xnew) .lt. 0.d0) .and. ((yn-ynew) .lt. 0.d0)) then
!~ 					thetai = pi + atan((yn-ynew)/(xn-xnew))
!~ 				else if (((xn-xnew) .gt. 0.d0) .and. ((yn-ynew) .lt. 0.d0)) then
!~ 					thetai = 2*pi + atan((yn-ynew)/(xn-xnew))
!~ 				else if ((xn-xnew) .eq. 0.d0) then
!~ 					if ((yn-ynew) .gt. 0.d0) then
!~ 						thetai = pi/2.d0
!~ 					else if ((yn-ynew) .lt. 0.d0) then
!~ 						thetai = 3*(pi/2.d0)
!~ 					else
!~ 						!print*, 'An error occurred! 3'
!~ 					endif
!~ 				else
!~ 					!print*, 'An error occurred! 4'
!~ 				endif
!~ 				if (((thetap .ge. 0.d0) .and. (thetap .lt. (pi/2.d0))) .and. ((thetai .gt. (3*(pi/2.d0))) .and. (thetai .lt. (2*pi)))) then
!~ 					thetai = thetai - (2*pi)
!~ 				else if (((thetap .gt. (3*(pi/2.d0))) .and. (thetap .lt. (2*pi))) .and. ((thetai .ge. 0.d0) .and. (thetai .lt. (pi/2.d0)))) then
!~ 					thetap = thetap - (2*pi)
!~ 				endif
!~ 				dtheta_sq = (thetai-thetap)**2
!~ 				!if (dtheta_sq .gt. 0.1d0) print*, dtheta_sq, thetai, thetap, i
!~ 				theta_avg = theta_avg + dtheta_sq
!~ 				counts_a = counts_a + 1
!~ 			endif
!~ 			counts_v = counts_v + 1
!~ 		else if ((rdist1 .gt. cutoff)) then
!~ 			v_t = 0.d0
!~ 			v_avg = v_avg + v_t
!~ 			counts_v = counts_v + 1
!~ 			!print*, 'The particle did not move!'
!~ 		endif
!~ 		i = i + 1
!~ 	enddo
	
!~ 	v_avg = v_avg/counts_v
!~ 	print*, 'v at the beginning', v_avg, counts_v
!~ 	theta_avg = 0.d0
!~ 	counts_a = 0
!~ 	cutoff = 49.d0
!~ 	i = 0
	
!~ 	do while (counts_a < 1000)
!~ 		xnew = x_arr(i)
!~ 		ynew = y_arr(i)
!~ 		xp = x_arr(i-1)
!~ 		yp = y_arr(i-1)
!~ 		xn = x_arr(i+1)
!~ 		yn = y_arr(i+1)
!~ 		rdist1 = radius-dsqrt(xnew**2+ynew**2)
!~ 		if ((rdist1 .gt. cutoff) .and. (((xnew .ne. xp) .or. (ynew .ne. yp)) .and. ((xnew .ne. xn) .or. (ynew .ne. yn)))) then
!~ 			dx = xnew-xp
!~ 			dy = ynew-yp
!~ 			d_r = dsqrt(dx**2+dy**2)
!~ 			!v_t = d_r/dt
!~ 			!v_avg = v_avg + v_t
!~ 			if (i .lt. 136349) then
!~ 			if (i .lt. 757041) then
!~ 			if (i .lt. 487966) then
!~ 				if ((dx .gt. 0.d0) .and. (dy .ge. 0.d0)) then
!~ 					thetap = atan(dy/dx)
!~ 				else if ((dx .lt. 0.d0) .and. (dy .ge. 0.d0)) then
!~ 					thetap = pi + atan(dy/dx)
!~ 				else if ((dx .lt. 0.d0) .and. (dy .lt. 0.d0)) then
!~ 					thetap = pi + atan(dy/dx)
!~ 				else if ((dx .gt. 0.d0) .and. (dy .lt. 0.d0)) then
!~ 					thetap = 2*pi + atan(dy/dx)
!~ 				else if (dx .eq. 0.d0) then
!~ 					if (dy .gt. 0.d0) then
!~ 						thetap = pi/2.d0
!~ 					else if (dy .lt. 0.d0) then
!~ 						thetap = 3*(pi/2.d0)
!~ 					else
!~ 						!print*, 'An error occurred! 1'
!~ 					endif
!~ 				else
!~ 					!print*, 'An error occurred! 2'
!~ 				endif
!~ 				if (((xn-xnew) .gt. 0.d0) .and. ((yn-ynew) .ge. 0.d0)) then
!~ 					thetai = atan((yn-ynew)/(xn-xnew))
!~ 				else if (((xn-xnew) .lt. 0.d0) .and. ((yn-ynew) .ge. 0.d0)) then
!~ 					thetai = pi + atan((yn-ynew)/(xn-xnew))
!~ 				else if (((xn-xnew) .lt. 0.d0) .and. ((yn-ynew) .lt. 0.d0)) then
!~ 					thetai = pi + atan((yn-ynew)/(xn-xnew))
!~ 				else if (((xn-xnew) .gt. 0.d0) .and. ((yn-ynew) .lt. 0.d0)) then
!~ 					thetai = 2*pi + atan((yn-ynew)/(xn-xnew))
!~ 				else if ((xn-xnew) .eq. 0.d0) then
!~ 					if ((yn-ynew) .gt. 0.d0) then
!~ 						thetai = pi/2.d0
!~ 					else if ((yn-ynew) .lt. 0.d0) then
!~ 						thetai = 3*(pi/2.d0)
!~ 					else
!~ 						!print*, 'An error occurred! 3'
!~ 					endif
!~ 				else
!~ 					!print*, 'An error occurred! 4'
!~ 				endif
!~ 				if (((thetap .ge. 0.d0) .and. (thetap .lt. (pi/2.d0))) .and. ((thetai .gt. (3*(pi/2.d0))) .and. (thetai .lt. (2*pi)))) then
!~ 					thetai = thetai - (2*pi)
!~ 				else if (((thetap .gt. (3*(pi/2.d0))) .and. (thetap .lt. (2*pi))) .and. ((thetai .ge. 0.d0) .and. (thetai .lt. (pi/2.d0)))) then
!~ 					thetap = thetap - (2*pi)
!~ 				endif
!~ 				dtheta_sq = (thetai-thetap)**2
!~ 				!if (dtheta_sq .gt. 0.1d0) print*, dtheta_sq, thetai, thetap, i
!~ 				theta_avg = theta_avg + dtheta_sq
!~ 				counts_a = counts_a + 1
!~ 			endif
!~ 			!counts_v = counts_v + 1
!~ 		else if ((rdist1 .gt. cutoff)) then
			!v_t = 0.d0
			!v_avg = v_avg + v_t
			!counts_v = counts_v + 1
			!print*, 'The particle did not move!'
!~ 		endif
!~ 		i = i + 1
!~ 	enddo
	
!~ 	theta_avg = theta_avg/counts_a
!~ 	Dr_ext = theta_avg/(2*dt)
!~ 	print*, 'Dr at the beginning', Dr_ext, counts_a
	
!~ 	v_avg = 0.d0
!~ 	theta_avg = 0.d0
!~ 	counts_v = 0
!~ 	counts_a = 0
!~ 	cutoff = 49.d0
	
!~ 	do i = 350000,350100
!~ 		xnew = x_arr(i)
!~ 		ynew = y_arr(i)
!~ 		xp = x_arr(i-1)
!~ 		yp = y_arr(i-1)
!~ 		xn = x_arr(i+1)
!~ 		yn = y_arr(i+1)
!~ 		rdist1 = radius-dsqrt(xnew**2+ynew**2)
!~ 		if ((rdist1 .gt. cutoff) .and. (((xnew .ne. xp) .or. (ynew .ne. yp)) .and. ((xnew .ne. xn) .or. (ynew .ne. yn)))) then
!~ 			dx = xnew-xp
!~ 			dy = ynew-yp
!~ 			d_r = dsqrt(dx**2+dy**2)
!~ 			v_t = d_r/dt
!~ 			v_avg = v_avg + v_t
!~ 			!if (i .lt. 136349) then
!~ 			if (i .lt. 757041) then
!~ 			!if (i .lt. 487966) then
!~ 				if ((dx .gt. 0.d0) .and. (dy .ge. 0.d0)) then
!~ 					thetap = atan(dy/dx)
!~ 				else if ((dx .lt. 0.d0) .and. (dy .ge. 0.d0)) then
!~ 					thetap = pi + atan(dy/dx)
!~ 				else if ((dx .lt. 0.d0) .and. (dy .lt. 0.d0)) then
!~ 					thetap = pi + atan(dy/dx)
!~ 				else if ((dx .gt. 0.d0) .and. (dy .lt. 0.d0)) then
!~ 					thetap = 2*pi + atan(dy/dx)
!~ 				else if (dx .eq. 0.d0) then
!~ 					if (dy .gt. 0.d0) then
!~ 						thetap = pi/2.d0
!~ 					else if (dy .lt. 0.d0) then
!~ 						thetap = 3*(pi/2.d0)
!~ 					else
!~ 						!print*, 'An error occurred! 1'
!~ 					endif
!~ 				else
!~ 					!print*, 'An error occurred! 2'
!~ 				endif
!~ 				if (((xn-xnew) .gt. 0.d0) .and. ((yn-ynew) .ge. 0.d0)) then
!~ 					thetai = atan((yn-ynew)/(xn-xnew))
!~ 				else if (((xn-xnew) .lt. 0.d0) .and. ((yn-ynew) .ge. 0.d0)) then
!~ 					thetai = pi + atan((yn-ynew)/(xn-xnew))
!~ 				else if (((xn-xnew) .lt. 0.d0) .and. ((yn-ynew) .lt. 0.d0)) then
!~ 					thetai = pi + atan((yn-ynew)/(xn-xnew))
!~ 				else if (((xn-xnew) .gt. 0.d0) .and. ((yn-ynew) .lt. 0.d0)) then
!~ 					thetai = 2*pi + atan((yn-ynew)/(xn-xnew))
!~ 				else if ((xn-xnew) .eq. 0.d0) then
!~ 					if ((yn-ynew) .gt. 0.d0) then
!~ 						thetai = pi/2.d0
!~ 					else if ((yn-ynew) .lt. 0.d0) then
!~ 						thetai = 3*(pi/2.d0)
!~ 					else
!~ 						!print*, 'An error occurred! 3'
!~ 					endif
!~ 				else
!~ 					!print*, 'An error occurred! 4'
!~ 				endif
!~ 				if (((thetap .ge. 0.d0) .and. (thetap .lt. (pi/2.d0))) .and. ((thetai .gt. (3*(pi/2.d0))) .and. (thetai .lt. (2*pi)))) then
!~ 					thetai = thetai - (2*pi)
!~ 				else if (((thetap .gt. (3*(pi/2.d0))) .and. (thetap .lt. (2*pi))) .and. ((thetai .ge. 0.d0) .and. (thetai .lt. (pi/2.d0)))) then
!~ 					thetap = thetap - (2*pi)
!~ 				endif
!~ 				dtheta_sq = (thetai-thetap)**2
!~ 				!if (dtheta_sq .gt. 0.1d0) print*, dtheta_sq, thetai, thetap, i
!~ 				theta_avg = theta_avg + dtheta_sq
!~ 				counts_a = counts_a + 1
!~ 			endif
!~ 			counts_v = counts_v + 1
!~ 		else if ((rdist1 .gt. cutoff)) then
!~ 			v_t = 0.d0
!~ 			v_avg = v_avg + v_t
!~ 			counts_v = counts_v + 1
!~ 			!print*, 'The particle did not move!'
!~ 		endif
!~ 	enddo
	
!~ 	v_avg = v_avg/counts_v
!~ 	theta_avg = theta_avg/counts_a
!~ 	Dr_ext = theta_avg/(2*dt)
!~ 	print*, 'Parameters at the middle', v_avg, counts_v, Dr_ext, counts_a
	
!~ 	v_avg = 0.d0
!~ 	theta_avg = 0.d0
!~ 	counts_v = 0
!~ 	counts_a = 0
!~ 	cutoff = 49.d0
!~ 	i = 132349
!~ 	i = 755041
!~ 	i = 483966
	
	!do i = 756941,757041
!~ 	do while (counts_v<1000)
!~ 		xnew = x_arr(i)
!~ 		ynew = y_arr(i)
!~ 		xp = x_arr(i-1)
!~ 		yp = y_arr(i-1)
!~ 		xn = x_arr(i+1)
!~ 		yn = y_arr(i+1)
!~ 		rdist1 = radius-dsqrt(xnew**2+ynew**2)
!~ 		if ((rdist1 .gt. cutoff) .and. (((xnew .ne. xp) .or. (ynew .ne. yp)) .and. ((xnew .ne. xn) .or. (ynew .ne. yn)))) then
!~ 			dx = xnew-xp
!~ 			dy = ynew-yp
!~ 			d_r = dsqrt(dx**2+dy**2)
!~ 			v_t = d_r/dt
!~ 			v_avg = v_avg + v_t
!~ 			if (i .lt. 136349) then
!~ 			if (i .lt. 757041) then
!~ 			if (i .lt. 487966) then
!~ 				if ((dx .gt. 0.d0) .and. (dy .ge. 0.d0)) then
!~ 					thetap = atan(dy/dx)
!~ 				else if ((dx .lt. 0.d0) .and. (dy .ge. 0.d0)) then
!~ 					thetap = pi + atan(dy/dx)
!~ 				else if ((dx .lt. 0.d0) .and. (dy .lt. 0.d0)) then
!~ 					thetap = pi + atan(dy/dx)
!~ 				else if ((dx .gt. 0.d0) .and. (dy .lt. 0.d0)) then
!~ 					thetap = 2*pi + atan(dy/dx)
!~ 				else if (dx .eq. 0.d0) then
!~ 					if (dy .gt. 0.d0) then
!~ 						thetap = pi/2.d0
!~ 					else if (dy .lt. 0.d0) then
!~ 						thetap = 3*(pi/2.d0)
!~ 					else
!~ 						!print*, 'An error occurred! 1'
!~ 					endif
!~ 				else
!~ 					!print*, 'An error occurred! 2'
!~ 				endif
!~ 				if (((xn-xnew) .gt. 0.d0) .and. ((yn-ynew) .ge. 0.d0)) then
!~ 					thetai = atan((yn-ynew)/(xn-xnew))
!~ 				else if (((xn-xnew) .lt. 0.d0) .and. ((yn-ynew) .ge. 0.d0)) then
!~ 					thetai = pi + atan((yn-ynew)/(xn-xnew))
!~ 				else if (((xn-xnew) .lt. 0.d0) .and. ((yn-ynew) .lt. 0.d0)) then
!~ 					thetai = pi + atan((yn-ynew)/(xn-xnew))
!~ 				else if (((xn-xnew) .gt. 0.d0) .and. ((yn-ynew) .lt. 0.d0)) then
!~ 					thetai = 2*pi + atan((yn-ynew)/(xn-xnew))
!~ 				else if ((xn-xnew) .eq. 0.d0) then
!~ 					if ((yn-ynew) .gt. 0.d0) then
!~ 						thetai = pi/2.d0
!~ 					else if ((yn-ynew) .lt. 0.d0) then
!~ 						thetai = 3*(pi/2.d0)
!~ 					else
!~ 						!print*, 'An error occurred! 3'
!~ 					endif
!~ 				else
!~ 					!print*, 'An error occurred! 4'
!~ 				endif
!~ 				if (((thetap .ge. 0.d0) .and. (thetap .lt. (pi/2.d0))) .and. ((thetai .gt. (3*(pi/2.d0))) .and. (thetai .lt. (2*pi)))) then
!~ 					thetai = thetai - (2*pi)
!~ 				else if (((thetap .gt. (3*(pi/2.d0))) .and. (thetap .lt. (2*pi))) .and. ((thetai .ge. 0.d0) .and. (thetai .lt. (pi/2.d0)))) then
!~ 					thetap = thetap - (2*pi)
!~ 				endif
!~ 				dtheta_sq = (thetai-thetap)**2
!~ 				theta_avg = theta_avg + dtheta_sq
!~ 				counts_a = counts_a + 1
!~ 			endif
!~ 			counts_v = counts_v + 1
!~ 		else if ((rdist1 .gt. cutoff)) then
!~ 			v_t = 0.d0
!~ 			v_avg = v_avg + v_t
!~ 			counts_v = counts_v + 1
!~ 		endif
!~ 		i = i + 1
!~ 	enddo
	
!~ 	v_avg = v_avg/counts_v
!~ 	print*, 'v at the end', v_avg, counts_v
!~ 	theta_avg = 0.d0
!~ 	counts_a = 0
!~ 	cutoff = 49.d0
!~ 	i = 132349
!~ 	i = 753041
!~ 	i = 483966
	
!~ 	do while (counts_a<1000)
!~ 		xnew = x_arr(i)
!~ 		ynew = y_arr(i)
!~ 		xp = x_arr(i-1)
!~ 		yp = y_arr(i-1)
!~ 		xn = x_arr(i+1)
!~ 		yn = y_arr(i+1)
!~ 		rdist1 = radius-dsqrt(xnew**2+ynew**2)
!~ 		if ((rdist1 .gt. cutoff) .and. (((xnew .ne. xp) .or. (ynew .ne. yp)) .and. ((xnew .ne. xn) .or. (ynew .ne. yn)))) then
!~ 			dx = xnew-xp
!~ 			dy = ynew-yp
!~ 			d_r = dsqrt(dx**2+dy**2)
!~ 			!v_t = d_r/dt
!~ 			!v_avg = v_avg + v_t
!~ 			if (i .lt. 136349) then
!~ 			if (i .lt. 757041) then
!~ 			if (i .lt. 487966) then
!~ 				if ((dx .gt. 0.d0) .and. (dy .ge. 0.d0)) then
!~ 					thetap = atan(dy/dx)
!~ 				else if ((dx .lt. 0.d0) .and. (dy .ge. 0.d0)) then
!~ 					thetap = pi + atan(dy/dx)
!~ 				else if ((dx .lt. 0.d0) .and. (dy .lt. 0.d0)) then
!~ 					thetap = pi + atan(dy/dx)
!~ 				else if ((dx .gt. 0.d0) .and. (dy .lt. 0.d0)) then
!~ 					thetap = 2*pi + atan(dy/dx)
!~ 				else if (dx .eq. 0.d0) then
!~ 					if (dy .gt. 0.d0) then
!~ 						thetap = pi/2.d0
!~ 					else if (dy .lt. 0.d0) then
!~ 						thetap = 3*(pi/2.d0)
!~ 					else
!~ 						!print*, 'An error occurred! 1'
!~ 					endif
!~ 				else
!~ 					!print*, 'An error occurred! 2'
!~ 				endif
!~ 				if (((xn-xnew) .gt. 0.d0) .and. ((yn-ynew) .ge. 0.d0)) then
!~ 					thetai = atan((yn-ynew)/(xn-xnew))
!~ 				else if (((xn-xnew) .lt. 0.d0) .and. ((yn-ynew) .ge. 0.d0)) then
!~ 					thetai = pi + atan((yn-ynew)/(xn-xnew))
!~ 				else if (((xn-xnew) .lt. 0.d0) .and. ((yn-ynew) .lt. 0.d0)) then
!~ 					thetai = pi + atan((yn-ynew)/(xn-xnew))
!~ 				else if (((xn-xnew) .gt. 0.d0) .and. ((yn-ynew) .lt. 0.d0)) then
!~ 					thetai = 2*pi + atan((yn-ynew)/(xn-xnew))
!~ 				else if ((xn-xnew) .eq. 0.d0) then
!~ 					if ((yn-ynew) .gt. 0.d0) then
!~ 						thetai = pi/2.d0
!~ 					else if ((yn-ynew) .lt. 0.d0) then
!~ 						thetai = 3*(pi/2.d0)
!~ 					else
!~ 						!print*, 'An error occurred! 3'
!~ 					endif
!~ 				else
!~ 					!print*, 'An error occurred! 4'
!~ 				endif
!~ 				if (((thetap .ge. 0.d0) .and. (thetap .lt. (pi/2.d0))) .and. ((thetai .gt. (3*(pi/2.d0))) .and. (thetai .lt. (2*pi)))) then
!~ 					thetai = thetai - (2*pi)
!~ 				else if (((thetap .gt. (3*(pi/2.d0))) .and. (thetap .lt. (2*pi))) .and. ((thetai .ge. 0.d0) .and. (thetai .lt. (pi/2.d0)))) then
!~ 					thetap = thetap - (2*pi)
!~ 				endif
!~ 				dtheta_sq = (thetai-thetap)**2
!~ 				theta_avg = theta_avg + dtheta_sq
!~ 				counts_a = counts_a + 1
!~ 			endif
!~ 			!counts_v = counts_v + 1
!~ 		else if ((rdist1 .gt. cutoff)) then
!~ 			!v_t = 0.d0
!~ 			!v_avg = v_avg + v_t
!~ 			!counts_v = counts_v + 1
!~ 		endif
!~ 		i = i + 1
!~ 	enddo
	
!~ 	theta_avg = theta_avg/counts_a
!~ 	Dr_ext = theta_avg/(2*dt)
!~ 	print*, 'Dr at the end', Dr_ext, counts_a
	
	tr = 0
	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
	
	l = 0
!~ 	do i = 0,136348
!~ 	do i = 0,757040
!~ 	do i = 0,487965
!~ 	do i = 0,278770
	do i = 0,271500
		xn = x_arr(i+1)
		yn = y_arr(i+1)
		xnew = x_arr(i)
		ynew = y_arr(i)
		dx = xn-xnew
		dy = yn-ynew
		rdistc = dsqrt((xnew-xc1)**2+(ynew-yc1)**2)
		rdistr = dsqrt(xnew**2+ynew**2)
		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! 5', i
			endif
		else
!~ 			print*, 'An error occurred! 6'
		endif
		if ((rdistc .le. rc) .and. (l .eq. 0)) then
			inphi = thetai
			l = 1
		else if ((rdistc .le. rc) .and. (l .eq. 1)) then
			inphi = thetai
		else if ((rdistc .gt. rc) .and. ((l .eq. 1) .or. (l .eq. 2))) then
			if (l .eq. 1) dummy = i
			l = 2
			if (rdistr .le. rr) then
				l = 0
!~ 				istarts(tr) = dummy
!~ 				istops(tr) = i
				mm=mod(inphi,2*pi)
				if (mm.lt.0) mm=mm+2*pi
				k = int(mm/grid_dtheta)
				trajang1(k) = trajang1(k) + 1.d0
				tr = tr+1
			else if (xnew .le. xtloc) then
				l = 0
!~ 				istarts(tr) = dummy
!~ 				istops(tr) = i
				mm=mod(inphi,2*pi)
				if (mm.lt.0) mm=mm+2*pi
				k = int(mm/grid_dtheta)
				trajang1(k) = trajang1(k) + 1.d0
				q1(k) = q1(k) + 1.d0
				tr = tr+1
				fpt = i-dummy
				write(9,*) fpt, k, 1
			endif
		endif
	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
	
	l = 0
!~ 	do i = 0,136348
!~ 	do i = 0,757040
!~ 	do i = 0,487965
!~ 	do i = 0,278770
	do i = 0,271500
		xn = x_arr(i+1)
		yn = y_arr(i+1)
		xnew = x_arr(i)
		ynew = y_arr(i)
		dx = xn-xnew
		dy = yn-ynew
		rdistc = dsqrt((xnew-xc2)**2+(ynew-yc2)**2)
		rdistr = dsqrt(xnew**2+ynew**2)
		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! 7'
			endif
		else
!~ 			print*, 'An error occurred! 8'
		endif
		if ((rdistc .le. rc) .and. (l .eq. 0)) then
			inphi = thetai
			l = 1
		else if ((rdistc .le. rc) .and. (l .eq. 1)) then
			inphi = thetai
		else if ((rdistc .gt. rc) .and. ((l .eq. 1) .or. (l .eq. 2))) then
			if (l .eq. 1) dummy = i
			l = 2
			if (rdistr .le. rr) then
				l = 0
!~ 				istarts(tr) = dummy
!~ 				istops(tr) = i
				mm=mod(inphi,2*pi)
				if (mm.lt.0) mm=mm+2*pi
				k = int(mm/grid_dtheta)
				trajang2(k) = trajang2(k) + 1.d0
				tr = tr+1
			else if (xnew .le. xtloc) then
				l = 0
!~ 				istarts(tr) = dummy
!~ 				istops(tr) = i
				mm=mod(inphi,2*pi)
				if (mm.lt.0) mm=mm+2*pi
				k = int(mm/grid_dtheta)
				trajang2(k) = trajang2(k) + 1.d0
				q2(k) = q2(k) + 1.d0
				tr = tr+1
				fpt = i-dummy
				write(9,*) fpt, k, 2
			endif
		endif
	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
	
	l = 0
!~ 	do i = 0,136348
!~ 	do i = 0,757040
!~ 	do i = 0,487965
!~ 	do i = 0,278770
	do i = 0,271500
		xn = x_arr(i+1)
		yn = y_arr(i+1)
		xnew = x_arr(i)
		ynew = y_arr(i)
		dx = xn-xnew
		dy = yn-ynew
		rdistc = dsqrt((xnew-xc3)**2+(ynew-yc3)**2)
		rdistr = dsqrt(xnew**2+ynew**2)
		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! 9'
			endif
		else
!~ 			print*, 'An error occurred! 0'
		endif
		if ((rdistc .le. rc) .and. (l .eq. 0)) then
			inphi = thetai
			l = 1
		else if ((rdistc .le. rc) .and. (l .eq. 1)) then
			inphi = thetai
		else if ((rdistc .gt. rc) .and. ((l .eq. 1) .or. (l .eq. 2))) then
			if (l .eq. 1) dummy = i
			l = 2
			if (rdistr .le. rr) then
				l = 0
!~ 				istarts(tr) = dummy
!~ 				istops(tr) = i
				mm=mod(inphi,2*pi)
				if (mm.lt.0) mm=mm+2*pi
				k = int(mm/grid_dtheta)
				trajang3(k) = trajang3(k) + 1.d0
				tr = tr+1
			else if (xnew .le. xtloc) then
				l = 0
!~ 				istarts(tr) = dummy
!~ 				istops(tr) = i
				mm=mod(inphi,2*pi)
				if (mm.lt.0) mm=mm+2*pi
				k = int(mm/grid_dtheta)
				trajang3(k) = trajang3(k) + 1.d0
				q3(k) = q3(k) + 1.d0
				tr = tr+1
				fpt = i-dummy
				write(9,*) fpt, k, 3
			endif
		endif
	enddo
	
	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
	
	print*, 'There have been ',tr,' complete trajectories.'
!~ 	print*, q1
!~ 	print*, trajang1
!~ 	print*, q2
!~ 	print*, trajang2
!~ 	print*, q3
!~ 	print*, trajang3
	
	open(unit=10,file=file_name3,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_name4,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_name5,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_name2,status='unknown')
!~ 	do i = 0,282
!~ 		write(10,*) istarts(i),istops(i)
!~ 	enddo
!~ 	close(10)
	
!~ 	open(unit=10,file='final_traj.dat',status='unknown')
!~ 	do i = 706518,757041
!~ 		write(10,*) x_arr(i),y_arr(i)
!~ 	enddo
!~ 	close(10)

	close(9)

end program LangevinDynamics2d_active