program tracking_master
  use bmad
  use runge_kutta_mod

  implicit none

  type (lat_struct), target:: lat
  type (coord_struct), allocatable :: co(:)
  type (em_field_struct) field

  integer pgopen, istat1, istat2
  integer i, j, k,l

  integer lun
  integer nargs, iargc, ios
  integer end, nturns
  integer n
  integer track_state
  integer nmuons
  integer ix

  real(rp) efield, bfield, scale_efield_bfield/1./
  real(rp) Q(3), deltae, deltax

  character*120 line, lat_file, lat_file_name/''/
 
  namelist/input/lat_file_name, nturns, nmuons,efield, deltae, bfield, scale_efield_bfield, deltax

  OPEN (UNIT=5, FILE='input.dat', STATUS='old', IOSTAT=ios)
    READ (5, NML=input, IOSTAT=ios)
    WRITE(6,NML=input)
  print *, 'ios=', ios
   rewind(unit=5)
!    READ (5, NML=input)
  CLOSE(5)

 bmad_com%spin_tracking_on = .true.

nargs = cesr_iargc()
  if (nargs == 1) then
    call cesr_getarg(1, lat_file)
    !print *, 'Using ', trim(lat_file)
  else if(lat_file_name /= '')then
     lat_file = lat_file_name
  else
    lat_file = 'bmad.'
    print '(a,$)','  Lattice file name?  (Default = bmad.)  '
    read(5,'(a)') line
    call string_trim(line, line, ix)
    lat_file = line
    if (ix==0) lat_file = 'bmad.'
  endif
    print *, ' lat_file = ', lat_file

  bmad_com%auto_bookkeeper=.false.
  bmad_com%min_ds_adaptive_tracking = 0.00001

  call bmad_parser (lat_file, lat)
  lun=lunget()
  call reallocate_coord (co, lat%n_ele_max)
  do i=1,lat%n_ele_track
  if(index( lat%ele(i)%name,'FREE')/= 0)lat%ele(i)%value(custom_attribute1$) = efield * scale_efield_bfield
  if(index( lat%ele(i)%name,'FREE')/= 0)lat%ele(i)%value(custom_attribute2$) = bfield * scale_efield_bfield
  end do

 co(0)%vec=0
 call track_all(lat,co)
 co(0)%vec=0
 co(0)%vec(6) = deltae
 co(0)%vec(1) = deltax
 co(0)%spin = (/(0.,0.),(1.,0.)/) 
 do i=1,nturns
    call track_all(lat,co)
    do j=1,lat%n_ele_track
     call compute_spins(co(j),Q)
     write(11, '(2i10,es12.4,6es12.4, 3es12.4, es12.4)')i,j,(i-1)+lat%ele(j)%s/lat%ele(lat%n_ele_track)%s,co(j)%vec, Q, co(j)%t
     if(i == 1)write(12,'(a,5es12.4)')lat%ele(j)%name,lat%ele(j)%floor%r,lat%ele(i)%floor%theta,lat%ele(i)%floor%phi
    end do
     co(0) = co(lat%n_ele_track)
 end do
 
end