c  *********************************************************
      subroutine ACCEL(fpar,xyz,tt,e,b)
c
c  Compute fields in RF cavity
c
      include 'icommon.inc'
c
      save jsold,jcold,kcold,cavnum
	save tstart,gap,v0,v1,v2,v3
!
      real*8 xyz(3),tt,e(3),b(3)
      real*4 k,kc,fpar(15)
      integer*4 cavnum
      integer*2 idm(136)
      character*4 prevft,ftg(4)
      data jcold/-1/,kcold/-1/,jsold/-1/,tstart/0./
c
      do i=1,3
      e(i) = 0.
      b(i) = 0.
      end do
!
      model = fpar(1)
      omrf = fpar(2) * twopi * 1.e6  ! RF angular frequency
      emagrf = fpar(3) * 1.e6        ! Ez on axis at center of cavity
! ---------------------------------------------------------
!
!     Initialization
!
! ---------------------------------------------------------
      if( jsrg.ne.jsold ) then    ! only do this at start of region 
!  Get RF phase
        phase0 = fpar(4) * pi / 180. ! initial phase for cavity #1
        dphase = fpar(5) * pi / 180. ! phase shift for adjacent cavities
! phase shift over an absorber inside a cell
        dphasea = fpar(6) * pi / 180. 
! phase shift over an absorber across a cell boundary
        dphasec = fpar(7) * pi / 180.
!  ........................................................ 
        if( phasemodel .eq. 1 ) then
!  Take phase directly from input deck
          phase = phase0
!  ........................................................
        else if( phasemodel .eq. 2 ) then
!  Use multi-cavity phase shift scheme from input deck
!
!  Find out if the previous region was an rf cavity
          prevft = 'xxxx'
!
          do 50 j=j7rec-1,1,-1
          read(7,rec=j) js,idm,ftg
          if( js .lt. 1 ) go to 50    ! watch out for pseudo-regions
          do i=1,4
          if( ftg(i) .eq. 'ACCE' ) prevft = ftg(i)
          end do
          go to 55
50        continue
!
55        continue
!   Reset cavity number to 0 if new cell
          if( jcel.ne.jcold .or. kcel.ne.kcold ) cavnum = 0              
          cavnum = cavnum + 1
!
          if( cavnum .eq. 1 ) then     ! 1st cavity in string 
            if( kcel .lt. 2 ) then     ! 1st cell ( or no cells )
              phase = phase0
            else                       ! not the first cell
              phase = phase + dphasec
            end if
!                                
          else                     ! not 1st cavity in string
!                            
            if( prevft .eq. 'ACCE' ) then  ! previous region was cavity
              phase = phase + dphase
            else                   ! previous region was not a cavity
              phase = phase + dphasea
            end if
!
          end if              ! end of test on cavnum
!  ........................................................
!  Cavity phase from reference particle
!   constant(3) or sawtooth(4) velocity profile
        else if( phasemodel.eq.3 .or. phasemodel.eq.4 ) then
          phase = -1.57079633 - omrf * trefmean + phase0
!  ........................................................
!  Cavity properties from external file
        else if( phasemodel .eq. 5 ) then
          read(rfphase,*,iostat=ioc) id,fphase,fomrf,femagrf
!
          if( ioc .ne. 0 ) then
            iflag = -51
            go to 900
          end if
          if( id .ne. jsrg ) then
            iflag = -51
            go to 900
          end if
!
          if( fphase .ne. 0. ) phase = fphase
          if( fomrf .ne. 0. ) omrf = fomrf * twopi * 1.e6
          if( femagrf .ne. 0. ) emagrf = femagrf * 1.e6
!
          phase = phase * pi / 180. + phase0
        end if      ! end of test on phasemodel
! .........................................................  
        phaserf = 360./twopi * MOD( phase, twopi )
!       WRITE(2,1000) JSRG,JPAR,JCEL,KCEL,CAVNUM,
!     &   PREVFT,XYZ(3),PHASE*180./pi
!1000   FORMAT(5i5,2x,a4,2x,2f12.4)
!  ........................................................
        if( model .eq. 4 ) then
          call CAV4_SET(fpar) 
        end if    
!
        if( model.eq.5 .and. (jcel.ne.jcold .or.
     &    (jcel.eq.jcold .and. kcel.eq.0) ) ) then     
	    call MAF2D_SET(fpar)        ! 2D Mafia RF field
        end if
!
        if( model .eq. 6 ) then
	    dt = fpar(2)      ! time offset
	    gap = fpar(3)
	    v0 = fpar(4)
	    v1 = fpar(5)
	    v2 = fpar(6)
	    v3 = fpar(7)
!
          if( fpar(8) .eq. 1. ) then
	      tstart = tpref - 2*(tpref-trefmean) + dt
          end if
        end if                
!  ........................................................
        jsold = jsrg  ! flag initialization done
        jcold = jcel
        kcold = kcel
!
      end if                ! end of test on new jsrg
! ---------------------------------------------------------
!
!     Normal tracking starts here
!
! ---------------------------------------------------------
      rp = sqrt( xyz(1)**2 + xyz(2)**2 )
      if( rp .gt. 0. ) then
        cphi = xyz(1) / rp
        sphi = xyz(2) / rp
      else
        cphi = 1.
        sphi = 0.
      end if
!  ........................................................
      if ( model .eq. 1 ) then    ! transverse-independent gap
        mode = fpar(8)
        if( mode .eq. 0 ) then    ! time-independent
          e(3) = emagrf
        else                      ! time-dependent
          targ = omrf * tt + phase
          e(3) = emagrf * COS(targ)
        end if
!  ........................................................
      else if( model .eq. 2 ) then     ! TM01p cylindrical resonator
        mode = fpar(8)
!
        if( mode .eq. 0 ) then              ! TM010 mode
          kc = omrf / cc
          rcav = 2.405 / kc
          if( rp .gt. rcav ) then
            iflag = -39
            go to 900
          end if
          rarg = kc * rp
          targ = omrf * tt + phase
!
          e(3) = emagrf * BESSJ0(rarg) * COS(targ)
          bphi = -emagrf / cc * BESSJ1(rarg) * SIN(targ)
!
          b(1) = -bphi * sphi
          b(2) =  bphi * cphi
!
        else          ! TM011 mode
          arg = (omrf/cc)**2 - (pi/slen)**2 
!
          if( arg .lt. 0. ) then              ! non-physical
            iflag = -44
            go to 900
!
          else if( arg .eq. 0. ) then        ! TM011 pi mode
            rarg = rp / 2.
            zarg = pi * xyz(3) / slen
            targ = omrf * tt + phase
            sz = SIN(zarg)
            ct = COS(targ)
!
            e(3) = emagrf * sz * ct
            er = -pi*emagrf/slen * rarg * COS(zarg) * ct
            e(1) = er * cphi
            e(2) = er * sphi
            bphi = -emagrf*omrf/(cc*cc) * rarg * sz * SIN(targ)
            b(1) = -bphi * sphi
            b(2) =  bphi * cphi
!
          else       ! TM011 (non-pi) mode
            kc = SQRT( arg )
            rcav = 2.405 / kc
            if( rp .gt. rcav ) then
              iflag = -39
              go to 900
            end if
            rarg = kc * rp
            zarg = pi * xyz(3) / slen
            targ = omrf * tt + phase
            r1 = BESSJ1(rarg)
            sz = SIN(zarg)
            ct = COS(targ)
!
            e(3) = emagrf * BESSJ0(rarg) * sz * ct
            er = -pi*emagrf/kc/slen * r1 * COS(zarg) * ct
            e(1) = er * cphi
            e(2) = er * sphi
            bphi = -emagrf*omrf/(cc*cc*kc) * r1 * sz * SIN(targ)
            b(1) = -bphi * sphi
            b(2) =  bphi * cphi
          end if                      ! end of test on arg
        end if                        ! end of test on mode 
!  ........................................................
      else if( model .eq. 3 ) then    ! travelling wave cavity
! Ref: Handbuch der Physik, V. 44, 1959, p. 345
        betw = fpar(8)
        if( betw.le.0. .or. betw.ge.1. ) then
          iflag = -45
          go to 900
        end if
        root = SQRT( 1. - betw*betw )
        k = omrf / cc
        rarg = k * rp * root / betw
        targ = omrf*tt - k*xyz(3)/betw + phase
        fac = emagrf * BESSI1(rarg) * SIN(targ)
!
        er   = -1. / root * fac
        e(1) = er * cphi
        e(2) = er * sphi
        e(3) = emagrf * BESSI0(rarg)* COS(targ)
        bphi = -betw / root / cc * fac
        b(1) = -bphi * sphi
        b(2) =  bphi * cphi
!  ........................................................
      else if( model .eq. 4  ) then
!
        call CAV4_GET(xyz,tt,rp,sphi,cphi,e,b)
!
      else if( model .eq. 5  ) then
!
        call MAF2D_GET(xyz,tt,rp,sphi,cphi,e,b)
!
      else if( model .eq. 6  ) then   ! induction linac
	  tr = tp - tstart
	  e(3) = (v0+v1*tr+v2*tr**2+v3*tr**3) / gap
!
      end if            ! end of test on models
!  --------------------------------------------------------
!
      ezrf = e(3)
      bphirf = bphi
!
900   continue
      return
      end

c  *********************************************************
      subroutine BENTSOL(fpar,x,e,b)
c
c  Compute field of bent solenoid
c
      include 'icommon.inc'
c
      real*8 x(3),e(3),b(3)
      real*4 lams,lamd,lamq,fpar(15)
!
      do i=1,3
      e(i) = 0.
      b(i) = 0.
      end do
!
      model = fpar(1)
      bsol = fpar(2)      ! solenoid strength
      dip = fpar(3)       ! dipole strength
      clens = fpar(4)     ! length of solenoid
      pref = fpar(5)
      quad = fpar(6)      ! quad strength in dipole field
! ---------------------------------------------------------
      if( model .eq. 1 ) then   ! first order Bs and Bd
!        if( x(3).lt.0. .or. x(3).gt.clens ) go to 900
!
        b(1) = quad * x(2)
        b(2) = dip + quad*x(1)
        htrk = pcharge / pref * dip
        b(3) = bsol * ( 1. - htrk*x(1) )
! n.b. there is no 1st order dipole contribution to Bs
!  cf. Carey, p. 26
        hptrk = 0.        ! constant curvature
        ordtrk = 1
!  ........................................................
      else if( model .eq. 2 ) then
! 2nd order fields with hyperbolic s-dependence
!  Require solenoid curvature = dipole curvature
        elens = fpar(7)
        lams = fpar(8)
        clend = fpar(9)
        elend = fpar(10)
        lamd = fpar(11)
        clenq = fpar(12)
        elenq = fpar(13)
        lamq = fpar(14)
        bcen = fpar(15)
!
!        zlen = MAX( clens+2.*elens, clend+2.*elend, clenq+2.*elenq )
!        if( x(3).lt.0. .or. x(3).gt.zlen ) go to 900
!
!  First get contribution from the sector dipole
        z = x(3) - elend
        a10 = dip * FTANH(z,clend,lamd,0)
        a10p = dip * FTANH(z,clend,lamd,1)
        a10pp = dip * FTANH(z,clend,lamd,2)
!
        if( quad .ne. 0. ) then
          z = x(3) - elenq
          a11 = quad * FTANH(z,clenq,lamq,0)
          a11p = quad * FTANH(z,clenq,lamq,1)
          a11pp = quad * FTANH(z,clenq,lamq,2)
	  else  
	    a11 = 0.
          a11p = 0.
	    a11pp = 0.
        end if
