subroutine scheff1(dist,mt,wtsauv)
c         scheff version Lloyd Young
c         sce(1)=beam current in amperes.
c         sce(2)=1 program number
c         sce(3)=impact parameter,cm for point to point space-charge
c         calculation. if zero use mesh space-charge calculation.
c         sce(4)=radial mesh size in cm.
c         sce(5)=longitudinal mesh size in cm.
c         sce(6)=no. of radial mesh intervals (le 20)
c         sce(7)=no. of longitudinal mesh intervals (le 200)
c         sce(8)=no. of adjacent bunches
c         sce(9)=distance between bunches.  (if zero, pl=sce(5))
c         sce(10)=opt.  if opt.lt.0, use only one ring.  if)opt=0,
c                use 2x2 array. if opt.gt.0, number of rings is0,
c                determined by gaus depending on aspect ratio.
c                if opt even radial mesh intervals give equal volume rings.
c                if opt odd radial mesh entervals equal dr.
c                if opt.eq.4 or 5  assume metal wall at z=0. and calculate
c                image charge.
c         sce(11)=frm, remeshing factor (default=1.5)
c         sce(12)=rwall, radius of conducting wall. if zero, no wall
c--------------------------------------------------------------------------
      save
c
      include 'param_sz.h'
      include 'constcom.h'
      include 'coordcom.h'
      include 'fldcom.h'
      include 'misccom.h'
      include 'ncordscom.h'
      include 'sccom.h'
      include 'syscom.h'
      include 'ucom.h'
c
      parameter (maxdim=84000)
      common/cathcur/rcur
      real dxp(imaa),dyp(imaa),dzp(imaa)
      real rsce(nsce)
      logical fexist,fopen,endfile,qflag
      dimension qol(200)
c     4221=(20+1)*(200+1), 84000=20*(20+1)*200
      dimension rm(21), zm(201), ers(maxdim), ezs(maxdim), er(4221),
     1 ez(4221), aa(4000), ismax(201), iemax(201),rssq(20),zzs(201)
c--------------------------------------------------------------------------
c*
      data irout,izout/0,0/fopen,endfile,qflag/.false.,.false.,.true./
      if (ngood.le.1) return
       if(sce(3).ne.0.)go to 400
      if (dist) 50,10,50
   10 continue
      write(nnout,*) ' Space charge subroutine number ',iprog
      call rtrans(gt,bt)
       if(dist.eq.0.)then
       irout=0
       izout=0
       endif
   15 continue
      gmesh=gt
c          set up field tables
   16 continue
      beami=sce(1)
      rmax=sce(4)
      dzmax=sce(5)*gmesh
      nr=sce(6)
      nz=sce(7)
      nr1=nr+1
      nz1=nz+1
      im1=nr*nz
      im2=nr1*nz1
      im3=nr1*nz
      im4=nr*nr1
      nip=sce(8)
      opt=sce(10)
      if(amod(opt,2.).eq.0.)then
      dr=rmax*rmax/sce(6)
      else
      dr=rmax/sce(6)
      endif
      dz=dzmax/sce(7)
      pl=sce(9)*gmesh
      if(pl.le.0.)pl=dzmax
      frm=sce(11)
      if(frm.le.0.)frm=1.5
      rwall=sce(12)
      if(nr.gt.20 .or. nz.gt.200 .or. (nr*nr1*nz.gt.maxdim .and.
     *(maxdim/nr1**2*dz.lt.10*rwall .or. maxdim/nr1**2*dz.lt.20*rmax)))
     *then
      call appendparm
      stop' Abnormal stop too many mesh points in space charge grid '
      endif
c          load rm, zm, rssq, zs
      rm(1)=0.0
           do 20 i=2,nr1
      if(amod(opt,2.).eq.0.)then
      rm(i)=sqrt(float(i-1)*dr)
      else
           rm(i)=float(i-1)*dr
      endif
           rssq(i-1)=.5*(rm(i-1)**2+rm(i)**2)
  20       continue
      zs=.5*dz
           do 30 i=1,nz1
           zm(i)=float(i-1)*dz
           zzs(i)=zm(i)+zs
  30       continue
           hl=.5*dzmax
c---if there is a reference particle, set n1=2; otherwise, n1=1
      n1=1
c---  if (w0.gt.0.) n1=2
      xnp=npoints+1-n1
