c ***************************************************************************
c Copyright(C) M.Kotelyanskii, November, 1997
c This code is just an example, it can be freely distributed and modified
c as soon as this notice remains in the file.
c It has been tested on the IBM Risc6000 under AIX. Other systems may require
c changes, particularly the functions calls for the random numbers generation.
c Even though it works properly to the best of the author's knowledge,
c he does not assume any responsibility for any damage caused by it
c ***************************************************************************
c ***************************************************************************
c *****  This is example of the Monte Carlo method for 1-dimensional integral
c *****  of the function f(x) on the interval [a,b]                             
c ***************************************************************************

       implicit real*8(a-h,o-z)

c **** Specify parameters:

c **** Interval 
       parameter(a=0.d0, b=1.d0)
c **** Number of trials
       parameter(Ncycle=10000, Nprint=Ncycle/10)
c **** Random number generator's seed
       parameter(iseed=1234567)

c **** set counters
       fav  = 0.0d0
       f2av = 0.0d0

c *** set seed for the random numbers generator (notice, that the name srand
c *** may be different on different architectures )
       call srand(iseed)

c *** Start MC cycle (assuming rand() returns random real*1 from 0 to 1)

       do icycle = 1, Ncycle
           fx  = f( a+rand()*(b-a) )
           fav = fav + fx
           f2av= f2av + fx*fx

c **** print out the running average
           if(mod(icycle,Nprint).eq.1) then
             frunav = fav / icycle
             f2runav = f2av / icycle
             f2runav = sqrt((f2runav-frunav*frunav)/icycle)
             write(6,2) icycle, frunav*(b-a), f2runav*(b-a)
2            format(1x,i8, "cycles Integral = ",e12.5," +/- ", e12.5)
           endif 
       enddo

c *** Get the average and  the standard deviation

       fav  = fav  / Ncycle
       f2av = f2av / Ncycle
       f2av = sqrt((f2av-fav*fav)/Ncycle)

c *** Write out the result

       write(6,1) Ncycle, fav*(b-a), f2av*(b-a)
1      format(1x,i8, "cycles Integral = ",e12.5," +/- ", e12.5)

c *** That's all folks

        end

c **** Function we integrating 
       real*8 function f(x)
       real*8 x  
         f=x*x
       return
       end	
