program initialdist
       implicit real*8 (a-h,o-z) 
       integer, parameter :: Npt = 1000
       real*8, dimension(Npt) :: x,px,y,py,z,dg

       twopi = 4*asin(1.0)
       emass = 0.511e6
       sigdg = 2e3/emass
       sig = 0.2e-3
       gamma0 = 100e6/emass
       emitun= 0.5e-6/gamma0
       sigxp = emitun/sig
       do i = 1, Npt
         call random_number(r0)
         z(i) = 3e-3*r0-1.5e-3
         call gaussian(rg)
         dg(i) = sigdg*rg
         call gaussian(rg)
         x(i) = sig*rg
         call gaussian(rg)
         px(i) = sig*rg
         call gaussian(rg)
         y(i) = sig*rg
         call gaussian(rg)
         py(i) = sig*rg
         write(10,110)x(i),px(i),y(i),py(i),z(i),dg(i)
       enddo
110    format(6(1x,e15.7))

       end program


       subroutine gaussian(rg)
       real*8 :: r1,r2,rr,rg,phi
         twopi = 4*asin(1.0)
         call random_number(r1)
         call random_number(r2)
         phi = twopi*r1
         rr = sqrt(-2*log(r2))
         rg = rr*cos(phi)
       end subroutine