      program sp620                 ! sp620.f
      implicit none
      integer i, j

      integer iy, Nr
      parameter( Nr=200000 )
      real*8 r
      data iy / 1234567 /

      integer Nx, Nx1
      real*8 Xmin, Xmax
      parameter( Nx=500, Nx1=Nx+1, Xmin=0.0, Xmax=1.0 )
      real*8    dx, x(Nx1), Box(Nx), BoxP(Nx)
      data Box / Nx*0.0 /

      integer Ny, Ny1
      real*8 Ymin, Ymax
      parameter( Ny=100, Ny1=Ny+1, Ymin=0.5, Ymax=1.5 )
      real*8 s, dy, y(Ny1), Boy(Ny), BoyP(Ny)
      data Boy / Ny*0.0 /
      real*8 sum, sum2, average, std_Dev
      data sum, sum2 / 0.0, 0.0 /
      
      real*8 pai, a, z, gauss
      parameter( pai=3.14159 )
      
      dx  = ( Xmax - Xmin ) / Nx
      x(1)= 0
      do j=1, Nx
        x(j+1) = x(j) + dx
      enddo
      do i=1, Nr
        call ransu( iy, r )             ! sp600.f
        do j=1, Nx
          if( x(j).lt.r .and. r.le.x(j+1) ) then
            Box(j) = Box(j) + 1         ! Box:Number of ransu
            go to 100
          endif
        enddo
  100   continue
      enddo
      do j=1, Nx
        BoxP(j) = Box(j) / dx / Nr      ! BoxP:Probability
      enddo
      open(2, file = 'sp620.dat' )
      do j=1, Nx
        write(2,'(i4,1p4e12.4)') j, x(j), x(j)+dx/2, Box(j), BoxP(j)
      enddo
      close(2)
 
      dy  = ( Ymax - ____ ) / __
      y(1)= Ymin
      do j=1, Ny
        y(j+1) = ____ + __
      enddo
      do i=1, Nx
        s    = BoxP(i)
        do j = 1, Ny
          if( ____.lt._ .and. _.le.______ ) then
            Boy(j) = ______ + _         ! Boy:Number of frequency
            go to 200
          endif
        enddo
  200   continue
        sum  = sum  + s
        sum2 = sum2 + s**2
      enddo
      average = sum/Nx                          ! Average
      std_Dev = sqrt( ____/__ - _______**2 )    ! Standard deviation
      write(*,'(a9,1pe15.7,a11,1pe15.7)')   'Average =', average, 
     &                                    '  Std_DeV =', std_Dev
      do j = 1, Ny
        BoyP(j) = Boy(j) / __ / __      ! BoyP:Probability
      enddo
      open(3, file = 'sp621.dat' )
      do j = 1, Ny
        write(3,'(i4,1p4e12.4)') j, y(j), y(j)+dy/2, Boy(j), BoyP(j)
      enddo
      close(3)
      
      open(4, file = 'sp622.dat' )
      a = 1.0/(sqrt(2*pai)*std_Dev)
      do j = 1, Ny
        z     = y(j) + dy/2
        gauss = a*exp(-(z-average)**2/(2*std_Dev**2))   ! Gauss function
        write(4,*) z, gauss
      enddo
      close(4)
      
      stop
      end
