subroutine scheff2(dist,mt,wtsauv) ! new version of scheff c dist = distance travelled by light since last call to scheff c sce(1)=beam current in amperes. c if < 0, = -(no. of particles in the bunch) c sce(2)=2 program number c sce(3) = 0, use mesh space-charge calculation. c sce(3) > 0, use point to point space-charge calculation c sce(3) < 0, add images in the plane z = 0 to the point to point calc. c abs(sce(3))= nsig. Routine PVOL assigns a size to a particle which is c nsig*(rms bunch length)/ngood**.3333. Then if two c particles are close to one another the supercharge c is reduced in routine EH. c If nsig = 999, skip the call to PVOL 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.ge.2 radial mesh intervals give equal volume rings. c sce(11)=frm, remeshing factor (default=1.5) c sce(12)=rwall, radius of conducting wall. if zero, no wall 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 'syscom.h' include 'ucom.h' c parameter (maxdim=84000) logical fexist,fopen,endfil common/ehparm/xi,yi,zi,ex,ey,eez,hx,hy,hz,xj,yj,zj,gbxj,gbyj, . gbzj,gamma,qfac common/intype/intype common/out6/nbphase,nbremesh real xpipe(imaa),ypipe(imaa) real*4 dbgx(imaa),dbgy(imaa),dbgz(imaa) integer*4 ipipe(imaa) dimension qol(200) dimension rsce(nsce) 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,endfil/.false.,.false./ c if (ngood.le.1) return c skip if want point-by-point space charge calc. 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) 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=im1*nr1 nip=sce(8) inip = nip opt=sce(10) iopt = opt if(opt.ge.2)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, rs, zs rm(1)=0.0 rm(nr1)=rmax do 20 i=2,nr1 if(opt.ge.2)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) c rs(i-1)=sqrt(rssq(i-1)) 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. 1/epsilon_0 =11.294e6 cm mev/coul. q=beami/(freq*1.e6*xnp) c when beami < 0, it is -(no. of particles in the bunch) if (beami.lt.0.) then q = -1.602e-19*beami/npoints endif c printout during initialization if((ip.eq.1.and.dist.eq.0.) .or. (ip.eq.2.and.dist.eq.0.))then write(nnout,300) 300 format (' distribute charge over rings for space charge calc.') if (beami.gt.0.) write(nnout,310) beami 310 format (' beam current =',t40,1pg12.3,' amps') abeami = abs(beami) if (beami.lt.0.) write(nnout,320) abeami 320 format (' no. of electrons in bunch =',t40,1pg12.3) write(nnout,330) q,rmax,sce(5),nr,nz,inip,sce(9),iopt,frm,rwall 330 format ( . ' charge per superparticle =',t40,1pg12.3,' coulombs'/ . ' 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---- 1/epsilon_0 verification 22/02/94 c1=11.294e6*q/erest call inscgrid(fexist) if(.not.fopen)then open(unit=nschef,file='scgrid',access='sequential', * status='unknown',form='unformatted') fopen=.true. endif j=1 c if-construct change (jump into if-block not allowed) 04/94 if(endfil)then backspace(nschef) write(ndiag,*)'not using stored data rdz,dz,i,rsce(i),sce(i)', * rdz,dz,j,rsce(j),sce(j) else if(fexist)then read(nschef,err=34,end=33)rdz,rc1,rsce,ers,ezs if(abs(rdz-dz).gt.0.00001*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(rc1.ne.c1)then do 35 i=1,im4 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' endfil=.true. go to 32 34 continue write(ndiag,*)'error on reading file scgrid' endfil=.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 do 40 k=1,nr do 40 j=1,nz do 40 i=1,nr1 rp=rm(i) zp=zm(j+1) call gaus(rm(k),rm(k+1),zm(1),zm(2),opt,er1,ez1) l=l+1 ers(l)=c1*er1 ezs(l)=c1*ez1 40 continue write(nschef)dz,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) c3=dist/gt c4=c3 c evaluate ng, xbar, ybar, and zbar. eng=ngood+1-n1 xbar=0. ybar=0. xsq=0. ysq=0. zbar=0. do 90 np=n1,ngood c skip if particle not yet emitted nenp = cord(7,np) if (nenp.eq.0.and.intype.eq.10) go to 90 if (nenp.eq.0.and.intype.eq.11) go to 90 if (nenp.eq.0.and.intype.eq.15) go to 90 gstar=gt*(gam(np)-bt*cord(6,np)) bgstar=gt*(cord(6,np)-bt*gam(np)) c gstar=sqrt(1.+cord(2,np)**2+cord(4,np)**2+bgstar**2) zstar=gt*(cord(5,np)-zc) gam(np)=gstar cord(6,np)=bgstar c xstar=cord(1,np)+cord(2,np)*bt*zstar/gstar c ystar=cord(3,np)+cord(4,np)*bt*zstar/gstar xstar=cord(1,np) ystar=cord(3,np) xbar=xbar+xstar ybar=ybar+ystar xsq=xsq+xstar**2 ysq=ysq+ystar**2 c 80 zbar=zbar+zstar*(1.+bt*cord(6,np)/gstar) 80 zbar=zbar+zstar 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*(sandom(dum)-.5) c clear and load bins do 100 i=1,im1 aa(i)=0 100 continue ng=0 n108=0 do 110 np=2,ngood c skip if particle not yet emitted nenp = cord(7,np) if (cord(5,np).le.0.) go to 110 if (nenp.eq.0.and.intype.eq.10) go to 110 if (nenp.eq.0.and.intype.eq.11) go to 110 if (nenp.eq.0.and.intype.eq.15) go to 110 zstar=gt*(cord(5,np)-zc) c xstar=cord(1,np)+bt*zstar*cord(2,np)/gam(np) c ystar=cord(3,np)+bt*zstar*cord(4,np)/gam(np) xstar=cord(1,np) ystar=cord(3,np) rsq=repsq*(xstar-xbar)**2+epsq*(ystar-ybar)**2 if(opt.ge.2)then i=rsq/dr+1. else i=sqrt(rsq)/dr+1. endif if (i.gt.nr) go to 106 c z=zstar*(1.+bt*cord(6,np)/gam(np)) z=zstar if (abs(z).le.hl) go to 105 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(j1.lt.1)j1=nz if(j1.gt.nz)j1=1 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 k=k+i1-i aa(k)=aa(k)+b*cc k=(j1-1)*nr+i aa(k)=aa(k)+a*d k=k+i1-i aa(k)=aa(k)+b*d go to 110 106 continue if(irout.ne.0 .or. cord(5,np).lt.0.)go to 110 write(nnout,107) np,xstar,ystar,nenp 107 format(' warning--particle no.',i4,' is outside radial mesh'/ * ' xstar=',f8.5,', ystar=',f8.5,' in element ',i4) irout=1 go to 110 108 continue n108=n108+1 if(izout.ne.0 .or. cord(5,np).lt.0.)go to 110 write(nnout,109) np,zstar,cord(5,np),nenp 109 format(' warning--particle no.',i4,' is outside longitudinal mesh' * / ' zstar=',f10.5,', z=',f10.5, ' in element ',i4) 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)) 121,121,131 121 continue k=k-1 120 continue 131 continue ismax(j)=k 130 continue 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 c---------------------------------------------------------------------- c the following 'do nothing line' actually causes the compiler to insert c FINIT instruction and a FLDCW instruction to reset the 80287. if(ez(2).gt.1.)m1=m1 c---------------------------------------------------------------------- do 200 js=1,nz 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 l=(is-1)*im3 do 170 je=1,js k1=l+(js-je)*nr1 m1=(je-1)*nr1 iem=iemax(je) if(iem.le.1)goto 170 do 168 ie=1,iem c n=m1+ie c k=k1+ie c ------------------------------------------------------------------------- c The following instructions us the 80287 and it must be reset before c they work right. However, they rarely fail without the 80287 being reset! er(m1+ie)=er(m1+ie)+a1*ers(k1+ie) ez(m1+ie)=ez(m1+ie)-a1*ezs(k1+ie) c -------------------------------------------------------------------------- 168 continue 170 continue do 180 je=js1,nz1 k1=l+(je-js1)*nr1 m1=(je-1)*nr1 iem=iemax(je) if(iem.le.1) 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 evaluate and apply impulse do 240 np=n1,ngood c skip if particle not yet emitted nenp = cord(7,np) if (nenp.eq.0.and.intype.eq.10) go to 240 if (nenp.eq.0.and.intype.eq.11) go to 240 if (nenp.eq.0.and.intype.eq.15) go to 240 bgfstr=cord(6,np) if(cord(5,np).le.0.) then go to 235 else endif zstar=gt*(cord(5,np)-zc) c xstar=cord(1,np)+bt*zstar*cord(2,np)/gam(np) c ystar=cord(3,np)+bt*zstar*cord(4,np)/gam(np) xstar=cord(1,np) ystar=cord(3,np) r=sqrt(repsq*(xstar-xbar)**2+epsq*(ystar-ybar)**2) if(np.eq.1) go to 201 if(r.lt.0.000001) r=.000001 xor=(xstar-xbar)*xfac/r yor=(ystar-ybar)*yfac/r go to 203 201 continue xor=0. yor=0. c z=zstar*(1.+bt*cord(6,np)/gam(np)) 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(opt.ge.2.and.r.le.rmax)then rb=sqrt(r*r/dr) i=1.+rb a=(r-rm(i))/(rm(i+1)-rm(i)) else rb=r/dr i=1.0+rb a=rb-float(i-1) endif if(r.gt.rmax)then a=(r-rmax)/(rwall-rmax) i=nr1 endif b=1.0-a zb=(z+hl)/dz 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)then c the following two lines assume all the charge is within rmax. cbgr=c1*c3*qol(j)/(2.*pi*r**2) cbgz=c4*(d*b*ez(l)+c*b*ez(m)) else cbgr=c3*(d*(a*er(l+1)+b*er(l))+c*(a*er(m+1)+b*er(m))) cbgz=c4*(d*(a*ez(l+1)+b*ez(l))+c*(a*ez(m+1)+b*ez(m))) endif go to 230 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 cord(2,np)=cord(2,np)+cbgr*xor cord(4,np)=cord(4,np)+cbgr*yor bgfstr=cord(6,np)+cbgz 235 continue gfstar=sqrt(1.+cord(2,np)**2+cord(4,np)**2+bgfstr**2) cord(6,np)=gt*(bgfstr+bt*gfstar) gam(np)=gt*(gfstar+bt*bgfstr) c gam(np)=sqrt(1.+cord(2,np)**2+cord(4,np)**2+cord(6,np)**2) 240 continue return c c point-by-point calculation of the space charge K.T. McDonald c To invoke this option, set parmeter sce(3) nonzero c If sce(3) < 0, calculate image fields due to a metallic surface c at z = 0 400 continue c initialization if (dist.eq.0.) then c q = no. of electrons represented by each simulated particle xnp = npoints q = beami/(freq*1.e6*xnp)/1.602e-19 if (beami.lt.0.) then q = -beami/xnp endif do 405 i = 1,ngood ipipe(i) = 0 405 continue pvfac = abs(sce(3)) write(nnout,*) ' Space charge subroutine number ',iprog if (ip.eq.1 .or. ip.eq.2) then write(nnout,410) q 410 format (' point to point space charge calculation'/ . ' q = ',1pg11.3,' electrons per superparticle') abeami = abs(beami) if (beami.gt.0.) write(nnout,420) abeami 420 format (' beam current = ',1pg12.3,' amps') if (beami.le.0.) write(nnout,430) abeami 430 format (' no. of electrons in bunch =',1pg12.3) if (sce(3).lt.0.) write(nnout,440) 440 format (' include images in the plane z = 0') write(nnout,450) pvfac 450 format (' particle size factor =',t40,1pg12.3) endif pvfac = pvfac**2 return endif c----r0 = e / ((m0c**2/e)*(4*pi*epsilon0)) in cm r0 = 2.818e-13 c----integration factor = r0 * wavelength * dphi(radians) /2 pi rfac = r0*dist c----qrfac includes the integration factor qrfac = q*rfac if (ngood.lt.2) return c c setup for image calc. in pipe walls if (npipe.eq.0) go to 520 i = 1 500 continue if (i.gt.ngood) then go to 520 else c find which pipe the particle is in zi = cord(5,i) jmin = 1 if (ipipe(i).gt.0) jmin = ipipe(i) do 510 j = jmin,npipe if (zi.ge.zlpipe(j).and.zi.le.zhpipe(j)) go to 515 510 continue j = 0 515 continue ipipe(i) = j if (j.gt.0) then xi = cord(1,i) yi = cord(3,i) risq = xi**2 + yi**2 rpsq = rpipe(j)**2 if (risq.ge.rpsq) then c particle is lost if outside the pipe c dum is a dummy argument because in pardyn with have to send wtp in swap call swap(i,dum) i = i - 1 elseif (risq.eq.0.) then c skip if particle is on axis ipipe(i) = 0 else c store position of image charge rnew = rpsq/risq xpipe(i) = xi*rnew ypipe(i) = yi*rnew endif endif i = i + 1 endif go to 500 520 continue c c c find volume of a typical superparticle. If a second particle c lies within the volume of the first, reduce the effective c charge in subroutine eh xsqsiz = 0. ysqsiz = 0. zsqsiz = 0. xyzvol = 1. if (abs(sce(3)).ne.999.) call pvoldp c calc fields due to rest of bunch, including image fields c but don't apply space charge fields on the reference particle do 1300 i = 2,ngood nei = cord(7,i) c skip if particle i has not yet been emitted if (nei.eq.0.and.intype.eq.10) go to 1300 if (nei.eq.0.and.intype.eq.11) go to 1300 if (nei.eq.0.and.intype.eq.15) go to 1300 xi = cord(1,i) yi = cord(3,i) zi = cord(5,i) ex = 0. ey = 0. c name -ez- already in use eez = 0. hx = 0. hy = 0. hz = 0. do 1200 j = 1,ngood nej = cord(7,j) c skip if particle j has not yet been emitted if (nej.eq.0.and.intype.eq.10) go to 1200 if (nej.eq.0.and.intype.eq.11) go to 1200 if (nej.eq.0.and.intype.eq.15) go to 1200 xj = cord(1,j) yj = cord(3,j) zj = cord(5,j) gamma = gam(j) gbxj = cord(2,j) gbyj = cord(4,j) gbzj = cord(6,j) c skip to image charge calc. if i = j if (i.eq.j) go to 1250 c calc. field at i of moving charge j qfac = qrfac call eh c c now do the image fields in the plane z = 0 1250 continue c skip unless sce(3) < 0, and particle is in 1st element if (sce(3).ge.0.0 .or. nej.ne.1) go to 1260 c skip image calculation once particle is more than 1 cm from cathode if (zj.gt.1.) go to 1260 zj = -zj gbzj = -gbzj c image charge has opposite sign qfac = -qrfac c for image of particle i, only use 1 electron charge, and not the c supercharge q if (i.eq.j) qfac = -rfac call eh c restore coords zj = -zj gbzj = -gbzj c c calc. image fields in beam pipe, if any 1260 continue if (ipipe(j).eq.0) go to 1200 xj = xpipe(j) yj = ypipe(j) c don't change magnitude of transverse velocity, thus keeping beta < 1 c (if revise this, must recalc. gamma also) gbxj = -gbxj gbyj = -gbyj c image charge has opposite sign qfac = -qrfac c for image of particle i, only use 1 electron charge, and not the c supercharge q if (i.eq.j) qfac = -rfac call eh 1200 continue c calculate and store the impulse on particle i gamma = gam(i) bxi = cord(2,i)/gamma byi = cord(4,i)/gamma bzi = cord(6,i)/gamma dbgx(i) = ex + byi*hz - bzi*hy dbgy(i) = ey + bzi*hx - bxi*hz dbgz(i) =eez + bxi*hy - byi*hx 1300 continue c apply the impulses do 1400 i = 2,ngood cord(2,i) = cord(2,i) + dbgx(i) cord(4,i) = cord(4,i) + dbgy(i) cord(6,i) = cord(6,i) + dbgz(i) gam(i) = sqrt(1 + cord(2,i)**2 + cord(4,i)**2 + cord(6,i)**2) 1400 continue return end c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*