source: PSPA/parmelaPSPA/trunk/scheff3.f @ 404

      subroutine scheff3(dist,mt,wtsauv)
c
c     point to point space charge calculations by Bourat
c     ref : thesis Bourat 1988 - orsay
c
c----
c         dist = distance travelled by light since last call to scheff
c         sce(1)=beam current in amperes.
c         sce(2)=3 program number
c         sce(3)=opt particle size factor
c         sce(4)=no. of adjacent bunches
c         sce(5)=distance between bunches.
c
c-------------------------------------------------------------------------
      save
c
      include 'param_sz.h'
      include 'constcom.h'
      include 'coordcom.h'
      include 'fldcom.h'
      include 'misccom.h'
      include 'ncordscom.h'
      include 'pipecom.h'
      include 'psizescom.h'
      include 'sccom.h'
      include 'ucom.h'
c
      common/intype/intype
c----
      data cte/1.758804786/     ! e/m0 in cgs units
c
      if (dist) 20,10,20
c      intialization
  10   continue
       beami=sce(1)
       opt=abs(sce(3))
       if (opt.le.0.)opt=1.
       nip=sce(4)
       pl=sce(5)
       if(pl.le.0)pl=wavel
       xnp=float(npoints-1)
       qmacro=beami/(freq*xnp)
       rmacro=opt*(qmacro**.33333333333)*1.e-8
       facq=qmacro*cte
c
       write(nnout,*) ' Space charge subroutine number ',iprog
       if (ip.eq.1 .or. ip.eq.2 .or. ip.eq.4) then
          qmacroprint=qmacro*1.e+03
          write (nnout,410) qmacroprint
  410     format (' point to point space charge calculation'/
     .    '  qmacro = ',1pg11.3,' nC per superparticle')
          abeami = abs(beami)
          if (beami.gt.0.) write (nnout,420) abeami
  420     format (' beam current = ',1pg12.3,' amps')
          write (nnout,450) opt
  450     format (' particle size factor =',t40,1pg12.3)
       endif
c
       return
c
c      space charge calculation
c
  20   continue
       do 30 j=2,ngood
c        calcul of impuls of j
       xj=cord(1,j)
       gbxj=cord(2,j)
       yj=cord(3,j)
       gbyj=cord(4,j)
       zj=cord(5,j)
       if(zj.le.0)goto 30   ! jg (zj.le.z0)
       gbzj=cord(6,j)
       gamj=gam(j)
       bxj=gbxj/gamj
       byj=gbyj/gamj
       bzj=gbzj/gamj
c
       do 40 i=2,ngood
c       field create by i
       if(i.eq.j)goto 40
       xi=cord(1,i)
       gbxi=cord(2,i)
       yi=cord(3,i)
       gbyi=cord(4,i)
       gbzi=cord(6,i)
       gami=gam(i)
       bxi=gbxi/gami
       byi=gbyi/gami
       bzi=gbzi/gami
c
       do 50 ni=-nip,nip
c        adjacent bunches
       zi=cord(5,i)+ni*pl
c
       xij=xj-xi
       yij=yj-yi
       zij=zj-zi
       rijbi=xij*bxi+yij*byi+zij*bzi
       gi2=gami*gami/(gami+1.)
       grb=gi2*rijbi
c
       xpij=xij+bxi*grb
       ypij=yij+byi*grb
       zpij=zij+bzi*grb
       rpij=sqrt(xpij*xpij+ypij*ypij+zpij*zpij)
       if(rpij.le.rmacro)rpij=rmacro
       rp3=rpij**3
c
       epxin=xpij/rp3
       epyin=ypij/rp3
       epzin=zpij/rp3
c
       biep=bxi*epxin+byi*epyin+bzi*epzin
       bjep=bxj*epxin+byj*epyin+bzj*epzin
       bjbi=bxj*bxi+byj*byi+bzj*bzi
       bjbi1=1.-bjbi
       ggi=gami/(gami+1.)
       facbi=bjep-ggi*biep
c
       fdp=facq*dist
c
       dgbxj=fdp*gami*(epxin*bjbi1+bxi*facbi)
       dgbyj=fdp*gami*(epyin*bjbi1+byi*facbi)
       dgbzj=fdp*gami*(epzin*bjbi1+bzi*facbi)
c
       cord(2,j)=cord(2,j)+dgbxj
       cord(4,j)=cord(4,j)+dgbyj
       cord(6,j)=cord(6,j)+dgbzj
       gam(j)=sqrt(1.+cord(2,j)**2+cord(4,j)**2+cord(6,j)**2)
c
  50   continue
  40   continue
  30   continue
c
       return
       end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*