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,dummy,fpt
	real*8 :: x_arr (0:20000000)
	real*8 :: y_arr (0:20000000)
	real*8 :: theta_arr (0:20000000)
	integer :: istarts(0:34614)
	integer :: istops(0:34614)
	integer :: angular_slice(0:34614)
	integer :: T_reach(0:34614)
	

!~ 	file_name1 = 'trajectory_trial_400m_s.dat'
!~ 	file_name1 = 'trajectory_3tr_mypar.dat'
!~ 	file_name1 = 'trajectory_trial_parv001Dr005.dat'
!~ 	file_name1 = 'trajectory_trasl_parv001Dr005.dat'
	file_name1 = 'trajectory_final_par_v001Dr005_correct.dat'
!~ 	file_name1 = 'trajectory_passive.dat'
	file_name2 = 'start_stop_sim.dat'
!~ 	file_name3 = 'comm_left_sim_trial_parv001Dr005.dat'
!~ 	file_name4 = 'comm_up_sim_trial_parv001Dr005.dat'
!~ 	file_name5 = 'comm_right_sim_trial_parv001Dr005.dat'
!~ 	file_name3 = 'comm_left_sim_trasl_parv001Dr005.dat'
!~ 	file_name4 = 'comm_up_sim_trasl_parv001Dr005.dat'
!~ 	file_name5 = 'comm_right_sim_trasl_parv001Dr005.dat'
	file_name3 = 'comm_left_sim_final_par_v001Dr005_correct.dat'
	file_name4 = 'comm_up_sim_final_par_v001Dr005_correct.dat'
	file_name5 = 'comm_right_sim_final_par_v001Dr005_correct.dat'
	file_name6 = 'fpt_sim_final_par_v001Dr005_correct.dat'

	open(unit=10,file=file_name1,status='unknown')
	do i =0,20000000
!~ 		read(10,*) x_arr(i), y_arr(i), theta_arr(i)
		read(10,*) x_arr(i), y_arr(i)
!~ 		print*, x_arr(i)
	enddo
	close(10)
	
	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 = 0.05d0																			!0.05 is 5 times the hard radius of the particle, to be sure that these points are not influenced by the potential
	
!~ 	do i = 1,1999998
!~ 		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. 1000000) 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!'
!~ 					endif
!~ 				else
!~ 					print*, 'An error occurred!'
!~ 				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!'
!~ 					endif
!~ 				else
!~ 					print*, 'An error occurred!'
!~ 				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
!~ 		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
	
	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
	
	do i = 0,34614
		istarts(i) = 0
		istops(i) = 0
		angular_slice(i) = 0
		T_reach(i) = 0
	enddo
	
	l = 0
	do i = 0,19999999
		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!'
			endif
		else
			print*, 'An error occurred!'
		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)
!~ 				angular_slice(tr) = k
				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)
!~ 				angular_slice(tr) = k
!~ 				T_reach(tr) = 1
				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,19999999
		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!'
			endif
		else
			print*, 'An error occurred!'
		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)
!~ 				angular_slice(tr) = k
				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)
!~ 				angular_slice(tr) = k
!~ 				T_reach(tr) = 1
				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,19999999
		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!'
			endif
		else
			print*, 'An error occurred!'
		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)
!~ 				angular_slice(tr) = k
				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)
!~ 				angular_slice(tr) = k
!~ 				T_reach(tr) = 1
				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,34614
!~ 		write(10,*) istarts(i),istops(i),angular_slice(i),T_reach(i)
!~ 	enddo
!~ 	close(10)

	close(9)
	

end program LangevinDynamics2d_active
