program CAVEL
!
! Examine electron distribution inside RF cavity and external B field
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
!
      version = '1.17'
!
	call SETUP
	call SIMULATE
!	call DISPLAY
!
      end
! *********************************************************
      BLOCK DATA
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
!
! cont
      data bndname/'outaut.txt'/
      data cc/2.99792458d8/
      data mel/0.51099892d0/        ! [MeV]
      data qel/299.792458d0/
! parm
      data cscale/1d0/,launch/1/,magdim/2/,rfnrm/1d0/
      data sigp/0./,ntemp/1/,nphi/1/,offlaun/1d-6/
      data rtuple/.false./,secem/0/,rskip/1/,ascrf/0/,nfix/10/
      data step/1d-10/,epsstep/1d-5/,maxsteps/2000/
      data rnseed/-1/,robth/0d0/,robphi/0d0/
      data seemax/160e-6/,sedelm/2.0/,seeps0/150e-6/,sefang/0.45/
      data stepmin/1d-12/,stepmax/1d-9/
      data zoff/0d0/,otherside/.false./
!
      end
!  *********************************************************
      subroutine BILINEAR(u,v,jx,jy,fg,mx,my,f)
!
! Make bilinear interpolation on grid
!
      implicit double precision (a-h,o-z)
      dimension fg(mx,my)
!
      f = (1.-u)*(1.-v)*fg(jx,jy)   &
            + u*(1.-v)*fg(jx+1,jy)  &
            + u*v*fg(jx+1,jy+1)     &
            + (1.-u)*v*fg(jx,jy+1)
!
      end
! *********************************************************
      subroutine CAVITY_BOUNDARY
!
! Extract cavity boundary from SuperFish file OUTAUT.TXT
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
!
      dimension tbx(2000),tby(2000)
      logical boundary,trigger,tit
      character(80) line
      character(16) chx,chy
!
      open(unit=11,file=bndname,status='old',iostat=ioc)
      if( ioc .ne. 0 ) then
        write(*,*) 'Cant open file ',bndname
        write(2,*) 'Cant open file ',bndname
        stop
      end if
      ncbpt = 1         ! save room for point with x=0
      boundary = .false.
      trigger = .false.
      tit = .false.
!
100   continue
      read(11,'(a80)',iostat=ioc) line
      if( ioc .ne. 0 ) then
        write(*,*) 'Error reading SuperFish boundary'
        write(2,*) 'Error reading SuperFish boundary'
        stop
      end if
      if( line(1:13) .eq. 'Region 1 done' ) go to 200
      if( tit ) then
        write(2,'(a80)') line
        tit = .false.
      end if
      if( line(1:19) .eq. 'Problem description' ) tit = .true.
      if( line(1:11) .eq. 'XMAX      =' ) then
        chx = line(13:20)
        xmax = FIXFMT( chx )
      end if
      if( line(1:11) .eq. 'YMAX      =' ) then
        chy = line(13:20)
        ymax = FIXFMT( chy )
      end if
!
      if( boundary ) then
        chx = line(16:33)
        chy = line(34:50)
        tmpx = FIXFMT( chx )
        tmpy = FIXFMT( chy )
        if( (tmpx.ne.0d0.and.tmpy.ne.0d0) .or. &
         (ascrf.ne.0.and.tmpy.ne.0d0) ) then      
          ncbpt = ncbpt + 1
          tbx(ncbpt) = tmpx
          tby(ncbpt) = tmpy
        end if
        go to 100
      end if
      if( trigger ) boundary = .true.     ! skip column labels
      if( line(1:20) .eq. 'Region 1 mesh points' ) then
        trigger = .true.
      end if
      go to 100
!
200   continue
      if( ncbpt .le. 0 ) then
        write(*,*) 'Error finding boundary points'
        write(2,*) 'Error finding boundary points'
        stop
      end if
      close(11)
!
      n = ncbpt
! If SuperFish file is L-R symmetric, create 2nd half of boundary
! Convert cm to m
      if( ascrf .eq. 0 ) then       ! symmetric       
        do i=2,n
        cbx(i) = -tbx(n+2-i) * 1d-2
        cby(i) = tby(n+2-i) * 1d-2
        cbx(n+i) = tbx(i) * 1d-2
        cby(n+i) = tby(i) * 1d-2
        end do
! Add in missing point at top
        cbx(n+1) = 0d0
        cby(n+1) = ymax * 1d-2
        ncbpt = 2*ncbpt 
!
      else                    ! not symmetric
        do i=2,n
        cbx(i) = tbx(i) * 1d-2
        cby(i) = tby(i) * 1d-2
        end do
      end if
!
! Add in horizontal end points
      cbx(1) = cbx(2)
      cby(1) = 0d0
      ncbpt = ncbpt + 1
      cbx(ncbpt) = cbx(ncbpt-1)
      cby(ncbpt) = 0d0
!
      xmax = xmax * 1d-2
      ymax = ymax * 1d-2
      write(2,'(a,f12.4,4x,a,f12.4)') 'xmax = ',xmax,' ymax = ',ymax
!
      write(2,'(a)')'Cavity boundary'
      do i=1,ncbpt
      write(2,251) i,cbx(i),cby(i)
251   format(i5,2f12.6)
      end do
!
! Determine the order the user entered the boundary points
      th1 = ATAN2(cby(ncbpt/10),cbx(ncbpt/10))
      th2 = ATAN2(cby(ncbpt/2),cbx(ncbpt/2))
      if( th2 .gt. th1 ) then
        ccw = .true.
      else
        ccw = .false.
      end if
!
      end
! *********************************************************
      subroutine CAVITY_FIELDS
!
! Get cavity RF fields from SuperFish Parmela file  name.T7
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
!
      open(unit=12,file=rfname,status='old',iostat=ioc)
      if( ioc .ne. 0 ) then
        write(*,*) 'Cant open file',rfname
        write(2,*) 'Cant open file',rfname
        stop
      end if
!
      read(12,*,iostat=ioc) zmin,zmax,nzint
      if( ioc .ne. 0 ) then
        write(2,*)'Error reading Superfish z data'
        stop
      end if
      if( nzint.le.0 .or. nzint.gt.2000 ) then
        write(2,*)'Error in Superfish: nzint'
        stop
      end if
!
      read(12,*,iostat=ioc) freq
      if( ioc .ne. 0 ) then
        write(2,*)'Error reading Superfish frequency'
        stop
      end if
      freq = freq * 1d6                         ! [Hz]
!
      read(12,*,iostat=ioc) rmin,rmax,nrint
      if( ioc .ne. 0 ) then
        write(2,*)'Error reading Superfish r data'
        stop
      end if
      if( nrint.le.0 .or. nrint.gt.2000 ) then
        write(2,*)'Error in Superfish: nrint'
        stop
      end if
!
      dzcrf = (zmax-zmin)*1d-2 / nzint          ! [m]
      drcrf = (rmax-rmin)*1d-2 / nrint          ! [m]
      nzcrf = nzint + 1
      nrcrf = nrint + 1
!
      do i=1,nzcrf
      zcrf(i) = zmin*1d-2 + (i-1)*dzcrf
      end do
!
      do i=1,nrcrf
      rcrf(i) = (i-1)*drcrf
      end do
      nrec = (nzint+1) * (nrint+1)
