program Fortran_Figure6_5

implicit none

save

 integer, parameter :: long = selected_real_kind(15,50)
 integer, parameter :: ranmax=1000000
 integer, parameter :: ranmax2=500000
 integer, parameter :: maxmolA=400
 integer, parameter :: maxmolB=2000

 real(kind=long), parameter :: DA=1.0_long
 real(kind=long), parameter :: DB=1.0_long
 real(kind=long), parameter :: dt=0.0004_long 
 real(kind=long), parameter :: k1=0.03_long ! A+A->0
 real(kind=long), parameter :: k2=0.3_long ! A+B->0
 real(kind=long), parameter :: k3=0.0001_long ! A->0 
 real(kind=long), parameter :: k4=0.0001_long ! B->0
 real(kind=long), parameter :: k5=0.0000001_long ! 0->A
 real(kind=long), parameter :: k6=0.000001_long ! 0->B
 real(kind=long), parameter :: L=1000.0_long
 real(kind=long), parameter :: h=25.0_long 
 real(kind=long), parameter :: bindingradius1=0.04_long
 real(kind=long), parameter :: bindingradius2=0.04_long
 real(kind=long), parameter :: probrem1 = 0.095_long
 real(kind=long), parameter :: probrem2 = 0.595_long
 real(kind=long), dimension(2,ranmax2) :: RanNum 
 real(kind=long), dimension(ranmax) :: RanNumNorm 
 real(kind=long), dimension(ranmax) :: RanNumUnif 
 logical, dimension(maxmolA) :: Amolecule 
 logical, dimension(maxmolB) :: Bmolecule 
 integer, dimension(40) :: Ahistogram,Bhistogram 
 real(kind=long), dimension(maxmolA,3) :: Apositions ! positions of A molecules
 real(kind=long), dimension(maxmolB,3) :: Bpositions ! positions of B molecules 
 real(kind=long) :: time,pom1,pom2,dist
 integer :: ii,jj,zz,irn,iru,tt,oldtt,qq,maxusedA,maxusedB,ip,pomIn
 
 integer size, seed(2)
 data seed /4, 4/ 
      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,ranmax2
    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

 oldtt=0
 time=0.0_long

 maxusedA=0
 maxusedB=0
 
 do ii=1,maxmolA
    Amolecule(ii)=.FALSE.
    Apositions(ii,1)=0.0_long
    Apositions(ii,2)=0.0_long
    Apositions(ii,3)=0.0_long
 end do
 do ii=1,maxmolB
    Bmolecule(ii)=.FALSE.
    Bpositions(ii,1)=0.0_long
    Bpositions(ii,2)=0.0_long
    Bpositions(ii,3)=0.0_long
 end do
        
 qq=0
    
 do while (time<12000.0_long)
     
    time=time+dt
    
    ! diffusion of A molecules
    
    do ii=1,maxusedA
    
       if (Amolecule(ii)) then
          
          irn=irn+1
          if (irn>ranmax) then
             call random_number(RanNum)
             do zz=1,ranmax2
                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 
    
          Apositions(ii,1)=Apositions(ii,1)+sqrt(2.0_long*DA*dt)*RanNumNorm(irn)
	  if (Apositions(ii,1)<0.0_long) then
	     Apositions(ii,1)=-Apositions(ii,1)
          end if
	  if (Apositions(ii,1)>L) then
	     Apositions(ii,1)=2*L-Apositions(ii,1)
          end if
	  
          irn=irn+1
          if (irn>ranmax) then
             call random_number(RanNum)
             do zz=1,ranmax2
                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 
    
          Apositions(ii,2)=Apositions(ii,2)+sqrt(2.0_long*DA*dt)*RanNumNorm(irn)
	  if (Apositions(ii,2)<0.0_long) then
	     Apositions(ii,2)=-Apositions(ii,2)
          end if
	  if (Apositions(ii,2)>h) then
	     Apositions(ii,2)=2*h-Apositions(ii,2)
          end if

          irn=irn+1
          if (irn>ranmax) then
             call random_number(RanNum)
             do zz=1,ranmax2
                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     
          
	  Apositions(ii,3)=Apositions(ii,3)+sqrt(2.0_long*DA*dt)*RanNumNorm(irn)
	  if (Apositions(ii,3)<0.0_long) then
	     Apositions(ii,3)=-Apositions(ii,3)
          end if
	  if (Apositions(ii,3)>h) then
	     Apositions(ii,3)=2*h-Apositions(ii,3)
          end if
	  
       end if    
    end do
    
    ! diffusion of B molecules
    
    do ii=1,maxusedB
    
       if (Bmolecule(ii)) then
          
          irn=irn+1
          if (irn>ranmax) then
             call random_number(RanNum)
             do zz=1,ranmax2
                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 
    
          Bpositions(ii,1)=Bpositions(ii,1)+sqrt(2.0_long*DB*dt)*RanNumNorm(irn)
	  if (Bpositions(ii,1)<0.0_long) then
	     Bpositions(ii,1)=-Bpositions(ii,1)
          end if
	  if (Bpositions(ii,1)>L) then
	     Bpositions(ii,1)=2*L-Bpositions(ii,1)
          end if
	  
          irn=irn+1
          if (irn>ranmax) then
             call random_number(RanNum)
             do zz=1,ranmax2
                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 
    
          Bpositions(ii,2)=Bpositions(ii,2)+sqrt(2.0_long*DB*dt)*RanNumNorm(irn)
	  if (Bpositions(ii,2)<0.0_long) then
	     Bpositions(ii,2)=-Bpositions(ii,2)
          end if
	  if (Bpositions(ii,2)>h) then
	     Bpositions(ii,2)=2*h-Bpositions(ii,2)
          end if

          irn=irn+1
          if (irn>ranmax) then
             call random_number(RanNum)
             do zz=1,ranmax2
                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     
          
	  Bpositions(ii,3)=Bpositions(ii,3)+sqrt(2.0_long*DB*dt)*RanNumNorm(irn)
	  if (Bpositions(ii,3)<0.0_long) then
	     Bpositions(ii,3)=-Bpositions(ii,3)
          end if
	  if (Bpositions(ii,3)>h) then
	     Bpositions(ii,3)=2*h-Bpositions(ii,3)
          end if
	  
       end if    
    end do

    ! creation of A
        
    iru=iru+1
    if (iru>ranmax) then
       call random_number(RanNumUnif)
       iru=1
    end if
    
    if (RanNumUnif(iru)<k5*h*h*L*dt) then
       ii=1
       do while ((Amolecule(ii)).AND.(ii<maxmolA))
          ii=ii+1
       end do
       if (Amolecule(ii)) then
          write(*,*) 'Error: number of A molecules greater than maxmolA (increase maxmolA)'
	  STOP
       else
          Amolecule(ii)=.TRUE.
          iru=iru+1
          if (iru>ranmax) then
             call random_number(RanNumUnif)
             iru=1
          end if
	  Apositions(ii,1)=RanNumUnif(iru)*L
          iru=iru+1
          if (iru>ranmax) then
             call random_number(RanNumUnif)
             iru=1
          end if
	  Apositions(ii,2)=RanNumUnif(iru)*h
          iru=iru+1
          if (iru>ranmax) then
             call random_number(RanNumUnif)
             iru=1
          end if
	  Apositions(ii,3)=RanNumUnif(iru)*h
	  if (ii>maxusedA) then
	     maxusedA=ii
	  end if   
       end if  	  	  
    end if
    
    ! creation of B
        
    iru=iru+1
    if (iru>ranmax) then
       call random_number(RanNumUnif)
       iru=1
    end if
    
    if (RanNumUnif(iru)<k6*h*h*L*dt*0.4_long) then
       ii=1
       do while ((Bmolecule(ii)).AND.(ii<maxmolB))
          ii=ii+1
       end do
       if (Bmolecule(ii)) then
          write(*,*) 'Error: number of B molecules greater than maxmolB (increase maxmolB)'
	  STOP
       else
          Bmolecule(ii)=.TRUE.
          iru=iru+1
          if (iru>ranmax) then
             call random_number(RanNumUnif)
             iru=1
          end if
	  Bpositions(ii,1)=0.6_long*L+RanNumUnif(iru)*L*0.4_long
          iru=iru+1
          if (iru>ranmax) then
             call random_number(RanNumUnif)
             iru=1
          end if
	  Bpositions(ii,2)=RanNumUnif(iru)*h
          iru=iru+1
          if (iru>ranmax) then
             call random_number(RanNumUnif)
             iru=1
          end if
	  Bpositions(ii,3)=RanNumUnif(iru)*h
	  if (ii>maxusedB) then
	     maxusedB=ii
	  end if   
       end if  	  	  
    end if
    
    ! bimolecular reaction A + A -> 0
            
    pomIn=maxusedA	
	
    do ii=1,pomIn-1
    
       if (Amolecule(ii)) then
       
          jj=ii
	  
          do while ((jj<maxusedA).AND.(Amolecule(ii)))

             jj=jj+1
	     	     
             if ((Amolecule(jj)).AND.(abs(Apositions(ii,1)-Apositions(jj,1))<bindingradius1)) then
       
	        dist=(Apositions(ii,1)-Apositions(jj,1))*(Apositions(ii,1)-Apositions(jj,1)) 
                dist=dist+(Apositions(ii,2)-Apositions(jj,2))*(Apositions(ii,2)-Apositions(jj,2)) 
                dist=dist+(Apositions(ii,3)-Apositions(jj,3))*(Apositions(ii,3)-Apositions(jj,3)) 
                dist=sqrt(dist)
	        if (dist<bindingradius1) then
                   iru=iru+1
                   if (iru>ranmax) then
                      call random_number(RanNumUnif)
                      iru=1
                   end if
	           if (RanNumUnif(iru)<probrem1) then
  	              Amolecule(ii)=.FALSE.
  	              Amolecule(jj)=.FALSE.
	              if (jj.EQ.maxusedA) then
	                 maxusedA=maxusedA-1
                      end if		
	              if (ii.EQ.maxusedA) then
	                 maxusedA=maxusedA-1
                      end if		
		   end if
                end if
             end if
          end do
       end if	  
    end do
              
    ! bimolecular reaction A + B -> 0
            
    pomIn=maxusedA	
	
    do ii=1,pomIn
    
       if (Amolecule(ii)) then
       
          jj=0
	  
          do while ((jj<maxusedB).AND.(Amolecule(ii)))

             jj=jj+1
	     	     
             if ((Bmolecule(jj)).AND.(abs(Apositions(ii,1)-Bpositions(jj,1))<bindingradius2)) then
       
	        dist=(Apositions(ii,1)-Bpositions(jj,1))*(Apositions(ii,1)-Bpositions(jj,1)) 
                dist=dist+(Apositions(ii,2)-Bpositions(jj,2))*(Apositions(ii,2)-Bpositions(jj,2)) 
                dist=dist+(Apositions(ii,3)-Bpositions(jj,3))*(Apositions(ii,3)-Bpositions(jj,3)) 
                dist=sqrt(dist)
	        if (dist<bindingradius2) then
                   iru=iru+1
                   if (iru>ranmax) then
                      call random_number(RanNumUnif)
                      iru=1
                   end if
	           if (RanNumUnif(iru)<probrem2) then
  	              Amolecule(ii)=.FALSE.
  	              Bmolecule(jj)=.FALSE.
	              if (jj.EQ.maxusedB) then
	                 maxusedB=maxusedB-1
                      end if		
	              if (ii.EQ.maxusedA) then
	                 maxusedA=maxusedA-1
                      end if		
		   end if
                end if
             end if
          end do
       end if	  
    end do
    
              
    ! degradation A -> 0
            
    pomIn=maxusedA		
    do ii=1,pomIn    
       if (Amolecule(ii)) then       
           iru=iru+1
           if (iru>ranmax) then
               call random_number(RanNumUnif)
               iru=1
           end if
	   if (RanNumUnif(iru)<k3*dt) then
  	      Amolecule(ii)=.FALSE.
	      if (ii.EQ.maxusedA) then
	         maxusedA=maxusedA-1
              end if		
           end if
       end if   	  
    end do
               
   ! degradation B -> 0
            
    pomIn=maxusedB		
    do ii=1,pomIn    
       if (Bmolecule(ii)) then       
           iru=iru+1
           if (iru>ranmax) then
               call random_number(RanNumUnif)
               iru=1
           end if
	   if (RanNumUnif(iru)<k4*dt) then
  	      Bmolecule(ii)=.FALSE.
	      if (ii.EQ.maxusedB) then
	         maxusedB=maxusedB-1
              end if		
           end if
       end if   	  
    end do
    
    tt=nint(time*0.1_long)
    if (tt>oldtt) then
       write(*,*) 'Time computed:',nint(time),'seconds (final time 12000 sec)'
       oldtt=tt
    end if
end do  
 
 do ii=1,40
    Ahistogram(ii)=0
    Bhistogram(ii)=0
 end do   
 
 do ii=1,maxusedA    
    if (Amolecule(ii)) then
       jj=nint(Apositions(ii,1)/25.0_long+0.5_long)
       if (jj.EQ.41) then
          jj=40
       end if
       Ahistogram(jj)=Ahistogram(jj)+1
    end if
 end do       
     
 do ii=1,maxusedB    
    if (Bmolecule(ii)) then
       jj=nint(Bpositions(ii,1)/25.0_long+0.5_long)
       if (jj.EQ.41) then
          jj=40
       end if
       Bhistogram(jj)=Bhistogram(jj)+1
    end if
 end do       
     
           
 open(unit=1, file='data_Figure6_5.dat', status='replace')

100 format(i5,x,i5)
 
 do ii=1,40
    write(1,100) Ahistogram(ii),Bhistogram(ii)
 end do

 close(unit=1)
 
end program Fortran_Figure6_5