!
        h = pcharge * a10 / pref
        hp = pcharge * a10p / pref
!
! sextupole constraints to make div B = 0
        if( h .ne. 0. ) then
          a12 = 2.*a10pp + a11*h - a11pp/h + a10p*hp/h
          a30 = -a12 - h*a11 - a10pp
        else
          a12 = 0.
          a30 = 0.
        end if
!
        b(1) = a11*x(2) + a12*x(1)*x(2) 
        b(2) = a10 + a11*x(1) + 0.5*a12*x(1)**2 + 0.5*a30*x(2)**2
        b(3) = a10p*x(2) + (a11p-h*a10p)*x(1)*x(2)
!
!  Now get contribution from the bent solenoid itself
        z = x(3) - elens
        bs = bcen + bsol * FTANH(z,clens,lams,0)
        bsp = bsol * FTANH(z,clens,lams,1)
        hx = h * x(1)
        b(1) = b(1) - 0.5 * ( bsp*x(1) - 0.5*bs*hp*x(1)**2
     &          - 1.5*bsp*hx*x(1) )
        b(2) = b(2) - 0.5 * ( bsp*x(2) - 2.*bsp*hx*x(2)
     &          - bs*hp*x(1)*x(2) )
        b(3) = b(3) + bs * ( 1. - hx + hx*hx )
!
        htrk = h      ! load common variables for DERIV
        hptrk = hp   
        ordtrk = 2.
!
      end if      ! end of test on model
c
900   continue
      return
      end

c  *********************************************************
      subroutine BLOOP(a,z,rho,bz,br)
!
!  Get field from a single circular current loop    (unnormalized)
!
c   cf. W.R. Smythe, Static and Dynamic Electricity, p.266-7
!
      real*4 k 
      data p2/1.5707963/
!
      f1 = (a + rho)**2 + z*z
      k = 4.*a*rho / f1
      k = sqrt(k)
      f2 = (a*a+rho**2+z*z) / ((a-rho)**2+z*z)
      f3 = (a*a-rho**2-z*z) / ((a-rho)**2+z*z)
      ele = ELLE(p2,k)
      elk = ELLF(p2,k)
      if( rho .ne. 0. ) then
        br = z / (rho*sqrt(f1)) * ( -elk + f2*ele )
      else
        br = 0.
      end if
      bz = 1. / sqrt(f1) * ( elk + f3*ele )
!
      return
      end

c  *********************************************************
      subroutine BSHEET(zl,clen,rho,a,bz,br)
!
!  Get field from a single solenoidal current sheet  (unnormalized)
!
!  Ref. M. Garrett, J. Appl. Phys. 34:2567, 1963
!
      real*4 k 
      data p2/1.5707963/
!
      z = zl
!
!  Left-end of sheet
      r1 = SQRT( (a + rho)**2 + z*z )
      k = 2. / r1 * SQRT( a*rho )
      c2 = -4.*a*rho / (a+rho)**2
      ele = ELLE(p2,k)
      elk = ELLF(p2,k)
      elp = ELLPI(p2,c2,k)
!
      f1 = z*a / ((a+rho)*r1)
      f2 = (a-rho) / (2.*a)
      bz1 = f1 * ( elk + f2*(elp-elk) )
!
      if( rho .ne. 0. ) then
        f1 = r1 / (4.*rho)
        br1 = f1 * ( 2.*(elk-ele) - k*k*elk )
      else
        br1 = 0.
      end if
!
      bz = -bz1
      br = -br1      
!
!  Right-end of sheet
      z = z + clen
      r1 = SQRT( (a + rho)**2 + z*z )
      k = 2. / r1 * SQRT( a*rho )
      c2 = -4.*a*rho / (a+rho)**2
      ele = ELLE(p2,k)
      elk = ELLF(p2,k)
      elp = ELLPI(p2,c2,k)
!
      f1 = z*a / ((a+rho)*r1)
      f2 = (a-rho) / (2.*a)
      bz1 = f1 * ( elk + f2*(elp-elk) )
!
      if( rho .ne. 0. ) then
        f1 = r1 / (4.*rho)
        br1 = f1 * ( 2.*(elk-ele) - k*k*elk )
      else
        br1 = 0.
      end if
!
      bz = bz + bz1       ! note change of sign
      br = br + br1
!      
      return
      end

c  *********************************************************
      subroutine CAV4_GET(xyz,tt,rp,sphi,cphi,e,b)
c
!  Interpolate model=4 cavity fields
!
      include 'icommon.inc'
c
      real*8 xyz(3),tt,e(3),b(3)
      integer*2 sym(2)
      data epss/1e-6/
!
      z = xyz(3)               ! relative z into cavity
c
      if( rp .gt. rdm4 ) then   ! lost it radially
        ifail = -37
        go to 900
      end if
c
c  Leave the RF field zero in the drift tubes
      if((z.lt.(zcavm4-gm4)/2.).or.(z.gt.(zcavm4+gm4)/2.)) go to 900
c
c  Find field components in the gap by interpolating
      zg = ABS(z - zcavm4/2.)    ! relative z from center of gap
      call HUNT(zcrfez,nzom4,zg+epss,izlo) ! find appropriate grid element
      call HUNT(rcrfez,nrom4,rp+epss,jrlo)
      if( (izlo.lt.1) .or. (izlo.gt.nzom4-1) .or. (jrlo.lt.1) .or.
     +  jrlo.gt.nrom4-1) then
        ifail = -38
        go to 900
      end if
c
c  Electric fields
      u = ( zg - zcrfez(izlo) ) / delzm4    ! bilinear interpolation
      v = ( rp - rcrfez(jrlo) ) / delrm4
      er = (1.-u)*(1.-v)*ercrf(jrlo,izlo) + u*(1.-v)*ercrf(jrlo,izlo+1)
     +  + u*v*ercrf(jrlo+1,izlo+1) + (1.-u)*v*ercrf(jrlo+1,izlo)
      ez = (1.-u)*(1.-v)*ezcrf(jrlo,izlo) + u*(1.-v)*ezcrf(jrlo,izlo+1)
     +  + u*v*ezcrf(jrlo+1,izlo+1) + (1.-u)*v*ezcrf(jrlo+1,izlo)
c
c  Magnetic fields
c   B is defined over a grid del/2 smaller on all sides than E
c   We can interpolate for small rp by assuming B=0 on axis
c   We have to extrapolate B into the other border regions
c
      dr2 = delrm4 / 2.
      dz2 = delzm4 / 2.
      g2 = gm4 / 2.
!
      if( rp .lt. dr2 ) then
        bp = bpcrf(1,izlo) / dr2 * rp
c
      else if( rp .gt. rdm4-dr2) then
        slp = (bpcrf(nrom4-1,izlo) - bpcrf(nrom4-2,izlo)) / delrm4
        bp = slp * (rp - (rdm4-dr2)) + bpcrf(nrom4-1,izlo)
c
      else if( zg.lt.dz2 .and. rp.ge.dr2 .and. rp.le.rdm4-dr2 ) then
        slp = (bpcrf(jrlo,2) - bpcrf(jrlo,1)) / delzm4
        bp = -slp * (dz2-zg) + bpcrf(jrlo,1)
c
      else if( zg.gt.g2-dz2 .and. rp.ge.dr2 .and. rp.le.rdm4-dr2) then
        slp = (bpcrf(jrlo,nzom4-1) - bpcrf(jrlo,nzom4-2)) / delzm4
        bp = slp * (zg-(g2-dz2)) + bpcrf(jrlo,nzom4-1)
c
      else
        zb = zg - dz2    ! use shifted grid for B
        rb = rp - dr2
        u = ( zb - zcrfez(izlo) ) / delzm4
        v = ( rb - rcrfez(jrlo) ) / delrm4
        bp = (1.-u)*(1.-v)*bpcrf(jrlo,izlo)
     &    + u*(1.-v)*bpcrf(jrlo,izlo+1)
     &    + u*v*bpcrf(jrlo+1,izlo+1) + (1.-u)*v*bpcrf(jrlo+1,izlo)
      end if
c
c  Determine field symmetry, depending on which quadrant particle is in
      do i=1,2
      sym(i) = 1
      end do
c
      if( z .lt. zcavm4/2. )then    ! 2nd and 3rd quad
        sym(1) = -1
        sym(2) = -1
      end if
c
      arg = omrf * tt + phase
      carg = COS(arg)
      sarg = SIN(arg)
c
      e(1) = er * cphi * carg * sym(1)
      e(2) = er * sphi * carg * sym(2)
      e(3) = ez * carg
      b(1) = -bp * sphi * sarg
      b(2) = bp * cphi * sarg
!
900   continue  
      return
      end
c  *********************************************************
      subroutine CAV4_SET(fpar)
c
!  Set up model=4 cavity fields
!
      include 'icommon.inc'
!
      parameter ( nsrc = 15 )
      real*4 fpar(15),rsrc(nsrc),zsrc(nsrc)
c
      zcavm4 = fpar(7)     ! total length of cavity
      gm4 = fpar(8)        ! gap
      rdm4 = fpar(9)       ! radius of drift tube
      rn = fpar(10)        ! radius of nose piece
!
c  Generate field grid  -----------------------------------
      g2 = gm4 / 2.
      rcn = rdm4 + rn
      delzm4 = g2 / nzom4
      delrm4 = rdm4 / nrom4
      dr2 = delrm4 / 2.
      dz2 = delzm4 / 2.
!
c  Generate relative source point locations on nose
      th = 0.
      dth = pi / 2. / nsrc
c
      do i=1,nsrc
      rsrc(i) = rn * cos(th)
      zsrc(i) = rn * sin(th)
      th = th + dth
      end do
c
c  Calculate electric field components over grid of observation
c   points in the gap
c
      do 200 iz=1,nzom4
      zcrfez(iz) = (iz-1) * delzm4
c
      do 150 ir=1,nrom4
      rcrfez(ir) = (ir-1) * delrm4
      ez = 0.
      er = 0.
c
      do 100 is=1,nsrc
      dr = rcrfez(ir) - ( rcn - rsrc(is) )    ! upper right nose
      dz = zcrfez(iz) - ( g2 - zsrc(is) )
      r3 = ( dr*dr + dz*dz )**(1.5)
      ez = ez - dz / r3
      er = er - dr / r3
c
      dr = rcrfez(ir) - ( rcn - rsrc(is) )    ! upper left nose
      dz = zcrfez(iz) - ( -g2 + zsrc(is) )
      r3 = ( dr*dr + dz*dz )**(1.5)
      ez = ez + dz / r3
      er = er + dr / r3
c
      dr = rcrfez(ir) - ( -rcn + rsrc(is) )   ! lower left nose
      dz = zcrfez(iz) - ( -g2 + zsrc(is) )
      r3 = ( dr*dr + dz*dz )**(1.5)
      ez = ez + dz / r3
      er = er + dr / r3
c
      dr = rcrfez(ir) - ( -rcn + rsrc(is) )   ! lower right nose
      dz = zcrfez(iz) - ( g2 - zsrc(is) )
      r3 = ( dr*dr + dz*dz )**(1.5)
      ez = ez - dz / r3
      er = er - dr / r3
100   continue
c
      if( ir .eq. 1 .and. iz .eq. 1 ) escal = emagrf / ez
      ezcrf(ir,iz) = ez * escal
      ercrf(ir,iz) = er * escal
150   continue
c
200   continue
c
c  Calculate magnetic field from Maxwell equations
c
      do 300 iz=1,nzom4-1