!
      do ir=1,nrcrf
      do iz=1,nzcrf
      read(12,*,iostat=ioc) ezcrf(ir,iz),ercrf(ir,iz),  &
        dum,bpcrf(ir,iz)
      if( ioc .ne. 0 ) then
        write(2,*)'Error reading Superfish field: iz,ir',iz,ir
        stop
      end if
      ezcrf(ir,iz) = ezcrf(ir,iz) * rfnrm * 1d6        ! [V/m]
      ercrf(ir,iz) = ercrf(ir,iz) * rfnrm * 1d6        ! [V/m]
      bpcrf(ir,iz) = -bpcrf(ir,iz) * rfnrm * 4d-7*pi   ! [T]
      end do
      end do
!
!      write(2,'(3e12.4)') zcrf(2),rcrf(2),ezcrf(2,2)
      close(12)
!
      end
!  *********************************************************
      subroutine DERIV(t,y,dydt)
!
!  Compute values of derivatives dYdt for the variables Y at t 
!
      implicit double precision (a-h,o-z)
      dimension y(6),dydt(6)
      include 'cavel.inc'
!
! Use time stepping
!   Y = {x, y, z, Px, Py, Pz} 
!
      call FIELD( y,t )
      gamma = SQRT(y(4)**2+y(5)**2+y(6)**2 + mel**2) / mel
      f = mel * gamma
      g = cc / f
      h = qel * cc
!
      dydt(1) = y(4) * g      
      dydt(2) = y(5) * g     
      dydt(3) = y(6) * g     
!
      dydt(4) = h * (efld(1)/cc + y(5)*bfld(3)/f - y(6)*bfld(2)/f)
      dydt(5) = h * (efld(2)/cc - y(4)*bfld(3)/f + y(6)*bfld(1)/f)
      dydt(6) = h * (efld(3)/cc + y(4)*bfld(2)/f - y(5)*bfld(1)/f) 
!
      end
! *********************************************************
!	subroutine DISPLAY
!	end
! *********************************************************
      subroutine DIRECTION
!
! Find initial particle direction at this location
! Assume launch takes place in y-z plane
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
!
      ppn(1) = 0d0
!
      if( launch .eq. 1) then       ! along peak E field
! Assume electric field is normal to the surface everywhere
        ppn(2) = -ABS(ercrf(iremx,izemx)) / emx * psurf
        if( zlaun .lt. 0d0 ) then
          ppn(3) = ABS(ezcrf(iremx,izemx)) / emx * psurf
        else
          ppn(3) = -ABS(ezcrf(iremx,izemx)) / emx * psurf
        end if
!
! Move launch point several microns off the surface
        th = ATAN2(ppn(2),ppn(3))
        rlaun = rlaun + offlaun * SIN(th)
        zlaun = zlaun + offlaun * COS(th)
! ---------------------------------------------------------
      else if( launch .eq. 2 ) then     ! use specified launch direction
        ppn(2) = pylaunch
        ppn(3) = pzlaunch
! ---------------------------------------------------------
      else if( launch .eq. 3 ) then
! Use geometric normals at specified boundary points
        do i=1,nbplst
        kk = ibplst(i)
        call NORMAL(kk,th)
        pybplst(i) = psurf * SIN(th)
        pzbplst(i) = psurf * COS(th)
!
! Move launch point several microns off the surface
        ybplst(i) = ybplst(i) + offlaun * SIN(th)
        zbplst(i) = zbplst(i) + offlaun * COS(th)
        end do
!
      end if
!
      end
!  ********************************************************
      subroutine EMISSION(ek,th,ph)
!
! Generate kinetic energy and exit angles for secondary emission
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
      dimension cpf(100)
      GAUSS(a,b,c) = 0.3989422804d0/c*EXP(-0.5d0*((a-b)/c)**2)
!
! Polar emission angle ~cos(theta)
! Ref. Furman et al, PRSTAB 5:124404 (2002), p.4
      th = ASIN( RAN1(irndum) )
      ph = RANV(0d0,twopi)
!
! Assume kinetic energy comes from three sources
!  cf. Cimino et al, PRL 93:014801, Fig. 1
!  1. true secondary (ts), 2. elastic (e), 3. uniform background (u)
!
! Set weight of each source depending on parent energy
      if( ekin .le. 20e-6 ) then
        wtts = 0.1
        wte = 0.6
      else if( ekin .ge. 100e-6 ) then
        wtts = 0.6
        wte = 0.1
      else
        wtts = 0.1 + (ekin-20e-6)*6250.
        wte = 0.6 - (ekin-20e-6)*6250.
      end if
!
! Compute cumulative probability function
      de = ekin / 100
      sig = 2e-6
! Assume height of uniform background is 5% of the larger gaussian
      htu = 0.05*0.3989/sig*MAX(wtts,wte)
!
      do i=1,100
      e = i*de - de/2.
      p = wtts*GAUSS(e,5d-6,sig) + wte*GAUSS(e,ekin,sig) + htu
      if( i .eq. 1 ) then
        cpf(1) = p
      else
        cpf(i) = cpf(i-1) + p
      end if
      end do
      do i=1,100        ! normalize
      cpf(i) = cpf(i) / cpf(100)
      end do
!
! Use transformation method
      rn = RAN1(irndum)
      do i=1,100
      if( rn .le. cpf(i) ) exit
      end do
!
      ek = i*de
!
      end
!  ********************************************************
      subroutine FIELD(xx,tt)
!
! Get E and B at current location
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
!
      dimension brf(3),be(3),xx(3),ber(3)
      logical inside,rot
      save jzlo,jrlo,jr,jx,jy,jz
!
      rp = SQRT(xx(1)**2 + xx(2)**2)
      if( rp .gt. 0d0 ) then
        cphi = xx(1) / rp
        sphi = xx(2) / rp
      else
        cphi = 1d0
        sphi = 0d0
      end if
!
      if( robth.ne.0d0 .or. robphi.ne.0d0 ) then
        rot = .true.
        ct = COS(robth*rad)
        st = SIN(robth*rad)
        cp = COS(robphi*rad)
        sp = SIN(robphi*rad)
      else
        rot = .false.
      end if
!
! Get z in coordinates of cavity grid
      if( ascrf .eq. 0 ) then     ! axial symmetric case
        zp = ABS(xx(3))           ! relative z from center of gap
      else                        ! no axial symmetry
        zp = xx(3)
      end if
!
! Find lower-left corner of grid box that contains current position 
      call LOCATE3(zcrf,nzcrf,zp,jzlo) 
      call LOCATE3(rcrf,nrcrf,rp,jrlo)
      if( (jzlo.le.0) .or. (jzlo.ge.nzcrf) .or.  &
        (jrlo.le.0) .or. (jrlo.ge.nrcrf) ) then
        stop 'RF field interpolation error' 
      end if
!
! Find number of grid box corners that are inside the boundary
      ninside = 0
      if( INSIDE(zcrf(jzlo),rcrf(jrlo)) ) ninside = ninside + 1
      if( INSIDE(zcrf(jzlo+1),rcrf(jrlo)) ) ninside = ninside + 1
      if( INSIDE(zcrf(jzlo),rcrf(jrlo+1)) ) ninside = ninside + 1
      if( INSIDE(zcrf(jzlo+1),rcrf(jrlo+1)) ) ninside = ninside + 1
