subroutine twiss_propagate_cm_general(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(:)
  
  real(rp) s, s_eff,c, step/0.01/
  real(rp) field_i
  integer lun
  integer nbranch
  integer ix_ele
  integer i, turns/10/
  logical eflag, first/.true./

  lun = lunget()
  open(unit = lun, file = '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,'(16a12)')'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)'
  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
  call twiss_propagate_all(lat)

  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)
   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.)

   write(lun,'(i10,16es12.4,1x,a16)')i,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
   s= s + step
  end do

return
end