program electron_tracking

  use bmad
  use muon_mod
  use muon_interface
  use runge_kutta_mod
  use materials_mod

  implicit none

  integer ios,status,system,lun
  integer seed
  integer track_state
  integer ix_end
  integer i
  integer datetime_values(8)

  real(rp) s_start, delta_s

  type (lat_struct),target:: lat
  type (em_field_struct) field
  type (electron_struct) electron
  type (electron_struct),allocatable::electrons(:)
  type (coord_struct), allocatable :: co_elec(:)
  type (branch_struct), pointer :: branch 
  type (coord_struct) start_orb, co_start, co_end

  character*120 lat_file_name
  character*16 date, time, zone

  logical err_flag

  ! Input namelist structure ("input.dat")
  namelist /input/ lat_file_name


  ! Read the lattice file.
  ! The lattice file also defines beam energy



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


   lat_file_name = 'bmad_osc_damping_wigglers.lat'

    print *, ' lat_file = ', lat_file_name

  ! Construct the lattice from the file
  
  bmad_com%auto_bookkeeper=.false.
  bmad_com%min_ds_adaptive_tracking = 0.00001


 ! call initializeMaterials()
  call bmad_parser (lat_file_name, lat)
 
!Initialize random number generator that is used to create phase space distribution 
   call ran_seed_put(seed)
   call ran_seed_get(seed)
 print *,' seed from ran_seed_get =',seed

 branch => lat%branch(0)

!define twiss parameter at the start of line
!  branch%ele(0)%a%beta  = abs(twiss%betax )
!  branch%ele(0)%b%beta  = abs(twiss%betay)
!  branch%ele(0)%a%alpha = twiss%alphax
!  branch%ele(0)%b%alpha = twiss%alphay
!  branch%ele(0)%x%eta   = twiss%etax 
!  branch%ele(0)%x%etap  = twiss%etapx 
!  branch%ele(0)%y%eta   = twiss%etay 
!  branch%ele(0)%y%etap  = twiss%etapy 


 

! start_orb%vec    = 0
  start_orb%t      = 0
  start_orb%s      = 0
  start_orb%vec(1) = 0
  start_orb%vec(2) = 0
  start_orb%vec(3) = 0
  start_orb%vec(4) = 0
  start_orb%vec(5) = 0
  start_orb%vec(6) = 0

 call init_coord(electron%coord,start_orb%vec,branch%ele(0),electron$)

  co_start = electron%coord

  print *, 'start_orb',co_start%p0c, co_start%species, co_start%state

  ix_end = branch%n_ele_track

! call track_from_s_to_s(lat, branch%ele(0)%s,branch%ele(ix_end)%s,co_start,co_end,all_orb=co_elec,ix_branch=0, track_state = track_state) !track electron to end of lattice
!  print *, 'end_orb',co_end%p0c, co_end%species, co_end%state

  
   lun = lunget()
   
   open(unit=lun,file= 'electron_traj.dat', status = 'replace')
   
   s_start=0.
   delta_s = 0.01

do i=1,branch%n_ele_track
  call track1(co_start,branch%ele(i),lat%param, co_end)
  print '(i10,1x,a,1x,7es12.4)',i,branch%ele(i)%name, branch%ele(i)%s,co_end%vec
   co_start=co_end
end do

  do while(s_start+delta_s <= branch%ele(branch%n_ele_track)%s)
   call track_from_s_to_s(lat, s_start,s_start+delta_s,co_start, co_end, all_orb = co_elec ,ix_branch=0, track_state = track_state)
   write(lun,'(2es12.4,6es12.4)')s_start+delta_s,co_end%t,co_end%vec
   s_start = s_start+delta_s
   co_start=co_end
  end do

!   do i= 0, ix_end
!   if (i>2 .and. co_elec(i)%t /= 0.) write (lun,'(8es12.4)')co_elec(i)%t,co_elec(i)%s,co_elec(i)%vec(:)
!   end do
   close(lun)

! if (track_state /= moving_forward$) then
!      electron%coord%state = co_elec(track_state)%state
      
! endif

deallocate(co_elec)

end