!
      if( ninside .eq. 4 ) then    ! whole box inside the boundary, interpolate
        u = ( zp - zcrf(jzlo) ) / dzcrf    
        v = ( rp - rcrf(jrlo) ) / drcrf
!
        call BILINEAR(v,u,jrlo,jzlo,ezcrf,mxrcrf,mxzcrf,ez)
        call BILINEAR(v,u,jrlo,jzlo,ercrf,mxrcrf,mxzcrf,er)
        call BILINEAR(v,u,jrlo,jzlo,bpcrf,mxrcrf,mxzcrf,bp)
!
      else if( ninside .eq. 0 ) then     ! whole box outside boundary, set field 0
        ez = 0d0
        er = 0d0
        bp = 0d0
!
      else        ! box crosses boundary, use nearest boundary point
        call NEAREST_GRID(zp,rp,jzlo,jrlo)
        if( iflag .lt. 0 ) go to 900
        ez = ezcrf(jrlo,jzlo)
        er = ercrf(jrlo,jzlo)
        bp = bpcrf(jrlo,jzlo)
      end if
!
      arg = twopi * freq * tt + phi0
      carg = COS(arg)
      sarg = SIN(arg)
!
      efld(1) = er * cphi * carg
      efld(2) = er * sphi * carg
      efld(3) = ez * carg
      brf(1) = -bp * sphi * sarg
      brf(2) = bp * cphi * sarg
      brf(3) = 0d0
!
! ---------------------------------------------------------
      if( magdim .eq. 2 ) then            ! 2D external magnetic field
        if( rot ) then        ! get coordinates in magnet grid
          x = xx(1)*cp - xx(2)*sp
          y = xx(1)*ct*sp + xx(2)*ct*cp - (xx(3)+zoff)*st
          zp = xx(1)*st*sp + xx(2)*st*cp + (xx(3)+zoff)*ct
          rp = SQRT( x*x + y*y )
        else
          zp = xx(3) + zoff              
        end if
!
        call LOCATE3(zgr,nzgr,zp,jz)
        call LOCATE3(rgr,nrgr,rp,jr)
!
        u = (zp - zgr(jz)) / delzgr          
        v = (rp - rgr(jr)) / delrgr
!
        call BILINEAR(u,v,jz,jr,bzgr,mxzgr,mxrgr,bz)
        call BILINEAR(u,v,jz,jr,brgr,mxzgr,mxrgr,br)
!
        be(1) = br * cphi
        be(2) = br * sphi
        be(3) = bz
! .........................................................
      else              ! 3D external magnetic field
        if( rot ) then        ! get coordinates in magnet grid
          x = xx(1)*cp - xx(2)*sp
          y = xx(1)*ct*sp + xx(2)*ct*cp - (xx(3)+zoff)*st
          z = xx(1)*st*sp + xx(2)*st*cp + (xx(3)+zoff)*ct
        else
          x = xx(1)
          y = xx(2)
          z = xx(3) + zoff
        end if
!
        call LOCATE3(xgr,nxgr,x,jx)   ! returns lower grid corner
        call LOCATE3(ygr,nygr,y,jy)   
        call LOCATE3(sgr,nsgr,z,jz)   
!
        u = (x-xgr(jx)) / dxgr
        v = (y-ygr(jy)) / dygr
        w = (z-sgr(jz)) / dsgr
!
        call TRILINEAR(u,v,w,jx,jy,jz,bxgr,mxxg,mxyg,mxzg,be(1))
        call TRILINEAR(u,v,w,jx,jy,jz,bygr,mxxg,mxyg,mxzg,be(2))
        call TRILINEAR(u,v,w,jx,jy,jz,bsgr,mxxg,mxyg,mxzg,be(3))
!
      end if
! ---------------------------------------------------------
      if( rot ) then          ! rotate field components into cavity coordinates
        ber(1) = be(1)*cp + be(2)*ct*sp + be(3)*st*sp
        ber(2) = -be(1)*sp + be(2)*ct*cp + be(3)*st*cp
        ber(3) = -be(2)*st + be(3)*ct
!
      else
        do i=1,3
        ber(i) = be(i)
        end do
      end if
!
! Get total magnetic field
      bfld(1) = brf(1) + ber(1)
      bfld(2) = brf(2) + ber(2)
      bfld(3) = brf(3) + ber(3)
!
900   continue
      end
!  *********************************************************
      function FIXFMT(string)
!
! Fix format problems getting f90 to read the Superfish TXT file
!
      real*8 fixfmt
      integer dp,expo
      character(16) string
      character(1) ch
!
! Find length of the string without trailing blanks
      nch = LEN_TRIM(string)
!
! Parse the string looking for decimal point or exponent
      dp = 0
      expo = 0
      do i=1,nch
      ch = string(i:i)
      if( ch .eq. '.' ) dp = i
      if( ch .eq. 'E' ) expo = i
      end do
!
      if( expo .eq. 0 ) then        ! no exponent
        if( dp .gt. 0 ) then        ! has decimal point
          read(string,'(f10.4)') fixfmt
        else                        ! no decimal point
          read(string,'(f10.0)') fixfmt
        end if
!
      else                          ! has exponent
        read(string,'(e16.4)') fixfmt
      end if
!
      end
!  ********************************************************
      FUNCTION GASDEV(idum)
      implicit double precision (a-h,o-z)
      INTEGER idum
!CU    USES ran1
      INTEGER iset
      SAVE iset,gset
      DATA iset/0/
!
      if (iset.eq.0) then
1       v1=2.*ran1(idum)-1.
        v2=2.*ran1(idum)-1.
        rsq=v1**2+v2**2
        if(rsq.ge.1..or.rsq.eq.0.)goto 1
        fac=sqrt(-2.*log(rsq)/rsq)
        gset=v1*fac
        gasdev=v2*fac
        iset=1
      else
        gasdev=gset
        iset=0
      endif
      return
      END
!  ********************************************************
      FUNCTION gammln(xx)
      implicit double precision (a-h,o-z)
      INTEGER j
      real*8 cof(6)
      SAVE cof,stp
      DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,  &
        24.01409824083091d0,-1.231739572450155d0,  &
        .1208650973866179d-2,-.5395239384953d-5,2.5066282746310005d0/
!
      x=xx
      y=x
      tmp=x+5.5d0
      tmp=(x+0.5d0)*log(tmp)-tmp
      ser=1.000000000190015d0
      do 11 j=1,6
        y=y+1.d0
        ser=ser+cof(j)/y
11    continue
      gammln=tmp+log(stp*ser/x)
      return
      END
!  *********************************************************
      function GAUSS(xm,sd)
!
!  Return a random value from a gaussian distribution
!   with mean  xm  and standard deviation  sd
!  uses Numerical Recipe function  GASDEV
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
!
!  Get random value from gaussian distribution with mean  0
!   and standard deviation  1
      gauss = GASDEV(irndum)
!
!  Renormalize to requested mean and standard deviation
      gauss = sd*gauss + xm
!
      end
! *********************************************************
      subroutine HEADINGS(ifl)
!
! Write coulumn headings for the data output
!
      write(ifl,'(a)')'Trajectory summary'
      write(ifl,50)
