subroutine init_phase_space(ele,emitx,emity,emitz,co,i_dim)
use bmad
implicit none
!type(lat_struct) lat
type(coord_struct), allocatable:: co(:)
type(ele_struct) ele
real(rp) emitx, emity, emitz
real(rp) vec(1:6)
real(rp) phi
integer i,i_dim

!ele = lat%ele(0)


phi = 0.
  vec(5) = sqrt(emitz * ele%mode3%c%beta)*cos(phi)
  vec(6) = sqrt(emitz / ele%mode3%c%beta)*(ele%mode3%c%alpha *cos(phi+twopi/4.) - sin(phi))
  vec(1) = sqrt(emitx * ele%a%beta)*cos(phi)
  vec(2) = sqrt(emitx /ele%a%beta)*(ele%a%alpha*cos(phi+twopi/4.) -sin(phi+twopi/4.))
  vec(3) = sqrt(emity * ele%b%beta)*cos(phi)
  vec(4) = sqrt(emity / ele%b%beta)*(ele%b%alpha*cos(phi+twopi/4.) - sin(phi+twopi/4.))

do i=1, i_dim  !the set of basis vectors in x,x',y,y',z,pz space 
  co(i)%vec(i) = vec(i)
end do

phi = twopi/8.
  vec(5) = sqrt(emitz * ele%mode3%c%beta)*cos(phi)
  vec(6) = sqrt(emitz / ele%mode3%c%beta)*(ele%mode3%c%alpha *cos(phi+twopi/4.) - sin(phi))
  vec(1) = sqrt(emitx * ele%a%beta)*cos(phi)
  vec(2) = sqrt(emitx /ele%a%beta)*(ele%a%alpha*cos(phi+twopi/4.) -sin(phi+twopi/4.))
  vec(3) = sqrt(emity * ele%b%beta)*cos(phi)
  vec(4) = sqrt(emity / ele%b%beta)*(ele%b%alpha*cos(phi+twopi/4.) - sin(phi+twopi/4.))
  co(i_dim+1)%vec(1:6) = vec(1:6)

return
end subroutine