!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Traffic Simulation using   !
!  Nagel-Schreckenberg Model  !
!  with optional Velocity     !
!  Dependent Randomization    !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!  Consider a circlular, one lane road "dimn" units  !!
!!  long, where a single unit is about a car length   !!
!!  in size.  The road is populated with "numstart"   !!
!!  vehicles each about a car length in size.  Cars   !!
!!  are not allowed to collide with or pass each      !!
!!  other.  Cars start with an initial velocity and   !!
!!  run through a four step simulation:  1)  Cars     !!
!!  accelerate by one unit, 2)  Collision is checked  !!
!!  for and avoided, 3)  Cars are randomly slowed     !!
!!  down (this is the Monte Carlo part), 4)  Each     !!
!!  car advances to it's new position (time step).    !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!  compiled with gfortran -fdefault-real-8 nasch_final.f95  !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  ABW 2015  !
!!!!!!!!!!!!!!
program nasch_new
  use sparse_display  !  Module used for sparse matrix plotting.  See J. Burkardt, FSU.
  implicit none
  
  integer,parameter :: dimn = 10000  !  Length of road
  integer,parameter :: tsmax = 40000 !  Number of time steps
  integer,parameter :: numstart = 2500  !  Number of vehicles
  integer,parameter :: vmax = 5  !  Maximum velocity
  real(kind=8),parameter :: pslowm = 0.01  !  Probability of a moving vehicle slowing (VDR)
  real(kind=8),parameter :: pslows = 0.5  !  Probability of a standing vehicle slowing (VDR)
  real(kind=8),dimension(numstart) :: randvel,pslow
  integer,dimension(numstart) :: vel,d
  integer,dimension(tsmax,dimn) :: posn
  integer,dimension(dimn) :: rhold
  integer :: i,j,dcount,icount,newpos,kcount,tstep,jhold
  real(kind=8) :: prob
  
  call init_random_seed  !  Initialize random_seed to system clock
  
  !  Initialize tsmax x dimn posn matrix
  !  0 represents empty road, 1 represents a single car
  do i = 1,tsmax
     do j = 1,dimn
	posn(i,j) = 0
     enddo
  enddo
  
  !  Set numstart vehicles randomly using the Knuth-Fisher-Yates shuffle algorithm
  !  Comment this part out if you want to set numstart vehicles uniformly
  do i = 1,dimn
     rhold(i) = i
  enddo
  call kfyshuffle(rhold,numstart,dimn)
  do i = 1,numstart
     posn(1,rhold(i)) = 1
  enddo
  
  !  Set numstart vehicles uniformly
  !  Comment this part out if you want to set numstart vehicles randomly
!  kcount = 0
!  do j = 1,dimn
!     if (mod(j,dimn/numstart) == 0) then
!	kcount = kcount + 1
!	if (kcount < numstart+1) then
!	   posn(1,j) = 1
!	endif
!     endif
!  enddo

  !  Set starting velocity with random_number or set to zero
  do i = 1,numstart
     call random_number(randvel(i))
     randvel(i) = randvel(i)*vmax
     vel(i) = Int(randvel(i))
!     vel(i) = 0
  enddo

  do tstep = 1,tsmax-1  !  Time step from 1 to tsmax - 1
     !  Dynamic randomization setup for Step THREE
     do i = 1,numstart
	if (vel(i) > 0) then
	   pslow(i) = pslowm
	elseif (vel(i) == 0) then
	   pslow(i) = pslows
	endif
     enddo

     !  Begin Step ONE, Acceleration:  Accelerate every vehicle with vel < vmax by one unit
     do i = 1,numstart
	if (vel(i) < vmax) then
	   vel(i) = vel(i) + 1
	endif
     enddo

     !  Begin Step TWO, Braking:  Do not allow cars to collide (exist in same piece of road)
     !  			  If a car's vel > d(distance to next vehicle) set vel = d
     icount = 0
     dcount = 0
     do j = 1,dimn
	if (j == 1 .and. posn(tstep,j) == 1) then
	   icount = 1
	   jhold = 1
	endif
	if (j > 1 .and. posn(tstep,j) == 1 .and. icount < numstart) then
	   if (icount > 0) then
	      d(icount) = dcount
	      dcount = 0
	      icount = icount + 1
	   elseif (icount < 1) then
	      dcount = 0
	      jhold = j
	      icount = icount + 1
	   endif
	endif
	if (posn(tstep,j) == 0 .and. icount < numstart) then
	   dcount = dcount + 1
	endif
	if (icount == numstart) then
	   d(numstart) = (dimn - j) + (jhold - 1)
	   icount = 0
	endif
     enddo
     do i = 1,numstart
	if (vel(i) > d(i)) then
	   vel(i) = d(i)
	endif
     enddo

     !  Begin Step THREE, Slowing:  With probability pslow, cars with vel > 0 slow by one unit
     !			 	    Turn off VDR by changing pslow(i) below to pslows or pslowm
     do i = 1,numstart
	call random_number(prob)
	if (prob <= pslow(i) .and. vel(i) > 0) then
	   vel(i) = vel(i) - 1
	endif
     enddo

     !  Begin Step FOUR, Motion:  Move numstart vehicles from their current position, j, by vel
     !				  Each Step FOUR represents the end of a single time step
     icount = 0
     newpos = 0
     do j = 1,dimn
	if (posn(tstep,j) == 1) then
	   newpos = j + vel(icount+1)
	   if (newpos > dimn) then
	      newpos = newpos - dimn
	   endif   
	   posn(tstep+1,newpos) = 1
	   icount = icount + 1
	endif
     enddo

  enddo
!  call mat(posn,tsmax,dimn)
  call spy_ge(tsmax,dimn,posn,'vdr_dimn10000_st2500_vm5_p001-05')  !  sparse_display subroutine

contains
  
  !  Knuth-Fisher-Yates Shuffle:  Takes in array, A, number of vehicles, n, length of road, m.
  !				  Returns A(1:n) filled with unique, random integers.
  subroutine kfyshuffle(A,n,m)
    integer,intent(inout) :: A(:)
    integer,intent(in) :: n,m
    integer :: i,randpos,temp,icount
    real :: r
    
    icount = 0
    do i = m, 2, -1
       call random_number(r)
       randpos = Int(r*i) + 1
       temp = A(randpos)
       A(randpos) = A(i)
       A(i) = temp
       icount = icount + 1
       if (icount > n) then
	  exit
       endif
    end do
    do i = 1,n
       temp = A(m - i + 1)
       A(i) = temp
    enddo
  end subroutine kfyshuffle
  
  !  Neatly print integer matrices:  Takes in n(column) x m(row) matrix, A.
  subroutine mat(A,n,m)
    integer :: A(:,:)
    integer :: n,m,i,j
      
    do i = 1,n
       do j = 1,m
	  write(*,"(I1,$)") A(i,j)
       enddo
       write(*,*)
    enddo
  end subroutine mat
  
  !  Initialize random seed:  Initilizes random_seed to system clock. (used with random_number)
  subroutine init_random_seed()
    integer :: i,n,clock
    integer,dimension(:),allocatable :: seed
      
    call random_seed(size = n)
    allocate(seed(n))      
    call system_clock(count = clock)      
    seed = clock + 37*(/ (i-1,i=1,n) /)
    call random_seed(put = seed)        
    deallocate(seed)
  end subroutine init_random_seed

end program nasch_new
