!+
! Subroutine em_field_custom (ele, param, s_rel, time, orb, local_ref_frame, field, calc_dfield, err_flag)
!
! Routine for handling custom (user supplied) EM fields.
! This routine is called when ele%field_calc = custom$ or when ele is a custom element (ele%key = custom$)
! In order to be used, this stub file must be modified appropriately. See the Bmad manual for more details.
!
! Note: Unlike all other elements, "s_rel" and "orb" areguments for a patch element are with respect to 
! the exit reference frame of the element. See the Bmad manual for more details.
!
! Note: Fields should not have any unphysical discontinuities. 
! Discontinuities may cause Runge-Kutta integration to fail resulting in particles getting marked as "lost".
! The mode of failure here is that RK will try smaller and smaller steps to integrate through the 
! discontinuity until the step size gets lower than bmad_com%min_ds_adaptive_tracking. At this
! point the particle gets marked as lost.
!
! General rule: Your code may NOT modify any argument that is not listed as
! an output agument below.
!
! Input:
!   ele         -- Ele_struct: Custom element.
!   param       -- lat_param_struct: Lattice parameters.
!   s_rel       -- Real(rp): Longitudinal position relative to the start of the element.
!   time   -- Real(rp): Particle time.
!                 For absolute time tracking this is the absolute time.
!                 For relative time tracking this is relative to the reference particle entering the element.
!   orb        -- Coord_struct: Coords with respect to the reference particle.
!   local_ref_frame 
!               -- Logical, If True then take the 
!                     input coordinates and output fields as being with 
!                     respect to the frame of referene of the element. 
!   calc_dfield -- Logical, optional: If present and True then the field 
!                     derivative matrix is wanted by the calling program.
!
! Output:
!   field    -- Em_field_struct: Structure hoding the field values.
!   err_flag -- Logical, optional: Set true if there is an error. False otherwise.
!-
subroutine em_field_custom(ele, param, s_rel, time, orb, local_ref_frame, field, calc_dfield, err_flag)
	use bmad_interface, except_dummy => em_field_custom
	use ecloudmod
	implicit none

	! === SUBROUTINE INPUT ===
	type (ele_struct) :: ele
	type (lat_param_struct) param
	type (coord_struct), intent(in) :: orb
	real(rp), intent(in) :: s_rel, time
	logical local_ref_frame
	type (em_field_struct) :: field
	logical, optional :: calc_dfield, err_flag
	character(32) :: r_name = 'em_field_custom'
	! === ================ ===

	! === 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
	! === ================ ===

	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. !
	!==================================================================!
	if (present(calc_dfield) .and. calc_dfield) then
		print *, "Wants the derivative"
                ! Hasn't yet happened, so not yet implemented
	end if

	select case ( int(ele%value(calculation_index$)) )
			case (math_model$)
                call calculateEM(ele,orb%t,orb,field)
			case (qdens_math$)
				call calculateEMFromDiscrete_QDens(ele, orb%t, orb, ecloud_data, field, custom_err)
			case (qdens_snapshot$)
				call calculateEMFromDiscrete_SingleBunch(ele, orb%t, orb, ecloud_data, field, custom_err)
			case (file_driven$)
				call calculateEMFromDiscrete_Multibunch(ele, orb%t, orb, ecloud_data, field, custom_err)
	end select

	if (custom_err) then
		print *, " ","[!] CALCULATION ERROR"
		if (present(err_flag)) then
			err_flag = custom_err
		else
			print *," ", "[!] PANIC!!!!"
			err_flag = .true.
		end if
	end if
	!=================================================================================!
	! Now add the B field which we would get from a normal bend with given attributes !
	!=================================================================================!
	field%B(2) = ele%value(B_field$)
end subroutine