50    format(t3,'par',t13,'stp',t19,'flg',t25,'phi [deg]',  &
        t38,'x [m]',t52,'y [m]',t66,'z [m]',  &
        t80,'px [MeV/c]',t94,'py [MeV/c]',t108,'pz [MeV/c]',  &
        t122,'t [s]',t136,'Ek [MeV]',t150,'ang [deg]',  &
        t161,'sey',  &
        t170,'x0 [m]',t184,'y0 [m]',t198,'z0 [m]',  &
        t212,'px0 [MeV/c]',t226,'py0 [MeV/c]',t240,'pz0 [MeV/c]',  &
        t254,'t0 [s]',  &
        t268,'Bx [T]',t282,'By [T]',t296,'Bz [T]',  &
        t310,'Ex [V/m]',t324,'Ey [V/m]',t338,'Ez [V/m]',  &
        t352,'step [s]',t366,'diff [m]',t380,'cbp')
!
      end
! *********************************************************
      function INSIDE(z,r)
!
! Return  true  if point (z,r) is inside the cavity boundary
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
      logical inside
!
! Find closest boundary point (cbp) assuming azimuthal symmetry
      dmin =1d9
      do i=1,ncbpt
      d = SQRT( (cbx(i)-z)**2 + (cby(i)-r)**2 )
      if( d .lt. dmin ) then
        dmin = d
        jcbp = i
        kc = i
      end if
      end do
!
! Get angles from the cbp to adjacent boundary points and particle location
      if( kc.gt.1 .and. kc.lt.ncbpt  ) then
        thp = ATAN2(cby(kc+1)-cby(kc),cbx(kc+1)-cbx(kc))
        if( thp .lt. 0d0 ) thp = thp + twopi
        thm = ATAN2(cby(kc-1)-cby(kc),cbx(kc-1)-cbx(kc))
        if( thm .lt. 0d0 ) thm = thm + twopi
        th = ATAN2(r-cby(kc),z-cbx(kc))
        if( th .lt. 0d0 ) th = th + twopi
!
      else if( kc .eq. 1 ) then
        thp = ATAN2(cby(2)-cby(1),cbx(2)-cbx(1))
        if( thp .lt. 0d0 ) thp = thp + twopi
        thm = 3d0*piovr2
        th = ATAN2(r-cby(kc),z-cbx(kc))
        if( th .lt. 0d0 ) th = th + twopi
!
      else if( kc .eq. ncbpt ) then
        thp = 3d0*piovr2
        thm = ATAN2(cby(kc-1)-cby(kc),cbx(kc-1)-cbx(kc))
        if( thm .lt. 0d0 ) thm = thm + twopi
        th = ATAN2(r-cby(kc),z-cbx(kc))
        if( th .lt. 0d0 ) th = th + twopi
      end if
!
! Compare angles to see if the point is inside the cavity
      if( ccw ) then          ! point entered counter clockwise
        if( thp .eq. 0d0 ) then       ! upper flat surface
          if( th.le.thp .and. th.ge.thm ) then
            inside = .false.
          else
            inside = .true.
          end if
!
        else                  ! not on upper flat surface
          if( cbx(kc) .ge. 0d0 ) then           ! on right side of cavity
            if( th.gt.thp .and. th.lt.thm ) then
              inside = .true.
            else
              inside = .false.
            end if
!
          else                                  ! on left side of cavity
            if( th.le.thp .and. th.ge.thm ) then
              inside = .false.
            else
              inside = .true.
            end if
          end if
        end if
 ! --------------------------------------------------------
      else                    ! point entered clockwise
        if( thp .eq. 0d0 ) then       ! upper flat surface
          if( th.le.thm .and. th.ge.thp ) then
            inside = .false.
          else
            inside = .true.
          end if
!
        else                  ! not on upper flat surface
          if( cbx(kc) .ge. 0d0 ) then           ! on right side of cavity
            if( th.gt.thm .and. th.lt.thp ) then
              inside = .true.
            else
              inside = .false.
            end if
!
          else                                  ! on left side of cavity
            if( th.le.thm .and. th.ge.thp ) then
              inside = .false.
            else
              inside = .true.
            end if
          end if
        end if
      end if            ! end of test on ordering of points
!
      end
! *********************************************************
      subroutine LAUNCH_POINT
!
! Find launch point (zlaun,rlaun) for the electrons
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
!
      if( launch .eq. 1 ) then      ! use the peak electric field location
        emx = 0d0
        zemx = 0.
        remx = 0.
        izemx = 1
        iremx = 1
        nbplst = 1
!
        do iz=1,nzcrf
        do ir=1,nrcrf
        emag = SQRT(ezcrf(ir,iz)**2 + ercrf(ir,iz)**2)
        if( emag .gt. emx ) then
          emx = emag
          zemx = zcrf(iz)
          remx = rcrf(ir)
          izemx = iz
          iremx = ir
        end if
        end do
        end do
        write(2,90) emx,izemx,iremx,zemx,remx
90      format('Maximum E field: ',e12.4,3x,2i4,3x,2f12.4)
! Find nearest cavity boundary point
        dmin = 1e9
        imin = 1
        do ib=1,ncbpt
        d = SQRT((cbx(ib)-zemx)**2 + (cby(ib)-remx)**2)
        if( d .lt. dmin ) then
          dmin = d
          imin = ib
        end if
        end do
        rlaun = cby(imin)
        zlaun = cbx(imin)
        if( otherside .and. ascrf.eq.0 ) zlaun = -zlaun
! ---------------------------------------------------------
      else if( launch .eq. 2 ) then      ! user specified location
        nbplst = 1
        rlaun = rlaunch
        zlaun = zlaunch
! ---------------------------------------------------------
      else if( launch .eq. 3 ) then      ! user specified boundary points
        open(unit=8,file='cavel.lst',status='old',iostat=ioc)
        if( ioc .ne. 0 ) then
          write(*,*) 'Cant open file cavel.lst'
          stop
        end if
        read(8,*) nbplst
        read(8,*) (ibplst(i),i=1,nbplst)
        do i=1,nbplst
        kk = ibplst(i)
        ybplst(i) = cby(kk)
        zbplst(i) = cbx(kk)
        end do
        close(8)
!
      end if
!
      end
!  ********************************************************
      subroutine LOCATE3(xx,n,x,j)
!
! Locate position J of variable X in the array XX(1:N)
!
      implicit double precision (a-h,o-z)
      dimension xx(n)
!
      jl = 0
      jh = n + 1
!
10    continue
      if( jh-jl .gt. 1 ) then
        jm = (jl + jh) / 2
        if( x .gt. xx(jm) ) then
          jl = jm
        else
          jh = jm
        end if
        go to 10
      end if
!
      j = jl
!
      if( j .eq. 0 ) j = 1
      if( j .eq. n ) j = n - 1
!
      end
! *********************************************************
      subroutine NEAREST_GRID(z,r,jzlo,jrlo)
!
! Find nearest RF grid point to (z,r) that is inside cavity boundary
!
! On entry (jzlo,jrlo) contains LL corner of grid box containing (z,r)
! On exit (jzlo,jrlo) contains nearest RF grid point to (z,r)
!  that is inside the cavity
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
      logical inside
!
! Search the given grid box and its nearest neighbor boxes
      dmin =1d9
      ninside = 0