c
      do 250 ir=1,nrom4-1
      derdz = ercrf(ir,iz+1)-ercrf(ir,iz)+ercrf(ir+1,iz+1)
     &  -ercrf(ir+1,iz)
      derdz = derdz / 2. / delzm4
      dezdr = ezcrf(ir+1,iz)-ezcrf(ir,iz)+ezcrf(ir+1,iz+1)
     &  -ezcrf(ir,iz+1)
      dezdr = dezdr / 2. / delrm4
      bpcrf(ir,iz) =  (dezdr - derdz) / omrf
250   continue
c
300   continue
c
900   continue
      return
      end
c  *********************************************************
      subroutine COIL(fpar,xyz,e,b)
c
c  Compute field of circular coil segment
c
      include 'icommon.inc'
c
      save jsold,jcold
      real*8 xyz(3),e(3),b(3)
      real*4 fpar(15)
      character*10 fname
      data bnorm/2.e-7/   
      data jcold/-1/,jsold/-1/
c
      model = fpar(1)
c
c  Compute field from sum of coil segments
!
      do i=1,3
      e(i) = 0.
      b(i) = 0.
      end do
!
      rp = sqrt( xyz(1)**2 + xyz(2)**2 )
      if( xyz(1)*xyz(2) .ne. 0. ) then
        phi = ATAN2( xyz(2),xyz(1) )
      else
        phi = 0.
      end if
      sphi = sin(phi)
      cphi = cos(phi)
! ---------------------------------------------------------
      if( model .eq. 1 ) then   ! exact field from single loop
        a = fpar(3)              ! coil radius
        curr = fpar(4)           ! current
        z = xyz(3) - fpar(2)
!
        call BLOOP(a,z,rp,bz,br)
!
        b(1) = br * cphi* curr
        b(2) = br * sphi * curr
        b(3) = bz * curr
!
        do i=1,3
        b(i) = b(i) * bnorm
        end do
      end if
!  --------------------------------------------------------
!  Initialize for models 2 and 3
!
      if( (kbcell.eq.0.and.jsrg.ne.jsold) .or. 
     &      (kbcell.eq.1.and.jcel.ne.jcold) ) then
!
        if( model.eq.2 .or. model.eq.3 ) then   ! need data file
          mfile = fpar(2)
          write(fname,40) mfile
40        format('for0',i2,'.dat')
          open (
     &     unit=mfile,
     &     file=pathName(1:Len_Trim(pathName)) // fname,
     &     status='unknown',
     &     iostat=ioc)
          if( ioc .ne. 0 ) then
            iflag = -34
            go to 900
          end if
          write(2,*)' COIL SUMMARY for file',mfile
          read(mfile,*,iostat=ioc)ncoils
          if( ioc .ne. 0 ) then
            iflag = -35
            go to 900
          end if
          write(2,'(a,i5)') ' NCOILS = ',ncoils
!
          do i=1,ncoils
          read(mfile,*,iostat=ioc)idm,dzelgr(i),relgr(i),curelgr(i)
          if( ioc .ne. 0 ) then
            iflag = -35
            go to 900
          end if
          write(2,50) idm,dzelgr(i),relgr(i),curelgr(i)
50        format(' #,DZ,R,I: ',i5,3e12.4)
          end do
!
          close(mfile)
!         
        end if            ! end of test on model = 2 or 3
!  ........................................................
      if( model.eq.3 ) then              ! need grid set up
        delzgr = fpar(3)
        delrgr = fpar(4)
        zglen = fpar(5)
        rglen = fpar(6)
        dzcut = fpar(7)
        mfile = fpar(9)
        if( delzgr .le. 0. .or. delrgr .le. 0. ) then
          iflag = -41
          go to 900
        end if
        nzgr = ANINT( zglen / delzgr + 1. )
        nrgr = ANINT( rglen / delrgr + 1. )
        if( nzgr.gt.mxzgr .or. nrgr.gt.mxrgr .or. krad.ne.1) then
          iflag = -42
          go to 900
        end if
!  ........................................................ 
        do iz=1,nzgr
        if( MOD(iz,20).eq.1 .and. termout ) then
          write(*,80)' Computing grid:',iz,' /',nzgr,' done'
80        format(a,i5,a,i5,a)
        end if
        zgr(iz) = (iz-1) * delzgr              
!       
        do ir=1,nrgr
        rgr(ir) = (ir-1) * delrgr
!
! Check that the bozo isn't asking for a grid point on a coil
          do i=1,ncoils
          if( zgr(iz).eq.dzelgr(i) .and. rgr(ir).eq.relgr(i) ) then
            write(2,*) ' Grid point requested on a coil'
            if (termout) write(*,*) ' Grid point requested on a coil'
            jcold = jcel
            jsold = jsrg
            iflag = -1
            go to 900
          end if
          end do
!
        bzgr(iz,ir) = 0.
        brgr(iz,ir) = 0.
!
        do 100 ic=1,ncoils
        a = relgr(ic)
        z = zgr(iz) - dzelgr(ic)
        if( ABS(z) .gt. dzcut ) go to 100              
!
        call BLOOP( a,z,rgr(ir),bz,br )
!
        bzgr(iz,ir) = bzgr(iz,ir) + bz * curelgr(ic)
        brgr(iz,ir) = brgr(iz,ir) + br * curelgr(ic)
100     continue              ! sum on coils
!
        bzgr(iz,ir) = bzgr(iz,ir) * bnorm
        brgr(iz,ir) = brgr(iz,ir) * bnorm
        end do              ! r grid
!
        end do              ! z grid
!
          if( mfile .gt. 19 ) then    ! save field grid to file
            write(fname,40) mfile
            open (
     &       unit=mfile,
     &       file=pathName(1:Len_Trim(pathName)) // fname,
     &       status='unknown',
     &       iostat=ioc)
            if( ioc .ne. 0 ) then
              iflag = -50
              go to 900
            end if
            write(mfile,'(i8)') nzgr*nrgr
            do i=1,nzgr
            do j=1,nrgr
            write(mfile,'(2i6,4e12.4)') i,j,zgr(i),rgr(j),
     &        bzgr(i,j),brgr(i,j)
            end do
            end do
            close(mfile)
          end if
!
      end if           ! end of test on model = 3
!
      jcold = jcel
      jsold = jsrg
      end if                  ! end of test on kbcell ...
! End of initialization
!  --------------------------------------------------------
      if( model .eq. 2 ) then
! Exact field from sum of coils in data file
!
        do i=1,ncoils
        a = relgr(i)
        z = xyz(3) - dzelgr(i)
!
        call BLOOP(a,z,rp,bz,br)
!
        b(1) = b(1) + br * cphi * curelgr(i)
        b(2) = b(2) + br * sphi * curelgr(i)
        b(3) = b(3) + bz * curelgr(i)
        end do
c
        do i=1,3
        b(i) = b(i) * bnorm
        end do
!
      end if
!  --------------------------------------------------------
      if( model .eq. 3 ) then
! Interpolated grid field from sum of coils in data file
        z = xyz(3)               
        call HUNT(zgr,nzgr,z,jzlo)
!  Watch out for initial point exactly on left edge of grid
        if( jzlo .eq. 0 ) jzlo = 1                                          
        call HUNT(rgr,nrgr,rp,jrlo)
        if( jrlo .eq. 0 ) jrlo = 1                                          
        if( z .ge. zgr(nzgr) ) go to 900     ! B=0 beyond grid
        if( rp .ge. rgr(nrgr) ) go to 900
!            
        interp = fpar(8)      ! interpolation model
!
        if( interp .le. 1 ) then   ! bi-linear interpolation
         u = (z - zgr(jzlo)) / delzgr          
         v = (rp - rgr(jrlo)) / delrgr
         br = (1.-u)*(1.-v)*brgr(jzlo,jrlo)+ u*(1.-v)*brgr(jzlo+1,jrlo)
     &    + u*v*brgr(jzlo+1,jrlo+1) + (1.-u)*v*brgr(jzlo,jrlo+1)
         bz = (1.-u)*(1.-v)*bzgr(jzlo,jrlo)+ u*(1.-v)*bzgr(jzlo+1,jrlo)
     &    + u*v*bzgr(jzlo+1,jrlo+1) + (1.-u)*v*bzgr(jzlo,jrlo+1)
!
        else if( interp.eq.2 .or. interp.eq.3 ) then 
!  Polynomial interpolation
         call BIPOLYINT(jzlo,jrlo,zgr,rgr,bzgr,mxzgr,mxrgr,interp,
     &     z,rp,bz,df)     
         call BIPOLYINT(jzlo,jrlo,zgr,rgr,brgr,mxzgr,mxrgr,interp,
     &     z,rp,br,df)     
!
         end if        ! end of test on interpolation
!
        b(1) = br * cphi
        b(2) = br * sphi
        b(3) = bz
      end if          ! end of test on model
c
900   continue
      return
      end

c  *********************************************************
      subroutine DEFLECTION(fpar,xyz,tt,e,b)
c
c  Compute fields for RF deflection cavity
c
      include 'icommon.inc'
c
      real*8 xyz(3),tt,e(3),b(3),x(3)
      real*4 fpar(15)
c
      do i=1,3
      e(i) = 0.
      b(i) = 0.
      end do
!
      model = fpar(1)
      om = fpar(2) * 2.*pi * 1.e6  ! RF angular frequency
      emagrf = fpar(3) * 1.e6        ! electric field strength
      phase = fpar(4) * pi / 180.  ! initial phase
      a = fpar(5)                                                                       ! width
! 
      if( ABS( xp(1) ) .gt. a ) then
        iflag = -47
        go to 900
      end if
!
      wl = 2. * pi * cc / om
      d = slen
      targ = om * tt + phase
! .........................................................
      if ( model .eq. 1 ) then             ! TE011 mode
        arg = (2./wl)**2 - (1./d)**2
        if( arg .gt. 0. ) then
          h = 1. / SQRT( arg )
        else
          iflag = -46
          go to 900
        end if
!
        if( ABS( xp(2) ) .gt. h ) then
          iflag = -47
          go to 900
        end if
!
!  Translate to coordinate system with origin at corner of box
        x(2) = xyz(2) + h / 2.
!
        harg = pi * x(2) / h
        zarg = pi * xyz(3) / d
!
        e(1) = -emagrf * SIN(harg) * SIN(zarg) * SIN(targ)
        b(2) = -emagrf/cc*wl/(2.*d) * SIN(harg)* COS(zarg) * COS(targ)
        b(3) = emagrf/cc*wl/(2.*h) * COS(harg)* SIN(zarg) * COS(targ)
!
! .........................................................
      else if ( model .eq. 2 ) then       ! TM210 mode
        arg = (2./wl)**2 - (2./a)**2
        if( arg .gt. 0. ) then
          h = 1. / SQRT( arg )
        else
          iflag = -46
          go to 900
        end if
!
        if( ABS( xp(2) ) .gt. h ) then
          iflag = -47
          go to 900
        end if
!
!  Translate to coordinate system with origin at corner of box
        x(1) = xyz(1) + a / 2.
        x(2) = xyz(2) + h / 2.
!
        warg = 2. * pi * x(1) / a
        harg = pi * x(2) / h
!
        e(3) = emagrf * SIN(warg) * SIN(harg) * COS(targ)
        b(1) = -emagrf/cc * wl/2./h * SIN(warg)* COS(harg) * SIN(targ)
        b(2) = emagrf/cc *wl/a * COS(warg)* SIN(harg) * SIN(targ)
      end if
!
900   continue
      return
      end

c  *********************************************************
      subroutine DIPOLE(fpar,x,e,b)
c
c  Compute field of sector dipole
c
      include 'icommon.inc'
c
      real*8 x(3),e(3),b(3)
      real*4 lamd,lamq,lams,fpar(15)
