SUBROUTINE fields(ele,x,y,s,t,field)

USE bmad
USE parameters_bmad

! Field maps
USE fieldMaps
USE nrtype
USE nr

IMPLICIT NONE

TYPE(ele_struct) ele
TYPE(em_field_struct) field
REAL(rp) x, y, s, t
LOGICAL  initialize_kicker/.true./, initialize_quad/.true./

INTEGER i,j
REAL(RP) :: pulseTiming
REAL(RP) :: kick_tStart, kick_tPeak, kick_tEnd

! Move these to the input file?
REAL(RP) :: dtRise  ! E989 kicker rise time
REAL(RP) :: dtFall  ! E989 kicker fall time
REAL(RP) :: tauRC  = 5.e-6  ! quad time constant for charging/scraping

! Do these things only once
IF (index(ele%name,'KICK')/= 0 .and. initialize_kicker) THEN
  ! Kickers
  IF (kickerFieldType==1) call initializeKickerMultipoleCoeffs(KMC)
  IF (kickerFieldType==2)then

      print *,' ele%descrip = ' , ele%descrip
      kicker%filename = ele%descrip

      call initializeKickerMaps()
  endif
  IF (kickerCircuit==821) call initializeE821KickerParams()
!write(44, '(9a12)')' kick','x','y','s','t','field%B(1)','field%B(2)','field%B(3)','PulseTim'
  initialize_kicker = .false.
ENDIF

IF (initialize_quad)then
  ! Quads
  IF (quadFieldType==2)   call initializeQuadMaps()
  call initializeQuadMultipoleCoeffs(QMC)
  ! Do only once
  initialize_quad = .false.
ENDIF

! Get the storage B-field of the ring
field%B = [ 0.0_rp, ele%value(B_field$), 0.0_rp ]