!
      do inc=2,15
        do iz=-inc,inc
        jz = jzlo + iz
        if( (jz.ge.1) .and. (jz.le.nzcrf) ) then
          do ir=-inc,inc
          jr = jrlo + ir
          if( (jr.ge.1) .and. (jr.le.nrcrf) ) then
            if( INSIDE(zcrf(jz),rcrf(jr)) ) then
              ninside = ninside + 1
              d = SQRT( (zcrf(jz)-z)**2 + (rcrf(jr)-r)**2 )
              if( d .lt. dmin ) then
                dmin = d
                jzn = jz
                jrn = jr
              end if
            end if
          end if
          end do
        end if
        end do
        if( ninside .gt. 0 ) go to 100
      end do            ! end loop on inc
!
      iflag = -11
      go to 900
!      stop 'Cant find nearest grid point'
!
100   continue
      jzlo = jzn
      jrlo = jrn
!
900   continue
      end
!  *********************************************************
      subroutine NORMAL(kk,th)
!
! Find angle TH of the inwards normal at boundary point KK
!
! Return TH in range  [-pi,0]
! Return result independent of cw or ccw point entry
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
!
      if( kk.gt.1 .and. kk.lt.ncbpt ) then
        dx = cbx(kk+1) - cbx(kk-1)
        dy = cby(kk+1) - cby(kk-1)
      else if( kk .eq. 1 ) then
        dx = cbx(2) - cbx(1)
        dy = cby(2) - cby(1)
      else if( kk .eq. ncbpt ) then
        dx = cbx(ncbpt) - cbx(ncbpt-1)
        dy = cby(ncbpt) - cby(ncbpt-1)
      end if
!
      if( ccw ) then
        dx = -dx
        dy = -dy
      end if
!
      if( dx .ne. 0d0 ) then
        slp = dy / dx
        if( slp .ne. 0d0 ) then
          alp = ATAN(-1d0/slp)
          if( cbx(kk) .lt. 0d0 ) then
            th = alp
          else
            th = alp - pi
          end if
        else            ! slope = 0
!          if( dx .gt. 0d0 ) then
            th = -piovr2
!          else
!            th = piovr2
!          end if
        end if
!
      else              ! dx = 0
        if( dy .gt. 0d0 ) then
          th = 0d0
        else
          th = -pi
        end if
      end if
!         
      end
!  *********************************************************
      function ON_BOUNDARY()
!
! Determine if electron path crosses cavity boundary
! Use smaller fixed steps near boundary
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
      logical on_boundary,inside,in
!
      on_boundary = .false.
!
      rp = SQRT(xp(1)**2 + xp(2)**2)
!
! Did this step cross the boundary?
! Get closest boundary point to present location
      in = INSIDE(xp(3),rp)       ! loads jcbp
      diff = SQRT(xp(3)**2 + rp**2) - SQRT(cbx(jcbp)**2 + cby(jcbp)**2)
!
      if( .not.in ) then
        if( near_boundary ) then
          on_boundary = .true.
        else
          do i=1,3
          xp(i) = xp1(i)
          end do
          tp = tp1
          hfix = hdid/nfix
          near_boundary = .true.
          iflag = 3
        end if
      end if
!
      end
!  *********************************************************
      subroutine OUTPUT(ifl)
!
! Write out particle information
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
!
      phi = phi0 * 180d0 / pi
!
      write(ifl,50) nparc,kstp,iflag,phi,xp,pp,tp,ekin,sang,  &
        sey,xp0,pp0,tp0,bfld,efld,hdid,diff,jcbp
50    format(a7,i8,i5,f11.3,3x,8e14.6,f11.3,f6.1,3x,15e14.6,i6)
!
      end
!  ********************************************************
      FUNCTION poidev(xm,idum)
      implicit double precision (a-h,o-z)
      INTEGER idum
      PARAMETER (PI=3.141592654)
!CU    USES gammln,ran1
      SAVE alxm,g,oldm,sq
      DATA oldm /-1./
!
      if (xm.lt.12.)then
        if (xm.ne.oldm) then
          oldm=xm
          g=exp(-xm)
        endif
        em=-1
        t=1.
2       em=em+1.
        t=t*ran1(idum)
        if (t.gt.g) goto 2
      else
        if (xm.ne.oldm) then
          oldm=xm
          sq=sqrt(2.*xm)
          alxm=log(xm)
          g=xm*alxm-gammln(xm+1.)
        endif
1       y=tan(PI*ran1(idum))
        em=sq*y+xm
        if (em.lt.0.) goto 1
! ---------------------------------------------------------
!  RCF modified 28 Sep 99
!   prevent large value of em (e.g. 559) from causing
!   a floating point underflow in the calculation of t, below
        if( em .gt. 30. ) go to 1
! ---------------------------------------------------------              
        em=int(em)
        t=0.9*(1.+y**2)*exp(em*alxm-gammln(em+1.)-g)
        if (ran1(idum).gt.t) goto 1
      endif
      poidev=em
      return
      END
!  *********************************************************
      function RANV(xlo,xhi)
!
!  Return uniform distribution of variable over specified range
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
!
      rn = RAN1(irndum)
      RANV = (xhi-xlo)*rn + xlo
!
      end
!  ********************************************************
      FUNCTION RAN1(idum)
      implicit double precision (a-h,o-z)
      INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
      PARAMETER (IA=16807,IM=2147483647,AM=1.d0/IM,IQ=127773,  &
        IR=2836,NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=2.3d-16,RNMX=1.d0-EPS)
      INTEGER j,k,iv(NTAB),iy
      SAVE iv,iy
      DATA iv /NTAB*0/, iy /0/
!
      if (idum.le.0.or.iy.eq.0) then
        idum=max(-idum,1)
        do 11 j=NTAB+8,1,-1
          k=idum/IQ
          idum=IA*(idum-k*IQ)-IR*k
          if (idum.lt.0) idum=idum+IM
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0) idum=idum+IM
      j=1+iy/NDIV
      iy=iv(j)
      iv(j)=idum
      ran1=min(AM*iy,RNMX)
      return
      END
!  *********************************************************
      SUBROUTINE RKCK(y,dydx,n,x,h,yout,yerr,derivs)
!
! Take one 5th order Runge Kutta step
!
      implicit double precision (a-h,o-z)
      INTEGER n,NMAX
      dimension dydx(n),y(n),yerr(n),yout(n)
      EXTERNAL derivs
      PARAMETER (NMAX=6)
!    USES derivs
      INTEGER i
      dimension ak2(NMAX),ak3(NMAX),ak4(NMAX),ak5(NMAX),ak6(NMAX)
      dimension ytemp(NMAX)
!
      PARAMETER (A2=.2d0,A3=.3d0,A4=.6d0,A5=1.d0,A6=.875d0)
      parameter (B21=.2d0,B31=3.d0/40.d0,B32=9.d0/40.d0)
      parameter (B41=.3d0,B42=-.9d0,B43=1.2d0,B51=-11.d0/54.d0)
      parameter (B52=2.5d0,B53=-70.d0/27.d0,B54=35.d0/27.d0)
      parameter (B61=1631.d0/55296.d0,B62=175.d0/512.d0)
      parameter (B63=575.d0/13824.d0,B64=44275.d0/110592.d0)
      parameter (B65=253.d0/4096.d0,C1=37.d0/378.d0)
      parameter (C3=250.d0/621.d0,C4=125.d0/594.d0,C6=512.d0/1771.d0)
      parameter (DC1=C1-2825.d0/27648.d0,DC3=C3-18575.d0/48384.d0)
      parameter (DC4=C4-13525.d0/55296.d0,DC5=-277.d0/14336.d0)
      parameter (DC6=C6-.25d0)