!
      do i=1,3
      e(i) = 0.
      b(i) = 0.
      end do
!
      model = fpar(1)
      d0  = fpar(2)            ! dipole strength
      ordtrk = fpar(3)
      pref = fpar(4)           ! reference momentum
!
!  --------------------------------------------------------
      if( model .eq. 1 ) then          ! constant field
        b(2) = d0
        h = pcharge / pref * d0
        hp = 0.        ! constant curvature
!  --------------------------------------------------------
!
      else if( model .eq. 2 ) then     ! TANH(s) model
!
        clend = fpar(5)            ! length
        elend = fpar(6)            ! dipole end length
        lamd = fpar(7)             ! dipole attenuation length
        g0 = fpar(8)               ! quadrupole component strength
        s0 = fpar(9)               ! setupole component strength
        clenq = fpar(10)
        elenq = fpar(11)
        lamq = fpar(12)
        clens = fpar(13)
        elens = fpar(14)
        lams = fpar(15)
!
        z = x(3) - elend 
        a10 = d0 * FTANH(z,clend,lamd,0)
        a10p = d0 * FTANH(z,clend,lamd,1)
        a10pp = d0 * FTANH(z,clend,lamd,2)
        h = pcharge / pref * a10
        hp = pcharge / pref * a10p
!
        if( g0 .ne. 0. ) then
          z = x(3) - elenq 
          a11 = g0 * FTANH(z,clenq,lamq,0)
          a11p = g0 * FTANH(z,clenq,lamq,1)
          a11pp = g0 * FTANH(z,clenq,lamq,2)
        else
          a11 = 0.
          a11p = 0.
          a11pp = 0.
        end if
!              
        if( s0 .ne. 0. ) then
          z = x(3) - elens 
          a12 = s0 * FTANH(z,clens,lams,0)
          a12p = s0 * FTANH(z,clens,lams,1)
          a12pp = s0 * FTANH(z,clens,lams,2)
        else
	    a12 = 0.
	    a12p = 0.
	    a12pp = 0.
        end if
! .........................................................             
        if( ordtrk .eq. 1 ) then
          b(1) = g0 * x(2)
          b(2) = d0 + g0*x(1)
          b(3) = 0.
! n.b. there is no 1st order dipole contribution to Bs
!  cf. Carey, p. 26
          h = pcharge / pref * d0
          hp = 0.        ! constant curvature
!  ........................................................
        else if( ordtrk .eq. 2 ) then
!
! Sextupole constraints to make div B = 0
          a12 = 2.*a10pp + a11*h - a11pp/h + a10p*hp/h
          a30 = -a12 - h*a11 - a10pp
!
          b(1) = a11*x(2) + a12*x(1)*x(2) 
          b(2) = a10 + a11*x(1) + 0.5*a12*x(1)**2 + 0.5*a30*x(2)**2
          b(3) = a10p*x(2) + (a11p - h*a10p)*x(1)*x(2)
!  ........................................................
        else if( ordtrk .eq. 3 .or. ordtrk .eq. 4 ) then
!
          a10ppp = d0 * FTANH(z,clend,lamd,3)
!  Octupole constraint
          a13 = -a12pp/h + 4.*a11pp - 6.*a10pp*h + 2.*a12*h
     &      - 2.*a11*h**2              
          a30 = -a12 - h*a11 - a10pp
          a31 = -a13 - a11pp + 2.*h*a10pp + h*h*a11 + hp*a10p - h*a12
          a30p = -a12p - a11p*h - a10ppp
!
          b(1) = a11*x(2) + a12*x(1)*x(2) +0.5*a13*x(1)**2*x(2)
     &     + a31/6.*x(2)**3
          b(2) = a10 + a11*x(1) + 0.5*a12*x(1)**2 + 0.5*a30*x(2)**2
     &     + 0.5*a31*x(1)*x(2)**2 + a13/6.*x(1)**3
          b(3) = a10p*x(2) + (a11p - h*a10p)*x(1)*x(2)
     &     + (0.5*a12p - a11p*h + a10p*h*h)*x(1)**2*x(2)
     &     + a30p/6.*x(2)**3
!
        end if              ! end test on order
! .........................................................
      end if                  ! end test on model
!  --------------------------------------------------------
      htrk = h
      hptrk = hp
!
900   continue
      return
      end

c  *********************************************************
      subroutine FOFO(fpar,xyz,e,b)
c
c  Get field from solenoidal FOFO lattice
c
      include 'icommon.inc'
c
      real*8 xyz(3),e(3),b(3)
      real*4 fpar(15)
c
      do i=1,3
      b(i) = 0.
      e(i) = 0.
      end do
!
      model = fpar(1)
      bmag = fpar(2) 
      bcen = fpar(3)
      period = fpar(4)
!      zlen = fpar(6)      ! total length of FOFO
!      if( xyz(3).lt.0. .or. xyz(3).gt.zlen ) go to 900
!             
      rp = sqrt( xyz(1)**2 + xyz(2)**2 )
      if( xyz(1)*xyz(2) .ne. 0. ) then
        phi = ATAN2( xyz(2),xyz(1) )
      else
        phi = 0.
      end if
!
      d = xyz(3) - fpar(5)
c  ........................................................
      if( model .eq. 1 ) then             ! linear ramp
        d = MOD( d+period, period )
        p4 = period / 4.
        g = bmag / p4
!
        if( d .le. p4 ) then
          b(3) = bcen + g * d
          br = -0.5 * g * rp
        else if( (p4.lt.d) .and. (d.le.2*p4) )then
          b(3) = bcen + g * ( 2.*p4 - d )
          br = 0.5 * g * rp
        else if( (2*p4.lt.d) .and. (d.le.3*p4) )then
          b(3) = bcen - g * ( d - 2.*p4 )
          br = 0.5 * g * rp
        else 
          b(3) = bcen - g * ( period - d )
          br = -0.5 * g * rp
        end if
!  ........................................................
      else if( model .eq. 2 ) then    ! sinusoidal ramp with Bessel r
        xk = 2. * pi / period
        b(3) = bcen + bmag * SIN( xk*d ) * BESSI0( xk*rp )
        br= -bmag * COS( xk*d ) * BESSI1( xk*rp )
      end if
c
      b(1) = br * cos(phi)
      b(2) = br * sin(phi)
c
900   continue
      return
      end

c  *********************************************************
      function FTANH(z,clen,lamb,order)
!
!  Calculates difference of TANH functions and derivatives
!   for z dependence of field strengths
!
      real*4 lamb
      integer*4 order
!
      d = z / lamb
      g = (z - clen) / lamb
      scd = 1. / COSH(d)
      scg = 1. / COSH(g)
      thd = TANH(d)
      thg = TANH(g)
!
      if( order .eq. 0 ) then
        ftanh = 0.5 * ( thd - thg )
!
      else if( order .eq. 1 ) then
        ftanh = 0.5/lamb * ( scd**2 - scg**2 )
!
      else if( order .eq. 2 ) then
        ftanh = 1./lamb**2 * ( -thd * scd**2 + thg * scg**2 )
!
      else if( order .eq. 3 ) then
        ftanh = 1./lamb**3 * ( -scd**4 + 2.*thd**2 * scd**2
     &      + scg**4 - 2.*thg**2 * scg**2 )
!
      else if( order .eq. 4 ) then
        ftanh = 4./lamb**4 * ( 2.*scd**4 * thd
     &    - thd**3 * scd**2 - 2.*scg**4 * thg + thg**3 * scg**2 )
!
      else if( order .eq. 5 ) then
        ftanh = 4./lamb**5 * ( -11.* scd**4 * thd**2
     &    + 2.*scd**6 + 2.*thd**4 * scd**2
     &    + 11.*scg**4 * thg**2 - 2.*scg**6 - 2.*thg**4 * scg**2 )
!
      else if( order .eq. 6 ) then
        ftanh = 8./lamb**6 * ( 26.*scd**4 * thd**3
     &      - 17.*thd * scd**6 - 2.*scd**2 * thd**5
     &      - 26.*thg**3 * scg**4 + 17.*scg**6 * thg
     &      + 2.*thg**5 * scg**2 )
!
      else if( order .eq. 7 ) then
        ftanh = 8./lamb**7 * ( -114.* scd**4 * thd**4
     &    + 180.*scd**6 * thd**2 - 17.* scd**8
     &    + 4.*scd**2 * thd**6 + 114.*scg**4 * thg**4
     &    - 180.*thg**2 * scg**6 + 17.*scg**8 - 4.*thg**6 * scg**2 )
!
      end if
!
      return
      end
c  *********************************************************
      subroutine MAF2D_GET(xyz,tt,rp,sphi,cphi,e,b)
c
!  Interpolate 2D Mafia cavity fields
!
      include 'icommon.inc'
c
      real*8 xyz(3),tt,e(3),b(3)
      integer*2 sym(2)
      data epss/1.e-6/
!
      if( rp .gt. rccrf ) then   ! lost it radially
        ifail = -57
        go to 900
      end if
!
!      if( xyz(3) .le. z2crf ) then
!	  zp = xyz(3)
!      else
!	  zp = 2.*z2crf - xyz(3)
!      end if
!
      zp = ABS(xyz(3) - z2crf)    ! relative z from center of gap
!
! For all 3 field components 
!  The radial grid is assumed to start at r=0 and extend to Rmax > rcrrf
!  The axial grid is assumed to start at z=0 and extend to Zmax > z2crf
! We assume the basic grid is defined in the 2nd quadrant
! .........................................................
! Ez component
!   
      call HUNT(zcrfez,nzcrfez,zp+epss,izlo) 
      call HUNT(rcrfez,nrcrfez,rp+epss,jrlo)
      if( (izlo.lt.1) .or. (izlo.gt.nzcrfez)
     &  .or. (jrlo.lt.1) .or. jrlo.gt.nrcrfez ) then
	    stop 'Ez grid error in MAF2D_GET'
      end if
!
      dzcrf = zcrfez(izlo+1) - zcrfez(izlo)
      drcrf = rcrfez(jrlo+1) - rcrfez(jrlo)
!
      u = ( zp - zcrfez(izlo) ) / dzcrf    ! bilinear interpolation
      v = ( rp - rcrfez(jrlo) ) / drcrf
      ez = (1.-u)*(1.-v)*ezcrf(jrlo,izlo) + u*(1.-v)*ezcrf(jrlo,izlo+1)
     &  + u*v*ezcrf(jrlo+1,izlo+1) + (1.-u)*v*ezcrf(jrlo+1,izlo)
! .........................................................
! Er component
!   
      call HUNT(zcrfer,nzcrfer,zp+epss,izlo) 
      call HUNT(rcrfer,nrcrfer,rp+epss,jrlo)
      if( (izlo.lt.1) .or. (izlo.gt.nzcrfer)
     &  .or. (jrlo.lt.1) .or. jrlo.gt.nrcrfer ) then
	    stop 'Er grid error in MAF2D_GET'
      end if
!
      dzcrf = zcrfer(izlo+1) - zcrfer(izlo)
      drcrf = rcrfer(jrlo+1) - rcrfer(jrlo)
!
      u = ( zp - zcrfer(izlo) ) / dzcrf    ! bilinear interpolation
      v = ( rp - rcrfer(jrlo) ) / drcrf
      er = (1.-u)*(1.-v)*ercrf(jrlo,izlo) + u*(1.-v)*ercrf(jrlo,izlo+1)
     &  + u*v*ercrf(jrlo+1,izlo+1) + (1.-u)*v*ercrf(jrlo+1,izlo)
