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, ffield
type(all_pointer_struct) a_ptr
REAL(rp) x, y, s, t, kicker_peak_field
real(rp) time_offset, time_units
LOGICAL  initialize_kicker/.true./
logical err_flag/.false./
logical found/.false./
logical first_344/.true./

INTEGER i,j
integer kicker_id

REAL(RP) :: pulseTiming, kick_width
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
real(rp) tempBy


character*120 string

! 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()
!  print '()','kickerCircuit == 
!write(44, '(9a12)')' kick','x','y','s','t','field%B(1)','field%B(2)','field%B(3)','PulseTim'
  initialize_kicker = .false.
ENDIF


! Get the storage B-field of the ring
field%B = [ B_radial, ele%value(B_field$) + ele%value(DB_field$), 0.0_rp ]
! if(index(ele%name,'KIC')/=0)return !temp debug 20141005
! Add kicker B-field

if(index(ele%name,'FOCUS_ELE')/= 0)then
   call focus_ele(ele,x,y,t,ffield)
   field%B=field%B+ffield%B
   field%E=field%E+ffield%E
endif

! Add quad E-field?
!IF (index(ele%name,'QUAD')/=0)  call quad_fields(ele,x,y,s,t,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 = value_of_attribute(ele, 'DTRISE_AND_FALL')  !obtain rise time and fall time here
  kick_width = value_of_attribute(ele, 'KICK_WIDTH')
  dtFall = dtRise                       !since the number of custom_attributes is limited, right now dtFall = dtRise

  ! Get/Calculate the kicker start time
  kick_tStart=0.

    call pointer_to_attribute(ele,'KICKER_ID', .false., a_ptr,err_flag)
    if(err_flag)then
      print *,' subroutine fields: line 151: KICKER_ID not found'
      stop
    endif
    kicker_id = a_ptr%r
  
  IF (kicker_id == 1) 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 - kick_width/2. - dtRise ! 20ns rise time
      ELSE ! Square
        kick_tStart = kick_tPeak - kick_width/2.
      ENDIF
    ENDIF
  ELSEIF (kicker_id == 2) 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 - kick_width/2. - dtRise ! 20ns rise time
      ELSE ! Square
        kick_tStart = kick_tPeak - kick_width/2.
      ENDIF
    ENDIF
  ELSEIF (kicker_id == 3) 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) ! 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 - kick_width/2. - dtRise ! 20ns rise time
      ELSE ! Square
        kick_tStart = kick_tPeak - kick_width/2.
      ENDIF
    ENDIF
  ELSEIF (kicker_id == 4) THEN
    IF (kicker_tStart(4) .ne. -1.) THEN
      ! User-defined kicker start time
      kick_tStart = kicker_tStart(4)
    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 - kick_width/2. - dtRise ! 20ns rise time
      ELSEIF (kickerCircuit==1)then ! Square
        kick_tStart = kick_tPeak - 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 + kick_width/2.
    pulseTiming = kickE989Timing(t,kick_tStart,dtRise,kick_width,dtFall)