c          load ers and ezs
c     mesh dimensions are in cm. ers and ezs are in 1/cm.
c     c1, c2 and c3 are in cm., and c4 is in MeV-cm.
c     q=coulombs/point.   e/epsilon =11.294e6 cm MeV/Coul.
       if(qflag)then
      q=beami/(freq*1.e6*xnp)
      do 21 i=2,npoints
      if(weight(i).eq.0.)weight(i)=q
   21 continue
       qflag=.false.
       endif
c---printout during initialization
      if((ip.eq.1.and.dist.eq.0.) .or. (ip.eq.2.and.dist.eq.0.))then
        inip=sce(8)
        iopt=sce(10)
        write(nnout,301)
  301   format (' distribute charge over rings for space charge calc.')
        write(nnout,310) beami
  310   format (' beam current =',t40,1pg12.3,' amps')
        write(nnout,330) rmax,sce(5),nr,nz,inip,sce(9),iopt,frm,rwall
  330   format (
     .  ' radial size of mesh =',t40,g12.3,' cm'/
     .  ' z length of mesh =',t40,g12.3,' cm'/
     .  ' no. of radial mesh points =',t40,i8/
     .  ' no. of z mesh points =',t40,i8/
     .  ' no. of "identical pulses" =',t40,i8/
     .  ' pulse length =',t40,g12.3/
     .  ' OPT =',t40,i8/
     .  ' refresh when energy grows by factor',t40,g12.3/
     .  ' wall radius for images =',t40,g12.3)
      endif
c---
      c1=11.294e6/erest
c-- the above line has replaced the following line because the weight
c-- is the charge of the point
c     c1=11.294e6*q/erest
      inquire(file='scgrid',exist=fexist)
      if(.not.fopen)then
           open(unit=nschef,file='scgrid',access='sequential',
     *          status='unknown',form='unformatted')
           fopen=.true.
      endif
         j=1
      if(fexist)then
      if(endfile)go to 32
         read(nschef,err=34,end=33)rdz,rgmesh,rdzmax,rc1,rsce,ers,ezs
         if(abs(rdz-dz).gt.0.1*dz)go to 32
         do 31 i=4,12
         j=i
         if(sce(i).ne.rsce(i))go to 32
  31     continue
c stored data good
         if(dz.ne.rdz)then
         dz=rdz
         gmesh=rgmesh
         dzmax=rdzmax
           do 530 i=1,nz1
           zm(i)=float(i-1)*dz
           zzs(i)=zm(i)+zs
 530       continue
           hl=.5*dzmax
         endif
         if(rc1.ne.c1)then
            do 35 i=1,maxdim
            ers(i)=ers(i)*c1/rc1
            ezs(i)=ezs(i)*c1/rc1
  35        continue
         endif
         write(ndiag,*)' using stored scheff mesh data'
         go to 41
  33     continue
         write(ndiag,*)' end of file'
         endfile=.true.
         go to 32
  34     continue
         write(ndiag,*)' error on reading file pascgrid'
         endfile=.true.
c stored data no good
  32  continue
      backspace(nschef)
      write(ndiag,*)' not using stored data rdz,dz,i,rsce(i),sce(i)',
     *  rdz, dz,j,rsce(j),sce(j)
      endif
      l=0
      zm1=zm(1)
      zm2=zm(2)
      d2zm=(zm2-zm1)*2
      if(d2zm.lt.rmax)d2zm=rmax
c           do 40 k=1,nr the order of these two do loops are reversed
c                        to put the nearest ers and ezs in first so that
c                        those faraway can be neglected.
                do 40 j=1,nz
c*****          zp=zm(j+1)
                zp=zzs(j)
                if(zp.gt.d2zm)then
                topt=0
                elseif(zp.gt.2*d2zm)then
                topt=-1
                else
                topt=opt
                endif
            do 40 k=1,nr
                if(l.ge.maxdim-nr1)go to 42
                do 40 i=1,nr1
                     rp=rm(i)
                call gaus(rm(k),rm(k+1),zm1,zm2,topt,er1,ez1)
                     l=l+1
                     ers(l)=c1*er1
                     ezs(l)=c1*ez1
   40                continue
   42 continue
      write(nschef)dz,gmesh,dzmax,c1,sce,ers,ezs
      endfile(nschef)
   41 continue
      if (dist.eq.0.)return
      go to 60
c           evaluate and apply space charge effects.
   50 continue
      call rtrans(gt,bt)
      if(gmesh.le.0.)go to 15
      if(gt/gmesh.gt.frm)then
      gmesh=gmesh*frm
      go to 16
      endif
   60 continue
      zc=cord(5,1)