! .........................................................
! Bphi component
!   
      call HUNT(zcrfbp,nzcrfbp,zp+epss,izlo) 
      call HUNT(rcrfbp,nrcrfbp,rp+epss,jrlo)
      if( (izlo.lt.1) .or. (izlo.gt.nzcrfbp)
     &  .or. (jrlo.lt.1) .or. jrlo.gt.nrcrfbp ) then
	    stop 'Bphi grid error in SYMRF_GET'
      end if
!
      dzcrf = zcrfbp(izlo+1) - zcrfbp(izlo)
      drcrf = rcrfbp(jrlo+1) - rcrfbp(jrlo)
!
      u = ( zp - zcrfbp(izlo) ) / dzcrf    ! bilinear interpolation
      v = ( rp - rcrfbp(jrlo) ) / drcrf
      bp = (1.-u)*(1.-v)*bpcrf(jrlo,izlo) + u*(1.-v)*bpcrf(jrlo,izlo+1)
     &  + u*v*bpcrf(jrlo+1,izlo+1) + (1.-u)*v*bpcrf(jrlo+1,izlo)
! .........................................................
!
!  Determine field symmetry, depending on which quadrant particle is in
      do i=1,2
      sym(i) = 1
      end do
!
      if( xyz(3) .gt. z2crf ) then
        sym(1) = -1
        sym(2) = -1
      end if
c
      arg = omrf * tt + phase
      carg = COS(arg)
      sarg = SIN(arg)
c
      e(1) = er * cphi * carg * sym(1)
      e(2) = er * sphi * carg * sym(2)
      e(3) = ez * carg
      b(1) = -bp * sphi * sarg
      b(2) = bp * cphi * sarg
!
900   continue  
      return
      end
c  *********************************************************
      subroutine MAF2D_SET(fpar)
!
! Load arrays from user-supplied azimuthally symmetric RF field
!
      include 'icommon.inc'
!
      real*4 fpar(15)
	character*80 title
	character*10 fname
!
	fnrm = fpar(3)
	if( fnrm .eq. 0. ) fnrm = 1.
      mfile = fpar(5)
	rccrf = fpar(6)
	z2crf = fpar(7)
	nrez = fpar(8)
	nzez = fpar(9)
	nrer = fpar(10)
	nzer = fpar(11)
	nrbp = fpar(12)
	nzbp = fpar(13)
!
      write(fname,105) mfile
105   format('for0',i2,'.dat')
      open (
     & unit=mfile,
     & file=pathName(1:Len_Trim(pathName)) // fname,
     & status='unknown',
     & iostat=ioc)
      if( ioc .ne. 0 ) then
        iflag = -54
        go to 900
      end if
	write(2,107) mfile,j7rec
107   format(' Opened user RF file: ',i3,'  for region: ',i5)
!
      read(mfile,'(a80)',iostat=ioc) title
	write(2,'(a80)') title
! .........................................................
      read(mfile,*,iostat=ioc) nrcrfez,nzcrfez
      if( ioc .ne. 0 ) then
        iflag = -55
        go to 900
      end if
!
      if( nrcrfez.ne.nrez .or. nzcrfez.ne.nzez .or.
     &   nrcrfez.le.0 .or. nzcrfez.le.0 ) then
	  iflag = -56
	  go to 900
	end if
!
! 2D-Mafia produces an Ez grid starting with r=0 and finite z
	do i=1,nrcrfez       
      read(mfile,*,iostat=ioc) id,rcrfez(i)
      if( ioc .ne. 0 ) then
        iflag = -55
        go to 900
      end if
	end do
!
! Add a new mesh point at z=0
      zcrfez(1) = 0.
	nzcrfez = nzcrfez + 1
!
	do i=2,nzcrfez       
      read(mfile,*,iostat=ioc) id,zcrfez(i)
      if( ioc .ne. 0 ) then
        iflag = -55
        go to 900
      end if
	end do
!
	do i=2,nzcrfez
      do j=1,nrcrfez       
      read(mfile,*,iostat=ioc) id,id,ezcrf(j,i)
      if( ioc .ne. 0 ) then
        iflag = -55
        go to 900
      end if
	ezcrf(j,i) = ezcrf(j,i) * fnrm
	end do
	end do
!
! Extrapolate Ez(r,0) values
      do i=1,nrcrfez
	slp = (ezcrf(i,3)-ezcrf(i,2))/(zcrfez(3)-zcrfez(2))
	ezcrf(i,1) = ezcrf(i,2) - slp*zcrfez(2)
	end do
!
! .........................................................
      read(mfile,*,iostat=ioc) nrcrfer,nzcrfer
      if( ioc .ne. 0 ) then
        iflag = -55
        go to 900
      end if
!
      if( nrcrfer.ne.nrer .or. nzcrfer.ne.nzer .or.
     &   nrcrfer.le.0 .or. nzcrfer.le.0 ) then
	  iflag = -56
	  go to 900
	end if
!
! 2D-Mafia produces an Er grid starting with z=0 and finite r
! Add a new mesh point at r=0
      rcrfer(1) = 0.
	nrcrfer = nrcrfer + 1
!
	do i=2,nrcrfer       
      read(mfile,*,iostat=ioc) id,rcrfer(i)
      if( ioc .ne. 0 ) then
        iflag = -55
        go to 900
      end if
	end do
!
	do i=1,nzcrfer       
      read(mfile,*,iostat=ioc) id,zcrfer(i)
      if( ioc .ne. 0 ) then
        iflag = -55
        go to 900
      end if
	end do
!
	do i=1,nzcrfer
      ercrf(1,i) = 0.      ! set Er(0,z) = 0
      do j=2,nrcrfer       
      read(mfile,*,iostat=ioc) id,id,ercrf(j,i)
      if( ioc .ne. 0 ) then
        iflag = -55
        go to 900
      end if
	ercrf(j,i) = ercrf(j,i) * fnrm
	end do
	end do
!
! .........................................................
      read(mfile,*,iostat=ioc) nrcrfbp,nzcrfbp
      if( ioc .ne. 0 ) then
        iflag = -55
        go to 900
      end if
!
      if( nrcrfbp.ne.nrbp .or. nzcrfbp.ne.nzbp .or.
     &   nrcrfbp.le.0 .or. nzcrfbp.le.0 ) then
	  iflag = -56
	  go to 900
	end if
!
! 2D-Mafia produces an Bphi grid starting with finite r and z
! Add a new mesh point at r=0
      rcrfbp(1) = 0.
	nrcrfbp = nrcrfbp + 1
!
	do i=2,nrcrfbp       
      read(mfile,*,iostat=ioc) id,rcrfbp(i)
      if( ioc .ne. 0 ) then
        iflag = -55
        go to 900
      end if
	end do
!
! Add a new mesh point at z=0
      zcrfbp(1) = 0.
	nzcrfbp = nzcrfbp + 1
!
	do i=2,nzcrfbp       
      read(mfile,*,iostat=ioc) id,zcrfbp(i)
      if( ioc .ne. 0 ) then
        iflag = -55
        go to 900
      end if
	end do
!
	do i=2,nzcrfbp
      do j=2,nrcrfbp       
      read(mfile,*,iostat=ioc) id,id,bpcrf(j,i)
      if( ioc .ne. 0 ) then
        iflag = -55
        go to 900
      end if
	bpcrf(j,i) = bpcrf(j,i) * fnrm
	end do
	end do
!
! Set Bphi(0,z) = 0
      do i=1,nzcrfbp
	bpcrf(1,i) = 0.
	end do
!
! Extrapolate Bphi(r,0) values
      do i=1,nrcrfbp
	slp = (bpcrf(i,3)-bpcrf(i,2))/(zcrfbp(3)-zcrfbp(2))
	bpcrf(i,1) = bpcrf(i,2) - slp*zcrfbp(2)
	end do
!
      close(mfile)
! .........................................................
!
900   continue
	return
	end
c  *********************************************************
      subroutine QUADRUPOLE(fpar,x,e,b)
c
c  Get field of magnetic quadrupole at position  x
c
      include 'icommon.inc'
c
!      real*8 a11,a11p,a11pp,a12,a12p,a13,a31
      real*8 x(3),e(3),b(3)
      real*4 lamb,fpar(15)
c
      do i=1,3
      b(i) = 0.
      e(i) = 0.
      end do
!
      model = fpar(1)
      g0 = fpar(2)    ! gradient strength [T/m]
!      clen = fpar(3)  ! length
!
      if( model .eq. 1) then       ! constant gradient
!        if( x(3).lt.0. .or. x(3).gt.clen ) go to 900
!
        b(1) = g0 * x(2)
        b(2) = g0 * x(1)
!
      else if( model .eq. 2 ) then        ! TANH(s) model
        clen = fpar(3)  ! length
        order = fpar(4)
        elen = fpar(5)        ! end length
        lamb = fpar(6)        ! end attenuation length
        s0 = fpar(7)          ! setupole component strength
        o0 = fpar(8)          ! octupole component strength
!
!        if( x(3).lt.0. .or. x(3).gt.clen+2*elen ) go to 900
!
        a11 = 0.
        a11p = 0.
        a11pp = 0.
        a12 = 0.
        a12p = 0.
        a13 = 0.
!
        s = x(3) - elen
        fnc = FTANH(s,clen,lamb,0)
        d1fnc = FTANH(s,clen,lamb,1)
!             
        if( order .eq. 1 ) then
          a11 = g0 * fnc
!
        else if( order .eq. 2 ) then
          a11 = g0 * fnc
          a11p = g0 * d1fnc
          a12 = s0 * fnc
!
        else if( order .eq. 3 ) then
          a11 = g0 * fnc
          a11p = g0 * d1fnc
          a12 = s0 * fnc
          a12p = s0 * d1fnc
          a13 = o0 * fnc
        end if              ! end test on order
!
        a31 = -a13 - a11pp 
!     
        b(1) = a11*x(2) + a12*x(1)*x(2) + 0.5*a13*x(1)**2*x(2)
     &    + a31/6.*x(2)**3
        b(2) = a11*x(1) + 0.5*a12*x(1)**2 - 0.5*a12*x(2)**2
     &    + 0.5*a31*x(1)*x(2)**2 + a13/6.*x(1)**3
        b(3) = a11p *x(1)*x(2) + 0.5*a12p *x(1)**2*x(2)
     &    - a12p/6.*x(2)**3
      end if                            ! end test on model
!
900   continue
      return
      end
c  *********************************************************
      subroutine ROD(fpar,xyz,e,b)
c
c  Field from current-carrying axial rod
c
      include 'icommon.inc'
c
      real*8 xyz(3),rp,e(3),b(3)
      real*4 fpar(15)
c
      do i=1,3
      b(i) = 0.
      e(i) = 0.
      end do
!
      model = fpar(1)
      brr = fpar(2)    ! Bphi at outer radius
      rr = fpar(3)     ! radius of rod [m]
      clen = fpar(4)   ! length of rod [m]
      rp = SQRT( xyz(1)**2 + xyz(2)**2 )
      if( xyz(1)*xyz(2) .ne. 0. ) then
        phi = ATAN2(xyz(2),xyz(1))
      else
        phi = 0.
      end if
!
      if( rp .le. rr ) then
        bm = brr * rp / rr
      else
        bm = brr * rr / rp
      end if
c
      if( model .eq. 1 ) then        ! central field only
!        if( xyz(3).lt.0. .or. xyz(3).gt.clen ) go to 900 
        b(1) = bm * SIN(phi)    
        b(2) = -bm * COS(phi)
      end if
c  ........................................................
      if( (model .eq. 2) .or. (model .eq. 3) ) then
c  Correct for fringe field at ends
        z = xyz(3)                ! distance into rod
        zend = fpar(5)            ! end length
