! **************************************************** !
! ************** abs-sphere-unif.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 2D (circular) / 3D (spherical) channel of
! radius R.
! The whole circumference is completely absorbent.
! 
! The Brownian walkers start their path being uniformly
! distributed inside the sphere / circle. 
! The simulation for each of them ends when the 
! absorbent shell 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, generally set to 1.
!
! tau = 1 / (15 * D0) * R**2
! **************************************************** !

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

program abssphereunif

! 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
! Circle / Sphere radius
real(kind=dp), parameter :: R = 1.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 &
                = (1.0_dp) / (15.0_dp * D0) * R**2


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

! General counter
integer :: i
! Initial position of the particle
real(kind=dp) :: x0, y0, z0
! Current position of the particle
real(kind=dp) :: x, y, z
! 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

  ! Initially, the particles should be uniformly 
  ! distributed inside the whole volume of the sphere.
  ! For a cube in cartesian coordinates, this would be
  ! straightforward, but for spherical coordinates
  ! the strategy needs to be modified, as otherwise the
  ! points will be inside a cube of edge R, and some
  ! points will remain outside the sphere.

  ! The calculation can be made using the sampling
  ! technique, which involves a cubic root and
  ! a transformation from spherical to cartesian
  ! coordinates.
  ! This can be computationally expensive and not
  ! readily understood, and some prefer to reject
  ! the points outside the sphere and get a new
  ! position until the point  remains inside. The latter
  ! is the approach we are taking here.
  do
    x0 = urand(-R, R)
    y0 = urand(-R, R)
    z0 = urand(-R, R)

    ! If the point lies inside the sphere of radius
    ! R, then we can continue.
    if( sqrt(x0**2 + y0**2 + z0**2) < R ) then
      exit
    end if
  end do

  ! The particle's current position is set to the
  ! initial position.
  x = x0
  y = y0
  z = z0

  !!! The random walk starts. Follow the particle until 
  !!! it is removed by the shell at radius R.
  do
    ! Make a step in space
    x = x + nrand(MU, SIGMA)
    y = y + nrand(MU, SIGMA)
    z = z + nrand(MU, SIGMA)
    ! Make a step in time
    tau = tau + DT

    ! This code also works for a 2D problem (CIRCLE),
    ! if you remove the «z» position so it is not taken
    ! into account for the calculation of the distance.
    ! If it's a 2D problem, just uncomment the next line
    ! z = 0.0_dp

    ! Check for removal of the particle.
    ! As we are inside a sphere (or a circle), the
    ! distance is calculated from the origin using
    ! the Euclidean distance.
    if( sqrt(x**2 + y**2 + z**2) >= R ) then
      ! If it was removed from the channel, stop its 
      ! 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 abssphereunif