PROGRAM DOUBLE_SLIT_EXP
!*******************************************************************************
!	Simulation program for particles passing a double-slit
!	The following recursive equations are solved for time step i+1:
!	r[i+1]   = r[i] + vr*delta_t
!	phi[i+1] = phi[i] + v/r[i+1] * delta_t
!
!	Variables:
!	r[]:		radius of mass
!	phi[]:		polar angle of mass
!	v:		constant tangential velocity (due to constant torsion)
!	vr:		radial drift velocity
!	tmax:		maximum simulated time
!	delta_t:	time difference per iteration step
!	
!*******************************************************************************
	IMPLICIT NONE

	CHARACTER*80 text
	INTEGER nfmax, nscreen, nmax, nbar, nslit, nout, n, m, nglob, itmax, it
	parameter (nfmax=4)
	REAL*8 t, tmax, delta_t, vx0, vy0
	REAL*8 xd, yd, xs, q, y1, ys, rx0, yrange, ya(nfmax)
	REAL*8, ALLOCATABLE, DIMENSION (:) :: rx, ry, vx, vy, rx1, ry1, vx1, vy1, rbarx, rbary, screen
	REAL*8, ALLOCATABLE, DIMENSION (:,:) :: y
	REAL*8, PARAMETER                  :: TWOPI = 6.2831853071795864769252d0, eps = 1.d-10
	
	call input(text, nmax, nbar, nslit, nout, nscreen, tmax, delta_t, vx0, vy0, &
		rbarx, rbary, xd, yd, xs, ys, y1, q, rx0, yrange)
	
	ALLOCATE (rx(nmax), ry(nmax), vx(nmax),vy(nmax), rx1(nmax), ry1(nmax), &
		vx1(nmax), vy1(nmax), y(nmax,nfmax), screen(nscreen))
	screen = 0.
	
	! maximum number of iterations
	itmax = int(tmax/delta_t + 1.e-8)
	print *, 'itmax ', itmax

	! initial particle and velocity distribution
	if (nslit.eq.1) then
		do n=1,nmax
			rx(n)  = rx0
			ry(n)  = -(ys + y1/2.) + n*(2.*ys + y1)/dfloat(nmax+1)
			vx(n)  = vx0
			vy(n)  = vy0
			y(n,1)=rx(n)
			y(n,2)=ry(n)
			y(n,3)=vx(n)
			y(n,4)=vy(n)
		enddo
	else if (nslit.eq.2) then
		do n=1,nmax
			rx(n)  = rx0
			if (n <= nmax/2) then
				ry(n)  = -ys - y1/2. + n*ys/dfloat(nmax/2.+1)
			else
				ry(n)  = ry(n-nmax/2)+y1+ys
			endif
			vx(n)  = vx0
			vy(n)  = vy0
			y(n,1)=rx(n)
			y(n,2)=ry(n)
			y(n,3)=vx(n)
			y(n,4)=vy(n)
		enddo
	endif
	call output(0, 0.d0, y, nmax, rbarx, rbary, nbar)
	
	! loop over time steps
	do it=1,itmax
		t = it*delta_t
		do n=1,nmax
			!nglob = n
			ya(1:4)=y(n,1:4)
					rx1(n)=ya(1)
					ry1(n)=ya(2)
			call runge(4, delta_t, t, ya)
			call constraints(y, ya, xd, yd, screen, nmax, nfmax, nscreen, nslit, n)
			y(n,1:4)=ya(1:4)
		enddo
		if (mod(it,itmax/nout).eq.0) &
		call output(it, t, y, nmax, rbarx, rbary, nbar)
	enddo
	close (20)
	! write screen histogram for gnuplot
	call outhist(screen, yd, nscreen)

CONTAINS

