| 1 | subroutine scheff2(dist,mt,wtsauv) ! new version of scheff |
|---|
| 2 | c dist = distance travelled by light since last call to scheff |
|---|
| 3 | c sce(1)=beam current in amperes. |
|---|
| 4 | c if < 0, = -(no. of particles in the bunch) |
|---|
| 5 | c sce(2)=2 program number |
|---|
| 6 | c sce(3) = 0, use mesh space-charge calculation. |
|---|
| 7 | c sce(3) > 0, use point to point space-charge calculation |
|---|
| 8 | c sce(3) < 0, add images in the plane z = 0 to the point to point calc. |
|---|
| 9 | c abs(sce(3))= nsig. Routine PVOL assigns a size to a particle which is |
|---|
| 10 | c nsig*(rms bunch length)/ngood**.3333. Then if two |
|---|
| 11 | c particles are close to one another the supercharge |
|---|
| 12 | c is reduced in routine EH. |
|---|
| 13 | c If nsig = 999, skip the call to PVOL |
|---|
| 14 | c sce(4)=radial mesh size in cm. |
|---|
| 15 | c sce(5)=longitudinal mesh size in cm. |
|---|
| 16 | c sce(6)=no. of radial mesh intervals (le 20) |
|---|
| 17 | c sce(7)=no. of longitudinal mesh intervals (le 200) |
|---|
| 18 | c sce(8)=no. of adjacent bunches |
|---|
| 19 | c sce(9)=distance between bunches. (if zero, pl=sce(5)) |
|---|
| 20 | c sce(10)=opt. if opt.lt.0, use only one ring. if)opt=0, |
|---|
| 21 | c use 2x2 array. if opt.gt.0, number of rings is0, |
|---|
| 22 | c determined by gaus depending on aspect ratio. |
|---|
| 23 | c if opt.ge.2 radial mesh intervals give equal volume rings. |
|---|
| 24 | c sce(11)=frm, remeshing factor (default=1.5) |
|---|
| 25 | c sce(12)=rwall, radius of conducting wall. if zero, no wall |
|---|
| 26 | c |
|---|
| 27 | c----------------------------------------------------------------------- |
|---|
| 28 | save |
|---|
| 29 | c |
|---|
| 30 | include 'param_sz.h' |
|---|
| 31 | include 'constcom.h' |
|---|
| 32 | include 'coordcom.h' |
|---|
| 33 | include 'fldcom.h' |
|---|
| 34 | include 'misccom.h' |
|---|
| 35 | include 'ncordscom.h' |
|---|
| 36 | include 'pipecom.h' |
|---|
| 37 | include 'psizescom.h' |
|---|
| 38 | include 'sccom.h' |
|---|
| 39 | include 'syscom.h' |
|---|
| 40 | include 'ucom.h' |
|---|
| 41 | c |
|---|
| 42 | parameter (maxdim=84000) |
|---|
| 43 | logical fexist,fopen,endfil |
|---|
| 44 | common/ehparm/xi,yi,zi,ex,ey,eez,hx,hy,hz,xj,yj,zj,gbxj,gbyj, |
|---|
| 45 | . gbzj,gamma,qfac |
|---|
| 46 | common/intype/intype |
|---|
| 47 | common/out6/nbphase,nbremesh |
|---|
| 48 | real xpipe(imaa),ypipe(imaa) |
|---|
| 49 | real*4 dbgx(imaa),dbgy(imaa),dbgz(imaa) |
|---|
| 50 | integer*4 ipipe(imaa) |
|---|
| 51 | dimension qol(200) |
|---|
| 52 | dimension rsce(nsce) |
|---|
| 53 | c 4221=(20+1)*(200+1), 84000=20*(20+1)*200 |
|---|
| 54 | dimension rm(21), zm(201), ers(maxdim), ezs(maxdim), er(4221), |
|---|
| 55 | 1 ez(4221), aa(4000), ismax(201), iemax(201),rssq(20),zzs(201) |
|---|
| 56 | c-------------------------------------------------------------------------- |
|---|
| 57 | c* |
|---|
| 58 | data irout,izout/0,0/fopen,endfil/.false.,.false./ |
|---|
| 59 | c |
|---|
| 60 | if (ngood.le.1) return |
|---|
| 61 | c skip if want point-by-point space charge calc. |
|---|
| 62 | if (sce(3).ne.0.) go to 400 |
|---|
| 63 | if (dist) 50,10,50 |
|---|
| 64 | 10 continue |
|---|
| 65 | write(nnout,*) ' Space charge subroutine number ',iprog |
|---|
| 66 | call rtrans(gt,bt) |
|---|
| 67 | 15 continue |
|---|
| 68 | gmesh=gt |
|---|
| 69 | c set up field tables |
|---|
| 70 | 16 continue |
|---|
| 71 | beami=sce(1) |
|---|
| 72 | rmax=sce(4) |
|---|
| 73 | dzmax=sce(5)*gmesh |
|---|
| 74 | nr=sce(6) |
|---|
| 75 | nz=sce(7) |
|---|
| 76 | nr1=nr+1 |
|---|
| 77 | nz1=nz+1 |
|---|
| 78 | im1=nr*nz |
|---|
| 79 | im2=nr1*nz1 |
|---|
| 80 | im3=nr1*nz |
|---|
| 81 | im4=im1*nr1 |
|---|
| 82 | nip=sce(8) |
|---|
| 83 | inip = nip |
|---|
| 84 | opt=sce(10) |
|---|
| 85 | iopt = opt |
|---|
| 86 | if(opt.ge.2)then |
|---|
| 87 | dr=rmax*rmax/sce(6) |
|---|
| 88 | else |
|---|
| 89 | dr=rmax/sce(6) |
|---|
| 90 | endif |
|---|
| 91 | dz=dzmax/sce(7) |
|---|
| 92 | pl=sce(9)*gmesh |
|---|
| 93 | if(pl.le.0.)pl=dzmax |
|---|
| 94 | frm=sce(11) |
|---|
| 95 | if(frm.le.0.)frm=1.5 |
|---|
| 96 | rwall=sce(12) |
|---|
| 97 | if(nr.gt.20 .or. nz.gt.200 .or. (nr*nr1*nz.gt.maxdim .and. |
|---|
| 98 | *(maxdim/nr1**2*dz.lt.10*rwall .or. maxdim/nr1**2*dz.lt.20*rmax))) |
|---|
| 99 | *then |
|---|
| 100 | call appendparm |
|---|
| 101 | stop' Abnormal stop too many mesh points in space charge grid ' |
|---|
| 102 | endif |
|---|
| 103 | c load rm, zm, rs, zs |
|---|
| 104 | rm(1)=0.0 |
|---|
| 105 | rm(nr1)=rmax |
|---|
| 106 | do 20 i=2,nr1 |
|---|
| 107 | if(opt.ge.2)then |
|---|
| 108 | rm(i)=sqrt(float(i-1)*dr) |
|---|
| 109 | else |
|---|
| 110 | rm(i)=float(i-1)*dr |
|---|
| 111 | endif |
|---|
| 112 | rssq(i-1)=.5*(rm(i-1)**2+rm(i)**2) |
|---|
| 113 | c rs(i-1)=sqrt(rssq(i-1)) |
|---|
| 114 | 20 continue |
|---|
| 115 | zs=.5*dz |
|---|
| 116 | do 30 i=1,nz1 |
|---|
| 117 | zm(i)=float(i-1)*dz |
|---|
| 118 | zzs(i)=zm(i)+zs |
|---|
| 119 | 30 continue |
|---|
| 120 | hl=.5*dzmax |
|---|
| 121 | c---if there is a reference particle, set n1=2; otherwise, n1=1 |
|---|
| 122 | n1=1 |
|---|
| 123 | c--- if (w0.gt.0.) n1=2 |
|---|
| 124 | xnp=npoints+1-n1 |
|---|
| 125 | c load ers and ezs |
|---|
| 126 | c mesh dimensions are in cm. ers and ezs are in 1/cm. |
|---|
| 127 | c c1, c2 and c3 are in cm., and c4 is in mev-cm. |
|---|
| 128 | c q=coulombs/point. 1/epsilon_0 =11.294e6 cm mev/coul. |
|---|
| 129 | q=beami/(freq*1.e6*xnp) |
|---|
| 130 | c when beami < 0, it is -(no. of particles in the bunch) |
|---|
| 131 | if (beami.lt.0.) then |
|---|
| 132 | q = -1.602e-19*beami/npoints |
|---|
| 133 | endif |
|---|
| 134 | c printout during initialization |
|---|
| 135 | if((ip.eq.1.and.dist.eq.0.) .or. (ip.eq.2.and.dist.eq.0.))then |
|---|
| 136 | write(nnout,300) |
|---|
| 137 | 300 format (' distribute charge over rings for space charge calc.') |
|---|
| 138 | if (beami.gt.0.) write(nnout,310) beami |
|---|
| 139 | 310 format (' beam current =',t40,1pg12.3,' amps') |
|---|
| 140 | abeami = abs(beami) |
|---|
| 141 | if (beami.lt.0.) write(nnout,320) abeami |
|---|
| 142 | 320 format (' no. of electrons in bunch =',t40,1pg12.3) |
|---|
| 143 | write(nnout,330) q,rmax,sce(5),nr,nz,inip,sce(9),iopt,frm,rwall |
|---|
| 144 | 330 format ( |
|---|
| 145 | . ' charge per superparticle =',t40,1pg12.3,' coulombs'/ |
|---|
| 146 | . ' radial size of mesh =',t40,g12.3,' cm'/ |
|---|
| 147 | . ' z length of mesh =',t40,g12.3,' cm'/ |
|---|
| 148 | . ' no. of radial mesh points =',t40,i8/ |
|---|
| 149 | . ' no. of z mesh points =',t40,i8/ |
|---|
| 150 | . ' no. of "identical pulses" =',t40,i8/ |
|---|
| 151 | . ' pulse length =',t40,g12.3/ |
|---|
| 152 | . ' OPT =',t40,i8/ |
|---|
| 153 | . ' refresh when energy grows by factor',t40,g12.3/ |
|---|
| 154 | . ' wall radius for images =',t40,g12.3) |
|---|
| 155 | endif |
|---|
| 156 | c---- 1/epsilon_0 verification 22/02/94 |
|---|
| 157 | c1=11.294e6*q/erest |
|---|
| 158 | call inscgrid(fexist) |
|---|
| 159 | if(.not.fopen)then |
|---|
| 160 | open(unit=nschef,file='scgrid',access='sequential', |
|---|
| 161 | * status='unknown',form='unformatted') |
|---|
| 162 | fopen=.true. |
|---|
| 163 | endif |
|---|
| 164 | j=1 |
|---|
| 165 | c if-construct change (jump into if-block not allowed) 04/94 |
|---|
| 166 | if(endfil)then |
|---|
| 167 | backspace(nschef) |
|---|
| 168 | write(ndiag,*)'not using stored data rdz,dz,i,rsce(i),sce(i)', |
|---|
| 169 | * rdz,dz,j,rsce(j),sce(j) |
|---|
| 170 | else if(fexist)then |
|---|
| 171 | read(nschef,err=34,end=33)rdz,rc1,rsce,ers,ezs |
|---|
| 172 | if(abs(rdz-dz).gt.0.00001*dz)go to 32 |
|---|
| 173 | do 31 i=4,12 |
|---|
| 174 | j=i |
|---|
| 175 | if(sce(i).ne.rsce(i))go to 32 |
|---|
| 176 | 31 continue |
|---|
| 177 | c stored data good |
|---|
| 178 | if(rc1.ne.c1)then |
|---|
| 179 | do 35 i=1,im4 |
|---|
| 180 | ers(i)=ers(i)*c1/rc1 |
|---|
| 181 | ezs(i)=ezs(i)*c1/rc1 |
|---|
| 182 | 35 continue |
|---|
| 183 | endif |
|---|
| 184 | write(ndiag,*)'using stored scheff mesh data' |
|---|
| 185 | go to 41 |
|---|
| 186 | 33 continue |
|---|
| 187 | write(ndiag,*)'end of file' |
|---|
| 188 | endfil=.true. |
|---|
| 189 | go to 32 |
|---|
| 190 | 34 continue |
|---|
| 191 | write(ndiag,*)'error on reading file scgrid' |
|---|
| 192 | endfil=.true. |
|---|
| 193 | c stored data no good |
|---|
| 194 | 32 continue |
|---|
| 195 | backspace(nschef) |
|---|
| 196 | write(ndiag,*)'not using stored data rdz,dz,i,rsce(i),sce(i)', |
|---|
| 197 | * rdz,dz,j,rsce(j),sce(j) |
|---|
| 198 | endif |
|---|
| 199 | l=0 |
|---|
| 200 | do 40 k=1,nr |
|---|
| 201 | do 40 j=1,nz |
|---|
| 202 | do 40 i=1,nr1 |
|---|
| 203 | rp=rm(i) |
|---|
| 204 | zp=zm(j+1) |
|---|
| 205 | call gaus(rm(k),rm(k+1),zm(1),zm(2),opt,er1,ez1) |
|---|
| 206 | l=l+1 |
|---|
| 207 | ers(l)=c1*er1 |
|---|
| 208 | ezs(l)=c1*ez1 |
|---|
| 209 | 40 continue |
|---|
| 210 | write(nschef)dz,c1,sce,ers,ezs |
|---|
| 211 | endfile(nschef) |
|---|
| 212 | 41 continue |
|---|
| 213 | if (dist.eq.0.)return |
|---|
| 214 | go to 60 |
|---|
| 215 | c evaluate and apply space charge effects. |
|---|
| 216 | 50 continue |
|---|
| 217 | call rtrans(gt,bt) |
|---|
| 218 | if(gmesh.le.0.)go to 15 |
|---|
| 219 | if(gt/gmesh.gt.frm)then |
|---|
| 220 | gmesh=gmesh*frm |
|---|
| 221 | go to 16 |
|---|
| 222 | endif |
|---|
| 223 | 60 continue |
|---|
| 224 | zc=cord(5,1) |
|---|
| 225 | c3=dist/gt |
|---|
| 226 | c4=c3 |
|---|
| 227 | c evaluate ng, xbar, ybar, and zbar. |
|---|
| 228 | eng=ngood+1-n1 |
|---|
| 229 | xbar=0. |
|---|
| 230 | ybar=0. |
|---|
| 231 | xsq=0. |
|---|
| 232 | ysq=0. |
|---|
| 233 | zbar=0. |
|---|
| 234 | do 90 np=n1,ngood |
|---|
| 235 | c skip if particle not yet emitted |
|---|
| 236 | nenp = cord(7,np) |
|---|
| 237 | if (nenp.eq.0.and.intype.eq.10) go to 90 |
|---|
| 238 | if (nenp.eq.0.and.intype.eq.11) go to 90 |
|---|
| 239 | if (nenp.eq.0.and.intype.eq.15) go to 90 |
|---|
| 240 | gstar=gt*(gam(np)-bt*cord(6,np)) |
|---|
| 241 | bgstar=gt*(cord(6,np)-bt*gam(np)) |
|---|
| 242 | c gstar=sqrt(1.+cord(2,np)**2+cord(4,np)**2+bgstar**2) |
|---|
| 243 | zstar=gt*(cord(5,np)-zc) |
|---|
| 244 | gam(np)=gstar |
|---|
| 245 | cord(6,np)=bgstar |
|---|
| 246 | c xstar=cord(1,np)+cord(2,np)*bt*zstar/gstar |
|---|
| 247 | c ystar=cord(3,np)+cord(4,np)*bt*zstar/gstar |
|---|
| 248 | xstar=cord(1,np) |
|---|
| 249 | ystar=cord(3,np) |
|---|
| 250 | xbar=xbar+xstar |
|---|
| 251 | ybar=ybar+ystar |
|---|
| 252 | xsq=xsq+xstar**2 |
|---|
| 253 | ysq=ysq+ystar**2 |
|---|
| 254 | c 80 zbar=zbar+zstar*(1.+bt*cord(6,np)/gstar) |
|---|
| 255 | 80 zbar=zbar+zstar |
|---|
| 256 | 90 continue |
|---|
| 257 | xbar=xbar/eng |
|---|
| 258 | ybar=ybar/eng |
|---|
| 259 | xsq=xsq/eng |
|---|
| 260 | ysq=ysq/eng |
|---|
| 261 | c epsq=sqrt((xsq-xbar**2)/(ysq-ybar**2)) |
|---|
| 262 | epsq=1. |
|---|
| 263 | c repsq=1./epsq |
|---|
| 264 | repsq=1. |
|---|
| 265 | c xfac=2./(1.+epsq) |
|---|
| 266 | xfac=1. |
|---|
| 267 | c yfac=epsq*xfac |
|---|
| 268 | yfac=1. |
|---|
| 269 | zbar=zbar/eng |
|---|
| 270 | zc=zc+dzmax/nz*(sandom(dum)-.5) |
|---|
| 271 | c clear and load bins |
|---|
| 272 | do 100 i=1,im1 |
|---|
| 273 | aa(i)=0 |
|---|
| 274 | 100 continue |
|---|
| 275 | ng=0 |
|---|
| 276 | n108=0 |
|---|
| 277 | do 110 np=2,ngood |
|---|
| 278 | c skip if particle not yet emitted |
|---|
| 279 | nenp = cord(7,np) |
|---|
| 280 | if (cord(5,np).le.0.) go to 110 |
|---|
| 281 | if (nenp.eq.0.and.intype.eq.10) go to 110 |
|---|
| 282 | if (nenp.eq.0.and.intype.eq.11) go to 110 |
|---|
| 283 | if (nenp.eq.0.and.intype.eq.15) go to 110 |
|---|
| 284 | zstar=gt*(cord(5,np)-zc) |
|---|
| 285 | c xstar=cord(1,np)+bt*zstar*cord(2,np)/gam(np) |
|---|
| 286 | c ystar=cord(3,np)+bt*zstar*cord(4,np)/gam(np) |
|---|
| 287 | xstar=cord(1,np) |
|---|
| 288 | ystar=cord(3,np) |
|---|
| 289 | rsq=repsq*(xstar-xbar)**2+epsq*(ystar-ybar)**2 |
|---|
| 290 | if(opt.ge.2)then |
|---|
| 291 | i=rsq/dr+1. |
|---|
| 292 | else |
|---|
| 293 | i=sqrt(rsq)/dr+1. |
|---|
| 294 | endif |
|---|
| 295 | if (i.gt.nr) go to 106 |
|---|
| 296 | c z=zstar*(1.+bt*cord(6,np)/gam(np)) |
|---|
| 297 | z=zstar |
|---|
| 298 | if (abs(z).le.hl) go to 105 |
|---|
| 299 | if(nip.eq.0)go to 108 |
|---|
| 300 | if (pl.ne.dzmax) go to 108 |
|---|
| 301 | z=z-sign(pl,z) |
|---|
| 302 | c------distribute charge among adjacent bins. |
|---|
| 303 | 105 continue |
|---|
| 304 | ng=ng+1 |
|---|
| 305 | zz=z+hl |
|---|
| 306 | jm1=zz/dz+1. |
|---|
| 307 | i1=i+1 |
|---|
| 308 | if(rsq.lt.rssq(i)) i1=i-1 |
|---|
| 309 | if(i1.lt.1)i1=1 |
|---|
| 310 | if(i1.gt.nr) i1=nr |
|---|
| 311 | j1=jm1+1 |
|---|
| 312 | if(zz.lt.zzs(jm1)) j1=jm1-1 |
|---|
| 313 | if(j1.lt.1)j1=nz |
|---|
| 314 | if(j1.gt.nz)j1=1 |
|---|
| 315 | a=1. |
|---|
| 316 | if(i1.ne.i)a=(rsq-rssq(i1))/(rssq(i)-rssq(i1)) |
|---|
| 317 | b=1.-a |
|---|
| 318 | cc=1. |
|---|
| 319 | if(j1.ne.jm1) cc=(zz-zzs(j1))/(zzs(jm1)-zzs(j1)) |
|---|
| 320 | d=1.-cc |
|---|
| 321 | k=(jm1-1)*nr+i |
|---|
| 322 | aa(k)=aa(k)+a*cc |
|---|
| 323 | k=k+i1-i |
|---|
| 324 | aa(k)=aa(k)+b*cc |
|---|
| 325 | k=(j1-1)*nr+i |
|---|
| 326 | aa(k)=aa(k)+a*d |
|---|
| 327 | k=k+i1-i |
|---|
| 328 | aa(k)=aa(k)+b*d |
|---|
| 329 | go to 110 |
|---|
| 330 | 106 continue |
|---|
| 331 | if(irout.ne.0 .or. cord(5,np).lt.0.)go to 110 |
|---|
| 332 | write(nnout,107) np,xstar,ystar,nenp |
|---|
| 333 | 107 format(' warning--particle no.',i4,' is outside radial mesh'/ |
|---|
| 334 | * ' xstar=',f8.5,', ystar=',f8.5,' in element ',i4) |
|---|
| 335 | irout=1 |
|---|
| 336 | go to 110 |
|---|
| 337 | 108 continue |
|---|
| 338 | n108=n108+1 |
|---|
| 339 | if(izout.ne.0 .or. cord(5,np).lt.0.)go to 110 |
|---|
| 340 | write(nnout,109) np,zstar,cord(5,np),nenp |
|---|
| 341 | 109 format(' warning--particle no.',i4,' is outside longitudinal mesh' |
|---|
| 342 | * / ' zstar=',f10.5,', z=',f10.5, ' in element ',i4) |
|---|
| 343 | izout=1 |
|---|
| 344 | 110 continue |
|---|
| 345 | eng=ng |
|---|
| 346 | c---calculate no. of particles in each column |
|---|
| 347 | do 115 j=1,nz |
|---|
| 348 | qol(j)=0. |
|---|
| 349 | l=(j-1)*nr |
|---|
| 350 | do 114 i=1,nr |
|---|
| 351 | k=l+i |
|---|
| 352 | qol(j)=qol(j)+aa(k) |
|---|
| 353 | 114 continue |
|---|
| 354 | qol(j)=qol(j)/dz |
|---|
| 355 | 115 continue |
|---|
| 356 | c find ismax for each j |
|---|
| 357 | do 130 j=1,nz |
|---|
| 358 | l=(j-1)*nr |
|---|
| 359 | k=nr |
|---|
| 360 | do 120 i=1,nr |
|---|
| 361 | m=l+k |
|---|
| 362 | if (aa(m)) 121,121,131 |
|---|
| 363 | 121 continue |
|---|
| 364 | k=k-1 |
|---|
| 365 | 120 continue |
|---|
| 366 | 131 continue |
|---|
| 367 | ismax(j)=k |
|---|
| 368 | 130 continue |
|---|
| 369 | c find iemax for each j |
|---|
| 370 | iemax(1)=1+ismax(1) |
|---|
| 371 | do 140 j=2,nz |
|---|
| 372 | iemax(j)=1+max0(ismax(j-1),ismax(j)) |
|---|
| 373 | 140 continue |
|---|
| 374 | iemax(nz1)=1+ismax(nz) |
|---|
| 375 | c set er and ez to zero |
|---|
| 376 | do 150 i=1,im2 |
|---|
| 377 | er(i)=0.0 |
|---|
| 378 | ez(i)=0.0 |
|---|
| 379 | 150 continue |
|---|
| 380 | c sum up fields |
|---|
| 381 | c---------------------------------------------------------------------- |
|---|
| 382 | c the following 'do nothing line' actually causes the compiler to insert |
|---|
| 383 | c FINIT instruction and a FLDCW instruction to reset the 80287. |
|---|
| 384 | if(ez(2).gt.1.)m1=m1 |
|---|
| 385 | c---------------------------------------------------------------------- |
|---|
| 386 | do 200 js=1,nz |
|---|
| 387 | js1=js+1 |
|---|
| 388 | ism=ismax(js) |
|---|
| 389 | if(ism.eq.0) go to 200 |
|---|
| 390 | 160 continue |
|---|
| 391 | do 190 is=1,ism |
|---|
| 392 | l=(js-1)*nr+is |
|---|
| 393 | a1=aa(l) |
|---|
| 394 | if (a1.eq.0.) go to 190 |
|---|
| 395 | l=(is-1)*im3 |
|---|
| 396 | do 170 je=1,js |
|---|
| 397 | k1=l+(js-je)*nr1 |
|---|
| 398 | m1=(je-1)*nr1 |
|---|
| 399 | iem=iemax(je) |
|---|
| 400 | if(iem.le.1)goto 170 |
|---|
| 401 | do 168 ie=1,iem |
|---|
| 402 | c n=m1+ie |
|---|
| 403 | c k=k1+ie |
|---|
| 404 | c ------------------------------------------------------------------------- |
|---|
| 405 | c The following instructions us the 80287 and it must be reset before |
|---|
| 406 | c they work right. However, they rarely fail without the 80287 being reset! |
|---|
| 407 | er(m1+ie)=er(m1+ie)+a1*ers(k1+ie) |
|---|
| 408 | ez(m1+ie)=ez(m1+ie)-a1*ezs(k1+ie) |
|---|
| 409 | c -------------------------------------------------------------------------- |
|---|
| 410 | 168 continue |
|---|
| 411 | 170 continue |
|---|
| 412 | do 180 je=js1,nz1 |
|---|
| 413 | k1=l+(je-js1)*nr1 |
|---|
| 414 | m1=(je-1)*nr1 |
|---|
| 415 | iem=iemax(je) |
|---|
| 416 | if(iem.le.1) go to 180 |
|---|
| 417 | do 178 ie=1,iem |
|---|
| 418 | c n=m1+ie |
|---|
| 419 | c k=k1+ie |
|---|
| 420 | er(m1+ie)=er(m1+ie)+a1*ers(k1+ie) |
|---|
| 421 | ez(m1+ie)=ez(m1+ie)+a1*ezs(k1+ie) |
|---|
| 422 | 178 continue |
|---|
| 423 | 180 continue |
|---|
| 424 | 190 continue |
|---|
| 425 | 200 continue |
|---|
| 426 | c evaluate and apply impulse |
|---|
| 427 | do 240 np=n1,ngood |
|---|
| 428 | c skip if particle not yet emitted |
|---|
| 429 | nenp = cord(7,np) |
|---|
| 430 | if (nenp.eq.0.and.intype.eq.10) go to 240 |
|---|
| 431 | if (nenp.eq.0.and.intype.eq.11) go to 240 |
|---|
| 432 | if (nenp.eq.0.and.intype.eq.15) go to 240 |
|---|
| 433 | bgfstr=cord(6,np) |
|---|
| 434 | if(cord(5,np).le.0.) then |
|---|
| 435 | go to 235 |
|---|
| 436 | else |
|---|
| 437 | endif |
|---|
| 438 | zstar=gt*(cord(5,np)-zc) |
|---|
| 439 | c xstar=cord(1,np)+bt*zstar*cord(2,np)/gam(np) |
|---|
| 440 | c ystar=cord(3,np)+bt*zstar*cord(4,np)/gam(np) |
|---|
| 441 | xstar=cord(1,np) |
|---|
| 442 | ystar=cord(3,np) |
|---|
| 443 | r=sqrt(repsq*(xstar-xbar)**2+epsq*(ystar-ybar)**2) |
|---|
| 444 | if(np.eq.1) go to 201 |
|---|
| 445 | if(r.lt.0.000001) r=.000001 |
|---|
| 446 | xor=(xstar-xbar)*xfac/r |
|---|
| 447 | yor=(ystar-ybar)*yfac/r |
|---|
| 448 | go to 203 |
|---|
| 449 | 201 continue |
|---|
| 450 | xor=0. |
|---|
| 451 | yor=0. |
|---|
| 452 | c z=zstar*(1.+bt*cord(6,np)/gam(np)) |
|---|
| 453 | 203 continue |
|---|
| 454 | z=zstar |
|---|
| 455 | 205 continue |
|---|
| 456 | if(abs(z).le.hl) go to 210 |
|---|
| 457 | if (nip.le.0) go to 235 |
|---|
| 458 | if (pl.ne.dzmax) go to 235 |
|---|
| 459 | z=z-sign(pl,z) |
|---|
| 460 | go to 205 |
|---|
| 461 | c interpolate impulse within mesh. |
|---|
| 462 | c 210 if(r.gt.rmax) go to 220 |
|---|
| 463 | 210 continue |
|---|
| 464 | if(opt.ge.2.and.r.le.rmax)then |
|---|
| 465 | rb=sqrt(r*r/dr) |
|---|
| 466 | i=1.+rb |
|---|
| 467 | a=(r-rm(i))/(rm(i+1)-rm(i)) |
|---|
| 468 | else |
|---|
| 469 | rb=r/dr |
|---|
| 470 | i=1.0+rb |
|---|
| 471 | a=rb-float(i-1) |
|---|
| 472 | endif |
|---|
| 473 | if(r.gt.rmax)then |
|---|
| 474 | a=(r-rmax)/(rwall-rmax) |
|---|
| 475 | i=nr1 |
|---|
| 476 | endif |
|---|
| 477 | b=1.0-a |
|---|
| 478 | zb=(z+hl)/dz |
|---|
| 479 | j=1.0+zb |
|---|
| 480 | c=zb-float(j-1) |
|---|
| 481 | d=1.0-c |
|---|
| 482 | l=i+(j-1)*nr1 |
|---|
| 483 | m=l+nr1 |
|---|
| 484 | if(r.gt.rmax)then |
|---|
| 485 | c the following two lines assume all the charge is within rmax. |
|---|
| 486 | cbgr=c1*c3*qol(j)/(2.*pi*r**2) |
|---|
| 487 | cbgz=c4*(d*b*ez(l)+c*b*ez(m)) |
|---|
| 488 | else |
|---|
| 489 | cbgr=c3*(d*(a*er(l+1)+b*er(l))+c*(a*er(m+1)+b*er(m))) |
|---|
| 490 | cbgz=c4*(d*(a*ez(l+1)+b*ez(l))+c*(a*ez(m+1)+b*ez(m))) |
|---|
| 491 | endif |
|---|
| 492 | go to 230 |
|---|
| 493 | c estimate impulse based on point charge at xbar,ybar,zbar. |
|---|
| 494 | c220 z=(zstar-zbar)*(1.+bt*cord(6,np)/gam(np)) |
|---|
| 495 | c---the following lines are for a point charge at beam centroid. |
|---|
| 496 | c220 z=(zstar-zbar) |
|---|
| 497 | c d=sqrt(z**2+r**2) |
|---|
| 498 | c rod3=r/d**3 |
|---|
| 499 | c zod3=z/d**3 |
|---|
| 500 | c if(nip.eq.0) go to 228 |
|---|
| 501 | c include neighboring bunches. |
|---|
| 502 | c do 226 i=1,nip |
|---|
| 503 | c xi=i |
|---|
| 504 | c do 226 j=1,2 |
|---|
| 505 | c s=z+xi*pl |
|---|
| 506 | c d=sqrt(s**2+r**2) |
|---|
| 507 | c rod3=rod3+r/d**3 |
|---|
| 508 | c zod3=zod3+s/d**3 |
|---|
| 509 | c 226 xi=-xi |
|---|
| 510 | c evaluate impulse. |
|---|
| 511 | c 228 cbgr=eng*c1*c3*rod3/(4.*pi) |
|---|
| 512 | c cbgz=eng*c1*c4*zod3/(4.*pi) |
|---|
| 513 | c---following two lines are for line charge with density qol(j) |
|---|
| 514 | c 220 cbgz=0. |
|---|
| 515 | c cbgr=c1*c3*qol(j)/(2.*pi*r**2) |
|---|
| 516 | c apply impulse |
|---|
| 517 | 230 continue |
|---|
| 518 | cord(2,np)=cord(2,np)+cbgr*xor |
|---|
| 519 | cord(4,np)=cord(4,np)+cbgr*yor |
|---|
| 520 | bgfstr=cord(6,np)+cbgz |
|---|
| 521 | 235 continue |
|---|
| 522 | gfstar=sqrt(1.+cord(2,np)**2+cord(4,np)**2+bgfstr**2) |
|---|
| 523 | cord(6,np)=gt*(bgfstr+bt*gfstar) |
|---|
| 524 | gam(np)=gt*(gfstar+bt*bgfstr) |
|---|
| 525 | c gam(np)=sqrt(1.+cord(2,np)**2+cord(4,np)**2+cord(6,np)**2) |
|---|
| 526 | 240 continue |
|---|
| 527 | return |
|---|
| 528 | c |
|---|
| 529 | c point-by-point calculation of the space charge K.T. McDonald |
|---|
| 530 | c To invoke this option, set parmeter sce(3) nonzero |
|---|
| 531 | c If sce(3) < 0, calculate image fields due to a metallic surface |
|---|
| 532 | c at z = 0 |
|---|
| 533 | 400 continue |
|---|
| 534 | c initialization |
|---|
| 535 | if (dist.eq.0.) then |
|---|
| 536 | c q = no. of electrons represented by each simulated particle |
|---|
| 537 | xnp = npoints |
|---|
| 538 | q = beami/(freq*1.e6*xnp)/1.602e-19 |
|---|
| 539 | if (beami.lt.0.) then |
|---|
| 540 | q = -beami/xnp |
|---|
| 541 | endif |
|---|
| 542 | do 405 i = 1,ngood |
|---|
| 543 | ipipe(i) = 0 |
|---|
| 544 | 405 continue |
|---|
| 545 | pvfac = abs(sce(3)) |
|---|
| 546 | write(nnout,*) ' Space charge subroutine number ',iprog |
|---|
| 547 | if (ip.eq.1 .or. ip.eq.2) then |
|---|
| 548 | write(nnout,410) q |
|---|
| 549 | 410 format (' point to point space charge calculation'/ |
|---|
| 550 | . ' q = ',1pg11.3,' electrons per superparticle') |
|---|
| 551 | abeami = abs(beami) |
|---|
| 552 | if (beami.gt.0.) write(nnout,420) abeami |
|---|
| 553 | 420 format (' beam current = ',1pg12.3,' amps') |
|---|
| 554 | if (beami.le.0.) write(nnout,430) abeami |
|---|
| 555 | 430 format (' no. of electrons in bunch =',1pg12.3) |
|---|
| 556 | if (sce(3).lt.0.) write(nnout,440) |
|---|
| 557 | 440 format (' include images in the plane z = 0') |
|---|
| 558 | write(nnout,450) pvfac |
|---|
| 559 | 450 format (' particle size factor =',t40,1pg12.3) |
|---|
| 560 | endif |
|---|
| 561 | pvfac = pvfac**2 |
|---|
| 562 | return |
|---|
| 563 | endif |
|---|
| 564 | c----r0 = e / ((m0c**2/e)*(4*pi*epsilon0)) in cm |
|---|
| 565 | r0 = 2.818e-13 |
|---|
| 566 | c----integration factor = r0 * wavelength * dphi(radians) /2 pi |
|---|
| 567 | rfac = r0*dist |
|---|
| 568 | c----qrfac includes the integration factor |
|---|
| 569 | qrfac = q*rfac |
|---|
| 570 | if (ngood.lt.2) return |
|---|
| 571 | c |
|---|
| 572 | c setup for image calc. in pipe walls |
|---|
| 573 | if (npipe.eq.0) go to 520 |
|---|
| 574 | i = 1 |
|---|
| 575 | 500 continue |
|---|
| 576 | if (i.gt.ngood) then |
|---|
| 577 | go to 520 |
|---|
| 578 | else |
|---|
| 579 | c find which pipe the particle is in |
|---|
| 580 | zi = cord(5,i) |
|---|
| 581 | jmin = 1 |
|---|
| 582 | if (ipipe(i).gt.0) jmin = ipipe(i) |
|---|
| 583 | do 510 j = jmin,npipe |
|---|
| 584 | if (zi.ge.zlpipe(j).and.zi.le.zhpipe(j)) go to 515 |
|---|
| 585 | 510 continue |
|---|
| 586 | j = 0 |
|---|
| 587 | 515 continue |
|---|
| 588 | ipipe(i) = j |
|---|
| 589 | if (j.gt.0) then |
|---|
| 590 | xi = cord(1,i) |
|---|
| 591 | yi = cord(3,i) |
|---|
| 592 | risq = xi**2 + yi**2 |
|---|
| 593 | rpsq = rpipe(j)**2 |
|---|
| 594 | if (risq.ge.rpsq) then |
|---|
| 595 | c particle is lost if outside the pipe |
|---|
| 596 | c dum is a dummy argument because in pardyn with have to send wtp in swap |
|---|
| 597 | call swap(i,dum) |
|---|
| 598 | i = i - 1 |
|---|
| 599 | elseif (risq.eq.0.) then |
|---|
| 600 | c skip if particle is on axis |
|---|
| 601 | ipipe(i) = 0 |
|---|
| 602 | else |
|---|
| 603 | c store position of image charge |
|---|
| 604 | rnew = rpsq/risq |
|---|
| 605 | xpipe(i) = xi*rnew |
|---|
| 606 | ypipe(i) = yi*rnew |
|---|
| 607 | endif |
|---|
| 608 | endif |
|---|
| 609 | i = i + 1 |
|---|
| 610 | endif |
|---|
| 611 | go to 500 |
|---|
| 612 | 520 continue |
|---|
| 613 | c |
|---|
| 614 | c |
|---|
| 615 | c find volume of a typical superparticle. If a second particle |
|---|
| 616 | c lies within the volume of the first, reduce the effective |
|---|
| 617 | c charge in subroutine eh |
|---|
| 618 | xsqsiz = 0. |
|---|
| 619 | ysqsiz = 0. |
|---|
| 620 | zsqsiz = 0. |
|---|
| 621 | xyzvol = 1. |
|---|
| 622 | if (abs(sce(3)).ne.999.) call pvoldp |
|---|
| 623 | c calc fields due to rest of bunch, including image fields |
|---|
| 624 | c but don't apply space charge fields on the reference particle |
|---|
| 625 | do 1300 i = 2,ngood |
|---|
| 626 | nei = cord(7,i) |
|---|
| 627 | c skip if particle i has not yet been emitted |
|---|
| 628 | if (nei.eq.0.and.intype.eq.10) go to 1300 |
|---|
| 629 | if (nei.eq.0.and.intype.eq.11) go to 1300 |
|---|
| 630 | if (nei.eq.0.and.intype.eq.15) go to 1300 |
|---|
| 631 | xi = cord(1,i) |
|---|
| 632 | yi = cord(3,i) |
|---|
| 633 | zi = cord(5,i) |
|---|
| 634 | ex = 0. |
|---|
| 635 | ey = 0. |
|---|
| 636 | c name -ez- already in use |
|---|
| 637 | eez = 0. |
|---|
| 638 | hx = 0. |
|---|
| 639 | hy = 0. |
|---|
| 640 | hz = 0. |
|---|
| 641 | do 1200 j = 1,ngood |
|---|
| 642 | nej = cord(7,j) |
|---|
| 643 | c skip if particle j has not yet been emitted |
|---|
| 644 | if (nej.eq.0.and.intype.eq.10) go to 1200 |
|---|
| 645 | if (nej.eq.0.and.intype.eq.11) go to 1200 |
|---|
| 646 | if (nej.eq.0.and.intype.eq.15) go to 1200 |
|---|
| 647 | xj = cord(1,j) |
|---|
| 648 | yj = cord(3,j) |
|---|
| 649 | zj = cord(5,j) |
|---|
| 650 | gamma = gam(j) |
|---|
| 651 | gbxj = cord(2,j) |
|---|
| 652 | gbyj = cord(4,j) |
|---|
| 653 | gbzj = cord(6,j) |
|---|
| 654 | c skip to image charge calc. if i = j |
|---|
| 655 | if (i.eq.j) go to 1250 |
|---|
| 656 | c calc. field at i of moving charge j |
|---|
| 657 | qfac = qrfac |
|---|
| 658 | call eh |
|---|
| 659 | c |
|---|
| 660 | c now do the image fields in the plane z = 0 |
|---|
| 661 | 1250 continue |
|---|
| 662 | c skip unless sce(3) < 0, and particle is in 1st element |
|---|
| 663 | if (sce(3).ge.0.0 .or. nej.ne.1) go to 1260 |
|---|
| 664 | c skip image calculation once particle is more than 1 cm from cathode |
|---|
| 665 | if (zj.gt.1.) go to 1260 |
|---|
| 666 | zj = -zj |
|---|
| 667 | gbzj = -gbzj |
|---|
| 668 | c image charge has opposite sign |
|---|
| 669 | qfac = -qrfac |
|---|
| 670 | c for image of particle i, only use 1 electron charge, and not the |
|---|
| 671 | c supercharge q |
|---|
| 672 | if (i.eq.j) qfac = -rfac |
|---|
| 673 | call eh |
|---|
| 674 | c restore coords |
|---|
| 675 | zj = -zj |
|---|
| 676 | gbzj = -gbzj |
|---|
| 677 | c |
|---|
| 678 | c calc. image fields in beam pipe, if any |
|---|
| 679 | 1260 continue |
|---|
| 680 | if (ipipe(j).eq.0) go to 1200 |
|---|
| 681 | xj = xpipe(j) |
|---|
| 682 | yj = ypipe(j) |
|---|
| 683 | c don't change magnitude of transverse velocity, thus keeping beta < 1 |
|---|
| 684 | c (if revise this, must recalc. gamma also) |
|---|
| 685 | gbxj = -gbxj |
|---|
| 686 | gbyj = -gbyj |
|---|
| 687 | c image charge has opposite sign |
|---|
| 688 | qfac = -qrfac |
|---|
| 689 | c for image of particle i, only use 1 electron charge, and not the |
|---|
| 690 | c supercharge q |
|---|
| 691 | if (i.eq.j) qfac = -rfac |
|---|
| 692 | call eh |
|---|
| 693 | 1200 continue |
|---|
| 694 | c calculate and store the impulse on particle i |
|---|
| 695 | gamma = gam(i) |
|---|
| 696 | bxi = cord(2,i)/gamma |
|---|
| 697 | byi = cord(4,i)/gamma |
|---|
| 698 | bzi = cord(6,i)/gamma |
|---|
| 699 | dbgx(i) = ex + byi*hz - bzi*hy |
|---|
| 700 | dbgy(i) = ey + bzi*hx - bxi*hz |
|---|
| 701 | dbgz(i) =eez + bxi*hy - byi*hx |
|---|
| 702 | 1300 continue |
|---|
| 703 | c apply the impulses |
|---|
| 704 | do 1400 i = 2,ngood |
|---|
| 705 | cord(2,i) = cord(2,i) + dbgx(i) |
|---|
| 706 | cord(4,i) = cord(4,i) + dbgy(i) |
|---|
| 707 | cord(6,i) = cord(6,i) + dbgz(i) |
|---|
| 708 | gam(i) = sqrt(1 + cord(2,i)**2 + cord(4,i)**2 + cord(6,i)**2) |
|---|
| 709 | 1400 continue |
|---|
| 710 | return |
|---|
| 711 | end |
|---|
| 712 | c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* |
|---|
![(please configure the [header_logo] section in trac.ini)](/PSPA/chrome/site/pspa.gif)