c      c3=dist/gt
c      c4=c3
c          evaluate ng, xbar, ybar, and zbar.
c     eng=ngood+1-n1
      eng=0.
      xbar=0.
      ybar=0.
      xsq=0.
      ysq=0.
      zbar=0.
           do 90 np=n1,ngood
      dbgz=cord(6,np)
      dbgx=cord(2,np)
      dbgy=cord(4,np)
      dgam=sqrt(1.d0+dbgx*dbgx+dbgy*dbgy+dbgz*dbgz)
c           gam(np)=gt*(dgam-bt*cord(6,np))
            bgstar=gt*(dbgz-bt*dgam)
      zstar=gt*(cord(5,np)-zc)
      cord(6,np)=bgstar
      xstar=cord(1,np)
      ystar=cord(3,np)
           xbar=xbar+xstar*weight(np)
           ybar=ybar+ystar*weight(np)
           xsq=xsq+xstar**2*weight(np)
           ysq=ysq+ystar**2*weight(np)
           zbar=zbar+zstar*weight(np)
           eng=eng+weight(np)
   90      continue
      xbar=xbar/eng
      ybar=ybar/eng
      xsq=xsq/eng
      ysq=ysq/eng
c      epsq=sqrt((xsq-xbar**2)/(ysq-ybar**2))
      epsq=1.
c      repsq=1./epsq
      repsq=1.
c      xfac=2./(1.+epsq)
      xfac=1.
c      yfac=epsq*xfac
      yfac=1.
      zbar=zbar/eng
      zc=zc+dzmax/nz*(ranf()-.5)/gt
c         clear and load bins
           do 100 i=1,im1
           aa(i)=0.0
  100      continue
      ng=0
         do 110 np=2,ngood
         charge=1.
         zmstar=0.
         zstar=gt*(cord(5,np)-zc)
           if((opt.eq.4 .or. opt.eq.5) .and. cord(5,np).lt.hl)then
               if(rcur.ne.0. .and. cord(5,np).lt.rcur)then
               a=sqrt((rcur-cord(5,np))**2+cord(1,np)**2+cord(3,np)**2)
               mcharge=-rcur/a
               aprime=rcur**2/a
               xmstar=cord(1,np)*aprime/a
               ymstar=cord(3,np)*aprime/a
               zmstar=gt*(rcur-aprime/a*(rcur-cord(5,np))-zc)
               else
               mcharge=-1.
               zmstar=-gt*(cord(5,np)+zc)
               endif
           endif
         xstar=cord(1,np)
         ystar=cord(3,np)
         rsq=repsq*(xstar-xbar)**2+epsq*(ystar-ybar)**2
          if(rsq.gt.rmax**2)go to 106
          if((opt.eq.4 .or. opt.eq.5) .and. cord(5,np).lt.0.)go to 110
           if(amod(opt,2.).eq.0.)then
           i=rsq/dr+1.
           else
           i=sqrt(rsq)/dr+1.
           endif
         if (i.gt.nr) go to 106
         z=zstar
 104     continue
         if (abs(z).le.hl) go to 105
         if(charge.lt.0.)go to 110
         if(nip.eq.0)go to 108
         if (pl.ne.dzmax) go to 108
         z=z-sign(pl,z)
