subroutine  write_phase_space(lun, turn, beam)              
   use bmad
   use beam_mod

     implicit none

     type (beam_struct) beam
     integer lun
     integer turn

     integer units(0:10)/11*0/
     integer n/0/
     integer j

     write(lun,'(/,a1,1x,a7,i10,/)')'#','turn = ',turn
     if(all(units(0:n) /= lun))then
       n=n+1
       units(n)=lun
       write(lun,'(6a12)')'x','px','y','py','z','pz'
     endif

     do j=1,size(beam%bunch(1)%particle(:))
       if(beam%bunch(1)%particle(j)%state == alive$)write(lun,'(6es12.4)')beam%bunch(1)%particle(j)%vec(1:6)
     end do

      
   return
end