!
      do 11 i=1,n
        ytemp(i)=y(i)+B21*h*dydx(i)
11    continue
      call derivs(x+A2*h,ytemp,ak2)
      do 12 i=1,n
        ytemp(i)=y(i)+h*(B31*dydx(i)+B32*ak2(i))
12    continue
      call derivs(x+A3*h,ytemp,ak3)
      do 13 i=1,n
        ytemp(i)=y(i)+h*(B41*dydx(i)+B42*ak2(i)+B43*ak3(i))
13    continue
      call derivs(x+A4*h,ytemp,ak4)
      do 14 i=1,n
        ytemp(i)=y(i)+h*(B51*dydx(i)+B52*ak2(i)+B53*ak3(i)+B54*ak4(i))
14    continue
      call derivs(x+A5*h,ytemp,ak5)
      do 15 i=1,n
        ytemp(i)=y(i)+h*(B61*dydx(i)+B62*ak2(i)+B63*ak3(i)  &
         +B64*ak4(i)+B65*ak5(i))
15    continue
      call derivs(x+A6*h,ytemp,ak6)
      do 16 i=1,n
        yout(i)=y(i)+h*(C1*dydx(i)+C3*ak3(i)+C4*ak4(i)+C6*ak6(i))
16    continue
      do 17 i=1,n
        yerr(i)=h*(DC1*dydx(i)+DC3*ak3(i)+DC4*ak4(i)+DC5*ak5(i)  &
         +DC6*ak6(i))
17    continue
!
      END
!  *********************************************************
      SUBROUTINE RKQS(y,dydx,n,x,eps,yscal,derivs)
!      SUBROUTINE RKQS(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs)
!
! Adaptive stepsize controller for RKCK
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
      INTEGER n,NMAX
      dimension dydx(n),y(n),yscal(n)
      EXTERNAL derivs
      PARAMETER (NMAX=6)
!CU    USES derivs,rkck
      INTEGER i
      dimension yerr(NMAX),ytemp(NMAX)
      PARAMETER (SAFETY=0.9d0,PGROW=-.2d0,PSHRNK=-.25d0,ERRCON=1.89d-4)
!
      h=htry
1     call RKCK(y,dydx,n,x,h,ytemp,yerr,derivs)
      errmax=0.
      do 11 i=1,n
        errmax=max(errmax,abs(yerr(i)/yscal(i)))
11    continue
      errmax=errmax/eps
!
      if(errmax.gt.1.d0)then
        htemp=SAFETY*h*(errmax**PSHRNK)
        h=sign(max(abs(htemp),0.1d0*abs(h)),h)
        if( h .lt. hmin ) then
          h = hmin
          call RKCK(y,dydx,n,x,h,ytemp,yerr,derivs)
          hdid=h
          x=x+h
          do i=1,n
          y(i)=ytemp(i)
          end do
          hnext = h
          iflag = 1
          return
        end if
!
        xnew=x+h
        if(xnew.eq.x) then
          write(*,*) 'stepsize underflow in rkqs'
          read(*,'()')
        end if
        goto 1
      else
        if(errmax.gt.ERRCON)then
          hnext=SAFETY*h*(errmax**PGROW)
        else
          hnext=5.d0*h
        endif
        hdid=h
        x=x+h
        do 12 i=1,n
          y(i)=ytemp(i)
12      continue
        return
      endif
!
      END
! *********************************************************
      subroutine RZFIELD
!
! Read in 2D external magnetic field map
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
!
      character(80) title
!
      open(unit=13,file=magname,status='old',iostat=ioc)
      if( ioc .ne. 0 ) then
        write(*,*) 'Cant open file',magname
        write(2,*) 'Cant open file',magname
        stop
      end if
!
      read(13,'(a80)',iostat=ioc) title
      write(2,'(a80)') title
      read(13,*,iostat=ioc) nzgr,nrgr
      if( ioc .ne. 0 ) then
        write(2,*)'Error reading nzgr,nrgr'
        stop
      end if
      if( nzgr.gt.mxzgr .or. nrgr.gt.mxrgr ) then
        write(2,*)'Number of z or r grid points exceeds max'
        stop
      end if
!
      do i=1,nzgr
      do j=1,nrgr
      read(13,*,iostat=ioc)idm,jdm,zgr(i),rgr(j),bz,br
      if( ioc .ne. 0 ) then
        write(2,*)'Error reading field data file'
        stop
      end if
      bzgr(i,j) = bz * cscale
      brgr(i,j) = br * cscale
      end do
      end do
!
      delzgr = zgr(2) - zgr(1)
      delrgr = rgr(2) - rgr(1)
      close(13)
!
      end
! *********************************************************
      subroutine SECONDARY
!
! Driver routine for secondary electron emission
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
      parameter ( mxdau=8 )
      parameter ( mxgen=8 )
      parameter ( mxst=mxdau*mxgen )
      dimension xpst(3,mxst),ppst(3,mxst),tpst(mxst)
      dimension kdau(mxst),kgen(mxst)
      dimension ekd(mxdau),thd(mxdau),phd(mxdau)
      dimension pl(3,mxdau),pc(3,mxdau)
      character(2) nppc,cdau
      character(1) cgen
!
      ngen = 2          ! generation number
      kst = 0           ! stack pointer
      nsec = 0          ! total number of secondaries
!
! Get number of secondaries from a Poisson distribution
      ndau = POIDEV(sey,irndum)
      if( ndau .eq. 0 ) then
        go to 900
      else if( ndau .gt. mxdau ) then
        write(2,*)'*****     Limited number of daughters     *****'
        ndau = mxdau
      end if
!
100   continue
      if( ngen .gt. mxgen ) then
        go to 300
      end if
!
200   continue
      edau = 0d0        ! total secondary energy
      do i=1,ndau
      call EMISSION(ekd(i),thd(i),phd(i))
      edau = edau + ekd(i)
      end do
      do i=1,ndau       ! ensure conservation of energy
      ekd(i) = ekd(i) * ekin / edau
      end do
!
! Get momentum in local emission frame
      do id=1,ndau
      p = SQRT(ekd(id)**2 + 2d0*mel*ekd(id))
      cth = COS(thd(id))
      sth = SIN(thd(id))
      cph = COS(phd(id))
      sph = SIN(phd(id))
      pl(1,id) = p * sth * cph
      pl(2,id) = p * sth * sph
      pl(3,id) = p * cth
      end do
!
! Get theta = angle between z axis and surface normal at parent location
      call NORMAL(jcbp,th)
      sth = SIN(th)
      cth = COS(th)
! Get phi = azimuthal angle between y axis and parent location
      pxy = SQRT(pp(1)**2 + pp(2)**2)
      sph = pp(1) / pxy
      cph = pp(2) / pxy
!
! Rotate momentum to cavity frame
! Use pc = R(ph,z) R(th,x) pl
      do id=1,ndau
      pc(1,id) = pl(1,id)*cph + pl(2,id)*cth*sph + pl(3,id)*sth*sph
      pc(2,id) = -pl(1,id)*sph + pl(2,id)*cth*cph + pl(3,id)*sth*cph
      pc(3,id) = -pl(2,id)*sth + pl(3,id)*cth
      end do