! Add kicker B-field
IF (index(ele%name,'KICKER')/=0) THEN

  ! First, determine the timing of the pulse (easy/fast).  If the pulse magnitude 
  ! is very small, don't bother getting the field from the map or calculating it 
  ! from the multipole expansion.  The kicker should be at its peak when the center
  ! of the beam is at the center of the kicker.  Assume the kicker plates are s=1.76m
  ! long.  Remember, the beam is initially centered at t=0, and K2 has been split in
  ! in half for the phase-space marker.  

  dtRise = ele%value(dtRise_and_fall$)  !obtain rise time and fall time here
  dtFall = dtRise                       !since the number of custom_attributes is limited, right now dtFall = dtRise

  ! Get/Calculate the kicker start time
  kick_tStart=0.
  IF (index(ele%name,'1')/=0) THEN
    IF (kicker_tStart(1) .ne. -1.) THEN
      ! User-defined kicker start time
      kick_tStart = kicker_tStart(1)
    ELSE
      ! Automatically calculate kicker start time
      kick_tPeak = (ele%s - 1.76/2.)/(betaMagic * c_light) ! Center of K1
      IF (kickerCircuit==821) THEN
        kick_tStart = kick_tPeak - kickerE821_tStartToPeak
      ELSEIF (kickerCircuit==989) THEN
        kick_tStart = kick_tPeak - ele%value(kick_width$)/2. - dtRise ! 20ns rise time
      ELSE ! Square
        kick_tStart = kick_tPeak - ele%value(kick_width$)/2.
      ENDIF
    ENDIF
  ELSEIF (index(ele%name,'2A')/=0) THEN
    IF (kicker_tStart(2) .ne. -1.) THEN
      ! User-defined kicker start time
      kick_tStart = kicker_tStart(2)
    ELSE
      ! Automatically calculate kicker start time
      kick_tPeak = ele%s/(betaMagic * c_light) ! center of K2 (downstream end of K2A)
      IF (kickerCircuit==821) THEN
        kick_tStart = kick_tPeak - kickerE821_tStartToPeak
      ELSEIF (kickerCircuit==989) THEN
        kick_tStart = kick_tPeak - ele%value(kick_width$)/2. - dtRise ! 20ns rise time
      ELSE ! Square
        kick_tStart = kick_tPeak - ele%value(kick_width$)/2.
      ENDIF
    ENDIF
  ELSEIF (index(ele%name,'2B')/=0) THEN
    IF (kicker_tStart(2) .ne. -1.) THEN
      ! User-defined kicker start time
      kick_tStart = kicker_tStart(2)
    ELSE
      ! Automatically calculate kicker start time
      kick_tPeak = (ele%s - 1.76/2.)/(betaMagic * c_light) ! center of K2 (upstream end of K2B)
      IF (kickerCircuit==821) THEN
        kick_tStart = kick_tPeak - kickerE821_tStartToPeak
      ELSEIF (kickerCircuit==989) THEN
        kick_tStart = kick_tPeak - ele%value(kick_width$)/2. - dtRise ! 20ns rise time
      ELSE ! Square
        kick_tStart = kick_tPeak - ele%value(kick_width$)/2.
      ENDIF
    ENDIF
  ELSEIF (index(ele%name,'3')/=0) THEN
    IF (kicker_tStart(3) .ne. -1.) THEN
      ! User-defined kicker start time
      kick_tStart = kicker_tStart(3)
    ELSE
      ! Automatically calculate kicker start time
      kick_tPeak = (ele%s - 1.76/2.)/(betaMagic * c_light)
      IF (kickerCircuit==821) THEN
        kick_tStart = kick_tPeak - kickerE821_tStartToPeak
      ELSEIF (kickerCircuit==989) THEN
        kick_tStart = kick_tPeak - ele%value(kick_width$)/2. - dtRise ! 20ns rise time
      ELSE ! Square
        kick_tStart = kick_tPeak - ele%value(kick_width$)/2.
      ENDIF
    ENDIF
  ENDIF

  ! Now that we have the kicker pulse start time, let's calculate the amplitude of the kicker pulse at the particle's time (based on the start time)
  pulseTiming = 0.
  IF (kickerCircuit==821) THEN ! LCR
    IF (kick_tStart<t) pulseTiming = kickerE821_normalization*exp(-(t-kick_tStart)/kickerE821_tau)*sin(kickerE821_omega*(t-kick_tStart))
  ELSEIF (kickerCircuit==989) THEN ! Blumlein with 20 ns rise/fall time
    kick_tPeak = kick_tStart + dtRise + ele%value(kick_width$)/2.
    pulseTiming = kickE989Timing(t,kick_tStart,dtRise,ele%value(kick_width$),dtFall)
  ELSE
    IF (kick_tStart<t .and. t<(kick_tStart+ele%value(kick_width$))) pulseTiming = 1.
  ENDIF

  ! Lastly, calculate the total B-field
  IF (abs(pulseTiming)>1.e-6) THEN
    IF (kickerPlates==1) THEN
      ! Uniform kicker B-field
      field%B(2) = field%B(2) + pulseTiming * (-ele%value(kicker_field$))
    ELSEIF (kickerFieldType==1) THEN
      ! Multipole expansion
      field%B = field%B + pulseTiming * kickerBfield(ele,x,y)
    ELSE
      ! Mapped
      field%B(1) = field%B(1) - ele%value(kicker_field$) * pulseTiming * splin2( Kicker%x, Kicker%y, Kicker%Bx, Kicker%Bx2, 100*x, 100*y ) ! (x,y) in centimeters
      field%B(2) = field%B(2) - ele%value(kicker_field$) * pulseTiming * splin2( Kicker%x, Kicker%y, Kicker%By, Kicker%By2, 100*x, 100*y ) ! (x,y) in centimeters
    ENDIF
  ENDIF

!print *,'ele%alias=', ele%alias
!if(ele%alias == 'trajectory')write(44, '(9es12.4)')ele%value(kicker_field$),x,y,s,t,field%B,pulseTiming

