program  Fortran_Figure7_5c

implicit none

save

 integer, parameter :: long = selected_real_kind(15,50)
 integer, parameter :: ranpom=1000000
 integer, parameter :: ranpom2=500000
 integer, parameter :: nomesh=100
 real(kind=long), parameter :: DiffConst=1.0_long
 real(kind=long), parameter :: a1=0.5_long
 real(kind=long), parameter :: a2=0.0_long
 real(kind=long), parameter :: L=4.0_long
 real(kind=long), parameter :: radius=0.3_long
 
 real(kind=long), dimension(2,ranpom2) :: RanNum 
 real(kind=long), dimension(ranpom) :: RanNumNorm 
 real(kind=long), dimension(ranpom) :: RanNumUnif 
 
 real(kind=long), dimension(2) :: xpos,ypos,zpos
 real(kind=long), dimension(nomesh,nomesh) :: pdfxy
 
 
 integer :: irn,iru,ii,jj,zz,bouncross
 real(kind=long) :: alpha,time,dt,pom1,pom2,speedaver,dx,angle
 real(kind=long) :: fx,fy,fz,distancecube,r
 
 integer size, seed(2)
 data seed /10, 10/
      size = 2

 call random_seed(SIZE=size)
 call random_seed(PUT=seed(1:size))
 call random_number(RanNum)
 call random_number(RanNumUnif)
 
 do zz=1,ranpom2
    pom1=sqrt(2.0_long*log(1.0_long/RanNum(1,zz)))
    pom2=2.0_long*3.14159265358979_long*RanNum(2,zz)
    RanNumNorm(2*zz-1)=pom1*cos(pom2)
    RanNumNorm(2*zz)=pom1*sin(pom2)
 end do   
 
 irn=0
 iru=0
 
 bouncross=0
 
 alpha=sqrt(2.0_long*DiffConst)
 
 xpos(1)=0.0_long
 xpos(2)=1.0_long
 xpos(1)=0.0_long
 xpos(2)=1.0_long
 
 do ii=1,nomesh
    do jj=1,nomesh
       pdfxy(ii,jj)=0.0_long
    end do
 end do
        
 time=0.0_long
 dt=0.001_long
            
 do while (time<100000.0_long)

       time=time+dt
              
       ! interaction term between ions
       
       distancecube=(xpos(1)-xpos(2))*(xpos(1)-xpos(2))
       distancecube=distancecube+(ypos(1)-ypos(2))*(ypos(1)-ypos(2))
       distancecube=distancecube+(zpos(1)-zpos(2))*(zpos(1)-zpos(2))
       distancecube=sqrt(distancecube)*distancecube
       
       fx=a2*(xpos(1)-xpos(2))/distancecube
       fy=a2*(ypos(1)-ypos(2))/distancecube
       fz=a2*(zpos(1)-zpos(2))/distancecube
       	     
       ! evolve first molecule over one time step
       	  	      
       irn=irn+1
       if (irn>ranpom) then
           call random_number(RanNum)
	   do zz=1,ranpom2
              pom1=sqrt(2.0_long*log(1.0_long/RanNum(1,zz)))
              pom2=2.0_long*3.14159265358979_long*RanNum(2,zz)
	      RanNumNorm(2*zz-1)=pom1*cos(pom2)
	      RanNumNorm(2*zz)=pom1*sin(pom2)
	   end do   
           irn=1
       end if 
          	     
       xpos(1)=xpos(1)+(fx+a1)*dt+alpha*RanNumNorm(irn)*sqrt(dt)

       irn=irn+1
       if (irn>ranpom) then
           call random_number(RanNum)
	   do zz=1,ranpom2
              pom1=sqrt(2.0_long*log(1.0_long/RanNum(1,zz)))
              pom2=2.0_long*3.14159265358979_long*RanNum(2,zz)
	      RanNumNorm(2*zz-1)=pom1*cos(pom2)
	      RanNumNorm(2*zz)=pom1*sin(pom2)
	   end do   
           irn=1
       end if 
          	     
       ypos(1)=ypos(1)+fy*dt+alpha*RanNumNorm(irn)*sqrt(dt)

       irn=irn+1
       if (irn>ranpom) then
           call random_number(RanNum)
	   do zz=1,ranpom2
              pom1=sqrt(2.0_long*log(1.0_long/RanNum(1,zz)))
              pom2=2.0_long*3.14159265358979_long*RanNum(2,zz)
	      RanNumNorm(2*zz-1)=pom1*cos(pom2)
	      RanNumNorm(2*zz)=pom1*sin(pom2)
	   end do   
           irn=1
       end if 
          	     
       zpos(1)=zpos(1)+fz*dt+alpha*RanNumNorm(irn)*sqrt(dt)

       ! evolve first molecule over one time step
       	  	      
       irn=irn+1
       if (irn>ranpom) then
           call random_number(RanNum)
	   do zz=1,ranpom2
              pom1=sqrt(2.0_long*log(1.0_long/RanNum(1,zz)))
              pom2=2.0_long*3.14159265358979_long*RanNum(2,zz)
	      RanNumNorm(2*zz-1)=pom1*cos(pom2)
	      RanNumNorm(2*zz)=pom1*sin(pom2)
	   end do   
           irn=1
       end if 
          	     
       xpos(2)=xpos(2)+(-fx+a1)*dt+alpha*RanNumNorm(irn)*sqrt(dt)

       irn=irn+1
       if (irn>ranpom) then
           call random_number(RanNum)
	   do zz=1,ranpom2
              pom1=sqrt(2.0_long*log(1.0_long/RanNum(1,zz)))
              pom2=2.0_long*3.14159265358979_long*RanNum(2,zz)
	      RanNumNorm(2*zz-1)=pom1*cos(pom2)
	      RanNumNorm(2*zz)=pom1*sin(pom2)
	   end do   
           irn=1
       end if 
          	     
       ypos(2)=ypos(2)-fy*dt+alpha*RanNumNorm(irn)*sqrt(dt)

       irn=irn+1
       if (irn>ranpom) then
           call random_number(RanNum)
	   do zz=1,ranpom2
              pom1=sqrt(2.0_long*log(1.0_long/RanNum(1,zz)))
              pom2=2.0_long*3.14159265358979_long*RanNum(2,zz)
	      RanNumNorm(2*zz-1)=pom1*cos(pom2)
	      RanNumNorm(2*zz)=pom1*sin(pom2)
	   end do   
           irn=1
       end if 
          	     
       zpos(2)=zpos(2)-fz*dt+alpha*RanNumNorm(irn)*sqrt(dt)

       
       ! boundary conditions
              
       do jj=1,2 
              	  
          if (xpos(jj)>L) then
	      xpos(jj)=xpos(jj)-L
	      bouncross=bouncross+1
          end if
          if (xpos(jj)<0.0_long) then
	      xpos(jj)=xpos(jj)+L
	      bouncross=bouncross-1
          end if
	  r=sqrt(ypos(jj)*ypos(jj)+zpos(jj)*zpos(jj))
	  if (r>radius) then
	     if (r<2.0_long*radius) then 
	         ypos(jj)=ypos(jj)*(2.0_long*radius-r)/r
	         zpos(jj)=zpos(jj)*(2.0_long*radius-r)/r
	     else
	         ypos(jj)=0.0_long
		 zpos(jj)=0.0_long
             end if		 	 
	  end if   
       end do
              
       ! saving
       
       ii=nint(real(nomesh)*xpos(1)/L+0.5_long)
       jj=nint(real(nomesh)*xpos(2)/L+0.5_long)
       if (ii>nomesh) then
           ii=nomesh
       end if
       if (ii<1) then
           ii=1
       end if
       if (jj>nomesh) then
          jj=nomesh
       end if
       if (jj<1) then
          jj=1
       end if
       pdfxy(ii,jj)=pdfxy(ii,jj)+1.0_long    
 end do
 
 dx=L/real(nomesh)
  
 do ii=1,nomesh
    do jj=1,nomesh
       pdfxy(ii,jj)=pdfxy(ii,jj)*dt/time/dx/dx
    end do
 end do
 
 speedaver=real(bouncross)*L/time/2.0_long 
           	   
 open(unit=1, file='data_Figure7_5_c.dat', status='replace')

100 format(ES15.6E3)
 
    do ii=1,nomesh
       do jj=1,nomesh
          write(1,100) pdfxy(ii,jj)
       end do
    end do

 close(unit=1)

 write(*,*) bouncross,speedaver
 
 end program Fortran_Figure7_5c
