subroutine twiss_propagate_cm(lat, nbranch)

  use bmad
  use parameters_bmad
  implicit none

  type (lat_struct) lat
  type (ele_struct) ele_at_s
  type (coord_struct) co, orb_at_s
  type (coord_struct), allocatable :: closed_orb(:)
  type(em_field_struct) field
  
  real(rp) s, s_eff,c, step/0.01/
  real(rp) field_i
  real(rp) s_ele
  integer lun
  integer nbranch
  integer ix_ele
  integer i, turns/10/
  logical eflag, first/.true./
  logical err_flag, local_ref_frame/.false./

  lun = lunget()
  open(unit = lun, file = trim(directory)//'/'//'twiss_file.dat')
  s=0
  write(lun,'(a,f12.4)')'Horizontal Tune = ', lat%a%tune/twopi
  write(lun,'(a,f12.4)')'Vertical Tune   = ', lat%b%tune/twopi
  write(lun,'(a,f12.4)')'Field index     = ', 1.-(lat%a%tune/twopi)**2
  write(lun,'(17a12)')'ix_ele','s','s_eff','beta x','alpha x','beta y', 'alpha y','eta_x','etap_x','etay','etap_y','beta x(ele)','beta y(ele)','x','y','x(ele)','y(ele)',&
      'bfield'
  call reallocate_coord (closed_orb, lat%branch(nbranch)%n_ele_track)
  call closed_orbit_calc(lat, closed_orb, i_dim=4, ix_branch=nbranch)
  c=lat%branch(nbranch)%ele(lat%branch(nbranch)%n_ele_track)%s
!do i=1, turns
  i=1
  do while(s <= lat%branch(nbranch)%ele(lat%branch(nbranch)%n_ele_track)%s)
   ix_ele= element_at_s(lat,s,choose_max=.true., ix_branch=nbranch, err_flag=eflag, s_eff=s_eff, position=co)
!
!   if(index(lat%branch(nbranch)%ele(ix_ele)%name,'FOCUS_ELE')/= 0)then
!     if(first) field_i = lat%branch(nbranch)%ele(ix_ele)%value(field_index$)
!     first=.false.
!    if(i == 0)then
!      lat%branch(nbranch)%ele(ix_ele)%value(field_index$) = field_i
!     else
!     lat%branch(nbranch)%ele(ix_ele)%value(field_index$) = 0.
!    endif
!   endif
   call twiss_and_track_at_s(lat,s,ele_at_s,orb=closed_orb,orb_at_s=orb_at_s, ix_branch=nbranch, use_last=.false.)
   s_ele = lat%branch(nbranch)%ele(ix_ele)%s_start
   call em_field_calc(lat%branch(nbranch)%ele(ix_ele), lat%param,s-s_ele,orb_at_s, local_ref_frame, field, err_flag=err_flag)
   write(lun,'(i10,16es12.4,1x,a16,1x,3es12.4)')ix_ele,s+(i-1)*c,s_eff, ele_at_s%a%beta, ele_at_s%a%alpha, ele_at_s%b%beta, ele_at_s%b%alpha, ele_at_s%x%eta,&
                                          ele_at_s%x%etap, ele_at_s%y%eta, ele_at_s%y%etap, lat%branch(nbranch)%ele(ix_ele)%a%beta,&
                                          lat%branch(nbranch)%ele(ix_ele)%b%beta,orb_at_s%vec(1),orb_at_s%vec(3), &
                                          closed_orb(ix_ele)%vec(1), closed_orb(ix_ele)%vec(3),lat%branch(nbranch)%ele(ix_ele)%name, field%b
   s= s + step
  end do
! lat%branch(nbranch)%ele(0) = lat%branch(nbranch)%ele(lat%branch(nbranch)%n_ele_track)
! lat%branch(nbranch)%ele(0)%s=0.
! call twiss_propagate_all(lat,ix_branch=nbranch)
! s=0
!end do ! turns

return
end