!+                       
! Subroutine track_all_int_efield (lat, orbit, ix_branch, track_state,  sumBetaCrossE, err_flag, orbit0)
!
! Subroutine to track through the lat and integrate (beta cross E) dl.
!
! Note: If x_limit (or y_limit) for an element is zero then track_all will take
!       x_limit (or y_limit) as infinite.
!
! Modules Needed:
!   use bmad
!
! Input:
!   lat         -- lat_struct: Lat to track through.
!   orbit(0)    -- Coord_struct: Coordinates at beginning of lat.
!   ix_branch   -- Integer, optional: Index of branch to track. Default is 0 (main branch).
!
! Output:
!   orbit(0:)    -- Coord_struct: Orbit array.
!   track_state  -- Integer, optional: Set to moving_forward$ if everything is OK.
!                     Otherwise: set to index of element where particle was lost.
!   err_flag     -- Logical, optional: Set true if particle lost or error. False otherwise
!   orbit0(0:)   -- Coord_struct, optional: Orbit array for branch 0. Used to fill in the
!                     orbit at lord elemenets. 
!                     Only needed when orbit(:) is not the orbit for branch 0.
!  sumBetaCrossE(:)  -- Vector 
!-

subroutine track_all_int_efield (lat, orbit , BetaCrossE_start_time,  sumBetaCrossE, ix_branch, track_state, err_flag, orbit0, verbose)

use bmad_interface
use muon_mod
use parameters_bmad

implicit none

type (lat_struct), target :: lat
type (coord_struct), allocatable, target :: orbit(:)
type (coord_struct), optional, allocatable, target :: orbit0(:)
type (coord_struct), pointer :: orb(:)
type (ele_struct), pointer :: lord, slave
type (branch_struct), pointer :: branch
type (ele_struct), pointer :: ele
type (track_struct) track
type (sumBetaCrossE_struct), allocatable :: sumBetaCrossE(:)

real(rp) sum(3), px, py, pz, vec(6),E(3), dl, dt, dt_sum, sum_pitch(3), B(3), sum_path, sum_pathx,  sum_bf(3), sum_BetaCrossB(3), betahat(3)
real(rp), save :: b0
real(rp) BetaCrossB(3)
real(rp) BetaCrossE_start_time

integer n, i, nn, ix_br
integer, optional :: ix_branch, track_state
integer, save :: lun

logical, optional :: err_flag
logical, optional :: verbose
logical betaCrossE_verbose/.true./
logical err
logical :: debug = .false.
logical first/.true./
logical writetrackdetails/.false./


character(12), parameter :: r_name = 'track_all'

! init

if (present(err_flag)) err_flag = .true.
ix_br = integer_option (0, ix_branch)
branch => lat%branch(ix_br)
if (present(track_state)) track_state = moving_forward$

if (.not. allocated(orbit)) call reallocate_coord (orbit, branch%n_ele_max)
if (size(orbit) < branch%n_ele_max) call reallocate_coord (orbit, branch%n_ele_max)

if (orbit(0)%state == not_set$) call init_coord(orbit(0), orbit(0)%vec, branch%ele(0), downstream_end$)
if (bmad_com%auto_bookkeeper) call control_bookkeeper (lat)

orbit(0)%ix_ele   = 0
orbit(0)%location = downstream_end$

if (orbit(0)%species /= photon$) then
  call convert_pc_to (branch%ele(0)%value(p0c$) * (1 + orbit(0)%vec(6)), orbit(0)%species, beta = orbit(0)%beta)
  orbit(0)%p0c = branch%ele(0)%value(p0c$)
endif

! track through elements.

if (present(err_flag)) err_flag = .false.
if(present(verbose))betaCrossE_verbose = verbose

if(.not.  allocated(sumBetaCrossE))allocate(sumBetaCrossE(0:branch%n_ele_max))
forall(i=0:branch%n_ele_max)sumBetaCrossE(i)%vec(1:3)=0
forall(i=0:branch%n_ele_max)sumBetaCrossE(i)%vec_pitch(1:3)=0
forall(i=0:branch%n_ele_max)sumBetaCrossE(i)%dt=0
forall(i=0:branch%n_ele_max)sumBetaCrossE(i)%path=0
forall(i=0:branch%n_ele_max)sumBetaCrossE(i)%pathx=0
forall(i=0:branch%n_ele_max)sumBetaCrossE(i)%B_field(1:3) = 0
forall(i=0:branch%n_ele_max)sumBetaCrossE(i)%BetaCrossB(1:3) = 0

