! Simulation of aeroplane problem
!
! N seats on an aeroplane and N people to fill them, each with a unique
! ticket number. Passengers proceed one at a time to their seats. The
! first person ignores their ticket number and chooses a seat at random.
! All subsequent people sit at their seat if it is available, otherwise
! they choose a seat at random. What is the probability that the last person
! can sit in their allocated seat.
!
! This simulation program is written in Fortran. To build it you'll need a
! fortran compiler. If you haven't got one get g95 from http://www.g95.org
! - it's free and is available for most OSs. If you want high performance then
! you'll find the Intel compiler (free for personal use on Linux) much better -
! just replace g95 with ifort below.
!
! To build the program type
!                           g95 -o aeroplane -O3 aeroplane.f90
!
! To run the program with 100 seats type at a command prompt
!                           aeroplane
!
! If you want to try a different number (e.g. 10) type
!                           aeroplane 10
!
! # seats    % availability  % std.dev
!     2         50.035        0.4830
!     3         49.931        0.4910
!     4         50.022        0.4979
!     5         50.026        0.5639
!     6         50.023        0.5229
!     7         50.030        0.5292
!     8         49.962        0.4597
!     9         50.021        0.5095
!    10         50.081        0.4409
!
! The % availability of all seats up to the given number is printed out along
! with the standard-deviation. If you just want the result for a particular
! number of seats use a negative number
!                           aeroplane -10
!
! # seats    % availability  % std.dev
!    10         50.020        0.4789
!
! Written by Simon Geard, 9/5/2007
! simon@whiteowl.co.uk
!
! Please feel free to adapt/use/change this program as you like.
!
!==============================================================================

module simulate

  public :: hasSeat, tidy

  private

  integer              :: nseats = 0
  logical, allocatable :: seatOccupied(:)

contains

  subroutine tidy
    if (allocated(seatOccupied)) deallocate(seatOccupied)
  end subroutine tidy

  logical function hasSeat(ns)
    ! Run the simulation once and return .true. if the final seat
    ! is available
    implicit none
    integer, intent(in) :: ns
    integer :: c
    integer :: i

    nseats = ns
    if (.not. allocated(seatOccupied)) then
       allocate(seatOccupied(ns))
    end if
    seatOccupied = spread(.false.,1,nseats)

    ! Random choice of 1st passenger
    c = selectSeat()
    if (c == 1) then
       hasSeat = .true.  ! Everyone else will be able to sit in their allocated seat
       return
    end if
    if (c == nseats) then
       hasSeat = .false. ! Last seat now occupied
       return
    end if
    seatOccupied(c) = .true.

    ! Other passengers
    do i=2,nseats-1
       if (seatOccupied(i)) then
          ! Seat occupied so choose another one at random
          c = selectSeat()
          if (c == nseats) then
             hasSeat = .false. ! Last seat now occupied
             return
          end if
          seatOccupied(c) = .true.
       else
          ! Seat available
          seatOccupied(i) = .true.
       end if
    end do
    hasSeat = .not. seatOccupied(nseats) ! Seat available if not occupied
  end function hasSeat

  integer function selectSeat()
    ! Choose a random unoccupied seat
    implicit none
    real :: h
    do
       call random_number(h)
       selectSeat = 1+nint((nseats-1)*h)
       if (.not. seatOccupied(selectSeat)) exit
    end do
  end function selectSeat

end module simulate



program aeroplane
  use simulate
  implicit none

  integer, parameter :: N = 10000 ! Number of trials in a simulation
  integer, parameter :: L = 100   ! Number simulations
  integer            :: c         ! Count of the number of times in a trial the seat
                                  ! was available
  integer            :: ns        ! Number of seats in the trial
  integer            :: nargs     ! Number of arguments on the command-line
  character(len=10)  :: arg       ! Holder for the first argument
  real*8             :: mean      ! Mean of results
  real*8             :: stdev     ! Standard-deviation of results
  real*8             :: results(L)! Results vector, 1 result for each simulation
  integer            :: i, j, k   ! Dummy loop counters

  ! Get the number of seats from the command line (default to 100)
  nargs = command_argument_count()
  if (nargs >= 1) then
     call get_command_argument(1,arg)
     read(arg,*) ns
  else
     ns = 100
  end if
  
  write(*,'(a)') '# seats    % availability  % std.dev'
  if (ns < 0) then
     k = -ns
     call runSimulation
  else
     do k=2,ns
        call runSimulation
     end do
  end if

contains
  subroutine runSimulation
     ! Do the L simulations
     do j=1,L
        c = 0
        !  Do the N trials
        do i=1,N
           if (hasSeat(k)) then
              c = c+1
           end if
        end do
        results(j) = dble(c)/N
     end do

     mean = sum(results)/L
     stdev = sqrt(dot_product(results,results)/L - mean**2)

     write(*,'(3X,i3,8X,f7.3,8X,f6.4)') k,mean*100,stdev*100

     call tidy

  end subroutine runSimulation
end program aeroplane
