program Fortran_Figure7_3

implicit none

save

 integer, parameter :: long = selected_real_kind(15,50)
 integer, parameter :: totalbac=1000  ! number of bacteria
 integer, parameter :: ranpom=1000000
 integer, parameter :: noboxes=100
 real(kind=long), parameter :: speed=0.02_long ! 20 micron/second
 real(kind=long), parameter :: lambda=1.0_long ! (unbiased) turning rate
 
 real(kind=long), dimension(ranpom) :: RanNumUnif 
 
 integer, dimension(noboxes) :: histsave 
 real(kind=long), dimension(totalbac) :: positions,velocities,internal
 
 integer :: ii,jj,iru,ss,tt
 real(kind=long) :: dt,time,signal,lambday,maxsignal,dersig
 
 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
 dersig=4.0_long
 
 do ii=1,totalbac
    positions(ii)=1.0_long 
    iru=iru+1
    if (iru>ranpom) then
        call random_number(RanNumUnif)
        iru=1
    end if    
    if (RanNumUnif(iru)<0.5_long) then
       velocities(ii)=speed
    else
       velocities(ii)=-speed
    end if
    internal(ii)=dersig*positions(ii)      
 end do
    
 do ss=1,noboxes
       histsave(ss)=0
 end do
      
 time=0.0_long
 dt=0.01_long
 maxsignal=3.0_long
     
 do tt=1,360000
    
       time=time+dt
       
       ! evolve every cell over one time step
       
       do ii=1,totalbac
          iru=iru+1
          if (iru>ranpom) then
             call random_number(RanNumUnif)
             iru=1
          end if 
	  signal=dersig*(maxsignal-abs(positions(ii)-maxsignal))
	  lambday=lambda+internal(ii)-signal
	  if (RanNumUnif(iru)<lambday*dt) then
	     velocities(ii)=-velocities(ii)
          end if
	  positions(ii)=positions(ii)+velocities(ii)*dt
	  internal(ii)=internal(ii)*exp(-dt)+(1.0_long-exp(-dt))*signal
	  if (positions(ii)<maxsignal) then
  	     internal(ii)=internal(ii)+velocities(ii)*dersig*(dt+exp(-dt)-1.0_long)
  	  else
	     internal(ii)=internal(ii)-velocities(ii)*dersig*(dt+exp(-dt)-1.0_long)
          end if
       end do
              	      	                    	      
    if (int(time*0.01_long)-int(time*0.01_long-dt*0.01_long)>0.5_long) then
        write(*,*) 'Time computed:',nint(time),'seconds (final time 3600 sec)'
    end if
        
 end do
    
 do jj=1,totalbac
    ss=nint(25.0_long*positions(jj))
    if ((ss<noboxes+1).AND.(ss>0)) then
	histsave(ss)=histsave(ss)+1
    end if	
 end do   	         
     
  
 open(unit=1, file='data_Figure7_3.dat', status='replace')

100 format(i6,1X,f6.2)
 
       do ss=1,noboxes
          write(1,100) histsave(ss),real(ss)/25.0_long
       end do   	    

 close(unit=1)

end program Fortran_Figure7_3