c------distribute charge among adjacent bins.
  105      continue
           ng=ng+1
           zz=z+hl
           jm1=zz/dz+1.
           i1=i+1
           if(rsq.lt.rssq(i)) i1=i-1
           if(i1.lt.1)i1=1
           if(i1.gt.nr) i1=nr
           j1=jm1+1
           if(zz.lt.zzs(jm1)) j1=jm1-1
           if(pl.eq.dzmax)then
           if(j1.lt.1)j1=nz
           if(j1.gt.nz)j1=1
           else
           if(j1.lt.1)j1=1
           if(j1.gt.nz)j1=nz
           endif
           a=1.
           if(i1.ne.i)a=(rsq-rssq(i1))/(rssq(i)-rssq(i1))
           b=1.-a
           cc=1.
           if(j1.ne.jm1) cc=(zz-zzs(j1))/(zzs(jm1)-zzs(j1))
           d=1.-cc
           k=(jm1-1)*nr+i
           aa(k)=aa(k)+a*cc*weight(np)*charge
           k=k+i1-i
           aa(k)=aa(k)+b*cc*weight(np)*charge
           k=(j1-1)*nr+i
           aa(k)=aa(k)+a*d*weight(np)*charge
           k=k+i1-i
           aa(k)=aa(k)+b*d*weight(np)*charge
           if((opt.eq.4 .or. opt.eq.5.).and.zmstar.ne.0)then
             z=zmstar
                if(rcur.ne.0.and.cord(5,np).lt.rcur)then
                  xstar=xmstar
                  ystar=ymstar
                endif
             charge=mcharge
             zmstar=0.
             go to 104
           endif
           go to 110
  106 continue  
      if(irout.ne.0 .or. cord(5,np).lt.0.)go to 110
       write(ndiag,*) ' In element number ',cord(7,np)
       write(ndiag,107)np,xstar,ystar
  107 format(' warning--particle no.',i4,' is outside radial mesh'/
     *   ' xstar=',f8.5,', ystar=',f8.5)
      irout=1
      go to 110
  108 continue
      if(izout.ne.0 .or. cord(5,np).lt.0.)go to 110
      write(ndiag,*) ' In element number ',cord(7,np)
      write(ndiag,109) np,zstar,cord(5,np)
  109 format(' warning--particle no.',i4,' is outside longitudinal mesh'
     * / ' zstar=',f10.5,', z=',f10.5)
      izout=1
  110      continue
      eng=ng
c---calculate no. of particles in each column
      do 115 j=1,nz
      qol(j)=0.
      l=(j-1)*nr
      do 114 i=1,nr
      k=l+i
      qol(j)=qol(j)+aa(k)
  114 continue
      qol(j)=qol(j)/dz
  115 continue
c          find ismax for each j
           do 130 j=1,nz
           l=(j-1)*nr
           k=nr
                do 120 i=1,nr
                m=l+k
                if (aa(m)) 130,120,130
  120      continue
              k=k-1
  130      ismax(j)=k
c          find iemax for each j
      iemax(1)=1+ismax(1)
           do 140 j=2,nz
           iemax(j)=1+max0(ismax(j-1),ismax(j))
  140      continue
      iemax(nz1)=1+ismax(nz)
c          set er and ez to zero
           do 150 i=1,im2
           er(i)=0.0
           ez(i)=0.0
  150      continue
c          sum up fields
           do 200 js=1,nz
c           js1=js+1
       ism=ismax(js)
           if(ism.eq.0) go to 200
  160           continue
                do 190 is=1,ism
                l=(js-1)*nr+is
                a1=aa(l)
                if (a1.eq.0.) go to 190
c               l=(is-1)*im3 need to modify this because ers and ezs are
c                            loaded in different order
                l=(is-1)*nr1
c*****          do 170 je=1,js
                     do 170 je=1,js-1
c                    k1=l+(js-je)*nr1 need to modify this because ers and
c                                     ezs are loaded in different order
                     k1=l+(js-je)*im4
                     m1=(je-1)*nr1
                     iem=iemax(je)
                     if(iem.le.1 .or. k1.ge.maxdim-nr1)goto 170
                           do 168 ie=1,iem
c                          n=m1+ie
c                          k=k1+ie
                          er(m1+ie)=er(m1+ie)+a1*ers(k1+ie)
                          ez(m1+ie)=ez(m1+ie)-a1*ezs(k1+ie)
  168                continue
  170                continue
c*****                do 180 je=js1,nz1
                      do 180 je=js,nz
c                     k1=l+(je-js1)*nr1  need to modify this because ers and
c                                        ezs are loaded in different order
c*****                k1=l+(je-js1)*im4
                     k1=l+(je-js)*im4
                     m1=(je-1)*nr1
                     iem=iemax(je)
                     if(iem.le.1 .or. k1.ge.maxdim-nr1) go to 180
                          do 178 ie=1,iem
c                          n=m1+ie
c                          k=k1+ie
                           er(m1+ie)=er(m1+ie)+a1*ers(k1+ie)
                           ez(m1+ie)=ez(m1+ie)+a1*ezs(k1+ie)
  178                      continue
  180                continue
  190           continue
  200       continue
c*****
                do 300 i=1,nr1
                   ez(im3+i)=ez(i)
                   er(im3+1)=er(i)
  300              continue
c*****
c---evaluate and apply space charge impulse.
  201 continue
      do 240 np=n1,ngood
      bgfstar=cord(6,np)