!        if( z.lt.0. .or. z.gt.clen+2*zend ) go to 900 
c
        b(1) = bm * sin(phi)    
        b(2) = -bm * cos(phi)
!
        if( z .lt. zend ) then
          b(1) = b(1) * z / zend
          b(2) = b(2) * z / zend
c
        else if( z .gt. clen-zend ) then
          f = (clen - z) / zend
          b(1) = b(1) * f
          b(2) = b(2) * f
        end if
c  ........................................................
        if( model .eq. 3 ) then
c  Correct field for minimal nonlinearity
c  Use A. Lennox, FNAL pbar note #269, Fig. 5
          x = rp / rr
          if( x .lt. 0.3 ) then
            a2 = 0.13753
            a1 = 0.66775
            a0 = -0.00173
          else if( x .gt. 0.8 ) then
            a2 = -1.00000
            a1 = 2.37000
            a0 = -0.65000
          else
            a2 = 0.
            a1 = 0.78571
            a0 = -0.02548
          end if
c
          f = (a2*x*x + a1*x + a0)
          b(1) = b(1) * f
          b(2) = b(2) * f
        end if
      end if      ! end of test on model=2,3
c
      if( model .eq. 4 ) then     ! dTANH(z)
	  elen = fpar(5)
	  xlam = fpar(6)
        z = xyz(3)
	  if( z.le.0. .or. z.gt.clen+2.*elen ) go to 900
	  z = z - elen
	  bphi = bm * FTANH(z,clen,xlam,0)
        b(1) = bphi * SIN(phi)    
        b(2) = -bphi * COS(phi)
	end if
!
900   continue
      return
      end
! *********************************************************
      function SC_F1(u)
! Radial field for space charge model 3
!
      include 'icommon.inc'
!
      real*8 aux1,aux2,aux3,r,s,s3,sx
!
      r = radsc
	s = slongsc
	s3 = szbsc
	sx = srbsc
!
      aux1=sqrt(u)
      aux2=1./sqrt(2.*sx*sx+2.*u*(s3*s3-sx*sx))
      aux3=r*r/(2.*sx*sx)+s*s*0.5/(u*(s3*s3-sx*sx)+sx*sx)
      sc_f1=aux1*aux2*exp(-aux3*u)
!
	return
      end
c *********************************************************
      function SC_F2(u)
! Long field sigma_x component for space charge model 3
!
      include 'icommon.inc'
!
      real*8 aux0,aux1,aux2,aux3,r,s,s3,sx
!
      r = radsc
	s = slongsc
	s3 = szbsc
	sx = srbsc
!
      aux0=0.5/(u*(s3*s3-sx*sx)+sx*sx)
      aux1=sqrt(u*aux0)
      aux2=u*r*r-2.*sx*sx
      aux3=r*r*0.5/sx/sx+s*s*aux0
      sc_f2=aux1*aux2*exp(-aux3*u)
!
      return
      end
c *********************************************************
      function SC_F3(u)
! Long field sigma_3 component for space charge model 3
!
      include 'icommon.inc'
!
      real*8 aux0,aux1,aux2,aux3,r,s,s3,sx
!
      r = radsc
	s = slongsc
	s3 = szbsc
	sx = srbsc
!
      aux0=0.5/(u*(sx*sx-s3*s3)+s3*s3)
      aux1=sqrt(u)*aux0
      aux2=(2*s*s*u-2.*s3*s3)
      aux3=s*s*0.5/s3/s3+r*r*aux0
      sc_f3=aux1*aux2*exp(-aux3*u)
!
      return
      end
c *********************************************************
      function SC_F4(u)
! Longitudinal field for space charge model 3
!
      include 'icommon.inc'
!
      real*8 aux0,aux1,aux3,r,s,s3,sx
!
      r = radsc
	s = slongsc
	s3 = szbsc
	sx = srbsc
!
      aux1=sqrt(u)
      aux0=0.5/(s3*s3+u*(sx*sx-s3*s3))
      aux3=s*s/(2.*s3*s3)+r*r*aux0
      sc_f4=aux1*aux0*exp(-aux3*u)
!
	return
      end
c *********************************************************
      function SC_H1(u)
! Radial field for space charge model 3
!
      include 'icommon.inc'
!
      real*8 aux1,aux3,r,s,s3,sx
!
      r = radsc
	s = slongsc
	s3 = szbsc
	sx = srbsc
!
      aux1=sqrt(u)
      aux2=1.
      aux3=u*(r**2+s**2)/2./sx/sx
      sc_h1=aux1*aux2*exp(-aux3)
!
	return
      end
c *********************************************************
      function SC_H2(u)
! Long field sigma component near 0 for space charge model 3
!
      include 'icommon.inc'
!
      real*8 aux1,aux3,r,s,s3,sx
!
      r = radsc
	s = slongsc
	s3 = szbsc
	sx = srbsc
!
      aux1=sqrt(u)
      aux2=1.
      aux3=u*(r**2+s**2)/2./sx/sx
      sc_h2=aux1*aux2*exp(-aux3)
!
	return
      end
c *********************************************************
      function SC_H3(u)
! Long field sigma component for space charge model 3
!
      include 'icommon.inc'
!
      real*8 aux1,aux3,r,s,s3,sx
!
      r = radsc
	s = slongsc
	s3 = szbsc
	sx = srbsc
!
      aux1=u**(3/2)
      aux2=2*(r**2+s**2)*aux1-3*(2*sx*sx)*sqrt(u)
       aux3=u*(r**2+s**2)/(2*sx*sx)
      sc_h3=aux2*exp(-aux3)
!
	return
      end
! *********************************************************
      subroutine SC_INIT(spcharge)
!
! Initialize space-charge for this region
!
!  If we want space charge, set  spcharge = true
!
      include 'icommon.inc'
!
      logical*4 spcharge
!
      if( (spacelev.eq.1 .and. ftag(1).eq.'ACCE')
     & .or. (spacelev.ge.2) ) then
        spcharge = .true.
      end if
!
      return
      end
! *********************************************************
      subroutine SC_KICK(loc)
!
! Give bunch a space charge kick
!
      include 'icommon.inc'
!
      real*8 stem(2,2),pm,tavg,term
      real*4 esc(3),bsc(3),ddsc(3)
!
      do i=1,2
      do j=1,2
      stem(i,j) = 0.
      end do
      end do
      stev = 0.
!
      do 520 ip=1,npart     ! accumulate time statistics
      call READ_EVT(ip)
      if( ipflg .ne. 0 ) go to 520
      stev = stev + 1.  
      stem(1,1) = stem(1,1) + tp
520   continue
!
      if( stev .gt. 0. ) tavg = stem(1,1) / stev
!
      stem(1,1) = 0.
!
      do 540 ip=1,npart     ! accumulate bunch statistics
      call READ_EVT(ip)
      if( ipflg .ne. 0 ) go to 540
! Transform z-plane values into spatial bunch
      dt = tp - tavg
      call VMAG(pp,pm)
      gamma = SQRT( 1. + (pm/pmass)**2 )
      term = cc * dt / pmass / gamma
      do i=1,3
      ddsc(i) = pp(i) * term
      xp(i) = xp(i) + ddsc(i)
      end do
      call WRITE_EVT(ip)
!
! Accumulate length and radius statistics
      stem(1,1) = stem(1,1) + ddsc(3)
      stem(2,1) = stem(2,1) + ddsc(3) * ddsc(3)
      r = SQRT( xp(1)**2 + xp(2)**2 )
      stem(1,2) = stem(1,2) + r
      stem(2,2) = stem(2,2) + r * r
540   continue
!
      if( stev .gt. 0. ) then
        zavg = stem(1,1) / stev
        szbsc = SQRT( stem(2,1)/stev - zavg*zavg )
        ravg = stem(1,2) / stev
        srbsc = SQRT( stem(2,2)/stev - ravg*ravg )
      end if
!
! Apply space charge kick to all the particles
      do 560 ip=1,npart 
      call READ_EVT(ip) 
      if( ipflg .ne. 0 ) go to 560
      call VMAG(pp,pm)
      gamma = SQRT( 1. + (pm/pmass)**2 )
      radsc = SQRT( xp(1)**2 + xp(2)**2 )
      dt = tp - tavg
      term = cc * dt / pmass / gamma
	slongsc = pp(3) * term
!
      call SPACE_CHARGE(loc,gamma,esc,bsc)
!
      do i=1,3
      xp(i) = xp(i) - pp(i)*term   ! back to z-plane values
      end do
!
      call PKICK(gamma,esc,bsc)     ! update momentum
!
      call WRITE_EVT(ip)
!
560   continue
!
      return
      end
c  *********************************************************
      subroutine SEXTUPOLE(fpar,x,e,b)
c
c  Get field of sextupole magnet
c
      include 'icommon.inc'
c
!      real*8 a12,a12p,a13
      real*8 x(3),e(3),b(3)
      real*4 lamb,fpar(15)
c
      do i=1,3
      b(i) = 0.
      e(i) = 0.
      end do
!
      model = fpar(1)
      s0 = fpar(2)      ! sextupole strength [T/m^2]
!
      if( model .eq. 1) then        ! constant strength
!        if( x(3).lt.0. .or. x(3).gt.clen ) go to 900
!
        b(1) = s0 * x(1) * x(2)
        b(2) = 0.5 * s0 * (x(1)**2 - x(2)**2)
!
      else if( model .eq. 2 ) then      ! TANH(s) model
        clen = fpar(3)    ! central length
        order = fpar(4)
        elen = fpar(5)        ! end length
        lamb = fpar(6)        ! end attenuation length
        o0 = fpar(7)          ! octupole component strength
!
!        if( x(3).lt.0. .or. x(3).gt.clen+2*elen ) go to 900
!
        a12 = 0.
        a12p = 0.
        a13 = 0.
!     
        s = x(3) - elen
        fnc = FTANH(s,clen,lamb,0)
        d1fnc = FTANH(s,clen,lamb,1)
!              
        if( order .eq. 2 ) then
          a12 = s0 * fnc
!
        else if( order .eq. 3 ) then
          a12 = s0 * fnc
          a12p = s0 * d1fnc
          a13 = o0 * fnc
        end if              ! end test on order
!
        b(1) = a12*x(1)*x(2) + 0.5*a13*x(1)**2*x(2)
     &    - a13/6.*x(2)**3
        b(2) = 0.5*a12*x(1)**2 - 0.5*a12*x(2)**2
     &    - 0.5*a13*x(1)*x(2)**2 + a13/6.*x(1)**3
        b(3) = 0.5*a12p *x(1)**2*x(2) - a12p/6.*x(2)**3
!
      end if                ! end test on model
c
900   continue
      return
      end

!  *********************************************************
      subroutine SHEET(fpar,xyz,e,b)
!
!  Get field from solenoidal current sheet
!
      include 'icommon.inc'
c
      save jsold,jcold
      real*8 xyz(3),e(3),b(3)
      real*4 fpar(15)
      character*10 fname
      data bnorm/4.e-7/
      data jcold/-1/,jsold/-1/
c
      model = fpar(1)
!
      do i=1,3
      e(i) = 0.
      b(i) = 0.
      end do
!
      rp = SQRT( xyz(1)**2 + xyz(2)**2 )
      if( xyz(1)*xyz(2) .ne. 0. ) then
        phi = ATAN2( xyz(2),xyz(1) )
      else
        phi = 0.
      end if
      sphi = sin(phi)
      cphi = cos(phi)
! ---------------------------------------------------------
      if( model .eq. 1 ) then    ! exact field from single sheet
        a = fpar(3)              ! radius of sheet
        clen = fpar(4)           ! length of sheet
        curr = fpar(5)           ! current
        z = xyz(3) - fpar(2)