if(first)then
   if(betaCrossE_verbose)then
   lun = lunget()
   open(unit=lun, file = trim(directory)//'/'//'Trajectory_and_E.dat')
!      open(unit=lun, file = 'Trajectory_and_E.dat')
   write(lun,'(a10,1x,a16,1x,a10,1x,8a12,12x,a12,12x)')' Element ','Name','step','s','x','xp','y','yp','z','pz','t','E'
   endif
   first=.false.
   b0=0.
   do n=1,branch%n_ele_track
      if(trim(branch%ele(n)%name) == 'AFREE')then
            b0 = branch%ele(n)%value(b_field$)
            exit
      endif    
    end do
   endif


do n = 1, branch%n_ele_track

   ele => branch%ele(n)
  track%n_pt=-1
  call track1 (orbit(n-1), ele, branch%param, orbit(n), track=track, err_flag = err)

  if (err .or. .not. particle_is_moving_forward(orbit(n))) then
    if (present(track_state)) track_state = n
    call set_orbit_to_zero (orbit, n+1, branch%n_ele_track)
    if (orbit(n)%location == upstream_end$) orbit(n)%vec = 0 ! But do not reset orbit(n)%state
    if (present(err_flag)) err_flag = .true.
    exit
  endif

  if (debug) then
    call out_io (s_blank$, r_name, ele%name, '\6es16.6\ ', r_array = orbit(n)%vec)
  endif

  if(no_energy_change .and. index(ele%name,'QUAD')/= 0)orbit(n)%vec(6) = orbit(n-1)%vec(6) 
  sum = 0
  sum_pitch = 0
  sum_path = 0
  sum_pathx = 0
  dt_sum = 0
  E(1:3) = 0
  sum_bf(1:3) = 0
  sum_BetaCrossB = 0
  
  if(betacrossE_verbose) then
     write(lun,'(a1,/)')'#'
     write(lun, '(i10,1x,a16,1x,i10,1x,8es12.4,3es12.4)')(n,ele%name, i, track%pt(i)%orb%s, track%pt(i)%orb%vec, track%pt(i)%orb%t, track%pt(i)%field%E, i=0,track%n_pt)
     write(lun,'(/,a1)')'#'
  endif

  do i=1,track%n_pt
    vec(1:6) = track%pt(i)%orb%vec(1:6)
    px = vec(2)
    py = vec(4)
    pz = sqrt((vec(6)+1)**2 - vec(2)**2-vec(4)**2)
    betahat(1:3)=(/px,py,pz/)/sqrt(px**2+py**2+pz**2)
    E = track%pt(i)%field%E
    B = track%pt(i)%field%B
if(writetrackdetails)write(94,'(i10,12es12.4)')i,track%pt(i)%orb%s,vec(1:6),E(1:3),B(1:3)
!    dl = sqrt( (vec(1)-track%pt(i-1)%orb%vec(1))**2+(vec(3)-track%pt(i-1)%orb%vec(3))**2 + (track%pt(i)%orb%s-track%pt(i-1)%orb%s)**2)
    dt = track%pt(i)%orb%t - track%pt(i-1)%orb%t
   if(track%pt(i)%orb%t < BetaCrossE_start_time)cycle  !do not start to integrate until start of data taking
    sum = sum + (/py*E(3) - pz*E(2), pz*E(1)-px*E(3), px*E(2)-py*E(1)/)* track%pt(i)%orb%beta * dt
    sum_pitch = sum_pitch + ((px*B(1)+py*B(2)+B(3)*pz)* track%pt(i)%orb%beta)*(/px,py,pz/)*track%pt(i)%orb%beta * dt  !beta dot B X \vec beta
    BetaCrossB =  (/betahat(2)*B(3) - betahat(3)*B(2), betahat(3)*B(1)-betahat(1)*B(3), betahat(1)*B(2)-betahat(2)*B(1)/)   ! hat beta X B
    sum_BetaCrossB = sum_BetaCrossB +  (sqrt(dot_product(B,B))-sqrt(dot_product(BetaCrossB,BetaCrossB)))*dt
    sum_path = sum_path + dt * track%pt(i)%orb%beta * c_light
    sum_pathx = sum_pathx + vec(1)*dt
    dt_sum = dt_sum + dt
    sum_bf = sum_bf + (B - (/0.0_rp,b0, 0.0_rp/)) * dt
 !  print '(10e12.4,1x,es12.4,1x,3es12.4)', px,py,pz,B(1:3),((px*B(1)+py*B(2)+B(3)*pz)* track%pt(i)%orb%beta)*(/px,py,pz/)*track%pt(i)%orb%beta *dt, dt_sum, sum_pitch(1:3)
 enddo
  sumBetaCrossE(n)%vec = sumBetaCrossE(n-1)%vec+sum
  sumBetaCrossE(n)%dt = sumBetaCrossE(n-1)%dt + dt_sum
  sumBetaCrossE(n)%vec_pitch = sumBetaCrossE(n-1)%vec_pitch + sum_pitch
  sumBetaCrossE(n)%pathx = sumBetaCrossE(n-1)%pathx + sum_pathx
  sumBetaCrossE(n)%path = sumBetaCrossE(n-1)%path + sum_path
  sumBetaCrossE(n)%B_field = sumBetaCrossE(n-1)%B_field+ sum_bf
  sumBetaCrossE(n)%BetaCrossB = sumBetaCrossE(n-1)%BetaCrossB + sum_BetaCrossB
!  write(lun,'(i10,1x,a16,1x,8es12.4,3es12.4)')n,ele%name,ele%s,orbit(n)%vec(1:6),orbit(n)%t, E(1:3)
  if(betacrossE_verbose) write(lun,'(a)')'#'
  !  print '(6es12.4))',sumBetaCrossE(n)%vec_pitch(1:3), sum_pitch(1:3)  
 ! check for lost particles


enddo

! Fill in orbits for lord elements.

if (ix_br == 0) then
  orb => orbit
elseif (present(orbit0)) then
  orb => orbit0
else
  return   ! Nothing can be done
endif

do n = lat%n_ele_track+1, lat%n_ele_max
  lord => lat%ele(n) 
  if (lord%lord_status /= super_lord$) cycle
  slave => pointer_to_slave(lord, lord%n_slave)
  if (slave%ix_branch /= ix_br) cycle
  orb(lord%ix_ele) = orbit(slave%ix_ele)
enddo

end subroutine