! **************************************************** !
! ****************** box-muller.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 performing
! a Box-Müller transformation over uniformly
! distributed numbers to make them follow a Gaussian
! distribution.
! **************************************************** !

! **************************************************** !
! To compile with GFortran, you must include the
! helpers.f90 module.
!         gfortran helpers.f90 box_muller.f90
! Then you can run the program:
!         ./a.out
!
! After running, a couple of files are generated:
!  - «random_uniform.dat», uniformly distributed numbers
!  - «random_normal.dat», normally distributed numbers
! Their histograms can be plotted with any software
! of your choice. To do this using gnuplot, execute:
!         gnuplot -p -e 'width=0.125 ; bin(x, width) = \
!          width*floor(x/width) ; plot \ 
!          "random_uniform.dat" u (bin($1,width)):(1.0)\
!          smooth freq with boxes notitle'
! and
!         gnuplot -p -e 'width=0.5 ; bin(x, width) = \
!          width*floor(x/width) ; plot \ 
!          "random_normal.dat" u (bin($1,width)):(1.0)\
!          smooth freq with boxes notitle'
!
! -p means that the plot must persist until the
!    user closes it explicitly.
! -e is used to specify the instructions
!    for using  gnuplot directly in the commandline.
! **************************************************** !

program box_muller

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

! Every data type must be explicitly stated.
implicit none

! How many numbers will be generated.
integer(kind=i32), parameter :: n = 10000

! A unit number for the file where the 
! uniformly distributed numbers will be saved.
! Fortran >= 2003 will
! assign this automatically in the open call.
integer(kind=i32) :: UNIF_UNIT
! A unit number for the file where the
! normally distributed numbers will be saved.
! Fortran >= 2003 will
! assign this automatically in the open call.
integer(kind=i32) :: NORMAL_UNIT

! Normally distributed numbers
real(kind=dp) :: x = 1.0_dp
real(kind=dp) :: y = 1.0_dp

! Variables to hold the uniformly distributed numbers
real(kind=dp) :: u1, u2

! A counter
integer :: i

! Open a file to write the uniformly distributed numbers
open(newunit=UNIF_UNIT, file='random_uniform.dat')
! Open a file to write the normally distributed numbers
open(newunit=NORMAL_UNIT, file='random_normal.dat')

! The cycle is executed only n/2 times because in every
! loop, two numbers of each kind are obtained. You must
! make sure that «n» is an even number  to get the 
! desired quantity of numbers.
do i=1, n/2
  ! Get the two uniformly distributed numbers
  call random_number(u1)
  call random_number(u2)

  ! Apply the Box-Müller method
  x = sqrt( -2.0_dp * log(u2) ) * cos(2.0_dp * PI * u1)
  y = sqrt( -2.0_dp * log(u2) ) * sin(2.0_dp * PI * u1)

  ! Write the uniformly distributed numbers to the files
  write(UNIF_UNIT, *) u1
  write(UNIF_UNIT, *) u2

  ! Write the normally distributed numbers to the files
  write(NORMAL_UNIT, *) x
  write(NORMAL_UNIT, *) y
end do

! Close the files
close(UNIF_UNIT)
close(NORMAL_UNIT)

print *, nstr(n) // ' uniformly distributed numbers, ' &
      // 'and ' // nstr(n) // ' normally distributed' &
      // ' numbers were generated!'

end program box_muller