program LangevinDynamics2d_active
use TPS2d_active_def_enhanced
implicit none

	integer i,counts_v,counts_a,l,tr,k,m,xint,yint,j,dummy
	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,fpt,xr,yr,r1,r2,r3,r4,x1,x2,y1,y2
	real*8 :: x_arr (0:757041)
	real*8 :: y_arr (0:757041)
	character (LEN=75) ::firstline
	integer :: istarts(0:376)
	integer :: istops(0:376)
	integer :: T_reach(0:376)

	file_name1 = '2_traj.txt'
	file_name2 = 'start_stop_trajs_2.dat'
	file_name3 = 'tpd1_exp2.dat'
	file_name4 = 'tpd2_exp2.dat'
	file_name5 = 'tpd3_exp2.dat'

	open(unit=10,file=file_name1,status='unknown')
	read(10,*) firstline
	do i = 0,757041
		read(10,*) dumm, x_arr(i), y_arr(i)
	enddo
	close(10)
	
	do i = 0,757041
		x_arr(i) = x_arr(i) - 616.d0
		y_arr(i) = y_arr(i) - 476.d0
	enddo
	
	tr = 0
	
	do i = 0,376
		istarts(i) = 0
		istops(i) = 0
		T_reach(i) = 0
	enddo
	reactiveprobdensity1(:,:)=0.d0
	reactiveprobdensity2(:,:)=0.d0
	reactiveprobdensity3(:,:)=0.d0
	
	l = 0
	do i = 0,757040
		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
				tr = tr+1
			else if (xnew .le. xtloc) then
				l = 0
				istarts(tr) = dummy
				istops(tr) = i
				T_reach(tr) = 1
				do m = dummy,i
					xr = x_arr(m)/467.d0
					yr = y_arr(m)/467.d0
					xint = Nint(xr/grid_dx)
					yint = Nint(yr/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
				tr = tr+1
			endif
		endif
	enddo
	
	l = 0
	do i = 0,757040
		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
				tr = tr+1
			else if (xnew .le. xtloc) then
				l = 0
				istarts(tr) = dummy
				istops(tr) = i
				T_reach(tr) = 1
				do m = dummy,i
					xr = x_arr(m)/467.d0
					yr = y_arr(m)/467.d0
					xint = Nint(xr/grid_dx)
					yint = Nint(yr/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
				tr = tr+1
			endif
		endif
	enddo
	
	l = 0
	do i = 0,757040
		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
				tr = tr+1
			else if (xnew .le. xtloc) then
				l = 0
				istarts(tr) = dummy
				istops(tr) = i
				T_reach(tr) = 1
				do m = dummy,i
					xr = x_arr(m)/467.d0
					yr = y_arr(m)/467.d0
					xint = Nint(xr/grid_dx)
					yint = Nint(yr/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
				tr = tr+1
			endif
		endif
	enddo
	
	print*, 'There have been ',tr,' complete trajectories.'
	
	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_name3,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_name4,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_name5,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)

	open(unit=10,file=file_name2,status='unknown')
	do i = 0,251
		write(10,*) istarts(i),istops(i),T_reach(i)
	enddo
	close(10)

end program LangevinDynamics2d_active