! Add quad E-field?
ELSEIF (index(ele%name,'QUAD')/=0) THEN

  IF (quadFieldType==1) THEN
    ! Analytic (multipole expansion)
    field%E = quadEfield(ele,x,y)

  ELSE
    ! Mapped
    IF (index(ele%name,'1')/=0 .and. t<300.e-9) THEN ! save time
      IF (index(ele%name,'LONG')/=0 .and. quadPlates==989) THEN
        field%E(1) = splin2( Q1L%x, Q1L%y, Q1L%Ex, Q1L%Ex2, 100*x, 100*y ) ! (x,y) in centimeters, E-field in V/m
        field%E(2) = splin2( Q1L%x, Q1L%y, Q1L%Ey, Q1L%Ey2, 100*x, 100*y ) ! (x,y) in centimeters, E-field in V/m
      ELSE
        field%E(1) = splin2( Q1S%x, Q1S%y, Q1S%Ex, Q1S%Ex2, 100*x, 100*y ) ! (x,y) in centimeters, E-field in V/m
        field%E(2) = splin2( Q1S%x, Q1S%y, Q1S%Ey, Q1S%Ey2, 100*x, 100*y ) ! (x,y) in centimeters, E-field in V/m
      ENDIF
      ! Maps were generated for n=0.180
      field%E = field%E * ele%value(field_index$)/0.180
    ELSE ! Use the storage multipole field for t>300ns
      field%E = quadEfield(ele,x,y)
    ENDIF
  ENDIF

  ! Now we have the raw field -- let's incorporate the timing of the circuit
  IF (quadCircuit==0) THEN
    IF (index(ele%name,'1')/=0) THEN ! .or. index(ele%name,'3')/=0) THEN
      ! Off during injection, ramped up to full potential with a time constant of tauRC
      IF (t<100.e-9) THEN
        field%E(:) = 0.
      ELSE 
        field%E = field%E * ( 1-exp(-(t-100.e-9)/tauRC) )
      ENDIF
    ENDIF
  ELSEIF (quadCircuit==2) THEN
    ! TODO: Scraping -- linearly interpolate between multipole expansion coefficients
  ENDIF

ENDIF