!
! Push daughters onto stack
      do id=1,ndau
      kst = kst + 1
      nsec = nsec + 1
      if( nsec .gt. 99 ) then
        write(2,*)'*****     Exceeded maximum secondaries     *****'
        go to 900
      end if
      xpst(1,kst) = xp(1)
      xpst(2,kst) = xp(2) + offlaun*sth
      xpst(3,kst) = xp(3) + offlaun*cth
      do j=1,3
      ppst(j,kst) = pc(j,id)
      end do
      tpst(kst) = tp
      kdau(kst) = nsec
      kgen(kst) = ngen        
      end do      ! end of loop on daughters
!
300   continue
      if( kst .eq. 0 ) go to 900          ! stack is empty, quit
!
! Pop particle off bottom of stack
      do i=1,3
      xp0(i) = xpst(i,kst)
      pp0(i) = ppst(i,kst)
      end do
      tp0 = tpst(kst)
      idau = kdau(kst)
      ngen = kgen(kst)
!
      kst = kst - 1
      write(nppc,'(i2)') npar
      write(cdau,'(i2)') idau
      write(cgen,'(i1)') ngen
      nparc = nppc // '.' // cgen // '.' // cdau
!
      call TRACK
!
      if( iflag .lt. 0 ) then       ! this track failed
        call OUTPUT(2)
        call OUTPUT(9)
        go to 300
      end if
!
! Get kinetic energy
      ekin = SQRT(pp(1)**2 + pp(2)**2 + pp(3)**2 + mel**2 ) - mel
!
      call SURFACE_ANGLE    ! find angle that track hits the surface
!
      call SEC_YIELD
!
      call OUTPUT(2)    ! put final value in log file
      if( .not.rtuple ) call OUTPUT(9)
!
      ndau = POIDEV(sey,irndum)
!      IF( NGEN .EQ. 7 ) NDAU = 0    ! <<<<<<<<<<<<<<<<<<<<<<<<<
      if( ndau .eq. 0 ) then
        go to 300
      else if( ndau .gt. mxdau ) then
        ndau = mxdau
        ngen = ngen + 1
        go to 100
      else
        ngen = ngen + 1
        go to 100
      end if

900   continue
      end
! *********************************************************
      subroutine SEC_YIELD
!
! Compute secondary yield
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
!
! Compute secondary emission yield for perpendicular incidence
!  Ref. Cimino et al, PRL 93:014801, 2004
      s = 1.35
      x = ekin / seemax
      dtrue = sedelm*s*x / (s-1.+x**s)
      t1 = SQRT(ekin)
      t2 = SQRT(ekin+seeps0)
      delas = (t1-t2)**2 / (t1+t2)**2
      seyp = dtrue + delas
!
! Correct yield for non-normal incidence 
!  Ref. Kirby, SLAC-PUB-8380, 2000, p. 18
      sey = seyp * EXP(sefang*(1.-COS(sang)))
!
      end
! *********************************************************
      subroutine SETUP
!
! Set up problem description
!
!  File useage     1 : input parameters
!                  2 : simulation LOG file
!                  9 : electron data file
!                 11 : SuperFish cavity boundary
!                 12 : SuperFish RF fields
!                 13 : external B grid
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
!
      namelist/CONT/ascrf,bndname,cscale,epsstep,  &
        launch,magdim,magname,maxsteps,nfix,nphi,ntemp,  &
        offlaun,otherside,philo,phistp,psurf,pylaunch,  &
        pzlaunch,rfname,rfnrm,  &
        rlaunch,rnseed,robphi,robth,rskip,rtuple,secem,  &
        sedelm,seemax,seeps0,sefang,  &
        sigp,step,stepmax,stepmin,zlaunch,zoff
!
      pi = 4d0 * ATAN2(1d0,1d0)
      twopi = 2d0 * pi
      piovr2 = pi / 2d0
      rad = pi / 180d0
      irndum = rnseed
      dum = RAN1(irndum)        ! initialize random number seed
!
! Command file
      open(unit=1,file='cavel.inp',status='old',iostat=ioc)
      if( ioc .ne. 0 ) then
        write(*,*) 'Cant open file cavel.inp'
        stop
      end if
!
! Log file
      open(unit=2,file='cavel.log',status='unknown',iostat=ioc)   
      if( ioc .ne. 0 ) then
        write(*,*) 'Cant open file cavel.log'
        stop
      end if
!
      write(*,'(a,a)')' Welcome to CAVEL - version ',version
      write(2,'(a,a)')' Welcome to CAVEL - version ',version
!
! Read in problem commands
      read(1,cont,iostat=ioc)
      if( ioc .ne. 0 ) then
        write(2,*) 'Cant read control variable namelist'
        stop
      end if
      write(2,cont)
!
! Output data file
      open(unit=9,file='cavel.dat',status='unknown',iostat=ioc)   
      if( ioc .ne. 0 ) then
        write(*,*) 'Cant open file cavel.dat'
        stop
      end if
!
      call HEADINGS(9)  ! write headings on data file
!
! Read in cavity boundary
      call CAVITY_BOUNDARY
!
! Read in cavity RF fields
      call CAVITY_FIELDS
!
! Read in external magnetic field
      if( magdim .eq. 2 ) then
        call RZFIELD
      else
        call XYZFIELD
      end if
!
	end
! *********************************************************
      subroutine SIMULATE
!
! Simulate electron motion in RF cavity
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
      character(2) nppc
!
      npar = 0
      hmin = stepmin
      hmax = stepmax
!
! Find launch point for the electrons (zlaun,rlaun)
      call LAUNCH_POINT
!
      call HEADINGS(2)  ! write data headings in log file
!
! Find initial particle direction at this location (ppn)
! Set displacement from surface at start of simulation
      call DIRECTION
!
! ---------------------------------------------------------
!
      do il=1,nbplst          ! loop on launch points
      if( launch .eq. 3 ) then
        rlaun = ybplst(il)
        zlaun = zbplst(il)
        ppn(2) = pybplst(il)
        ppn(3) = pzbplst(il)
      end if
!
      xp0(1) = 0d0
      xp0(2) = rlaun
      xp0(3) = zlaun
      tp0 = 0d0
!
      do ip=1,nphi            ! loop on initial phases
      phi0 = philo*rad + (ip-1)*phistp*rad
!
      do it=1,ntemp           ! loop on thermal noise at launch
      npar = npar + 1
      write(nppc,'(i2)') npar
      nparc = nppc 
      write(*,'(a,i4)') 'Processing particle ',npar
      htry = step
      hdid = step             ! for initial output
      diff = 0d0              ! for initial output
      iflag = 0
      jcbp = 0
      nsec = 0
      sang = 0d0
!
! Apply thermal smearing of initial direction (pp0)
      if( sigp .gt. 0. ) then
        call SMEAR
      else
        do i=1,3
        pp0(i) = ppn(i)
        end do
      end if
!
      call TRACK        ! track electron thru cavity
!
      if( iflag .lt. 0 ) then
        call OUTPUT(2)
        call OUTPUT(9)
!        go to 900
        exit
      end if
!
! Get kinetic energy
      ekin = SQRT(pp(1)**2 + pp(2)**2 + pp(3)**2 + mel**2 ) - mel
!
      call SURFACE_ANGLE    ! find angle that track hits the surface
!
      if( secem .gt. 0 ) then
        call SEC_YIELD
      else
        sey = 0.
      end if
