!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Monte Carlo Approximation of Pi  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!  Consider one-quarter of a unit circle inscribed within   !!
!!  a square of side length 1.  Random (x,y) coordinate      !!
!!  pairs are generated using FORTRAN's built in random      !!
!!  number generator.  The number of (x,y) inside the        !!
!!  circle, N_c, approximates the area of the quarter        !!
!!  circle (1/4 pi r^2).  The number of (x,y) inside the     !!
!!  square, N_s, approximates the area of the square (r^2).  !!
!!  The ratio of the area of the quarter circle to the unit  !!
!!  square is then pi/4.  Finally, Pi_approx = 4*N_c/N_s.    !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!  compiled with gfortran -fno-range-check -fdefault-real-8 monte_pi.f95  !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  ABW 2015  !
!!!!!!!!!!!!!!
program thon
  implicit none
  
  integer,parameter :: lgint = selected_int_kind(16)  ! Use for large integers (> 2 000 000 000)
  integer(kind=lgint),parameter :: itmax = 100000000000  ! Maximum # of iterations
  real(kind=8),parameter :: pi = 3.1415926535897932d0  ! Actual value of pi to 16 digits
  real(kind=8) :: x,y,answ,exerr,var,stddev,PIsum,PIav
  integer(kind=lgint) :: i,n

  call init_random_seed()  !  Initialize random_seed
  
  n = 0
  PIsum = 0.0d0
  var = 0.0d0
  do i = 1,itmax  !  Iterate itmax times

     call random_number(x)
     call random_number(y)
     if (x*x + y*y <= 1.0d0) then
	n = n + 1  !  Counts the (x,y) inside the circle
     endif
       answ = 4.0d0*(real(n)/real(i))   !  pi_approx = 4*N_c/N_s = 4*n/i
       PIsum = PIsum + answ
       PIav=PIsum/real(i)
       var = var + (PIav - answ)**2  	!  Variance = 1/N SUM(i,1,N)[pi_average - pi_approx]^2
       stddev = sqrt(var/real(i))  	!  Standard Deviation = sqrt(Variance) ~ 1/sqrt(N)
       exerr = abs(answ - pi)		!  Exact Error = | Pi_approx - Pi_exact |
       if (mod(i,1000000) == 0) then	!  Occasionally print out data
	  open (unit=1,file="montepi_piactual.dat",action="write",access="append")
	  write(1,"(I14,F24.16)") i,answ
	  close (1)
	  open (unit=2,file="montepi_stddev.dat",action="write",access="append")
	  write(2,"(I14,F24.16)") i,stddev
	  close (2)
	  open (unit=3,file="montepi_exerr.dat",action="write",access="append")
	  write(3,"(I14,F24.16)") i,exerr
	  close (3)
	  open (unit=4,file="montepi_xy.dat",action="write",access="append")
	  write(4,"(F24.16,F24.16)") x,y
	  close (4)
       endif
  enddo
  
contains

  !  Initialize random seed:  Initilizes random_seed to system clock. (used with random_number)
  subroutine init_random_seed()
    integer :: i,n,clock
    integer,dimension(:),allocatable :: seed
      
    call random_seed(size = n)
    allocate(seed(n))      
    call system_clock(count = clock)      
    seed = clock + 37*(/ (i-1,i=1,n) /)
    call random_seed(put = seed)  
    deallocate(seed)
  end subroutine init_random_seed

end program thon
