      program ransu_dosu_bunpu_baratsuki_bunpu  ! ex060204.f
      implicit none
      integer i, j

      integer iy, Nr
      parameter( Nr=100000 )
      real*8 r
      common /comransu/iy
      save   /comransu/
c     data iy/1234567/

      integer Nx, Nx1
      real*8 Xmin, Xmax
      parameter( Nx=200, 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 )
      
      iy = 1234567
      
      dx  = ( Xmax - Xmin ) / Nx
      x(1)= 0
      do j=1, Nx
        x(j+1) = x(j) + dx
      enddo
      do i=1, Nr
        call ransu( r )
        do j=1, Nx
          if( x(j).lt.r .and. r.le.x(j+1) ) then
            Box(j) = Box(j) + 1
            go to 100
          endif
        enddo
  100   continue
      enddo
      do j=1, Nx
        BoxP(j) = Box(j) / dx / Nr 
      enddo

      open(2, file = 'ex06020402.dat' )
      do j=1, Nx
        write(2,'(1p4e12.4)') x(j), x(j)+dx/2, Box(j), BoxP(j)
c       write(*,'(1p4e12.4)') x(j), x(j)+dx/2, Box(j), BoxP(j)
      enddo
      close(2)
 
      dy  = ( Ymax - Ymin ) / Ny
      y(1)= Ymin
      do j=1, Ny
        y(j+1) = y(j) + dy
      enddo
      do i=1, Nx
        s    = BoxP(i)
        do j = 1, Ny
          if( y(j).lt.s .and. s.le.y(j+1) ) then
            Boy(j) = Boy(j) + 1 
            go to 200
          endif
        enddo
  200   continue
        sum  = sum  + s
        sum2 = sum2 + s**2
      enddo
      average = sum/Nx
      std_Dev = sqrt( sum2/Nx - average**2 )
      write(*,'(a9,1pe15.7,a11,1pe15.7)')   'Average =', average, 
     &                                     '  Std_DeV =', std_Dev
      do j = 1, Ny
        BoyP(j) = Boy(j) / dy / Nx
      enddo
      open(3, file = 'ex06020403.dat' )
      do j = 1, Ny
        write(3,'(1p4e12.4)') y(j), y(j)+dy/2, Boy(j), BoyP(j)
c       write(*,'(1p4e12.4)') y(j), y(j)+dy/2, Boy(j), BoyP(j)
      enddo
      close(3)
      
      open(4, file = 'ex06020404.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))
        write(4,'(1p2e12.4)') z, gauss
      enddo
      close(4)
      
      stop
      end
