!********************************************************************
! This program Shows how to generate Gaussian Random Variables using
! the Box - Muller Algorithm. To make the histogram of the output
! data, I also introduce an easy subroutine for 1D histograms.
! The program shows the basic usage of vectors in Fortran.
!********************************************************************
!********************************************************************
! Author: Giuseppe Forte
! First release: 23/09/2017
! e-mail: giuseppe.forte@giuseppeforte.me
! Website: https://www.giuseppeforte.me
! License: Creative Commons Attribution-ShareAlike 4.0
! International (CC BY-SA 4.0)
! (https://creativecommons.org/licenses/by-sa/4.0/)
!********************************************************************
!********************************************************************
! The following module contains the subroutine to make the Histogram
! of a 1D Random variable.
! Remember to use the current module in the main program by
! "uploading" it with the instruction
!
! use Histogram
!
! written immediately after the instruction program, i.e.
!
! program Prog_Name
!---> use Histogram
! implicit none
! [main program body]
! end program Prog_Name
!********************************************************************
module Histogram
contains
subroutine MakeHist(Vec,Nbin,Mids,Counts,Density)
implicit none
real*8, dimension(:),intent(in):: Vec ! Vec is a vector
! variable, its dimension is not specified (otherwise
! I would have used the instruction dimension(N), with N the
! predetermined dimension of the vector).
! That means that a certain memory space will be pre-allocated
! by fortran, thus the memory is left free and
! definitively allocated only elsewhere. This procedure is
! called dynamical memory allocation, in contrast with a
! fixed allocation, like the vector Counts(Nbin) below, for example,
! whose dimension is fixed and given by Nbin. In particular,
! Counts is a set of numbers, Counts(1), Counts(2), ...,
! Counts(Nbin). Such dimension (Nbin) cannot be changed during
! the execution of the program. On the other hand,
! Vec is an array whose dimension is not specified in the
! declaration section and can be allocated and REALLOCATED many
! times, as the program develops. Here we shall be concerned
! with a single allocation of the memory. An example with
! multiple allocations will be discussed in another note
integer, intent(in):: Nbin ! Number of bins
integer, intent(out):: Counts(Nbin) ! counting the number
! of events within a given bin
real*8, intent(out):: Mids(Nbin), Density(Nbin)
! ! Density = the renormalized counts to make the integral of
! the density equal to one
! Minds = mid points within the range of the random variable
real*8 xmax,xmin,dx,Pos,Norm,inflim,suplim
integer k
logical,allocatable:: Cond(:) ! this line is completely
! equivalent to the statement
! logical, dimension(:):: Cond (Here we are defining a logical
! vector, whose dimension is not declared, again, dynamical
! memory allocation)
! In order to allocate a vector dynamically allocated, like
! the vecotor "Vec" (or Cond) we should write an instruction like
! allocate(Vec(SpecifyDimension)), however, in this case the
! vector Vec is AUTOMATICALLY allocated. Indeed, Vec is an input
! argoment of the subroutine, meaning that, once we pass the
! argument to the subroutine, Vec will be allocated in
! agreement with the dimension of the passed vector
xmax = maxval(Vec)
xmin = minval(Vec)
dx = (xmax-xmin)/dble(Nbin)
Pos = xmin
k = 1
open(unit=15, file="Histo.dat",status="unknown")
do 108 while(Pos.le.xmax)
suplim = Pos + dx
inflim = Pos
Cond = Vec.lt.suplim ! here, the logical vector "Cond"
! is AUTOMATICALLY allocated by comparing it with the vector
! Vec. In order to make the equivalence syntactically correct
! Cond must be charcaterized by the same dimension of the expression
! on the right, thus it has to have the same dimension of Vec
Cond = Cond.and.(Vec.ge.inflim)
Counts(k) = count(Cond) !counts the number of true conditions
! in cond (i.e. the mumber of random variables in the current
! interval)
Mids(k) = inflim + 0.5d0*dx
Pos = Pos + dx
k = k + 1
108 end do
Norm = sum(Counts)*dx
Density = dble(Counts)/Norm ! sum(Density)/Norm = 1!!!!
do k=1, Nbin
write(15,*)real(Mids(k)),Counts(k),real(Density(k))
end do
close(15)
return
end subroutine MakeHist
end module Histogram
!********************************************************************
! MAIN PROGRAM
!********************************************************************
program BoxMULLER
use Histogram
implicit none
integer i !dummy index
integer Ns, Nbin ! Ns =Numbers of samples
parameter(Ns = 10000, Nbin = 100)
! Nbin = number of bins in the histogram
real*8 x(Ns)! x is a vector of dimension Ns, i.e.
! it has got Ns components. We can access the generic
! component "i" (i=1,2,3,...,Ns) in the vactor x
! simply by writing x(i)
real*8 Mids(Nbin), Density(Nbin)
integer Counts(Nbin)
! The abscissa of the histogram of the RV x runs from
! the minimum of x to the maximum of x, divided in
! Nbin intervals. The mid point of such intervals is stored
! in the vector Mids. The number of RVs falling in a given
! interval is stored in the vector Couts. Finally the
! normalized histogram is stored in the vector Density
real*8 u1,u2 !convenience variables
real*8 Mu, Sigma ! Mu = mean of the target Gaussian
! Sigma = Standard Deviation of the target
! Gaussian
parameter(Mu = 1.2d0, Sigma = 0.1d0)
real*8 TwoPi ! TwoPi = 2 * Greek Pi = 2 * 3.141593... =
! = 6.283185...
parameter(TwoPi = 8.0d0*atan(1.0d0))
! Initialization of the vector "x"
x = 0.0d0 ! this instruction puts ALL the components
! of the vector "x" equal to 0
!the above instruction is "equivalent" to the set
! of instructions
! do i=1, Ns
! x(i) = 0.0d0
! end do
! Actually, the instruction is not exactly the same. indeed, if we
! use the cycle do i=1...end do, Fortran will execute the instructions
! sequentially on the same processor. If we use the instruction x=0.0d0
! Fortran will run the instruction AUTOMATICALLY on ALL the AVAILABLE
! processors: in other words x=0.0d0 will be parallelized (not bad,
! considering that we are not making explicit use of MPI. This is a
! feauture which has been introduced only in the modern versions
! of Fortran.)
! Generation of the Gaussian random variables with mean Mu
! and Standard Deviation Sigma (Variance = Sigma^2)
do 100 i=1, Ns
call random_number(u1)
call random_number(u2)
! Implementation of the Box-Muller Algorithm
x(i) = Mu + Sigma*sqrt(-2.0d0*log(u1))*cos(TwoPi*u2)!
! Box-Muller
100 end do
! Construction of the histogram (frequencies) and
! associated probability density
call MakeHist(x,Nbin,Mids,Counts,Density)
! the subroutine MakeHist is contained in the module Histogram
! Such subroutine takes in input the vector "x", i.e. the
! "observed" random values (here we are simulating such values)
! and the number of bins (Nbin). The output are the vectors Mids,
! Counts and Density, TOGETHER with a data file "Histo.dat"
! containing three columns, specifically
!
! Mids(k), Counts(k), Density(k) (k=1,2,...,Nbin)
!
end program BoxMULLER