cbug      if(cord(5,np).le.-dist*bt) go to 235
cbug      if(cord(5,np).lt.0.and.(opt.eq.4. .or. opt.eq.5))go to 235
cbug        if(cord(5,np).le.dist*bt)then
cbug            if(cord(5,np).gt.0)then
cbug               c3=.5*cord(5,np)/(bt*gt)+.5*dist/gt
cbug            else
cbug               c3=.5*(dist*bt+cord(5,np))**2/(bt*gt*(dist*bt))
cbug            endif
cbug        else
cbug            c3=dist/gt
cbug        endif
c goto 235 no space charge kicks until particle is emitted from the cathode
      if(cord(5,np).le.0) go to 235
c turn-off screening at dist*bt, skip screening for metallic cathodes
      if(cord(5,np).le.dist*bt.and.(opt.ne.4. .or. opt.ne.5))then
        if(cord(5,np).gt.0)then
           c3=.5*cord(5,np)/(bt*gt)+.5*dist/gt
        else
           c3=.5*(dist*bt+cord(5,np))**2/(bt*gt*(dist*bt))
        endif
      else
        c3=dist/gt
      endif
      zstar=gt*(cord(5,np)-zc)
      xstar=cord(1,np)
      ystar=cord(3,np)
           r=sqrt(repsq*(xstar-xbar)**2+epsq*(ystar-ybar)**2)
           if(np.eq.1) go to 202
           if(r.lt.0.000001) r=.000001
           xor=(xstar-xbar)*xfac/r
           yor=(ystar-ybar)*yfac/r
           go to 203
  202      continue
           xor=0.
           yor=0.
  203      continue
           z=zstar
  205      continue
           if(abs(z).le.hl) go to 210
           if (nip.le.0) go to 235
           if (pl.ne.dzmax) go to 235
           z=z-sign(pl,z)
           go to 205
c          interpolate impulse within mesh.
c 210      if(r.gt.rmax) go to 220
 210       continue
           if(amod(opt,2.).eq.0..and.r.le.rmax)then
           rb=(r*r/dr)
           i=1.+rb
           a=(r-rm(i))/(rm(i+1)-rm(i))
           elseif(r.le.rmax)then
           rb=r/dr
           i=1.0+rb
           a=rb-float(i-1)
           endif
           if(r.gt.rmax.and.rwall.ne.0)then
              a=(r-rmax)/(rwall-rmax)
              i=nr1
c---if rwall is zero arbitraily asume an effective rwall of 2*rmax to
c   calculate a and i
           elseif(r.gt.rmax)then
              a=(r-rmax)/rmax
              i=nr1
           endif
           b=1.0-a
c******    zb=(z+hl)/dz
            zb=(z+hl-.5*dz)/dz
            if(zb.lt.0.)zb=zb+float(nz)
c*****
           j=1.0+zb
           c=zb-float(j-1)
           d=1.0-c
           l=i+(j-1)*nr1
           m=l+nr1
           if(r.gt.rmax.and.r.lt.rwall)then
c the following two lines assume all the charge is within rmax.
               cbgr=c3*(d*er(l)+c*er(m))*rmax/r
               cbgz=c3*(d*b*ez(l)+c*b*ez(m))
            elseif(r.lt.rwall .or. rwall.eq.0.)then
               cbgr=c3*(d*(a*er(l+1)+b*er(l))+c*(a*er(m+1)+b*er(m)))
               cbgz=c3*(d*(a*ez(l+1)+b*ez(l))+c*(a*ez(m+1)+b*ez(m)))
           endif
c---The following makes a correction to the kinetic energy for the
c--- space charge potential energy.
      if(cord(5,np).lt.dist*bt.and.cord(5,np).gt.0..and.rwall.ne.0.
     *  .and.cord(6,np).gt.0..and.rm(iemax(j)).gt.0.)then
        iem=min0(iemax(j),iemax(j+1))
        cbgz=cbgz-(log(rwall)-log(rm(iem)))*(d*er(iem+(j-1)*nr1)+
     *       c*er(iem+j*nr1))
        cbgz=cbgz-.5*(d*(a*er(l+1)+b*er(l))+c*(a*er(m+1)+b*er(m))
     *      +d*er(l+1)+c*er(m+1))*(rm(i+1)-r)
        do 215 k=i+2,iem
        cbgz=cbgz-.5*(d*(er(k+(j-1)*nr1)+er(k-1+(j-1)*nr1))
     *      +c*(er(k+j*nr1)+er(k-1+j*nr1)))*(rm(k)-rm(k-1))
 215    continue
      endif