!  ELSEIF (kickerCircuit==20150328 .or. kickerCircuit == 20150331 .or. kickerCircuit == 20150507 .or. &
!                      kickerCircuit == 20171113 .or. kickerCircuit == 20171118 .or. kickerCircuit == 20171119) THEN ! Blumlein with 20 ns rise/fall time
  ELSEIF (kickerCircuit>=20150328 ) THEN ! Blumlein with 20 ns rise/fall ti    kick_tPeak = kick_tStart + dtRise + kick_width/2.
    call pointer_to_attribute(ele,'KICKER_ID', .false., a_ptr,err_flag)
    if(err_flag)then
      print *,' subroutine fields: line 151: KICKER_ID not found'
      stop
    endif
    kicker_id = a_ptr%r
    pulseTiming = KickCustomTiming(kicker_id, t,kick_tStart,kickerCircuit, kickerPulseFile)
  ELSEIF (kickerCircuit == 1)then
    IF (kick_tStart<t .and. t<(kick_tStart+kick_width)) pulseTiming = 1.
  ENDIF


  ! Lastly, calculate the total B-field
  IF (abs(pulseTiming)>1.e-6) THEN
     kicker_peak_field = value_of_attribute(ele, 'KICKER_FIELD')
    IF (kickerPlates==1) THEN
      ! Uniform kicker B-field
      field%B(2) = field%B(2) + pulseTiming * (-kicker_peak_field)
    ELSEIF (kickerFieldType==1) THEN
      ! Multipole expansion
      field%B = field%B + pulseTiming * kickerBfield(ele,x,y)
    ELSE
      ! Mapped
      field%B(1) = field%B(1) - kicker_peak_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) - kicker_peak_field * pulseTiming * splin2( Kicker%x, Kicker%y, Kicker%By, Kicker%By2, 100*x, 100*y ) ! (x,y) in centimeters
      tempBy = - kicker_peak_field * pulseTiming * splin2( Kicker%x, Kicker%y, Kicker%By, Kicker%By2, 100*x, 100*y ) ! (x,y) in centimeters
     ENDIF
     if(ele%alias == 'trajectory')then
        if(first_344)write(344, '(10a12,2a11,a12)')' Peak Kicker B', 'x','y','s','t','field%B(1)','field%B(2)','field%B(3)','pulseTiming', 'kicker_peak_field','muon_number','kicker_id','kicker_peak_field'
        first_344=.false.
        write(344, '(10es12.4,2(1x,i10),es12.4)')value_of_attribute(ele,'KICKER_FIELD'), x,y,s,t,field%B,pulseTiming, pulseTiming*(-kicker_peak_field), &
              muon_number,kicker_id,kicker_peak_field
     endif
     
!write(6, '(10es12.4)')ele%value(kicker_field$),x,y,s,t,field%B,pulseTiming, tempBy
!   write(344, '(10es12.4,1x,i10)')value_of_attribute (ele, 'KICKER_FIELD'),x,y,s,t,field%B,pulseTiming, tempBy, kicker_id
  ENDIF

!print '(a,a,4es12.4)', ele%name, 'kick_tStart, dtRise, t, pulseTiming = ',kick_tStart, dtRise, t, pulseTiming
!print '(a,a,i10,4es12.4)', ele%name, 'kickercircuit,kick_width,ele%s,betaMagic, c_light =',kickercircuit,kick_width,ele%s,betaMagic, c_light

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

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 ] * (-kicker_peak_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 ] * (-kicker_peak_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 ] * (-kicker_peak_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 ] * (-kicker_peak_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 ] * (-kicker_peak_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 ] * (-kicker_peak_field)
      ENDIF
    ENDIF
  ENDIF

END FUNCTION kickerBfield




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

FUNCTION kickCustomTiming(kicker_id,t,tStart,kickerCircuit, kickerPulseFile)
  USE bmad
!  USE precision_def
  IMPLICIT NONE

 type (ele_struct) ele
  REAL(RP) :: kickCustomTiming
  REAL(RP), INTENT(IN) :: t, tStart !, dtFlat
!  REAL(RP), INTENT(INOUT) :: dtRise, dtFall
  REAL(RP) :: arg ! argument of tanh's
  REAL(RP) :: h127
  real(rp) time
  real(rp)  off/0./
  real(rp)  a127 /-2.84411/
  real(rp)  b127 /3.3351/
  real(rp)  c127 /0./
  real(rp)  d127 /1.27428/
  real(rp)  e127 /0.0970149/
  real(rp)  f127 /-3.92405/
  real(rp)  tconvert/1.e7/

  real(rp) a60/ -3.89347     /
  real(rp) b60/2.62643      /
  real(rp) d60/ 1.10409      /
  real(rp) e60/ 0.264364     /
  real(rp) f60/ -6.75644     /
  real(rp) h60, c60/0./
  real(rp), save :: pulse_time(10000), pulse_height(4,10000)
  real(rp) time_units, time_offset

  logical first/.true./
  integer kickerCircuit
  integer lun,n, ios
  integer, save :: npoints, max_pulse(4)
  integer kicker_id
  integer, save :: nwords

  character*60 file_name
  character*60 kickerPulseFile
  character*120 string

! convert to units of 100ns
!  if(first .and. (kickercircuit == 20150507 .or. kickercircuit == 20171113 .or. kickerCircuit == 20171118 .or. kickerCircuit == 20171119))then
  if(first .and. (kickercircuit >= 20150507))then
   print *, ' Use Kick Custom Timing: kickerCircuit = ', kickerCircuit
   lun=lunget()
   n=kickerCircuit