!*******************************************************************************
SUBROUTINE runge(nf, h, t, ya)
!*******************************************************************************
	IMPLICIT none
	INTEGER nf, nfmax, i
	REAL*8 h, t, ya(nf)
	!external dgl
	parameter (nfmax=10)
	REAL*8 k1(nfmax), k2(nfmax), k3(nfmax), k4(nfmax), yhilf(nfmax)
	if (nf.gt.nfmax) stop '"nf bigger than nfmax in RUNGE"'

	call dgl(nf, t, ya, k1)
	do i=1,nf
		yhilf(i) = ya(i) + h*k1(i)/2.
	enddo
	
	call dgl(nf, t+h/2., yhilf, k2)
	do i=1,nf
		yhilf(i) = ya(i) + h*k2(i)/2.
	enddo
	
	call dgl(nf, t+h/2., yhilf, k3)
	do i=1,nf
		yhilf(i) = ya(i) + h*k3(i)/2.
	enddo
	
	call dgl(nf, t+h, yhilf, k4)
	do i=1,nf
		ya(i) = ya(i) + h*(k1(i)/6. + k2(i)/3. + k3(i)/3. + k4(i)/6.)
	enddo
	return
END SUBROUTINE runge

!*******************************************************************************
SUBROUTINE dgl(nf, t, y, f)
! right-hand side of diff. eq.
! input:
! y(1) = rx
! y(2) = ry
! y(3) = vx
! y(4) = vy
! output:
! f(1..4) = RHS
!*******************************************************************************
	IMPLICIT none
	INTEGER nf, nd
	REAL*8 t, y(nf), f(nf)
	REAL*8 Fx, Fy, xnorm3
 
	! check for non-divergence
	!print*,t,y(1),y(2)
	do nd=1,nf
		if (dabs(y(nd)).gt.1.d100) then
			print*, 't =', t, '    y(', nd, ') = ', y(nd)
			stop '****** Divegence of variables!'
		endif
	enddo

	! compute right-hand side of DEQ
	Fx = 0.
	Fy = 0.
	do m=1,nbar
		xnorm3 = ((y(1)-rbarx(m))**2 + (y(2)-rbary(m))**2)**1.5
		if (abs(xnorm3).lt.eps) xnorm3=eps
		Fx = Fx + q*(y(1)-rbarx(m))/xnorm3
		Fy = Fy + q*(y(2)-rbary(m))/xnorm3
	enddo
	! Interaction between particles, gives no reasonable results!
	!do m=1,nmax
	!	if(m.lt.nglob) then
	!		xnorm3 = ((rx1(nglob)-rx1(m))**2 + (ry1(nglob)-ry1(m))**2)**1.5
	!		if (abs(xnorm3).lt.eps) xnorm3=eps
	!		Fx = Fx + q*(rx1(nglob)-rx1(m))/xnorm3
	!		Fy = Fy + q*(ry1(nglob)-ry1(m))/xnorm3
	!	endif
	!enddo

	f(3) = Fx
	f(4) = Fy
	f(1) = y(3)
	f(2) = y(4)
END SUBROUTINE dgl

!*******************************************************************************
SUBROUTINE constraints(y, ya, xd, yd, screen, nmax, nf, nscreen, nslit, n)
! Repel particle if it has diffused through the border
!*******************************************************************************
	IMPLICIT none
	INTEGER nmax, nf, nscreen, nslit, n, index
	REAL*8 y(nmax, nf), ya(nf), xd, yd, screen(nscreen)

	!check for zero transition of x coordinate
	if (y(n,1)<=0. .and. ya(1)>0.) then
		!invert x and vx
		if (nslit==1) then
			!print *, '1', n, y(n,1), ya(1)
			if (dabs(y(n,2)) > 2.*ys+y1) then
				ya(1)=-ya(1)
				ya(3)=-ya(3)
			endif
		else !nslit==2
			!print *, '2', n, y(n,1), ya(1)
			if (dabs(y(n,2)) > ys+y1/2. .or. dabs(y(n,2)) < y1/2.) then
				ya(1)=-ya(1)
				ya(3)=-ya(3)
			endif
		endif
	endif

	!Check for screen passage of particle
	if (y(n,1)<=xd .and. ya(1)>xd) then
		index = (ya(2)+yd*.5)/yd*nscreen + 1
		!print *, 'index, y = ', index, ya(2)
		if (index.lt.1 .or. index.gt.nscreen) then
			print *, '*** Warning: screen index = ', index
		else
			screen(index) = screen(index)+1.
		endif
	endif
	END SUBROUTINE constraints
		