c          estimate impulse based on point charge at xbar,ybar,zbar.
c220       z=(zstar-zbar)*(1.+bt*cord(6,np)/gam(np))
c---the following lines are for a point charge at beam centroid.
c220       z=(zstar-zbar)
c          d=sqrt(z**2+r**2)
c          rod3=r/d**3
c          zod3=z/d**3
c          if(nip.eq.0) go to 228
c          include neighboring bunches.
c          do 226 i=1,nip
c                    xi=i
c                    do 226 j=1,2
c                    s=z+xi*pl
c                    d=sqrt(s**2+r**2)
c                    rod3=rod3+r/d**3
c                    zod3=zod3+s/d**3
c 226                xi=-xi
c          evaluate impulse.
c 228      cbgr=eng*c1*c3*rod3/(4.*pi)
c          cbgz=eng*c1*c4*zod3/(4.*pi)
c---following two lines are for line charge with density qol(j)
c 220 cbgz=0.
c     cbgr=c1*c3*qol(j)/(2.*pi*r**2)
c          apply impulse
  230      continue
           dbgx=cord(2,np)+cbgr*xor
           cord(2,np)=dbgx
           dbgy=cord(4,np)+cbgr*yor
           cord(4,np)=dbgy
           bgfstar=cord(6,np)+cbgz
  231      continue
           gfstar=sqrt(1.d0+dbgx*dbgx+dbgy*dbgy+bgfstar*bgfstar)
           cord(6,np)=gt*(bgfstar+bt*gfstar)
           gam(np)=gt*(gfstar+bt*bgfstar)
           go to 240
  235      continue
           dbgx=cord(2,np)
           dbgy=cord(4,np)
           bgfstar=cord(6,np)
           go to 231
ctest du bug
cbug      write(89,*) cbgr*xor*gam(np)/dist,
cbug     *            cbgr*yor*gam(np)/dist,
cbug     *            cbgz/dist,
cbug     *            cord(1,np),cord(3,np),cord(5,np)
c
  240 continue
      return
c
c  very simple inflexible 3d space charge calculator
c  uses point-to-point interaction with finite radius
c  as suggested by G. Boicourt
c  relevant scheff parameters are:
c  sce(1) = beam current, amps
c  sce(3) = impact parameter, cm (and flag for point-to-point)
c
  400 continue
      if (dist.eq.0.) then
        write(nnout,*) ' Space charge subroutine number ',iprog
        write(nnout,*) ' Point to point interaction with finite radius'
        write(nnout,*) ' Impact parameter =',sce(3),' cm'        
        return
      endif
      call rtrans(gt,bt)
      a3 = sce(3)**2
      points = npoints
      q=sce(1)/(freq*1.e6*points)
      con=8.984e5*q*dist/(gt*erest)
      do 410 i=1,ngood
      dxp(i) = 0.
      dyp(i) = 0.
      dzp(i) = 0.
  410 continue
      do 450 np=1,ngood
c
c  only do calculation once for each pair of points
c
cdir$ IVDEP
      do 430 mp=np+1,ngood
      x = cord(1,np) - cord(1,mp)
      y = cord(3,np) - cord(3,mp)
      z =(cord(5,np) - cord(5,mp))*gt
      d = x*x + y*y + z*z + a3
      d3i= 1./(d*sqrt(d))
      dxp(np) = dxp(np) + x * d3i
      dyp(np) = dyp(np) + y * d3i
      dzp(np) = dzp(np) + z * d3i
      dxp(mp) = dxp(mp) - x * d3i
      dyp(mp) = dyp(mp) - y * d3i
      dzp(mp) = dzp(mp) - z * d3i
  430 continue
c
c  increment xprime, yprime and energy
c
      rr=amax1(0.,cord(5,np))/(cord(5,np)+1.e-20)
      cord(2,np) = cord(2,np) + con * dxp(np) * rr
      cord(4,np) = cord(4,np) + con * dyp(np) * rr
      bgfstar=gt*(cord(6,np)-bt*gam(np))+con*dzp(np)*rr
      gfstar=sqrt(1.+cord(2,np)**2+cord(4,np)**2+bgfstar**2)
      cord(6,np)=gt*(bgfstar+bt*gfstar)
      gam(np)=gt*(gfstar+bt*bgfstar)
  450 continue
      return
      end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*