! **************************************************** !
! ******************* abs-abs.f90 ******************** !
! **************************************************** !
!
!    Leonardo Dagdug, Ivan Pompa-García, Jason Peña
!
! From the book:
! **************************************************** !
!      Diffusion Under Confinement:
!               A Journey Through Counterintuition    
! **************************************************** !

! **************************************************** !
! This file contains the source code for the simulation 
! of a 1D channel of length L.
! It has two absorbent points at x=0, and x=L.
! 
! The Brownian walkers start their path at x=x0, and 
! the simulation for each of them ends when an 
! absorbent point is reached. Then, the Mean First 
! Passage Time (MFPT or <tau>) is obtained.
!
! The theoretical value is used to calculate the 
! relative error of the simulation. D0 is the bulk 
! diffusion constant, usually set to 1.
!
! tau = x0 * (L-x0) / (2.0 * D0)
! **************************************************** !

! **************************************************** !
! To compile with GFortran, you must include the
! helpers.f90 module.
!         gfortran helpers.f90 abs-abs.f90
! Then you can run the program:
!         ./a.out
! **************************************************** !

program absabs

! Load the helpers module, which contains functions,
! constants, and more...
use helpers

! Mandatory declaration of data type variables
! and constants.
implicit none

! **************************************************** !
!              Declaration of constants
! **************************************************** !

!!! Simulation parameters
! Seed of the PRNG
integer, parameter :: RSEED = 0
! Number of random walkers (particles)
integer, parameter :: NRW = 2500
! Channel length
real(kind=dp), parameter :: L = 1.0_dp
! Initial position of the particles is fixed
real(kind=dp), parameter :: x0 = 1 / 3.0_dp
! Temporal step size
real(kind=dp), parameter :: DT = 1.0e-6_dp

!!! Physical parameters of the system
! Diffusion constant (bulk)
real(kind=dp), parameter :: D0 = 1.0_dp
! Thermal energy inverse (1/k_b T)
real(kind=dp), parameter :: BETA = 1.0_dp
! Standard deviation for Brownian motion
real(kind=dp), parameter :: SIGMA &
                              = dsqrt(2.0_dp * D0 * DT)
! Mean distribution for Brownian motion
real(kind=dp), parameter :: MU = 0.0_dp

! Theoretical value
real(kind=dp), parameter :: THEOV &
                          = x0 * (L-x0) / (2.0_dp * D0)


! **************************************************** !
!               Declaration of variables
! **************************************************** !

! General counter
integer :: i
! Current position of the particle
real(kind=dp) :: x
! Current particle passage time
real(kind=dp) :: tau
! Sum of passage times
real(kind=dp) :: ttau = 0.0_dp
! MFPT (<tau>)
real(kind=dp) :: mfpt


! **************************************************** !
!           Main simulation program
! **************************************************** !

! Set the seed of the PRNG to achieve repeatable results
call setseed(RSEED)

! The loop iterates over each (i) particle. 
do i=1, NRW
  !!! Initializing the needed variables
  ! The passage time is set to zero
  tau = 0.0_dp
  ! The particle current position is set to x0 (initial)
  x = x0

  !!! The random walk starts. Follow the particle until 
  !!! it is removed by an absorbent point (x=0,x=L).
  do
    ! Make a step in space
    x = x + nrand(MU, SIGMA)
    ! Make a step in time
    tau = tau + DT

    ! Check for the removal of the particle
    if(x <= 0.0 .or. x >= L) then
      ! If it was removed from the channel, stop the 
      ! simulation.
      exit
    end if
  end do ! End of i-particle's random walk.

  ! We need to add the passage time of the i-particle to
  ! the total time.
  ttau = ttau + tau

  ! You can print out the MFPT up to this point.
  print *, &
          '[' // nstr(i) // ' particles simulated | ' &
          // nstr(NRW-i) // ' to go]: ' &
          // TNL // TAB &
          // '<tau-sim> = ' // nstr(ttau / i) &
          // TNL // TAB &
          // '<tau-theo> = ' // nstr(THEOV) &
          // TNL // TAB &
          // 'Error = ' &
          // nstr( abs( (ttau/i) - THEOV ) / THEOV &
                   * 100, 1) &
          // '%'
end do ! End of the particle's loop.

! At this point, every particle's simulation is done.
! Obtain the <tau>
mfpt = ttau / NRW

! Now we print out the MFPT, and show the theoretical value.
print *
print *, '=== Final result of simulation ==='
print *, TAB // '<tau-sim> = ' // nstr(mfpt)
print *, TAB // '<tau-theo> = ' // nstr(THEOV)
! Calculates and prints the percentage error.
print *, TAB // 'Error = ' &
          // nstr( abs( mfpt - THEOV ) &
                  / THEOV * 100, 1) &
          // '%'
! As a reminder, print out the number of walkers.
print *, TAB // 'NRW = ' // nstr(NRW)
! The time step.
print *, TAB // 'dt = ' // nstr(DT, 6)
! And the PRNG seed.
print *, TAB // 'PRNG seed = ' // nstr(RSEED)

end program absabs