      program sp710                 ! sp710.f neutrino interaction
      implicit none

      integer iy
      real*8 r
      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( iy, r )             ! -> sp600.f
        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                        ! Averagee
      std_Dev = sqrt( sum2/Nr - average**2 )    ! Standard deviation
      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='sp710.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