!   if(kickerCircuit == 20150507 ) file_name='pulse_127cm.dat'
!   if(n == 20171113)   file_name='kicker1data.txt'
!   if(kickerCircuit == 20171118)   file_name='kicker1data_halftail.txt'
!   if(kickerCircuit == 20171119)   file_name='kicker1data_05_tail.txt'
   if(kickerCircuit >= 20150507)file_name = kickerPulseFile
  print *,' Kick Custom TIming filename = ',file_name
   open(unit=lun, file=file_name, status='old')  !interpolate between points in table
   print *,' Open file = ',file_name
   n=0
   ios=0
   time_units  = 1.
   time_offset = 0.
   do while(ios==0)
      n=n+1
      read(lun,'(a)', IOSTAT = ios)string
      if(index(string,'#')/=0)cycle
      found = .false.
    if(index(string,'time_units')/=0)call get_value_from_string(string,'time_units',time_units, 1.0d0, found)
    if(index(string,'time_offset')/=0)call get_value_from_string(string,'time_offset',time_offset, 0.0d0,found)
    if(found)cycle
    call number_of_words(string, nwords)
    read(string,*)pulse_time(n), (pulse_height(i,n),i=1,nwords-1)
    pulse_time(n) = pulse_time(n)*time_units-time_offset
    npoints=n
    if(n > 9999)then
      print *,' Number of points > 10,000: fields.f90, line 440'
      stop
   endif
   end do

   do i=1,nwords-1
    max_pulse(i) = maxloc(pulse_height(i,:),1)
   end do
   if(kickerCircuit >= 20171113 .or. index(file_name,'KICKER1')/=0)then
     pulse_time(:) = (pulse_time(:)-200)/100.
     pulse_height(:,:) = pulse_height(:,:)/25.
   endif
!   do n=1,npoints
!    print '(i10,2es12.4)',n,pulse_time(n), pulse_height(n)
!   end do
   if(ios > 0)then
     print *,' Read Error KICK CUSTOM TIMING'
     stop
   endif
   first=.false.
  endif

  off = tstart * tconvert
  time = t*tconvert
if(kickerCircuit == 20150328)then
  h127 = b127*tanh(a127*(-d127+0.5))+c127 +b127*tanh(-f127*(e127+0.5))
  kickCustomTiming = (b127*tanh(a127*(time-d127-off))+c127 +b127*tanh(-f127*(time+e127-off)))/h127
endif

if(kickerCircuit == 20150331)then
  h60 = b60*tanh(a60*(-d60+0.3))+c60 +b60*tanh(-f60*(e60+0.3))
  kickCustomTiming = (b60*tanh(a60*(time-d60-off))+c60 +b60*tanh(-f60*(time+e60-off)))/h60
endif

!if(kickerCircuit == 20150507 .or. kickerCircuit ==20171113 .or. kickerCircuit == 20171118)then !interpolate between points in table
if(kickerCircuit >= 20150507)then !interpolate between points in table
!  h127 = b127*tanh(a127*(-d127+0.5))+c127 +b127*tanh(-f127*(e127+0.5))  !for normalization
    
    kickCustomTiming=0
    time = time - off

    if(kicker_id > nwords - 1) then
     print *,'kickCustomTiming:  there are more kickers than pulse shapes'
     stop
    endif
    do n=1,npoints
    if( time > pulse_time(n) .and. time < pulse_time(n+1))then
      kickCustomTiming = (pulse_height(kicker_id, n) + (pulse_height(kicker_id, n+1)-pulse_height(kicker_id,n)) /(pulse_time(n+1)-pulse_time(n)) * &
                               (time-pulse_time(n)))/pulse_height(kicker_id,max_pulse(kicker_id))
!      print '(a,es12.4,1x,i10,3es12.4)','time,n,pulse_time, pulse_height, kickCustomTiming ' ,time,n,pulse_time(n), pulse_height(n), kickCustomTiming  
      exit
    endif
  end do
endif ! interpolation 

  RETURN
END FUNCTION kickCustomTiming

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


   subroutine number_of_words(string_in, nwords)
    use bmad
    implicit none
    character*(*) string_in
    character*200 string
    integer nwords, word_length
    string = string_in
    nwords = 0
    call string_trim(string, string,word_length)
    do while(word_length > 0)
       nwords = nwords + 1
       call string_trim(string(word_length+1:), string,word_length)
    end do
  return
  end subroutine number_of_words


END SUBROUTINE fields