program Fortran_Figure7_6_a

implicit none

save

 integer, parameter :: long = selected_real_kind(15,50)
 integer, parameter :: totalmol=1000  ! number of molecules
 integer, parameter :: ranpom=1000000
 integer, parameter :: noboxes=40
 integer, parameter :: nosavings=4

 real(kind=long), parameter :: DiffConst=0.0001_long ! diffusion constant
 real(kind=long), parameter :: tmlom=0.02_long ! 2*maximallengtofmove
 
 real(kind=long), dimension(totalmol) :: position
 real(kind=long), dimension(noboxes,nosavings) :: distribution
 real(kind=long), dimension(ranpom) :: RanNumUnif
   
 integer :: in,iru,jj,qq,ss,noiterations,noburnin,ll
 real(kind=long) :: proposedposition,pomprob
 real(kind=long) :: signalatpos,signalatproposedpos
 
 integer size, seed(2)
 data seed /6, 5/ 
      size = 2

 call random_seed(SIZE=size)
 call random_seed(PUT=seed(1:size))
 call random_number(RanNumUnif)
  
 iru=0
     
 noiterations=10000
 noburnin=1000
    
 do ll=1,nosavings
    
    noiterations=noiterations*10
    noburnin=noburnin*10
    
    write(*,*) ll,noiterations,noburnin
 
    do jj=1,noboxes
       distribution(jj,ll)=0.0_long
    end do
    
    do jj=1,totalmol
       iru=iru+1
       if (iru>ranpom) then
          call random_number(RanNumUnif)
          iru=1
       end if    
       position(jj)=RanNumUnif(iru)
    end do
 
 !==============
 !burn in
 
 do in=1,noburnin

    ! choose randomly a molecule
       
    iru=iru+1
    if (iru>ranpom) then
       call random_number(RanNumUnif)
       iru=1
    end if    
    qq=int(real(totalmol)*RanNumUnif(iru))+1
    
    ! move randomly the chosen molecule
    
    iru=iru+1
    if (iru>ranpom) then
       call random_number(RanNumUnif)
       iru=1
    end if
    proposedposition=position(qq)+(RanNumUnif(iru)-0.5_long)*tmlom
    
    ! acceptance criterion
    
    if (position(qq)<0.5) then
       signalatpos=0.001_long*position(qq)
    else
       signalatpos=0.001_long-0.001_long*position(qq)
    end if      
    if (proposedposition<0.5) then
       signalatproposedpos=0.001_long*proposedposition
    else
       signalatproposedpos=0.001_long-0.001_long*proposedposition
    end if   
    pomprob=exp((signalatproposedpos-signalatpos)/DiffConst)
    if ((proposedposition<0.0).OR.(proposedposition>1.0)) then
       pomprob=0.0_long
    end if
    
    iru=iru+1
    if (iru>ranpom) then
       call random_number(RanNumUnif)
       iru=1
    end if
    if (RanNumUnif(iru)<pomprob) then
       position(qq)=proposedposition 
    end if
    
 end do
 
 ! ===============
 ! sampling
           
 do in=1,noiterations

    ! choose randomly a molecule
       
    iru=iru+1
    if (iru>ranpom) then
       call random_number(RanNumUnif)
       iru=1
    end if    
    qq=int(real(totalmol)*RanNumUnif(iru))+1
    
       
    ! move randomly the chosen molecule
    
    iru=iru+1
    if (iru>ranpom) then
       call random_number(RanNumUnif)
       iru=1
    end if
    proposedposition=position(qq)+(RanNumUnif(iru)-0.5_long)*tmlom
    
    ! acceptance criterion
    
    if (position(qq)<0.5) then
       signalatpos=0.001_long*position(qq)
    else
       signalatpos=0.001_long-0.001_long*position(qq)
    end if      
    if (proposedposition<0.5) then
       signalatproposedpos=0.001_long*proposedposition
    else
       signalatproposedpos=0.001_long-0.001_long*proposedposition
    end if   
    pomprob=exp((signalatproposedpos-signalatpos)/DiffConst)
    if ((proposedposition<0.0).OR.(proposedposition>1.0)) then
       pomprob=0.0_long
    end if
    
    iru=iru+1
    if (iru>ranpom) then
       call random_number(RanNumUnif)
       iru=1
    end if
    if (RanNumUnif(iru)<pomprob) then
       position(qq)=proposedposition 
    end if
    
    ! add to distribution
    
    do jj=1,totalmol
       ss=nint(real(noboxes)*position(jj)+0.5000_long)
       if (ss>noboxes) then
	   ss=noboxes
       end if	   
       if (ss<1) then
	   ss=1
       end if	   
       distribution(ss,ll)=distribution(ss,ll)+1.0_long
    end do   	         
 end do
 
 do ss=1,noboxes
    distribution(ss,ll)=distribution(ss,ll)/real(noiterations)
 end do   	         
     
 end do
 
 open(unit=1, file='data_Figure7_6_a.dat', status='replace')

100 format(f6.2)

    do ll=1,nosavings   
       do ss=1,noboxes
          write(1,100) distribution(ss,ll)
       end do
    end do      	    

 close(unit=1)

end program Fortran_Figure7_6_a
