C----ROB SALGADO (cps615g)
C 1b.f 09/27/92 (09/28/92 PROBLEM 1.1b)
C
C USE MONTE-CARLO TO CALCULATE THE VOLUME OF A d-SPHERE
C
      program dspheremonte
      integer iseed,n,r,hit,dim,fac
      real radius,x,l,monterad(1000),sum,mean,sumdevsq,truevol
      real pi

      pi=3.1415926535

      open (unit=3, file='d-sphere')
      write(3,*) 'd-Sphere Monte-Carlo [ROB SALGADO]'
      print*, 'd-Sphere Monte-Carlo [ROB SALGADO]'
      
      iseed=0

      print*, 'Dimension dim=?'
      read*, dim
      print*, 'Radius radius=?'
      read*, radius
C-----------------------------COMPUTE TRUE D-VOLUME
      fac=1
      do 5 i=dim,1,-2
         fac=fac*i
 5    continue
      print*, fac
      if (2*int(dim/2).eq.dim) goto 8
C-------DIM IS ODD
      truevol=2**((dim+1)/2)*pi**((dim-1)/2)*radius**dim/fac
      goto 9
C-------DIM IS EVEN
 8    truevol=(2*pi)**(dim/2)*radius**dim/fac
 9    continue
C-------------------------------------------------|



      print*, 'Number of samples per run?'
      read*, n
 10   print*, 'Number of runs?'
      read*, r
      if (r.gt.1000) goto 10

C------------------------------------------------------r RUNS
      sum=0.0
      write(3,*) 'Trial     Radius'
      do 200 j=1,r,1
         hit=0
C---------------------------------n SAMPLES
      do 100 i=1,n,1
         l=0.0
         do 50 k=1,dim,1
            x=rand(iseed)
            l=l+x*x
 50      continue
      if (l.gt.1.0) goto 100
         hit=hit+1
 100  continue
C------------------------------------------|
      monterad(j)=(2*radius)**dim*hit/n
      write(3,*) j, monterad(j)
      sum=sum+monterad(j)
 200  continue
C-----------------------------------------------------------|

C------------------------------STATISTICS
      mean=sum/r
      write(3,*) '----'
      write(3,*) r,' Runs, each with',n,' samples'
      write(3,*) 'Volume enclosed by a',dim,' dimensional sphere'
      write(3,*) '           of radius',radius
      sumdevsq=0.0
      do 300 j=1,r,1
         sumdevsq=sumdevsq+(mean-monterad(j))*(mean-monterad(j))
 300  continue
      write(3,*) 'Mean volume=',mean, '  Std Dev=', sqrt(sumdevsq/n)
      write(3,*) 'True volume=',truevol
      write(3,*) '% Error=', 100.0*(mean-truevol)/truevol
      end








