!+
! Subroutine track1_custom (start_orb, ele, param, end_orb, err_flag, finished, track)
!
! Dummy routine for custom tracking. 
! This routine needs to be replaced for a custom calculation.
! If not replaced and this routine is called, this routine will generate an error message.
!
! This routine is potentially called twice by track1. 
! The entry_pt argument indicates from which point in track1 this routine is being called.
! 
! The radiation_included argument should be set to True if this routine takes into account radiation 
! damping and/or excitation when bmad_com%radiation_damping_on and/or bmad_com%radiation_fluctuations_on is True.
! If not, the track1 routine will use track1_radiation to include the radiation effects.
! Note: If this routine calles symp_lie_bmad, the symp_lie_bmad routine does take into account radiation effects.
! 
! General rule: Your code may NOT modify any argument that is not listed as an output agument below.
!
! Modules Needed:
!   use bmad
!
! Input:
!   start_orb  -- coord_struct: Starting position.
!   ele        -- ele_struct: Element.
!   param      -- lat_param_struct: Lattice parameters.
!   entry_pt   -- integer: Flag indicating from which point in track1 this routine is called.
!                   Possibilities are: entry_pt1$, entry_pt2$.
!
! Output:
!   end_orb     -- coord_struct: End position.
!   track       -- track_struct, optional: Structure holding the track information if the 
!                    tracking method does tracking step-by-step.
!   err_flag    -- logical: Set true if there is an error. False otherwise.
!   finished    -- logical: When set True, track1 will halt processing and return to its calling routine.
!   radiation_included
!               -- logical: Should be set True if radiation damping/excitation is included in the tracking.
!-

subroutine track1_custom (start_orb, ele, param, end_orb,  err_flag, finished, track)

use bmad_interface, except_dummy => track1_custom

implicit none

type (coord_struct) :: start_orb
type (coord_struct) :: end_orb
type (ele_struct) :: ele
type (lat_param_struct) :: param
type (track_struct), optional :: track

logical err_flag, finished
logical first/.true./
character(32) :: r_name = 'track1_custom'

real(rp) x_norm, y_norm, r, rel_p
real(rp) kx, ky
real(rp) coef
real(rp) zfact
real(rp) z,sigz
real(rp) bbi_length, alpha, d, sigx,sigy

integer nslice
!

if(index(ele%name,'ECLOUD') == 0)then
  print '(a,a)',' Custom tracking called for ',ele%name
  print '(a)',' Custom tracking available only for ECLOUD element'
  stop
endif
 
finished = .false.
end_orb = start_orb

if(index(ele%name,'ECLOUD_BBI')/= 0)then
! z dependence
z=start_orb%vec(5)
sigz = ele%value(sig_z$)
bbi_length = ele%value(custom_attribute1$)
d = ele%value(custom_attribute2$)
alpha = z/bbi_length/sigz
x_norm = start_orb%vec(1)/(ele%value(sig_x$) * exp(alpha-alpha**2)) ! the idea is that the width of the cloud shrinks along the bunch
y_norm = start_orb%vec(3)/(ele%value(sig_y$) * exp(alpha-alpha**2)) !  and positive z means the particle is ahead of the reference.
!But on further inspection these functions don't make much sense
sigx = ele%value(sig_x$) / (d*exp(-alpha-alpha**2)+1)
sigy = ele%value(sig_y$) / (d*exp(-alpha-alpha**2)+1)
x_norm = start_orb%vec(1)/sigx ! the idea is that the width of the cloud shrinks along the bunch
y_norm = start_orb%vec(3)/sigy !  and positive z means the particle is ahead of the reference
r= ele%value(sig_y$)/ele%value(sig_x$)

call bbi_kick (x_norm, y_norm, r, kx, ky)

 nslice = ele%value(n_slice$)
 rel_p = 1 + start_orb%vec(6)
 coef = ele%value(bbi_const$) / (nslice * rel_p)
! this is where the dependence on longitudinal dependence appears
! zfact = 1-start_orb%vec(5)/ele%value(sig_z$) or alternatively
 zfact = 1.
 end_orb%vec(2) = end_orb%vec(2) + kx * coef * zfact
 end_orb%vec(4) = end_orb%vec(4) + ky * coef * zfact
! if(first) write(19,'(6a12)')'x_norm','y_norm','zfact','delta px','delta py','end_orb%t'
! if((abs(x_norm) > 0.00001) .or. (abs(y_norm) > 0.00001))write(19,'(6es12.4)')x_norm,y_norm,zfact,kx*coef*zfact,ky*coef*zfact, end_orb%t
endif
first=.false.

if(index(ele%name,'ECLOUD_MAP')/= 0)call ecloud_ele(start_orb, ele, end_orb)

 end_orb%s = ele%s 

err_flag = .true.

! Remember to also set end_orb%t


end subroutine