      program sp720                 ! Gauss ransu
      implicit none

      integer iy
      real*8 r1, r2
      data iy/1234567/

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

      real*8 x1, x2, EX, STD
      parameter( EX=0d0, STD=10d0 )

      real*8 f, t, pai
      parameter( pai=3.14159d0 )
      f(t) = 1d0/(sqrt(_*___)*___) * exp( -(t-__)**2/2/___**2 )

      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/2

        call ransu( iy, r1 )             ! -> sp600.f
        call ransu( iy, r2 )             ! -> sp600.f
        x1 = __ + ___*( cos(_*___*__)*sqrt( -2*log(__) ) )
        x2 = __ + ___*( sin(_*___*__)*sqrt( -2*log(__) ) )
!       write(*,*) i, r1, r2, x1, x2

        do j=1,Nx
          if( x(j).lt.x1 .and. x1.le.x(j+1) ) then
            Box(j) = Box(j) + 1
            go to 100
          endif
        enddo
 100    continue
        sum = sum + x1
        sum2= sum2+ x1**2

        do j=1,Nx
          if( x(j).lt.x2 .and. x2.le.x(j+1) ) then
            Box(j) = Box(j) + 1
            go to 200
          endif
        enddo
 200    continue
        sum = sum + x2
        sum2= sum2+ x2**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='sp720.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