!
        call BSHEET(-z,clen,rp,a,bz,br)
!
        b(1) = br * cphi* curr
        b(2) = br * sphi * curr
        b(3) = bz * curr
!
        do i=1,3
        b(i) = b(i) * bnorm
        end do
      end if
!  --------------------------------------------------------
!  Initialize for models 2 and 3
!
      if( (kbcell.eq.0.and.jsrg.ne.jsold) .or. 
     &      (kbcell.eq.1.and.jcel.ne.jcold) ) then
!
        if( model.eq.2 .or. model.eq.3 ) then   ! need data file
          mfile = fpar(2)
          write(fname,40) mfile
40        format('for0',i2,'.dat')
          open (
     &     unit=mfile,
     &     file=pathName(1:Len_Trim(pathName)) // fname,
     &     status='unknown',
     &     iostat=ioc)
          if( ioc .ne. 0 ) then
            iflag = -34
            go to 900
          end if
!
          write(2,*)' SHEET SUMMARY for file',mfile
          read(mfile,*,iostat=ioc)nsheets
          if( ioc .ne. 0 ) then
            iflag = -35
            go to 900
          end if
          write(2,'(a,i5)') ' NSHEETS = ',nsheets
!
          do i=1,nsheets
          read(mfile,*,iostat=ioc)idm,dzelgr(i),zlelgr(i),
     &      relgr(i),curelgr(i)
          if( ioc .ne. 0 ) then
            iflag = -35
            go to 900
          end if
          write(2,50) idm,dzelgr(i),zlelgr(i),
     &      relgr(i),curelgr(i)
50        format(' #,DZ,L,R,I: ',i5,4e12.4)
          end do
!
          close(mfile)
!         
        end if            ! end of test on model = 2 or 3
!  ........................................................
        if( model.eq.3 ) then              ! need grid set up
          delzgr = fpar(3)
          delrgr = fpar(4)
          zglen = fpar(5)
          rglen = fpar(6)
          dzcut = fpar(7)
          mfile = fpar(9)
          if( delzgr .le. 0. .or. delrgr .le. 0. ) then
            iflag = -41
            go to 900
          end if
          nzgr = ANINT( zglen / delzgr + 1. )
          nrgr = ANINT( rglen / delrgr + 1. )
          if( nzgr.gt.mxzgr .or. nrgr.gt.mxrgr .or. krad.ne.1) then
            iflag = -42
            go to 900
          end if
!  ........................................................ 
          do iz=1,nzgr
          if( MOD(iz,20) .eq. 1 .and. termout ) then
            write(*,80)' Computing grid:',iz,' /',nzgr,' done'
80          format(a,i5,a,i5,a)
          end if
          zgr(iz) = (iz-1) * delzgr              
!       
          do ir=1,nrgr
          rgr(ir) = (ir-1) * delrgr
!
! Check that the bozo isn't asking for a grid point on a sheet
          do i=1,nsheets
          zend = dzelgr(i) + zlelgr(i)
          if( zgr(iz).ge.dzelgr(i) .and. zgr(iz).le.zend
     &       .and. rgr(ir).eq.relgr(i) ) then
            write(2,*) ' Grid point requested on a current sheet'
            if (termout) 
     &       write(*,*) ' Grid point requested on a current sheet'
            jcold = jcel
            jsold = jsrg
            iflag = -1
            go to 900
          end if
          end do
!
          bzgr(iz,ir) = 0.
          brgr(iz,ir) = 0.
!
          do 100 ic=1,nsheets
          clen = zlelgr(ic)
          a = relgr(ic)
          z = dzelgr(ic) - zgr(iz) 
          if( ABS(z) .gt. dzcut ) go to 100              
!
          call BSHEET( z,clen,rgr(ir),a,bz,br )
!
          bzgr(iz,ir) = bzgr(iz,ir) + bz * curelgr(ic)
          brgr(iz,ir) = brgr(iz,ir) + br * curelgr(ic)
100       continue              ! sum on sheets
!
          bzgr(iz,ir) = bzgr(iz,ir) * bnorm
          brgr(iz,ir) = brgr(iz,ir) * bnorm
          end do              ! r grid
!
          end do              ! z grid
!
          if( mfile .gt. 19 ) then    ! save field grid to file
            write(fname,40) mfile
            open (
     &       unit=mfile,
     &       file=pathName(1:Len_Trim(pathName)) // fname,
     &       status='unknown',
     &       iostat=ioc)
            if( ioc .ne. 0 ) then
              iflag = -50
              go to 900
            end if
            write(mfile,'(i8)') nzgr*nrgr
            do i=1,nzgr
            do j=1,nrgr
            write(mfile,'(2i6,4e12.4)') i,j,zgr(i),rgr(j),
     &        bzgr(i,j),brgr(i,j)
            end do
            end do
            close(mfile)
          end if
!
        end if           ! end of test on model = 3
!
        jcold = jcel
        jsold = jsrg
      end if                  ! end of test on kbcell ...
! End of initialization
!
!  --------------------------------------------------------
      if( model .eq. 2 ) then
! Exact field from sum of sheets in data file
!
        do i=1,nsheets
        clen = zlelgr(i)
        a = relgr(i)
        z = xyz(3) - dzelgr(i)
!
        call BSHEET( -z,clen,rp,a,bz,br )
!
        b(1) = b(1) + br * cphi* curelgr(i)
        b(2) = b(2) + br * sphi * curelgr(i)
        b(3) = b(3) + bz * curelgr(i)
        end do
c
        do i=1,3
        b(i) = b(i) * bnorm
        end do
!
      end if
!  --------------------------------------------------------
      if( model .eq. 3 ) then
! Interpolated grid field from sum of sheets in data file
        z = xyz(3)               
        call HUNT(zgr,nzgr,z,jzlo)
!  Watch out for initial point exactly on left edge of grid
        if( jzlo .eq. 0 ) jzlo = 1                                          
        call HUNT(rgr,nrgr,rp,jrlo)
        if( jrlo .eq. 0 ) jrlo = 1                                          
        if( z .ge. zgr(nzgr) ) go to 900     ! B=0 beyond grid
        if( rp .ge. rgr(nrgr) ) go to 900
!            
        interp = fpar(8)      ! interpolation model
!
        if( interp .le. 1 ) then   ! bi-linear interpolation
         u = (z - zgr(jzlo)) / delzgr          
         v = (rp - rgr(jrlo)) / delrgr
         br = (1.-u)*(1.-v)*brgr(jzlo,jrlo)+ u*(1.-v)*brgr(jzlo+1,jrlo)
     &    + u*v*brgr(jzlo+1,jrlo+1) + (1.-u)*v*brgr(jzlo,jrlo+1)
         bz = (1.-u)*(1.-v)*bzgr(jzlo,jrlo)+ u*(1.-v)*bzgr(jzlo+1,jrlo)
     &    + u*v*bzgr(jzlo+1,jrlo+1) + (1.-u)*v*bzgr(jzlo,jrlo+1)
!
        else if( interp.eq.2 .or. interp.eq.3 ) then  
! Polynomial interpolation
!
         call BIPOLYINT(jzlo,jrlo,zgr,rgr,bzgr,mxzgr,mxrgr,interp,
     &     z,rp,bz,df)     
         call BIPOLYINT(jzlo,jrlo,zgr,rgr,brgr,mxzgr,mxrgr,interp,
     &     z,rp,br,df)     
!
        end if        ! end of test on interpolation
!
        b(1) = br * cphi
        b(2) = br * sphi
        b(3) = bz
      end if          ! end of test on model
c
900   continue
      return
      end

c  *********************************************************
      subroutine SOLENOID(fpar,xyz,e,b)
c
c  Get field from solenoid
c
      include 'icommon.inc'
c
      real*8 xyz(3),e(3),b(3)
      real*4 mu,lamb,fpar(15)
      data bnorm/2.e-7/
      data mu/12.566371e-7/
!
      do i=1,3
      e(i) = 0.
      b(i) = 0.
      end do
!
      model = fpar(1)
      b0 = fpar(2)              ! field at center [T]
      clen = fpar(3)            ! central length [m]
      elen = fpar(4)            ! end length [m]
!
!      if( xyz(3).lt.0. .or. xyz(3).gt.clen+2*elen ) go to 900
!
      rp = sqrt( xyz(1)**2 + xyz(2)**2 )
!
      if( xyz(1)*xyz(2) .eq. 0.) then
        phi = 0.
        sphi = 0.
        cphi = 1.
      else
        phi = ATAN2( xyz(2),xyz(1) )
        sphi = sin(phi)
        cphi = cos(phi)
      end if
!
c  ........................................................
      if( model .eq. 1) then    ! linear ends + constant center Bz
        bcen = fpar(5)
        z = xyz(3)
! 
        if( elen .gt. 0. )then
          g = b0 / elen
        else
          g = 0.
        end if
!
        if( z .lt. elen ) then    ! entrance fringe
          br = -0.5 * g * rp
          b(1) = br * cphi
          b(2) = br * sphi
          b(3) = g * z + bcen
        else if( z .gt. (elen+clen) ) then    ! exit fringe
          br = 0.5 * g * rp
          b(1) = br * cphi
          b(2) = br * sphi
          b(3) = b0 - g * (z - elen - clen) + bcen
        else    ! central region
          b(1) = 0.
          b(2) = 0.
          b(3) = b0 + bcen
        end if
c  ........................................................
      else if( model .eq. 2 ) then    ! hyperbolic tangent Bz
        order = fpar(5)
        lamb = fpar(6)
        bcen = fpar(7)
        z = xyz(3) - elen
        bz0 = bcen + b0 * FTANH(z,clen,lamb,0)
        dbz1 = b0 * FTANH(z,clen,lamb,1)
        br = -rp /2. * dbz1
        bz = bz0
!
        if( order .gt. 0 ) then
          dbz2 = b0 * FTANH(z,clen,lamb,2)
          dbz3 = b0 * FTANH(z,clen,lamb,3)
          br = br + rp**3 / 16. * dbz3
          bz = bz - rp**2 / 4. * dbz2
        end if
!
        if( order .gt. 1 ) then
          dbz4 = b0 * FTANH(z,clen,lamb,4)
          dbz5 = b0 * FTANH(z,clen,lamb,5)
          br = br - rp**5 / 384. * dbz5
          bz = bz + rp**4 / 64. * dbz4
        end if
!
        if( order .gt. 2 ) then
          dbz6 = b0 * FTANH(z,clen,lamb,6)
          dbz7 = b0 * FTANH(z,clen,lamb,7)
          br = br + rp**7 / 18432. * dbz7
          bz = bz - rp**6 / 2304. * dbz6
        end if
!
        b(1) = br * cphi
        b(2) = br * sphi
        b(3) = bz
!  ........................................................
      else if( model .eq. 3 ) then              ! sum of circular coils
!
        ncoils = fpar(5)
        a = fpar(6)
!
! Determine current that gives B0 in the center of solenoid
!
        if( ncoils .eq. 1 ) then
          cur = 2. * a * b0 / mu
!
        else
          dzc = clen / (ncoils-1)
          off = -clen/2.- dzc
!
          do i=1,ncoils               
          z = off + dzc*i
          bz = a*a / (a*a + z*z)**(1.5)
          b(3) = b(3) + bz
          end do
c
          b(3) = b(3) * mu / 2.
          cur = b0 / b(3) 
          b(3) = 0.
        end if
!
!  Compute field at particle location
!
        if( ncoils .eq. 1 ) then
          z = xyz(3) - elen
!       
          call BLOOP(a,z,rp,bz,br)
