!+
! Subroutine track1_custom (start_orb, ele, param, end_orb, track, err_flag)
!
! Dummy routine for custom tracking. 
! If called, this routine will generate an error message and quit.
! This routine needs to be replaced for a custom calculation.
!
! Note: This routine is not to be confused with track1_custom2.
! See the Bmad manual for more details.
!
! 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.
!
! 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.
!-

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

use bmad_interface, except_dummy => track1_custom
use ecloudmod

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

! === CUSTOM PARAMETER ===
integer, parameter :: calculation_index$ = custom_attribute1$
integer, parameter :: math_model$ = 1
integer, parameter :: qdens_math$ = 2
integer, parameter :: qdens_snapshot$ = 3
integer, parameter :: file_driven$ = 4
! === ================ ===

! Need to store field values
type (em_field_struct) :: field

real(rp) :: xkick, ykick, zkick

character(32) :: r_name = 'track1_custom'
	
	type (DiscreteEcloud), save :: ecloud_data
	logical, save :: init_needed = .true.
	logical :: custom_err = .false.
	integer lun1, lun2, lun3

	!===========================================!
	! First create the data structure if needed !
	!===========================================!
	if (init_needed) then
		print *, "[*] ECLOUD TRACKING: Performing first time setup"
		print *, "   [*] INPUT VALUES:"
		print *, "  (1): ", ele%value(custom_attribute1$)
		print *, "  (2): ", ele%value(custom_attribute2$)
		print *, "  (3): ", ele%value(custom_attribute3$)
		print *, "  (4): ", ele%value(custom_attribute4$)
		print *, "  (5): ", ele%value(custom_attribute5$)
		select case ( int(ele%value(calculation_index$)) )
			case (math_model$)
				print *, " ","[*] MATH MODEL"
				! ...
			case (qdens_math$)
				! Calculate by reading qdens from file then using trends
				print *, " ","[*] QDENS MATH"
				lun1 = lunget()
				open(unit=lun1,file='snapshot/qdens.data')
				call calculateDiscreteFromQdensFile(lun1,ecloud_data)
			case (qdens_snapshot$)
				! Take snapshot of bunch and scale with qdens from file
				print *, " ","[*] QDENS SNAPSHOT"
				lun1 = lunget()
				open(unit=lun1,file='snapshot/snapshot.data')
				lun2 = lunget()
				open(unit=lun2,file='snapshot/qdens.data')
				lun3 = lunget()
				open(unit=lun3,file='snapshot/bunchhead.data')
				call calculateDiscreteFromSnapshotsFiles(lun1,lun2,lun3,ecloud_data)
			case (file_driven$)
				! Calculate by reading EM from file
				print *, " ","[*] FULL TRAIN FILE"
				lun1 = lunget()
				open(unit=lun1,file='snapshot/snapshot.data')
				lun2 = lunget()
				open(unit=lun2,file='snapshot/qdens.data')
				lun3 = lunget()
				open(unit=lun3,file='snapshot/bunchhead.data')
				!call calculateDiscreteFromSnapshotsFiles_Full(lun1,ecloud_data)
		end select
		init_needed = .false.
		!call performTests(ele, ecloud_data)
	end if

	!==================================================================!
	! Now calculate the Ex and Ey given the position of the particles. !
	! This is either done using a Discrete structure, or the formulas. !
	!==================================================================!
	select case ( int(ele%value(calculation_index$)) )
			case (math_model$)
                call calculateEM(ele,start_orb%t,start_orb,field)
			case (qdens_math$)
				call calculateEMFromDiscrete_QDens(ele, start_orb%t, start_orb, ecloud_data, field, custom_err)
			case (qdens_snapshot$)
				call calculateEMFromDiscrete_SingleBunch(ele, start_orb%t, start_orb, ecloud_data, field, custom_err)
			case (file_driven$)
				call calculateEMFromDiscrete_Multibunch(ele, start_orb%t, start_orb, ecloud_data, field, custom_err)
	end select

	if (custom_err) then
		print *, " ","[!] CALCULATION ERROR"
		print *," ", "[!] PANIC!!!!"
		err_flag = .true.
	end if

! Now we have calculated the fields the particle will experience, we must now calcualte kicks
xkick = (field%E(1) / ele%value(E_tot$)) 
ykick = (field%E(2) / ele%value(E_tot$))
zkick = (field%E(3) / ele%value(E_tot$))

end_orb = start_orb
end_orb%vec(2) = end_orb%vec(2) + xkick
end_orb%vec(4) = end_orb%vec(4) + ykick
end_orb%vec(6) = end_orb%vec(6) + zkick ! should always be 0 added
end_orb%s = ele%s 

end subroutine