subroutine create_distribution(nmuons, epsx,epsy, epsdistr,tdistr,tlength,tsigma,pzdistr,pz,pzsigma, muons)

use parameters_bmad
use muon_mod
use random_mod

implicit none

type (muon_struct), allocatable :: muons(:), muons_raw(:)
real(rp) epsx,epsy, sigma_x,sigma_y, epsdistr,tdistr,tlength,tsigma,pzdistr,pzsigma,pz
integer nmuons

real(rp) epsx0,epsy0, sigma_xp,sigma_yp
real(rp) dtheta,theta
real(rp) s_target/-1000./ ! distance to target (m)
integer i

logical first/.true./

if(first)then
  epsx0=epsx
  epsy0=epsy
  first = .false.
endif

!================================================
! CREATE 6D PHASE-SPACE DISTRIBUTION (AT TARGET)
!================================================

! Assume for simplicity that transverse positions and momenta are uncorrelated at target, i.e. alpha_{x,y}=0
sigma_xp = epsx0/sigma_x
sigma_yp = epsy0/sigma_y

IF (epsdistr=='gaus') THEN
  ! Gaussian distribution
  call ran_gauss(muons(:)%gas)
  muons(:)%coord%vec(1) = muons(:)%gas * sigma_x
  call ran_gauss(muons(:)%gas)
  muons(:)%coord%vec(2) = muons(:)%gas * sigma_xp
  call ran_gauss(muons(:)%gas)
  muons(:)%coord%vec(3) = muons(:)%gas * sigma_y
  call ran_gauss(muons(:)%gas)
  muons(:)%coord%vec(4) = muons(:)%gas * sigma_yp
ELSEIF (epsdistr=='flat') THEN
  ! Uniform distribution (default).  The width of the flat distribution 
  ! is chosen so as to recover the user's input emittances numerically
  call ran_uniform(muons(:)%flat)
  muons(:)%coord%vec(1) = ((muons(:)%flat)-0.5) * 2*1.73421323168796*sigma_x
  call ran_uniform(muons(:)%flat)
  muons(:)%coord%vec(2) = ((muons(:)%flat)-0.5) * 2*1.73421323168796*sigma_xp
  call ran_uniform(muons(:)%flat)
  muons(:)%coord%vec(3) = ((muons(:)%flat)-0.5) * 2*1.73421323168796*sigma_y
  call ran_uniform(muons(:)%flat)
  muons(:)%coord%vec(4) = ((muons(:)%flat)-0.5) * 2*1.73421323168796*sigma_yp
ELSEIF (epsdistr == 'round')THEN
 dtheta = twopi/size(muons(:))
 theta=0.
 i=0
 print *,' epsx0 = ',epsx0,' epsy0= ',epsy0
 do while (theta <= twopi .and. i < nmuons)
  theta = theta + dtheta
  i=i+1
  muons(i)%coord%vec(1) = epsx0 * cos(theta)
  muons(i)%coord%vec(2) = 0.
  muons(i)%coord%vec(3) = epsy0 * sin(theta)
  muons(i)%coord%vec(4) = 0.
 end do
! print *,'line 144'
! do i=1,size(muons(:))
!  print '(i10,1x,6es12.4)',i,muons(i)%coord%vec
! end do
ELSE ! delta-function
  sigma_x  = 0.
  sigma_y  = 0.
  sigma_xp = 0.
  sigma_yp = 0.
  muons(:)%coord%vec(1) = 0.
  muons(:)%coord%vec(2) = 0.
  muons(:)%coord%vec(3) = 0.
  muons(:)%coord%vec(4) = 0.
ENDIF


! Generate a momentum at the target.  Assume that position and momentum are uncorrelated.
IF (pzdistr=="gaus" .or. pzdistr=='Gaus' .or. pzdistr=='GAUS') THEN
  DO i=1, nmuons
    ! Generate a random momentum
    call ran_gauss(muons(i)%gas)
    muons(i)%coord%vec(6) = muons(i)%gas * pzsigma
    ! Make sure the random momentum lies within [-pz +pz], i.e. "cut the tails off the Gaussian"
    DO WHILE ( abs(muons(i)%coord%vec(6))>pz )
      call ran_gauss(muons(i)%gas)
      muons(i)%coord%vec(6) = muons(i)%gas * pzsigma
    ENDDO
  ENDDO
ELSE
  ! Uniform distribution (default)
  call ran_uniform(muons(:)%flat)
  muons(:)%coord%vec(6) = ((muons(:)%flat)-0.5) * pz
ENDIF

  print *,'s_target/betaMagic/c_light =',s_target/betaMagic/c_light


 return
end subroutine create_distribution