CONTAINS
FUNCTION kickerBfield(ele,x,y)
  ! Analytic kicker B-field from multipole expansion.  Applies to kickerPlates={821,989}
  USE bmad
  USE precision_def
  IMPLICIT NONE

  TYPE(ele_struct) ele
  REAL(RP), INTENT(IN)   :: x,y
  REAL(RP), DIMENSION(3) :: kickerBfield

  INTEGER  :: n ! order of multipole
  REAL(RP) :: Bx, By, r, r0, theta
  REAL(RP) :: x_temp, Bx_inner, Bx_outer, By_inner, By_outer, frac_inner, frac_outer ! linear interpolation stuff

  ! Initialize (necessary for avoiding scoping issues for some reason)
  Bx=0.
  By=0.

  ! Calculate r, theta
  r = sqrt(x**2 + y**2)
  theta = atan2(y,x)

  IF (kickerPlates==821) THEN
    r0 = 0.045
    IF (abs(x)<0.05) THEN
      ! If the particle is between the kicker plates, use the multipole expansion
      DO n=0, size(KMC(1,:))-1
        Bx = Bx + (r/r0)**n * KMC(1,n) * sin(n*theta)
        By = By + (r/r0)**n * KMC(1,n) * cos(n*theta)
      ENDDO
      ! Assign the B-field
      kickerBfield = [ Bx, By, 0.0_rp ] * (-ele%value(kicker_field$))

    ELSEIF (0.05<abs(x) .and. abs(x)<0.05075) THEN
      ! Particle is inside the electrode -- linearly interpolate between field values at electrode surfaces
      IF (x>0.) x_temp =  0.05
      IF (x<0.) x_temp = -0.05
      r = sqrt(x_temp**2 + y**2)
      theta = atan2(y,x_temp)

      ! Get the field value at the inner electrode surface from the multipole expansion
      DO n=0, size(KMC(1,:))-1
        Bx = Bx + (r/r0)**n * KMC(1,n) * sin(n*theta)
        By = By + (r/r0)**n * KMC(1,n) * cos(n*theta)
      ENDDO

      ! Assign the field values at the inner/outer electrode surfaces
      Bx_inner =  Bx
      Bx_outer = -Bx
      By_inner =  By
      By_outer = -By

      ! Lineraly interpolate between field values at the electrode surfaces
      frac_inner = ( 0.05075-abs(x) )/( 0.05075-0.05000 )
      frac_outer = ( abs(x)-0.05000 )/( 0.05075-0.05000 )
      Bx = frac_inner*Bx_inner + frac_outer*Bx_outer
      By = frac_inner*By_inner + frac_outer*By_outer

      ! Assign the B-field
      kickerBfield = [ Bx, By, 0.0_rp ] * (-ele%value(kicker_field$))

    ELSE ! x>0.05, i.e. outside 
      ! Assume the electrode represents a plane of symmetry.  Calculate
      ! the field in the storage region, then flip both components of 
      ! the field to get the field outside the storage region
      IF (x>0.) x_temp =  0.05 - (x-0.05075)
      IF (x<0.) x_temp = -0.05 - (x+0.05075)
      r = sqrt(x_temp**2 + y**2)
      theta = atan2(y,x_temp)

      ! Get the field value from the multipole expansion
      DO n=0, size(KMC(1,:))-1
        Bx = Bx + (r/r0)**n * KMC(1,n) * sin(n*theta)
        By = By + (r/r0)**n * KMC(1,n) * cos(n*theta)
      ENDDO

      ! Flip the field
      Bx_outer = -Bx
      By_outer = -By
      Bx = Bx_outer
      By = By_outer

      ! Assign the B-field
      kickerBfield = [ Bx, By, 0.0_rp ] * (-ele%value(kicker_field$))

    ENDIF
  ELSE ! E989
    r0 = 0.04
    IF (r<0.045 .or. abs(x)<r*cos(asin(0.025/0.045))) THEN
      ! If the particle is between the kicker plates, use the multipole expansion
      DO n=0, size(KMC(2,:))-1
        Bx = Bx + (r/r0)**n * KMC(2,n) * sin(n*theta)
        By = By + (r/r0)**n * KMC(2,n) * cos(n*theta)
      ENDDO
      ! Assign the B-field
      kickerBfield = [ Bx, By, 0.0_rp ] * (-ele%value(kicker_field$))

    ELSEIF (abs(y)<0.025) THEN
      ! If the particle is outside r=45mm, yet it has |y|<25mm (i.e. the plate
      ! height), take the negative of the field calculated at the inner surface
      ! of the electrode. If the particle is inside the electrode, linearly interpolate.
      theta = asin(y/0.045)
      DO n=0, size(KMC(2,:))-1
        Bx = Bx + (0.045/r0)**n * KMC(2,n) * sin(n*theta)
        By = By + (0.045/r0)**n * KMC(2,n) * cos(n*theta)
      ENDDO

      IF (0.045<r .and. r<0.046) THEN
        ! If the particle is inside the electrode, linearly interpolate
        ! Assign the field values at the inner/outer electrode surfaces
        Bx_inner =  Bx
        Bx_outer = -Bx
        By_inner =  By
        By_outer = -By

        ! Lineraly interpolate between field values at the electrode surfaces
        frac_inner = ( 0.046-r )/( 0.046-0.045 )
        frac_outer = ( r-0.045 )/( 0.046-0.045 )
        Bx = frac_inner*Bx_inner + frac_outer*Bx_outer
        By = frac_inner*By_inner + frac_outer*By_outer
        IF (x<0.) Bx = -Bx

        ! Assign the B-field
        kickerBfield = [ Bx, By, 0.0_rp ] * (-ele%value(kicker_field$))
      ELSE
        ! If the particle is on the outside of the electrode, simply take the negative 
        ! of the B-field calculated at the inner electrode surface.
        Bx = -Bx
        By = -By
        IF (x<0.) Bx = -Bx

        ! Assign the B-field
        kickerBfield = [ Bx, By, 0.0_rp ] * (-ele%value(kicker_field$))
      ENDIF
    ENDIF
  ENDIF