!
          b(1) = br * cphi
          b(2) = br * sphi
          b(3) = bz
!
        else                            
          off = -elen + dzc
!
          do i=1,ncoils               
          z = xyz(3) + off - dzc*i
!
          call BLOOP( a,z,rp,bz,br )
!
          b(1) = b(1) + br * cphi
          b(2) = b(2) + br * sphi 
          b(3) = b(3) + bz 
          end do
        end if
c
        do i=1,3
        b(i) = b(i) * cur * bnorm
        end do
!  ........................................................
      else if( model .eq. 4 ) then    ! exact field from single sheet
        a = fpar(5)              ! radius of sheet
        z = xyz(3) - elen
!
! Determine current density that gives B0 in center of solenoid
        cur = b0/mu/clen * SQRT(clen*clen+4.*a*a) 
!
!  Compute field at particle location
!
        call BSHEET(-z,clen,rp,a,bz,br)
!
        b(1) = br * cphi* cur * 2.*bnorm
        b(2) = br * sphi * cur * 2.*bnorm
        b(3) = bz * cur * 2.*bnorm
!
      end if      ! end of test on models
!  --------------------------------------------------------
900   continue
      return
      end
! *********************************************************
      subroutine SPACE_CHARGE(loc,gamma,e,b)
!
!  Compute space-charge field at given location in a bunch
!
!     SZBSC,SRBSC: sigma length and radius of bunch
!     SLONGSC,RADSC: axial and radial position of this particle
!
      include 'icommon.inc'
!
      external sc_f1,sc_f2,sc_f3,sc_f4,sc_h1,sc_h2,sc_h3
      real*4 zero(20)
      real*4 e(3),b(3)    
!
      data r2pi/2.5066283/     ! SQRT(2 pi)
      data (zero(n),n=1,20)/2.4048255577,5.5200781103,8.6537279129,
     +  11.7915344391,14.9309177086,18.0710639679,21.2116366299,
     +  24.3524715308,27.4934791320,30.6346064684,33.7758202136,
     +  36.9170983537,40.0584257646,43.1997917132,46.3411883717,
     +  49.4826098974,52.6240518411,55.7655107550,58.9069839261,
     +  62.0484691902/
      data  epsi_0/8.8542d-12/,q/1.6022d-19/
!
      gamma2 = gamma**2
      beta = SQRT( 1. - (1./gamma2) )
      if( xp(1)*xp(2) .ne. 0. ) then
        phi = ATAN2( xp(2),xp(1) )
      else
        phi = 0.
      end if
      cphi = COS( phi )
      sphi = SIN( phi )
!  spcon = 1 / (4 pi eps0) * e * (# mu/bunch)
      spcon = 1.438e-9 * parbunsc
! ---------------------------------------------------------
      if( spacelev .eq. 1 ) then      
! Static approximation for bunch in cylindrical can
! Uniform transverse and longitudinal distributions 
!   J.C. Gallardo, 24 April 1998
!
        r = radsc
        s = slongsc
        xlb = 4. * szbsc    ! total uniform length
        rb = 2. * srbsc     ! total uniform radius
!  XL1 is distance from back of bunch to start of region
        if( loc .eq. 1 ) then
          xl1 = 0.
        else
          xl1 = slen - xlb 
        end if
        xlc = slen        ! length of cavity
        rc = rhigh(1)     ! radius of cavity
!
        ak=0
        aak=0
c  
        do 20 k=1,20
        u1=k*pi/xlc
        zaux=cos(u1*xl1)-cos(u1*(xl1+xlb))
        yaux=zaux*sin(u1*s)/k
        yyaux=zaux*cos(u1*s)        
        an=0
        aan=0
!
        do 15 n=1,20
        xnum=zero(n)**2+(u1*rc)**2
        xaux=bessj1(zero(n))**2
        u2=zero(n)/rc
        a=bessj1(u2*rb)*bessj1(u2*r)
        aa=bessj1(u2*rb)*bessj0(u2*r)
        an=an+a/xaux/xnum
        aan=aan+aa/xaux/xnum/zero(n)
15      end do      ! end of loop on  n
!
        ak=ak+an*yaux
        aak=aak+aan*yyaux
20      end do      ! end of loop on  k
!
        fac = spcon / xlb
        er = fac * 16/pi/rb*ak 
        bphi = beta / cc * er
        e(1) = er * cphi
        e(2) = er * sphi
        e(3) = -fac * 16*rc/rb/xlc/gamma2*aak 
        b(1) = -bphi * sphi
        b(2) = bphi * cphi
        b(3) = 0.
! .........................................................
      else if( spacelev .eq. 2 ) then
! Radial space charge from bunch in free space
!  Gaussian transverse and longitudinal distributions
!  Ref. Weidemann, V.2, p. 333
!
        z = slongsc 
        u = z**2 / (2.*szbsc*szbsc)
        xlz = 1./szbsc/r2pi * EXP( -u )
!
	  er = spcon * xlz / srbsc / srbsc
	  e(1) = er * xp(1)
	  e(2) = er * xp(2)
        e(3) = 0.  
!
        bphi = beta * er / cc 
	  b(1) = -bphi * xp(2)
	  b(2) = bphi * xp(1)
	  b(3) = 0.
! .........................................................
      else if( spacelev .eq. 3 ) then
!  Radial and longitudinal space charge from bunch in free space
!  Gaussian transverse and longitudinal distributions
!  Ref. J. Gallardo, Mucool Note #38, rev. 1, 5 October 1999
!
        dl = 0.
        du = 1.
        r = radsc
	  s = slongsc
	  s3 = szbsc
	  sx = srbsc
!
        qe=q/epsi_0  ! C/....
        coeff10= qe/sqrt(pi)*4.*parbunsc
!          
        if( sx .ne. s3 ) then
          coeff1=coeff10/(2.*sx*sx)
          coeff2=coeff10/sqrt(2.*s3*s3)
          call qromb(sc_f1,dl,du,result1)
c          call qromb(sc_f2,dl,du,result2)
c          call qromb(sc_f3,dl,du,result3)
          call qromb(sc_f4,dl,du,result4)
          er = result1*r*coeff1
c         e1=-(result2*dsxds/(2.*sx*sx)**2+
c     +      result3*ds3ds/(sqrt(2.*s3*s3))**3)*coeff10
          e2=1/gamma**2*result4*coeff2*s
c         field_s(ir,is)=e1+e2
          e(3) = e2
!
        else
          coeff1=coeff10/(sqrt(2*sx*sx))**3
          coeff2=0.5*coeff1/(2*sx*sx)
c              dl=0.
c              dini=0. 
          call qromb(sc_h1,dl,du,result1)
          er = result1*r*coeff1
          call qromb(sc_h3,dl,du,result3)
!          e1=-(result3*dsxds)*coeff2
          e2=1/gamma**2*result1*coeff1*s
!          e(3) = e1+e2
          e(3) = e2
        endif
!
	  e(1) = er * cphi
	  e(2) = er * sphi
        bphi = beta * er / cc 
        b(1) = -bphi * sphi
        b(2) = bphi * cphi
	  b(3) = 0.
      end if    ! end of test on spacelev
! ---------------------------------------------------------
!
      return
      end
c  *********************************************************
      subroutine TAPERED_ROD(fpar,xyz,e,b)
c
c  Get field from tapered current carrying rod
c
      include 'icommon.inc'
c
      real*8 xyz(3),e(3),b(3)
      real*4 fpar(15)
c
      do i=1,3
      b(i) = 0.
      e(i) = 0.
      end do
!
      bc = fpar(2)    ! central field [T]
      rc = fpar(3)              ! central radius [m]
      zc = fpar(4)              ! central length
      r1 = fpar(5)              ! entrance radius [m]
      z1 = fpar(6)              ! entrance length
      r2 = fpar(7)              ! exit radius [m]
      z2 = fpar(8)              ! exit length
!
      z = xyz(3) 
!      if( z.lt.0. .or. z.gt.z1+zc+z2 ) go to 900
! 
      rp = sqrt( xyz(1)**2 + xyz(2)**2 )
!
      if( z .lt. z1 ) then    ! entrance taper
        rr = r1 + (rc-r1)/z1 * z
      else if( z .gt. (z1+zc) ) then    ! exit taper
        rr = rc + (r2-rc)/z2 * (z - z1 - zc)
      else    ! central region
        rr = rc
      end if
!
      if( rp .ge. rr ) then
        iflag = -40
        go to 900
      end if
!
      if( xyz(1)*xyz(2) .eq. 0.) then              ! particle on-axis              
        b(1) = 0.
        b(2) = 0.
        go to 900
!
      else
        phi = ATAN2( xyz(2),xyz(1) )
        sphi = sin(phi)
        cphi = cos(phi)
        bm = bc * (rc / rr) * (rp / rr)
        b(1) = bm * sphi
        b(2) = -bm * cphi
      end if
c  ........................................................
900   continue
      return
      end

c  *********************************************************
      subroutine TAPSOL(fpar,xyz,e,b)
c
c  Get field from tapered solenoid
c
      include 'icommon.inc'
c
      real*8 xyz(3),e(3),b(3)
      real*4 fpar(15)
c
      do i=1,3
      b(i) = 0.
      e(i) = 0.
      end do
!
      bc = fpar(2)              ! central field [T]
      rc = fpar(3)              ! central radius [m]
      zc = fpar(4)              ! central length
      b1 = fpar(5)              ! entrance field [T]
      r1 = fpar(6)              ! entrance radius [m]
      z1 = fpar(7)              ! entrance length
      b2 = fpar(8)              ! exit field [T]
      r2 = fpar(9)              ! exit radius [m]
      z2 = fpar(10)             ! exit length
!
!      if( xyz(3).lt.0. .or. xyz(3).gt.z1+zc+z2 ) go to 900
!
      if( z1 .gt. 0. ) then
        g1 = (bc-b1) / z1
      end if
      if( z2 .gt. 0. ) then
        g2 = (b2-bc) / z2
      end if
!
      z = xyz(3) 
      rp = sqrt( xyz(1)**2 + xyz(2)**2 )
!  ........................................................ 
      if( xyz(1)*xyz(2) .eq. 0.) then              ! particle on-axis
        b(1) = 0.
        b(2) = 0.
!
      if( z .lt. z1 ) then    ! entrance fringe
        b(3) = b1 + g1 * z
        rs = r1 + (rc-r1)/z1 * z
      else if( z .gt. (z1+zc) ) then    ! exit fringe
        b(3) = bc + g2 * (z - z1 - zc)
        rs = rc + (r2-rc)/z2 * (z - z1 - zc)
      else    ! central region
        b(3) = bc
        rs = rc
      end if
        if( rp .ge. rs ) then
          iflag = -31
        end if
        go to 900
!
      else
        phi = ATAN2( xyz(2),xyz(1) )
        sphi = sin(phi)
        cphi = cos(phi)
      end if
c  ........................................................
!  General case
      if( z .lt. z1 ) then    ! entrance fringe
        br = -0.5 * g1 * rp
        b(1) = br * cphi
        b(2) = br * sphi
        b(3) = b1 + g1 * z
        rs = r1 + (rc-r1)/z1 * z
c
      else if( z .gt. (z1+zc) ) then    ! exit fringe
        br = -0.5 * g2 * rp
        b(1) = br * cphi
        b(2) = br * sphi
        b(3) = bc + g2 * (z - z1 - zc)
        rs = rc + (r2-rc)/z2 * (z - z1 - zc)
c
      else    ! central region
        b(1) = 0.
        b(2) = 0.
        b(3) = bc
        rs = rc
      end if
c
      if( rp .ge. rs ) then
        iflag = -31
      end if

900   continue
      return
      end
