!+
! 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.
!
! Also see:
!   track1_preprocess
!   track1_postprocess
!
! 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, a custom version of track1_preprocess should be 
! constructed to set its radiation_included argument to 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.
!
! Note: If ele%spin_tracking_method = tracking$, then it is expected that this routine will also handle
! spin tracking. The alternative is when ele%spin_tracking_method = custom$ in which case track1_spin_custom will
! be called after this routine. If doing spin tracking here, bmad_com%spin_tracking_on should be checked
! to see if spin tracking is actually wanted.
!
! General rule: Your code may NOT modify any argument that is not listed as an output agument below.
!
! Input:
!   start_orb  -- coord_struct: Starting position.
!   ele        -- ele_struct: Element.
!   param      -- lat_param_struct: Lattice parameters.
!
! Output:
!   end_orb     -- coord_struct: End position.
!   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.
!   track       -- track_struct, optional: Structure holding the track information if the 
!                    tracking method does tracking step-by-step.
!-

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

use bmad, 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

real(rp) m11, m12, m21, m22
real(rp) L0, K_und, B, vbar, gamma_rel
real(rp) lambda, k1

logical err_flag, finished

character(*), parameter :: r_name = 'track1_custom'

L0 = ele%value(n_pole$)*ele%value(l_pole$) !undulator length
lambda = ele%value(l_pole$)*2 !undulator wavelength
B = ele%value(b_max$) !peak field (T)
gamma_rel = sqrt((start_orb%p0c*(1.+start_orb%vec(6)))**2+m_electron**2)/m_electron !relativistic gamma

K_und = e_charge*B*lambda/(twopi*m_electron*e_charge/c_light) !undualtor paramater
vbar = c_light*(1.-1./2./gamma_rel/gamma_rel*(1.+K_und*K_und)) !average longitudinal velocity (m/s)

k1 = 1./2.*((e_charge*B)/(vbar*gamma_rel*m_electron*e_charge/c_light**2))**2 !effective focusing strength

!transfer matrix
m11 = cos(sqrt(k1)*L0)
m12 = 1./sqrt(k1)*sin(sqrt(k1)*L0)
m21 = -sqrt(k1)*sin(sqrt(k1)*L0)
m22 = cos(sqrt(k1)*L0)

end_orb = start_orb

end_orb%vec(1) = m11*start_orb%vec(1) + m12*start_orb%vec(2)
end_orb%vec(2) = m21*start_orb%vec(1) + m22*start_orb%vec(2)
end_orb%vec(3) = m11*start_orb%vec(3) + m12*start_orb%vec(4)
end_orb%vec(4) = m21*start_orb%vec(3) + m22*start_orb%vec(4)


end_orb%s = start_orb%s + L0
end_orb%t = start_orb%t + L0/vbar

err_flag = .false.
finished = .false.


end subroutine