!*******************************************************************************
SUBROUTINE input(text, nmax, nbar, nslit, nout, nscreen, tmax, delta_t, vx0, vy0, &
		rbarx, rbary, xd, yd, xs, ys, y1, q, rx0, yrange)
!*******************************************************************************
	IMPLICIT NONE
	CHARACTER*(*) text
	INTEGER nmax, nbar, nslit, nout, nscreen
	REAL*8  tmax, delta_t,  vx0, vy0, xd, yd, xs, ys, y1, q, rx0, yrange
	REAL*8, DIMENSION (:), ALLOCATABLE :: rbarx, rbary
	INTEGER n
	
   	open(unit=10, file="input.txt", form="formatted")
	read (10,'(A)') text
	read (10,*) nmax
	read (10,*) tmax, delta_t
	read (10,*) nslit, nout, nscreen
	read (10,*) vx0, vy0
	read (10,*) xd, yd
	read (10,*) xs, ys, y1
	read (10,*) q, rx0
	read (10,*) yrange
	read (10,*) nbar

	ALLOCATE (rbarx(nbar), rbary(nbar))

	do n=1,nbar
		read (10,*) rbarx(n), rbary(n)
	enddo	
	
	close (10)

	if (nslit.eq.2 .and. mod(nmax,2).eq.1) nmax=nmax-1

	write (*,*)
	write (*,'(1X,(A))') text
	write (*,*) 'nmax          = ', nmax
	write (*,*) 'tmax, delta_t = ', tmax, delta_t
	write (*,*) 'nslit, nout,ns= ', nslit, nout, nscreen
	write (*,*) 'vx0, vy0      = ', vx0, vy0
	write (*,*) 'xd, yd        = ', xd, yd
	write (*,*) 'xs, ys, y1    = ', xs, ys, y1
	write (*,*) 'q, rx0        = ', q, rx0
	write (*,*) 'yrange        = ', yrange
	write (*,*) 'nbar          = ', nbar
	write (*,*)

END SUBROUTINE input

!*******************************************************************************
SUBROUTINE output(it, t, y, nmax, rbarx, rbary, nbar)
!*******************************************************************************
	IMPLICIT NONE
	INTEGER it, nmax, nbar
	REAL*8 t
	REAL*8, DIMENSION (1:nmax) :: rbarx, rbary
	REAL*8, DIMENSION (1:nmax,nfmax) :: y

	INTEGER n
	LOGICAL :: isopen = .false.
	save isopen

	! first call
	if (.not.isopen) then
		isopen = .true.
   		open(unit=20, file="output.txt", form="formatted")
		do n=1,nbar
			write (20,101) n, rbarx(n), rbary(n)
		enddo
		write (20,*)
		write (20,*)
	endif
	
	write (20,100) it, t, nmax
	do n=1,nmax
		! sort out masses with negative vx
		!if (y(n,3).ge. 0.) write (20,101) n, y(n,1:4)
		write (20,101) n, y(n,1:4)
	enddo
	write (20,*)
	write (20,*)

100	format (1h#,i5,e15.6,i5)
101	format (i5,4e15.6)

END SUBROUTINE output

!*******************************************************************************
SUBROUTINE outhist(screen, yd, nscreen)
! write screen histogram
!*******************************************************************************
	IMPLICIT NONE
	INTEGER nscreen
	REAL*8, DIMENSION (1:nscreen) :: screen
	REAL*8  :: yd, y

	open(unit=30, file="outscreen.txt", form="formatted")

	do n=1,nscreen
		y = yd*(n-1)/dfloat(nscreen-1) - yd*.5
		write (30,102) y, screen(n)
	enddo
	close(30)

102	format (4e15.6)

END SUBROUTINE outhist

END PROGRAM DOUBLE_SLIT_EXP