END FUNCTION kickerBfield


FUNCTION quadEfield(ele,x,y)
  ! Analytic quad E-field from multipole expansion.
  USE bmad
  USE parameters_bmad
  IMPLICIT NONE
  type (ele_struct) ele
  REAL(RP), DIMENSION(3) :: quadEfield
  REAL(RP), INTENT(IN)  :: x,y

  REAL(RP) :: Ex, Ey, r, theta, Er, Etheta
  REAL(RP) :: x_temp, Ex_inner, Ex_outer
  REAL(RP) :: frac_inner, frac_outer ! linear interpolation while "inside" the quad plate
  REAL(RP) :: r0 = 0.045
  INTEGER n

  ! Initialize (avoid scoping issues)
  Ex=0.
  Ey=0.
  Er=0.
  Etheta=0.

  ! Calculate r, theta
  r = sqrt(x**2 + y**2)
  theta = atan2(y,x)

  IF (quadPlates==989 .and. index(ele%name,'1_LONG')/=0) THEN
    ! E989-style Quad1_Long
    IF (0.0705<x) THEN
      ! Assume the plane of the electrode represents a plane of vertical
      ! symmetry.  Calculate the field INSIDE the storage region, then "flip"
      ! the field about the vertical plane of symmetry to get the field OUTSIDE
      ! the storage region.
      x_temp = 0.0700 - (x-0.0705)
      r = sqrt(x_temp**2 + y**2)
      theta = atan2(y,x_temp)

      ! The 20-pole term (n=10) really screws things up in the multiplole
      ! expansion for x>5cm.  Therefore, only use the terms up to nmax=8.
      DO n=1, 8
        Er     = Er     + (n/r0) * (r/r0)**(n-1) * ( QMC(0,1,n)*cos(n*theta) + QMC(0,2,n)*sin(n*theta) )
        Etheta = Etheta + (n/r0) * (r/r0)**(n-1) * (-QMC(0,1,n)*sin(n*theta) + QMC(0,2,n)*cos(n*theta) )
      ENDDO
      Ex = cos(theta)*Er - sin(theta)*Etheta
      Ey = sin(theta)*Er + cos(theta)*Etheta

      ! Flip Ex about the plane of the electrode to go from being "inside" the
      ! storage region to being "outside" Q1_Long_Outer.
      Ex_inner =  Ex
      Ex_outer = -Ex_inner
      Ex       =  Ex_outer

      ! Assign the electric field.  Multipole coefficients came from E821 quad
      ! NIM paper with effective field index neff=0.141
      quadEfield = [ Ex, Ey, 0.0_rp ] * ele%value(field_index$)/0.141

    ELSEIF (0.7000<x .and. x<0.0705) THEN
      ! Particle is INSIDE Q1_Long_Outer electrode -- linearly interoplate between field values at electrode surfaces
      x_temp = 0.070
      r = sqrt(x_temp**2 + y**2)
      theta = atan2(y,x_temp)

      ! The 20-pole term (n=10) really screws things up in the multiplole
      ! expansion for x>5cm.  Therefore, only use the terms up to nmax=8.
      DO n=1, 8
        Er     = Er     + (n/r0) * (r/r0)**(n-1) * ( QMC(0,1,n)*cos(n*theta) + QMC(0,2,n)*sin(n*theta) )
        Etheta = Etheta + (n/r0) * (r/r0)**(n-1) * (-QMC(0,1,n)*sin(n*theta) + QMC(0,2,n)*cos(n*theta) )
      ENDDO
      Ex = cos(theta)*Er - sin(theta)*Etheta
      Ey = sin(theta)*Er + cos(theta)*Etheta

      ! Assign the field values at the inner/outer electrode surfaces
      Ex_inner =  Ex
      Ex_outer = -Ex

      ! Linearly interpolate between horizontal field values at the electrode surfaces
      frac_inner = ( 0.0705-x )/( 0.0705-0.0700 )
      frac_outer = ( x-0.0700 )/( 0.0705-0.0700 )
      Ex = frac_inner*Ex_inner + frac_outer*Ex_outer

      ! Assign the electric field.  Multipole coefficients came from E821 quad
      ! NIM paper with effective field index neff=0.141
      quadEfield = [ Ex, Ey, 0.0_rp ] * ele%value(field_index$)/0.141

    ELSE
      ! Particle is in the storage region.  The 20-pole term (n=10) really
      ! screws things up in the multipole expansion for x>5cm.  Therefore, only
      ! use the terms up to nmax=8 during injection (i.e. t<100ns).
      IF (t<100.e-9) THEN
        DO n=1,8
          Er     = Er     + (n/r0) * (r/r0)**(n-1) * ( QMC(0,1,n)*cos(n*theta) + QMC(0,2,n)*sin(n*theta) )
          Etheta = Etheta + (n/r0) * (r/r0)**(n-1) * (-QMC(0,1,n)*sin(n*theta) + QMC(0,2,n)*cos(n*theta) )
        ENDDO
      ELSE
        DO n=1,14
          Er     = Er     + (n/r0) * (r/r0)**(n-1) * ( QMC(0,1,n)*cos(n*theta) + QMC(0,2,n)*sin(n*theta) )
          Etheta = Etheta + (n/r0) * (r/r0)**(n-1) * (-QMC(0,1,n)*sin(n*theta) + QMC(0,2,n)*cos(n*theta) )
        ENDDO
      ENDIF
      Ex = cos(theta)*Er - sin(theta)*Etheta
      Ey = sin(theta)*Er + cos(theta)*Etheta

      ! Assign the electric field.  Multipole coefficients came from E821 quad
      ! NIM paper with effective field index neff=0.141
      quadEfield = [ Ex, Ey, 0.0_rp ] * ele%value(field_index$)/0.141

    ENDIF
  ELSE
    ! E821-style quad
    IF (0.0505<x) THEN
      ! Assume the plane of the electrode represents a plane of vertical
      ! symmetry.  Calculate the field INSIDE the storage region, then "flip"
      ! the field about the vertical plane of symmetry to get the field OUTSIDE
      ! the storage region.
      x_temp = 0.0500 - (x-0.0505)
      r = sqrt(x_temp**2 + y**2)
      theta = atan2(y,x_temp)

      ! Calculate the field from the full 14-term multipole expansion
      DO n=1, 14
        Er     = Er     + (n/r0) * (r/r0)**(n-1) * ( QMC(0,1,n)*cos(n*theta) + QMC(0,2,n)*sin(n*theta) )
        Etheta = Etheta + (n/r0) * (r/r0)**(n-1) * (-QMC(0,1,n)*sin(n*theta) + QMC(0,2,n)*cos(n*theta) )
      ENDDO
      Ex = cos(theta)*Er - sin(theta)*Etheta
      Ey = sin(theta)*Er + cos(theta)*Etheta

      ! Flip Ex about the plane of the electrode to go from being "inside" the
      ! storage region to being "outside" Q1_Long_Outer.
      Ex_inner =  Ex
      Ex_outer = -Ex_inner
      Ex       =  Ex_outer

      ! Assign the electric field.  Multipole coefficients came from E821 quad
      ! NIM paper with effective field index neff=0.141
      quadEfield = [ Ex, Ey, 0.0_rp ] * ele%value(field_index$)/0.141

    ELSEIF (0.5000<x .and. x<0.0505) THEN
      ! Particle is INSIDE Q1_Outer electrode -- linearly interoplate between field values at electrode surfaces
      x_temp = 0.050
      r = sqrt(x_temp**2 + y**2)
      theta = atan2(y,x_temp)

      ! Calculate the field from the full 14-term multipole expansion
      DO n=1,14
        Er     = Er     + (n/r0) * (r/r0)**(n-1) * ( QMC(0,1,n)*cos(n*theta) + QMC(0,2,n)*sin(n*theta) )
        Etheta = Etheta + (n/r0) * (r/r0)**(n-1) * (-QMC(0,1,n)*sin(n*theta) + QMC(0,2,n)*cos(n*theta) )
      ENDDO
      Ex = cos(theta)*Er - sin(theta)*Etheta
      Ey = sin(theta)*Er + cos(theta)*Etheta

      ! Assign the field values at the inner/outer electrode surfaces
      Ex_inner =  Ex
      Ex_outer = -Ex

      ! Linearly interpolate between horizontal field values at the electrode surfaces
      frac_inner = ( 0.0505-x )/( 0.0505-0.0500 )
      frac_outer = ( x-0.0500 )/( 0.0505-0.0500 )
      Ex = frac_inner*Ex_inner + frac_outer*Ex_outer

      ! Assign the electric field.  Multipole coefficients came from E821 quad
      ! NIM paper with effective field index neff=0.141
      quadEfield = [ Ex, Ey, 0.0_rp ] * ele%value(field_index$)/0.141

    ELSE
      ! Particle is in the storage region.  Use the full 14-term multipole expansion.
      DO n=1, 14
        Er     = Er     + (n/r0) * (r/r0)**(n-1) * ( QMC(0,1,n)*cos(n*theta) + QMC(0,2,n)*sin(n*theta) )
        Etheta = Etheta + (n/r0) * (r/r0)**(n-1) * (-QMC(0,1,n)*sin(n*theta) + QMC(0,2,n)*cos(n*theta) )
      ENDDO
      Ex = cos(theta)*Er - sin(theta)*Etheta
      Ey = sin(theta)*Er + cos(theta)*Etheta

      ! Assign the electric field.  Multipole coefficients came from E821 quad
      ! NIM paper with effective field index neff=0.141
      quadEfield = [ Ex, Ey, 0.0_rp ] * ele%value(field_index$)/0.141

    ENDIF
  ENDIF

