C
C  program python2
C
C  illustration of simple Monte Carlo integration
C  computes integral from 0 to 1 of 1/(1+x^2) = pi/4
C
C  modified from python1
C  includes weighting
C
C  CWJ SDSU 2/20/2005
C
C  function called: RAN3  in ranlib.f
C

      program python2
      implicit none

      real pi

C............... variables for random

      real ran3
      real x
      integer seed
      real y

C.................... variables for integration

      real sum
      real sum2
      real f         ! the function
      real w,wsum,wsum2
      real vsum,vsum2

      integer nsamp
      integer i
C.................... START.......................

      pi = 2.*asin(1.)

      print*,' PYTHON 2 '
      print*,' Illustration of Monte Carlo quadrature '

      print*,' '
      print*,' Enter seed for random numbers '
      read*,seed
      if(seed.gt.0)seed=-seed
  101 continue
      print*,' Enter # of points '
      read*,nsamp
      if(nsamp.le.0)stop
      
C.............. start integration.........

      sum = 0.0
      do i = 1,nsamp
         x = ran3(seed)
         f = 1./(1+x*x)
         sum = sum+f
         sum2 = sum2+f*f 
      enddo

      sum = sum/nsamp
      sum2 = sum2/nsamp - sum**2
C................NEW weighting...............

      wsum = 0.0
      wsum2 = 0.0

      do i = 1,nsamp
        y = ran3(seed)
        x = 2.-sqrt(4.-3.*y)
        w = 1./3.*(4.-2*x)
        f = 1./(1+x*x)
        wsum = wsum+f/w
        wsum2 = wsum2+(f/w)**2
      enddo
      wsum = wsum/nsamp
      wsum2 = wsum2/nsamp - wsum**2
C..................... von Neummann

      call vonneumann(nsamp,seed,vsum,vsum2) 
C............... output.....................
      write(6,*)'      MC  est error   exact err'
      write(6,303)sum,sqrt(sum2/nsamp),abs(sum-pi/4.)
      write(6,303)wsum,sqrt(wsum2/nsamp),abs(wsum-pi/4.)
      write(6,303)vsum,sqrt(vsum2/nsamp),abs(vsum-pi/4.)
  303 format(3f10.5)    
      goto 101
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


      subroutine vonneumann(nsamp,seed,sum,sum2)

      implicit none


      real w,wprime

      real x

      real f 
      real ran3
      real sum,sum2
      integer nsamp
      integer i
      integer seed
      real dice
      logical accept
      integer nct

      sum = 0.0
      sum2 = 0.0

      nct = 0
      do i = 1,nsamp
        accept=.false.
        do while(.not.accept)
          x = ran3(seed)
          w = 6./5.*(1.-0.5*x*x)
          wprime = 6./5. 
          dice = ran3(seed)
          nct = nct+1
          if(w/wprime .gt. dice)accept = .true.

        enddo 

        f= 1./(1.+x*x)
        sum = sum+f/w
        sum2 = sum2+(f/w)**2

      enddo
      sum = sum/nsamp
      sum2 = sum2/nsamp-sum*sum
      write(6,*)' % acceptance von Neumann = ',
     & 100*float(nsamp)/float(nct)
      write(6,*)' '
      return
      end