c HW #6  Computing PI using Monte Carlo

        program calcpi

        implicit real (a-h, o-z)
        implicit integer (i-n)

        write(*,*) 'Choose n: '
        read(*,*) n 
       
        id = -5
 
        nhit = 0 

        do 50 i=1, n
           x = ran1(id)   
           y = ran1(id)  

c Calculate radius
c  NOTE: no need to waste computation time on square root
c        since we are comparing to 1.
c
           r = x**2.0 + y**2.0

           if (r.le. 1.0) then 
              nhit=nhit+1 
           endif

  50    continue

* Find the ratio of points in the circle to the total points used

        write(*,*) 
        write(*,*) n, nhit
        pi = 4.0 * (real(nhit)/real(n))
c
c Calcualte Actual Value of Pi and % Error
c
        pi_act = acos(-1.0) 

        rerror = abs((pi_act - pi)/ pi_act)
       
        write(*,100) n, pi, pi_act, rerror
 100    format(2x, i8, f15.8, f15.8, f15.8)
c
        stop
        end
         

C return a uniform random deviate between 0.0 and 1.0. 
C Call with negative idum to initiate.
C
      real FUNCTION ran1(idum)
      INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
      real ran1,AM,EPS,RNMX
      PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
     *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
      INTEGER j,k,iv(NTAB),iy
      SAVE iv,iy
      DATA iv /NTAB*0/, iy /0/
      if (idum.le.0.or.iy.eq.0) then
        idum=max(-idum,1)
        do 11 j=NTAB+8,1,-1
          k=idum/IQ
          idum=IA*(idum-k*IQ)-IR*k
          if (idum.lt.0) idum=idum+IM
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0) idum=idum+IM
      j=1+iy/NDIV
      iy=iv(j)
      iv(j)=idum
      ran1=min(AM*iy,RNMX)
      return
      END