      program ex060301            ! ex060301.f  neutrino interaction
c                                 ! same as sp710.f
      implicit none

      integer iy
      real*8 r
      common /comransu/iy
      save   /comransu/
      data iy/1234567/

      integer i, j
      integer Nx, Nx1, Nr
      parameter( Nx=150, Nx1=Nx+1, Nr=100000 )
      real*8 xmin, xmax, dx, x(Nx1), Box(Nx), BoxP(Nx)
      parameter( xmin=0, xmax=15000d0 )
      data Box / Nx*0d0 /

      real*8 s, lambda
      parameter( lambda = 1000d0 )       ! 1000 km

      real*8 f, t
      f(t) = 1d0/lambda * exp(-t/lambda)

      real*8 sum, sum2, average, std_Dev
      data sum, sum2 / 0d0, 0d0 /

      dx  = ( xmax - xmin ) / Nx
      x(1)= xmin
      do j=1,Nx
        x(j+1) = x(j) + dx
      enddo
      do i=1,Nr
        call ransu( r )
        s = - lambda * log(r)
!       write(*,*) i, r, s
        do j=1,Nx
          if( x(j).lt.s .and. s.le.x(j+1) ) then
            Box(j) = Box(j) + 1
            go to 100
          endif
        enddo
 100    continue
        sum = sum + s
        sum2= sum2+ s**2
      enddo

      average = sum / Nr
      std_Dev = sqrt( sum2/Nr - average**2 )
      write(*,'(a9,1pe15.7,a11,1pe15.7)')   'Average =', average, 
     &                                    '  Std_DeV =', std_Dev

      do j=1,Nx
        BoxP(j) = Box(j) / dx / Nr
      enddo
      open(1, file='ex060301.dat')
      do j=1,Nx
        t = x(j) + dx/2
        write(1,'(1p5e12.4)') x(j), t, Box(j), BoxP(j), f(t)
      enddo
      close(1)

      stop
      end

      subroutine ransu(r)
      implicit real*8 (a-h,o-z)
      common /comransu/iy
      save   /comransu/
 1000 iy = iy * 48828125
      if( iy .lt. 0 ) iy = ( iy + 2147483647 ) + 1
      r = dfloat( iy ) * 4.6566128752457967d-10
      if( r .le. 0d0 ) go to 1000
      return
      end