!
      if( .not.rtuple ) call OUTPUT(9)
      call OUTPUT(2)    ! put final value in log file
!
      if( secem .eq. 2 ) then        ! check for secondary trackss
        call SECONDARY
      end if
!
      end do            ! end of transverse loop
!
      end do            ! end of loop on phi
!
      end do            ! end of loop on launch points
!
900   continue
	end
! *********************************************************
      subroutine SMEAR
!
! Smear out momentum vector thermally
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
!
! Get smeared vector
      pp0(1) = ppn(1) + GAUSS(0d0,sigp/3d0)
      pp0(2) = ppn(2) + GAUSS(0d0,sigp/3d0)
      pp0(3) = ppn(3) + GAUSS(0d0,sigp/3d0)
!
! Dont let smearing reverse the direction from the boundary
      if( pp0(3)*ppn(3) .lt. 0d0 ) pp0(3) = 0d0
!
      end
! *********************************************************
      subroutine SURFACE_ANGLE
!
! Find the angle that the track hits the surface wrt normal
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
      dimension p(3)
!
! Find azimuthal angle of intersection point wrt y-z plane
      rt = SQRT( xp(1)**2 + xp(2)**2 )
      if( rt .gt. 0d0 ) then
        sphi = xp(1) / rt
        cphi = SQRT( 1d0 - sphi**2 )
      else
        sphi = 0d0
        cphi = 1d0
      end if
!
! Rotate pp into y-z plane. Cavity is azimuthally symmetric around z
      p(1) = cphi*pp(1) - sphi*pp(2)
      p(2) = sphi*pp(1) + cphi*pp(2)
      p(3) = pp(3)
!
      alpp = ATAN2(p(2),p(3))       ! angle of particle wrt z axis
!
! Find angle of normal vector
      call NORMAL(jcbp,alpn)                ! inward normal
      alpn = alpn + pi                      ! outward normal
!
      alyz = alpp - alpn
      if( ABS(alyz) .gt. pi ) alyz = twopi - ABS(alyz)
      alxz = ATAN2(p(1),p(3))
      sang = ACOS(1d0/SQRT(1d0+TAN(alyz)**2+TAN(alxz)**2)) / rad
!
      end
! *********************************************************
      subroutine TRACK
!
! Track an electron through the fields
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
      logical on_boundary
      parameter( nmax=6 )
!      parameter( tiny=1d-30 )
      dimension y(nmax),dydt(nmax),yscal(nmax)
      external deriv
!
      kstp = 0
      do i=1,3          ! define for initial output
      xp(i) = xp0(i)
      pp(i) = pp0(i)
      end do
      tp = tp0
      call FIELD(xp0,tp0)
      if( iflag .lt. 0 ) go to 900
      if( rtuple ) call OUTPUT(9)   ! initial value
!
      do i=1,3
      xp1(i) = xp0(i)
      pp1(i) = pp0(i)
      end do
      tp1 = tp0
      hnext = htry
      near_boundary = .false.
!
! ---------------------------------------------------------
100   continue
      kstp = kstp + 1
      if( kstp .gt. maxsteps ) then
        iflag = -10
        if( rtuple ) call OUTPUT(9)
        go to 900
      end if
      do i=1,3
      y(i) = xp1(i)
      y(i+3) = pp1(i)
      end do
      tt = tp1
!
      htry = hnext
      call DERIV(tt,y,dydt)
      if( iflag .lt. 0 ) go to 900
      do i=1,nmax
!      yscal(i) = ABS(y(i)) + ABS(htry*dydt(i))  + tiny
      yscal(i) = 1d0         ! [m,MeV]
      end do
!
! Take one step in time
      call RKQS(y,dydt,nmax,tt,epsstep,yscal,deriv)
      if( iflag .lt. 0 ) go to 900
!      call RKQS(y,dydt,nmax,tt,htry,epsstep,yscal,hdid,hnext,deriv)
!
      do i=1,3
      xp(i) = y(i)
      pp(i) = y(i+3)
      end do
      tp = tt
      ekin = SQRT(pp(1)**2 + pp(2)**2 + pp(3)**2 + mel**2 ) - mel
!
      if( ON_BOUNDARY() ) then      ! done tracking this electron?
        call SURFACE_ANGLE    ! find angle that track hits the surface
        if( rtuple ) call OUTPUT(9)
        go to 900
      else
        if( rtuple .and. MOD(kstp,rskip).eq.0 ) call OUTPUT(9)
      end if
!
      iflag = 0         ! can get warning from RKQS
      hnext = MIN(hnext,hmax)
      hnext = MAX(hnext,hmin)
      if( near_boundary ) then
        if( hfix .lt. hmin ) then
          hnext = hmin
          iflag = 2
        else
          hnext = hfix
        end if
      end if
!
! Remember particle values for this step
      do i=1,3
      xp1(i) = xp(i)
      pp1(i) = pp(i)
      end do
      tp1 = tp
!
      go to 100
! ---------------------------------------------------------
!
900   continue
      end
!  *********************************************************
      subroutine TRILINEAR(u,v,w,jx,jy,jz,fg,mx,my,mz,f)
!
!  Do tri-linear interpolation 
!
      implicit double precision (a-h,o-z)
      dimension fg(mx,my,mz)
!
      f = (1.-u)*(1.-v)*(1.-w) * fg(jx,jy,jz)  &
         +  (1.-u)*(1.-v)*w * fg(jx,jy,jz+1)   &
         +  (1.-u)*v*w * fg(jx,jy+1,jz+1)      &
         +  (1.-u)*v*(1.-w) * fg(jx,jy+1,jz)   &
         +   u*(1.-v)*(1.-w) * fg(jx+1,jy,jz)  &
         +   u*(1.-v)*w * fg(jx+1,jy,jz+1)     &
         +   u*v*w * fg(jx+1,jy+1,jz+1)        &
         +   u*v*(1.-w) * fg(jx+1,jy+1,jz)
!
      end
! *********************************************************
      subroutine XYZFIELD
!
! Read in 3D external magnetic field map
!
      implicit double precision (a-h,o-z)
      include 'cavel.inc'
!
      character(80) title
!
      open(unit=13,file=magname,status='old',iostat=ioc)
      if( ioc .ne. 0 ) then
        write(*,*) 'Cant open file',magname
        write(2,*) 'Cant open file',magname
        stop
      end if
!
      read(13,'(a80)') title
      read(13,*) nxgr,nygr,nsgr
      if( nxgr.gt.mxxg .or. nygr.gt.mxyg .or. nsgr.gt.mxzg ) then
        write(2,*)'Dimensions exceed max limit for 3D field'
        stop
      end if
!
      read(13,*) (xgr(j),j=1,nxgr)
      read(13,*) (ygr(j),j=1,nygr)
      read(13,*) (sgr(j),j=1,nsgr)
!
      nrecords = nxgr*nygr*nsgr
      dxgr = xgr(2) - xgr(1)
      dygr = ygr(2) - ygr(1)
      dsgr = sgr(2) - sgr(1)
!
      do i=1,nrecords
      read(13,*)ix,iy,iz,bx,by,bz
      bxgr(ix,iy,iz) = bx * cscale
      bygr(ix,iy,iz) = by * cscale
      bsgr(ix,iy,iz) = bz * cscale
      end do
!
      close(13)
!
      end