END FUNCTION quadEfield


FUNCTION kickE989Timing(t,tStart,dtRise,dtFlat,dtFall)
  USE bmad
  USE precision_def
  IMPLICIT NONE

  REAL(RP) :: kickE989Timing
  REAL(RP), INTENT(IN) :: t, tStart, dtFlat
  REAL(RP), INTENT(INOUT) :: dtRise, dtFall
  REAL(RP) :: arg ! argument of tanh's
  REAL(RP) :: tPeak

  ! Avoid infinities
  IF (dtRise<1.e-12) dtRise = 1.e-12 ! ~= 0 nanoseconds
  IF (dtFall<1.e-12) dtFall = 1.e-12 ! ~= 0 nanoseconds

  ! Calculate the time at which the pulse is at its peak (i.e. the middle of the "flat" part)
  tPeak = tStart + dtRise + dtFlat/2.

  ! Let f=0.95 be the fraction of kick attained within dtRise/dtFall.  Then ArcTanh[f=0.95] = 1.83178.
  ! This number is used in computing the kick, which is based on splicing together a couple of shifted,
  ! dilated, tanh(x) functions (one of which is flipped horizontally to form the trailing part of the pulse.)
  IF (t<tPeak) THEN
    ! Leading part of the pulse
    arg =  2 * 1.83178 * (t - (tStart + dtRise/2.)) / dtRise
    kickE989Timing = ( 1 + tanh(arg) )/2.
  ELSE
    ! Trailing part of the pulse.  There is a subtlety here in which, if dtRise /= dtFall, 
    ! the pulse will have a discontinuity if left uncorrected.  As a result, we have to 
    ! shift the trailing part of the pulse up by a small amount.
    arg =  2 * 1.83178 * (tPeak - (tStart + dtRise/2.)) / dtRise
    kickE989Timing = ( 1 + tanh(arg) )/2.
    arg = -2 * 1.83178 * (tPeak - (tStart + dtRise + dtFlat + dtFall/2.)) / dtRise
    kickE989Timing = kickE989Timing - ( 1 + tanh(arg) )/2. ! discontinuity at t=tPeak
    arg = -2 * 1.83178 * (t - (tStart + dtRise + dtFlat + dtFall/2.)) / dtRise
    kickE989Timing = kickE989Timing + ( 1 + tanh(arg) )/2. ! trailing part of pulse for t>tPeak
  ENDIF

  RETURN
