program etracking
  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
  real(rp) gamma_magic,beta_magic, p0c_magic, amu_measured/0.00116592091/, amu
  character*120 line, lat_file, lat_file_name/''/
 
  namelist/ering_input/lat_file_name, nturns, nmuons,efield, deltae, bfield, scale_efield_bfield, deltax

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

 bmad_com%spin_tracking_on = .true.

nargs = command_argument_count()
  if (nargs == 1) then
    call get_command_argument(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%spin_tracking_on = .true.
  bmad_com%auto_bookkeeper=.false.
  bmad_com%min_ds_adaptive_tracking = 0.0
  bmad_com%rel_tol_adaptive_tracking = 1d-10        ! Runge-Kutta tracking relative tolerance.
  bmad_com%abs_tol_adaptive_tracking = 1d-12       ! Runge-Kutta tracking absolute tolerance.

amu = amu_measured
gamma_magic = sqrt(1+1/amu)
beta_magic= sqrt(1-1/gamma_magic**2)
p0c_magic = gamma_magic*beta_magic*m_muon
print '(4(a15,es14.6),a22,es14.6)',' amu = ',amu,' gamma_magic = ', gamma_magic,' beta_magic = ',beta_magic,' p0c_magic = ',p0c_magic, ' B[T]( => R=7.112m) = ', p0c_magic/c_light/7.112

amu = anomalous_mag_moment_muon
gamma_magic = sqrt(1+1/amu)
beta_magic= sqrt(1-1/gamma_magic**2)
p0c_magic = gamma_magic*beta_magic*m_muon
print '(4(a15,es14.6),a22,es14.6)',' amu = ',amu,' gamma_magic = ', gamma_magic,' beta_magic = ',beta_magic,' p0c_magic = ',p0c_magic, ' B[T]( => R=7.112m) = ', p0c_magic/c_light/7.112

  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)then
    pointer_to_attribute (lat%ele(i), 'KICK_WIDTH', .true., err) %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

 call init_coord(co(0), co(0)%vec, lat%ele(0), downstream_end$)
 co(0)%vec=0
 co(0)%vec(6) = deltae
 co(0)%vec(1) = deltax
 co(0)%spin = [0, 0, -1]
 open(unit=11, file = 'phase_space_spin.dat')
 write(11, '(2a10,a14,6a14, 3a14, a14)')'turn','element','s (turns)','x','px','y','py','z','delta E', 'Sx','Sy','Sz', 'time'

 open(unit=12, file = 'floor_plan.dat')
 write(12,'(a14, a16,5a12)')'Element number','Element name','x','y','z','theta','phi'
do i=1,nturns
    call track_all(lat,co)
    do j=1,lat%n_ele_track
     Q = co(j)%spin
     write(11, '(2i10,es14.6,6es14.6, 3es14.6, es14.6)')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,'(5x,i5,4x,a16,5es12.4)')j,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
 close(unit=11)
 close(unit=12)
end