END FUNCTION kickE989Timing


SUBROUTINE initializeKickerMaps()

  IMPLICIT NONE
  
  ! Assign the appropriate filename
!  IF (kickerPlates==821) THEN
!    Kicker%filename = "KICKER_E821.dat"
!  ELSE
!    Kicker%filename = "KICKER_E989.dat"
!  ENDIF
  ! Read the field map from the file and create splines
  print *,' kicker%filename =', kicker%filename
  CALL readFieldMap2D( Kicker%filename, Kicker%ngx, Kicker%ngy, Kicker%x, Kicker%y, Kicker%Bx, Kicker%By, Kicker%Bx2, Kicker%By2 )
  CALL splie2( Kicker%x, Kicker%y, Kicker%Bx, Kicker%Bx2 )
  CALL splie2( Kicker%x, Kicker%y, Kicker%By, Kicker%By2 )
END SUBROUTINE initializeKickerMaps


SUBROUTINE initializeQuadMaps()
  IMPLICIT NONE

  ! Transversally, Q1 of E821 is the same as Q1_Short of E989
  Q1S%filename = "QUAD1_Short.dat"
  CALL readFieldMap2D( Q1S%filename, Q1S%ngx, Q1S%ngy, Q1S%x, Q1S%y, Q1S%Ex, Q1S%Ey, Q1S%Ex2, Q1S%Ey2 )
  CALL splie2( Q1S%x, Q1S%y, Q1S%Ex, Q1S%Ex2 )
  CALL splie2( Q1S%x, Q1S%y, Q1S%Ey, Q1S%Ey2 )
  IF (quadPlates==989) THEN
    ! QUAD1_Long (downstream) section
    Q1L%filename = "QUAD1_Long.dat"
    CALL readFieldMap2D( Q1L%filename, Q1L%ngx, Q1L%ngy, Q1L%x, Q1L%y, Q1L%Ex, Q1L%Ey, Q1L%Ex2, Q1L%Ey2 )
    CALL splie2( Q1L%x, Q1L%y, Q1L%Ex, Q1L%Ex2 )
    CALL splie2( Q1L%x, Q1L%y, Q1L%Ey, Q1L%Ey2 )
  ENDIF
END SUBROUTINE initializeQuadMaps


SUBROUTINE initializeE821KickerParams()
  USE precision_def
  USE bmad
  USE parameters_bmad
  IMPLICIT NONE

  ! PRINT *, "INITIALIZE E821 KICKER PARAMETERS"
  kickerE821_omega = sqrt( (kickerL*kickerC)**(-1) - (kickerR/(2.*kickerL))**2 )   ! angular frequency
  kickerE821_tau = 2*kickerL/kickerR                                               ! decay time constant
  kickerE821_tStartToPeak = atan(kickerE821_omega*kickerE821_tau)/kickerE821_omega ! time interval from tstart to tpeak
  kickerE821_normalization = 1/( exp(-kickerE821_tStartToPeak/kickerE821_tau) * sin(kickerE821_omega*kickerE821_tStartToPeak) ) ! normalization factor
END SUBROUTINE initializeE821KickerParams


END SUBROUTINE fields