| 1 | subroutine scheff1(dist,mt,wtsauv) |
|---|
| 2 | c scheff version Lloyd Young |
|---|
| 3 | c sce(1)=beam current in amperes. |
|---|
| 4 | c sce(2)=1 program number |
|---|
| 5 | c sce(3)=impact parameter,cm for point to point space-charge |
|---|
| 6 | c calculation. if zero use mesh space-charge calculation. |
|---|
| 7 | c sce(4)=radial mesh size in cm. |
|---|
| 8 | c sce(5)=longitudinal mesh size in cm. |
|---|
| 9 | c sce(6)=no. of radial mesh intervals (le 20) |
|---|
| 10 | c sce(7)=no. of longitudinal mesh intervals (le 200) |
|---|
| 11 | c sce(8)=no. of adjacent bunches |
|---|
| 12 | c sce(9)=distance between bunches. (if zero, pl=sce(5)) |
|---|
| 13 | c sce(10)=opt. if opt.lt.0, use only one ring. if)opt=0, |
|---|
| 14 | c use 2x2 array. if opt.gt.0, number of rings is0, |
|---|
| 15 | c determined by gaus depending on aspect ratio. |
|---|
| 16 | c if opt even radial mesh intervals give equal volume rings. |
|---|
| 17 | c if opt odd radial mesh entervals equal dr. |
|---|
| 18 | c if opt.eq.4 or 5 assume metal wall at z=0. and calculate |
|---|
| 19 | c image charge. |
|---|
| 20 | c sce(11)=frm, remeshing factor (default=1.5) |
|---|
| 21 | c sce(12)=rwall, radius of conducting wall. if zero, no wall |
|---|
| 22 | c-------------------------------------------------------------------------- |
|---|
| 23 | save |
|---|
| 24 | c |
|---|
| 25 | include 'param_sz.h' |
|---|
| 26 | include 'constcom.h' |
|---|
| 27 | include 'coordcom.h' |
|---|
| 28 | include 'fldcom.h' |
|---|
| 29 | include 'misccom.h' |
|---|
| 30 | include 'ncordscom.h' |
|---|
| 31 | include 'sccom.h' |
|---|
| 32 | include 'syscom.h' |
|---|
| 33 | include 'ucom.h' |
|---|
| 34 | c |
|---|
| 35 | parameter (maxdim=84000) |
|---|
| 36 | common/cathcur/rcur |
|---|
| 37 | real dxp(imaa),dyp(imaa),dzp(imaa) |
|---|
| 38 | real rsce(nsce) |
|---|
| 39 | logical fexist,fopen,endfile,qflag |
|---|
| 40 | dimension qol(200) |
|---|
| 41 | c 4221=(20+1)*(200+1), 84000=20*(20+1)*200 |
|---|
| 42 | dimension rm(21), zm(201), ers(maxdim), ezs(maxdim), er(4221), |
|---|
| 43 | 1 ez(4221), aa(4000), ismax(201), iemax(201),rssq(20),zzs(201) |
|---|
| 44 | c-------------------------------------------------------------------------- |
|---|
| 45 | c* |
|---|
| 46 | data irout,izout/0,0/fopen,endfile,qflag/.false.,.false.,.true./ |
|---|
| 47 | if (ngood.le.1) return |
|---|
| 48 | if(sce(3).ne.0.)go to 400 |
|---|
| 49 | if (dist) 50,10,50 |
|---|
| 50 | 10 continue |
|---|
| 51 | write(nnout,*) ' Space charge subroutine number ',iprog |
|---|
| 52 | call rtrans(gt,bt) |
|---|
| 53 | if(dist.eq.0.)then |
|---|
| 54 | irout=0 |
|---|
| 55 | izout=0 |
|---|
| 56 | endif |
|---|
| 57 | 15 continue |
|---|
| 58 | gmesh=gt |
|---|
| 59 | c set up field tables |
|---|
| 60 | 16 continue |
|---|
| 61 | beami=sce(1) |
|---|
| 62 | rmax=sce(4) |
|---|
| 63 | dzmax=sce(5)*gmesh |
|---|
| 64 | nr=sce(6) |
|---|
| 65 | nz=sce(7) |
|---|
| 66 | nr1=nr+1 |
|---|
| 67 | nz1=nz+1 |
|---|
| 68 | im1=nr*nz |
|---|
| 69 | im2=nr1*nz1 |
|---|
| 70 | im3=nr1*nz |
|---|
| 71 | im4=nr*nr1 |
|---|
| 72 | nip=sce(8) |
|---|
| 73 | opt=sce(10) |
|---|
| 74 | if(amod(opt,2.).eq.0.)then |
|---|
| 75 | dr=rmax*rmax/sce(6) |
|---|
| 76 | else |
|---|
| 77 | dr=rmax/sce(6) |
|---|
| 78 | endif |
|---|
| 79 | dz=dzmax/sce(7) |
|---|
| 80 | pl=sce(9)*gmesh |
|---|
| 81 | if(pl.le.0.)pl=dzmax |
|---|
| 82 | frm=sce(11) |
|---|
| 83 | if(frm.le.0.)frm=1.5 |
|---|
| 84 | rwall=sce(12) |
|---|
| 85 | if(nr.gt.20 .or. nz.gt.200 .or. (nr*nr1*nz.gt.maxdim .and. |
|---|
| 86 | *(maxdim/nr1**2*dz.lt.10*rwall .or. maxdim/nr1**2*dz.lt.20*rmax))) |
|---|
| 87 | *then |
|---|
| 88 | call appendparm |
|---|
| 89 | stop' Abnormal stop too many mesh points in space charge grid ' |
|---|
| 90 | endif |
|---|
| 91 | c load rm, zm, rssq, zs |
|---|
| 92 | rm(1)=0.0 |
|---|
| 93 | do 20 i=2,nr1 |
|---|
| 94 | if(amod(opt,2.).eq.0.)then |
|---|
| 95 | rm(i)=sqrt(float(i-1)*dr) |
|---|
| 96 | else |
|---|
| 97 | rm(i)=float(i-1)*dr |
|---|
| 98 | endif |
|---|
| 99 | rssq(i-1)=.5*(rm(i-1)**2+rm(i)**2) |
|---|
| 100 | 20 continue |
|---|
| 101 | zs=.5*dz |
|---|
| 102 | do 30 i=1,nz1 |
|---|
| 103 | zm(i)=float(i-1)*dz |
|---|
| 104 | zzs(i)=zm(i)+zs |
|---|
| 105 | 30 continue |
|---|
| 106 | hl=.5*dzmax |
|---|
| 107 | c---if there is a reference particle, set n1=2; otherwise, n1=1 |
|---|
| 108 | n1=1 |
|---|
| 109 | c--- if (w0.gt.0.) n1=2 |
|---|
| 110 | xnp=npoints+1-n1 |
|---|
| 111 | c load ers and ezs |
|---|
| 112 | c mesh dimensions are in cm. ers and ezs are in 1/cm. |
|---|
| 113 | c c1, c2 and c3 are in cm., and c4 is in MeV-cm. |
|---|
| 114 | c q=coulombs/point. e/epsilon =11.294e6 cm MeV/Coul. |
|---|
| 115 | if(qflag)then |
|---|
| 116 | q=beami/(freq*1.e6*xnp) |
|---|
| 117 | do 21 i=2,npoints |
|---|
| 118 | if(weight(i).eq.0.)weight(i)=q |
|---|
| 119 | 21 continue |
|---|
| 120 | qflag=.false. |
|---|
| 121 | endif |
|---|
| 122 | c---printout during initialization |
|---|
| 123 | if((ip.eq.1.and.dist.eq.0.) .or. (ip.eq.2.and.dist.eq.0.))then |
|---|
| 124 | inip=sce(8) |
|---|
| 125 | iopt=sce(10) |
|---|
| 126 | write(nnout,301) |
|---|
| 127 | 301 format (' distribute charge over rings for space charge calc.') |
|---|
| 128 | write(nnout,310) beami |
|---|
| 129 | 310 format (' beam current =',t40,1pg12.3,' amps') |
|---|
| 130 | write(nnout,330) rmax,sce(5),nr,nz,inip,sce(9),iopt,frm,rwall |
|---|
| 131 | 330 format ( |
|---|
| 132 | . ' radial size of mesh =',t40,g12.3,' cm'/ |
|---|
| 133 | . ' z length of mesh =',t40,g12.3,' cm'/ |
|---|
| 134 | . ' no. of radial mesh points =',t40,i8/ |
|---|
| 135 | . ' no. of z mesh points =',t40,i8/ |
|---|
| 136 | . ' no. of "identical pulses" =',t40,i8/ |
|---|
| 137 | . ' pulse length =',t40,g12.3/ |
|---|
| 138 | . ' OPT =',t40,i8/ |
|---|
| 139 | . ' refresh when energy grows by factor',t40,g12.3/ |
|---|
| 140 | . ' wall radius for images =',t40,g12.3) |
|---|
| 141 | endif |
|---|
| 142 | c--- |
|---|
| 143 | c1=11.294e6/erest |
|---|
| 144 | c-- the above line has replaced the following line because the weight |
|---|
| 145 | c-- is the charge of the point |
|---|
| 146 | c c1=11.294e6*q/erest |
|---|
| 147 | inquire(file='scgrid',exist=fexist) |
|---|
| 148 | if(.not.fopen)then |
|---|
| 149 | open(unit=nschef,file='scgrid',access='sequential', |
|---|
| 150 | * status='unknown',form='unformatted') |
|---|
| 151 | fopen=.true. |
|---|
| 152 | endif |
|---|
| 153 | j=1 |
|---|
| 154 | if(fexist)then |
|---|
| 155 | if(endfile)go to 32 |
|---|
| 156 | read(nschef,err=34,end=33)rdz,rgmesh,rdzmax,rc1,rsce,ers,ezs |
|---|
| 157 | if(abs(rdz-dz).gt.0.1*dz)go to 32 |
|---|
| 158 | do 31 i=4,12 |
|---|
| 159 | j=i |
|---|
| 160 | if(sce(i).ne.rsce(i))go to 32 |
|---|
| 161 | 31 continue |
|---|
| 162 | c stored data good |
|---|
| 163 | if(dz.ne.rdz)then |
|---|
| 164 | dz=rdz |
|---|
| 165 | gmesh=rgmesh |
|---|
| 166 | dzmax=rdzmax |
|---|
| 167 | do 530 i=1,nz1 |
|---|
| 168 | zm(i)=float(i-1)*dz |
|---|
| 169 | zzs(i)=zm(i)+zs |
|---|
| 170 | 530 continue |
|---|
| 171 | hl=.5*dzmax |
|---|
| 172 | endif |
|---|
| 173 | if(rc1.ne.c1)then |
|---|
| 174 | do 35 i=1,maxdim |
|---|
| 175 | ers(i)=ers(i)*c1/rc1 |
|---|
| 176 | ezs(i)=ezs(i)*c1/rc1 |
|---|
| 177 | 35 continue |
|---|
| 178 | endif |
|---|
| 179 | write(ndiag,*)' using stored scheff mesh data' |
|---|
| 180 | go to 41 |
|---|
| 181 | 33 continue |
|---|
| 182 | write(ndiag,*)' end of file' |
|---|
| 183 | endfile=.true. |
|---|
| 184 | go to 32 |
|---|
| 185 | 34 continue |
|---|
| 186 | write(ndiag,*)' error on reading file pascgrid' |
|---|
| 187 | endfile=.true. |
|---|
| 188 | c stored data no good |
|---|
| 189 | 32 continue |
|---|
| 190 | backspace(nschef) |
|---|
| 191 | write(ndiag,*)' not using stored data rdz,dz,i,rsce(i),sce(i)', |
|---|
| 192 | * rdz, dz,j,rsce(j),sce(j) |
|---|
| 193 | endif |
|---|
| 194 | l=0 |
|---|
| 195 | zm1=zm(1) |
|---|
| 196 | zm2=zm(2) |
|---|
| 197 | d2zm=(zm2-zm1)*2 |
|---|
| 198 | if(d2zm.lt.rmax)d2zm=rmax |
|---|
| 199 | c do 40 k=1,nr the order of these two do loops are reversed |
|---|
| 200 | c to put the nearest ers and ezs in first so that |
|---|
| 201 | c those faraway can be neglected. |
|---|
| 202 | do 40 j=1,nz |
|---|
| 203 | c***** zp=zm(j+1) |
|---|
| 204 | zp=zzs(j) |
|---|
| 205 | if(zp.gt.d2zm)then |
|---|
| 206 | topt=0 |
|---|
| 207 | elseif(zp.gt.2*d2zm)then |
|---|
| 208 | topt=-1 |
|---|
| 209 | else |
|---|
| 210 | topt=opt |
|---|
| 211 | endif |
|---|
| 212 | do 40 k=1,nr |
|---|
| 213 | if(l.ge.maxdim-nr1)go to 42 |
|---|
| 214 | do 40 i=1,nr1 |
|---|
| 215 | rp=rm(i) |
|---|
| 216 | call gaus(rm(k),rm(k+1),zm1,zm2,topt,er1,ez1) |
|---|
| 217 | l=l+1 |
|---|
| 218 | ers(l)=c1*er1 |
|---|
| 219 | ezs(l)=c1*ez1 |
|---|
| 220 | 40 continue |
|---|
| 221 | 42 continue |
|---|
| 222 | write(nschef)dz,gmesh,dzmax,c1,sce,ers,ezs |
|---|
| 223 | endfile(nschef) |
|---|
| 224 | 41 continue |
|---|
| 225 | if (dist.eq.0.)return |
|---|
| 226 | go to 60 |
|---|
| 227 | c evaluate and apply space charge effects. |
|---|
| 228 | 50 continue |
|---|
| 229 | call rtrans(gt,bt) |
|---|
| 230 | if(gmesh.le.0.)go to 15 |
|---|
| 231 | if(gt/gmesh.gt.frm)then |
|---|
| 232 | gmesh=gmesh*frm |
|---|
| 233 | go to 16 |
|---|
| 234 | endif |
|---|
| 235 | 60 continue |
|---|
| 236 | zc=cord(5,1) |
|---|
| 237 | c c3=dist/gt |
|---|
| 238 | c c4=c3 |
|---|
| 239 | c evaluate ng, xbar, ybar, and zbar. |
|---|
| 240 | c eng=ngood+1-n1 |
|---|
| 241 | eng=0. |
|---|
| 242 | xbar=0. |
|---|
| 243 | ybar=0. |
|---|
| 244 | xsq=0. |
|---|
| 245 | ysq=0. |
|---|
| 246 | zbar=0. |
|---|
| 247 | do 90 np=n1,ngood |
|---|
| 248 | dbgz=cord(6,np) |
|---|
| 249 | dbgx=cord(2,np) |
|---|
| 250 | dbgy=cord(4,np) |
|---|
| 251 | dgam=sqrt(1.d0+dbgx*dbgx+dbgy*dbgy+dbgz*dbgz) |
|---|
| 252 | c gam(np)=gt*(dgam-bt*cord(6,np)) |
|---|
| 253 | bgstar=gt*(dbgz-bt*dgam) |
|---|
| 254 | zstar=gt*(cord(5,np)-zc) |
|---|
| 255 | cord(6,np)=bgstar |
|---|
| 256 | xstar=cord(1,np) |
|---|
| 257 | ystar=cord(3,np) |
|---|
| 258 | xbar=xbar+xstar*weight(np) |
|---|
| 259 | ybar=ybar+ystar*weight(np) |
|---|
| 260 | xsq=xsq+xstar**2*weight(np) |
|---|
| 261 | ysq=ysq+ystar**2*weight(np) |
|---|
| 262 | zbar=zbar+zstar*weight(np) |
|---|
| 263 | eng=eng+weight(np) |
|---|
| 264 | 90 continue |
|---|
| 265 | xbar=xbar/eng |
|---|
| 266 | ybar=ybar/eng |
|---|
| 267 | xsq=xsq/eng |
|---|
| 268 | ysq=ysq/eng |
|---|
| 269 | c epsq=sqrt((xsq-xbar**2)/(ysq-ybar**2)) |
|---|
| 270 | epsq=1. |
|---|
| 271 | c repsq=1./epsq |
|---|
| 272 | repsq=1. |
|---|
| 273 | c xfac=2./(1.+epsq) |
|---|
| 274 | xfac=1. |
|---|
| 275 | c yfac=epsq*xfac |
|---|
| 276 | yfac=1. |
|---|
| 277 | zbar=zbar/eng |
|---|
| 278 | zc=zc+dzmax/nz*(ranf()-.5)/gt |
|---|
| 279 | c clear and load bins |
|---|
| 280 | do 100 i=1,im1 |
|---|
| 281 | aa(i)=0.0 |
|---|
| 282 | 100 continue |
|---|
| 283 | ng=0 |
|---|
| 284 | do 110 np=2,ngood |
|---|
| 285 | charge=1. |
|---|
| 286 | zmstar=0. |
|---|
| 287 | zstar=gt*(cord(5,np)-zc) |
|---|
| 288 | if((opt.eq.4 .or. opt.eq.5) .and. cord(5,np).lt.hl)then |
|---|
| 289 | if(rcur.ne.0. .and. cord(5,np).lt.rcur)then |
|---|
| 290 | a=sqrt((rcur-cord(5,np))**2+cord(1,np)**2+cord(3,np)**2) |
|---|
| 291 | mcharge=-rcur/a |
|---|
| 292 | aprime=rcur**2/a |
|---|
| 293 | xmstar=cord(1,np)*aprime/a |
|---|
| 294 | ymstar=cord(3,np)*aprime/a |
|---|
| 295 | zmstar=gt*(rcur-aprime/a*(rcur-cord(5,np))-zc) |
|---|
| 296 | else |
|---|
| 297 | mcharge=-1. |
|---|
| 298 | zmstar=-gt*(cord(5,np)+zc) |
|---|
| 299 | endif |
|---|
| 300 | endif |
|---|
| 301 | xstar=cord(1,np) |
|---|
| 302 | ystar=cord(3,np) |
|---|
| 303 | rsq=repsq*(xstar-xbar)**2+epsq*(ystar-ybar)**2 |
|---|
| 304 | if(rsq.gt.rmax**2)go to 106 |
|---|
| 305 | if((opt.eq.4 .or. opt.eq.5) .and. cord(5,np).lt.0.)go to 110 |
|---|
| 306 | if(amod(opt,2.).eq.0.)then |
|---|
| 307 | i=rsq/dr+1. |
|---|
| 308 | else |
|---|
| 309 | i=sqrt(rsq)/dr+1. |
|---|
| 310 | endif |
|---|
| 311 | if (i.gt.nr) go to 106 |
|---|
| 312 | z=zstar |
|---|
| 313 | 104 continue |
|---|
| 314 | if (abs(z).le.hl) go to 105 |
|---|
| 315 | if(charge.lt.0.)go to 110 |
|---|
| 316 | if(nip.eq.0)go to 108 |
|---|
| 317 | if (pl.ne.dzmax) go to 108 |
|---|
| 318 | z=z-sign(pl,z) |
|---|
| 319 | c------distribute charge among adjacent bins. |
|---|
| 320 | 105 continue |
|---|
| 321 | ng=ng+1 |
|---|
| 322 | zz=z+hl |
|---|
| 323 | jm1=zz/dz+1. |
|---|
| 324 | i1=i+1 |
|---|
| 325 | if(rsq.lt.rssq(i)) i1=i-1 |
|---|
| 326 | if(i1.lt.1)i1=1 |
|---|
| 327 | if(i1.gt.nr) i1=nr |
|---|
| 328 | j1=jm1+1 |
|---|
| 329 | if(zz.lt.zzs(jm1)) j1=jm1-1 |
|---|
| 330 | if(pl.eq.dzmax)then |
|---|
| 331 | if(j1.lt.1)j1=nz |
|---|
| 332 | if(j1.gt.nz)j1=1 |
|---|
| 333 | else |
|---|
| 334 | if(j1.lt.1)j1=1 |
|---|
| 335 | if(j1.gt.nz)j1=nz |
|---|
| 336 | endif |
|---|
| 337 | a=1. |
|---|
| 338 | if(i1.ne.i)a=(rsq-rssq(i1))/(rssq(i)-rssq(i1)) |
|---|
| 339 | b=1.-a |
|---|
| 340 | cc=1. |
|---|
| 341 | if(j1.ne.jm1) cc=(zz-zzs(j1))/(zzs(jm1)-zzs(j1)) |
|---|
| 342 | d=1.-cc |
|---|
| 343 | k=(jm1-1)*nr+i |
|---|
| 344 | aa(k)=aa(k)+a*cc*weight(np)*charge |
|---|
| 345 | k=k+i1-i |
|---|
| 346 | aa(k)=aa(k)+b*cc*weight(np)*charge |
|---|
| 347 | k=(j1-1)*nr+i |
|---|
| 348 | aa(k)=aa(k)+a*d*weight(np)*charge |
|---|
| 349 | k=k+i1-i |
|---|
| 350 | aa(k)=aa(k)+b*d*weight(np)*charge |
|---|
| 351 | if((opt.eq.4 .or. opt.eq.5.).and.zmstar.ne.0)then |
|---|
| 352 | z=zmstar |
|---|
| 353 | if(rcur.ne.0.and.cord(5,np).lt.rcur)then |
|---|
| 354 | xstar=xmstar |
|---|
| 355 | ystar=ymstar |
|---|
| 356 | endif |
|---|
| 357 | charge=mcharge |
|---|
| 358 | zmstar=0. |
|---|
| 359 | go to 104 |
|---|
| 360 | endif |
|---|
| 361 | go to 110 |
|---|
| 362 | 106 continue |
|---|
| 363 | if(irout.ne.0 .or. cord(5,np).lt.0.)go to 110 |
|---|
| 364 | write(ndiag,*) ' In element number ',cord(7,np) |
|---|
| 365 | write(ndiag,107)np,xstar,ystar |
|---|
| 366 | 107 format(' warning--particle no.',i4,' is outside radial mesh'/ |
|---|
| 367 | * ' xstar=',f8.5,', ystar=',f8.5) |
|---|
| 368 | irout=1 |
|---|
| 369 | go to 110 |
|---|
| 370 | 108 continue |
|---|
| 371 | if(izout.ne.0 .or. cord(5,np).lt.0.)go to 110 |
|---|
| 372 | write(ndiag,*) ' In element number ',cord(7,np) |
|---|
| 373 | write(ndiag,109) np,zstar,cord(5,np) |
|---|
| 374 | 109 format(' warning--particle no.',i4,' is outside longitudinal mesh' |
|---|
| 375 | * / ' zstar=',f10.5,', z=',f10.5) |
|---|
| 376 | izout=1 |
|---|
| 377 | 110 continue |
|---|
| 378 | eng=ng |
|---|
| 379 | c---calculate no. of particles in each column |
|---|
| 380 | do 115 j=1,nz |
|---|
| 381 | qol(j)=0. |
|---|
| 382 | l=(j-1)*nr |
|---|
| 383 | do 114 i=1,nr |
|---|
| 384 | k=l+i |
|---|
| 385 | qol(j)=qol(j)+aa(k) |
|---|
| 386 | 114 continue |
|---|
| 387 | qol(j)=qol(j)/dz |
|---|
| 388 | 115 continue |
|---|
| 389 | c find ismax for each j |
|---|
| 390 | do 130 j=1,nz |
|---|
| 391 | l=(j-1)*nr |
|---|
| 392 | k=nr |
|---|
| 393 | do 120 i=1,nr |
|---|
| 394 | m=l+k |
|---|
| 395 | if (aa(m)) 130,120,130 |
|---|
| 396 | 120 continue |
|---|
| 397 | k=k-1 |
|---|
| 398 | 130 ismax(j)=k |
|---|
| 399 | c find iemax for each j |
|---|
| 400 | iemax(1)=1+ismax(1) |
|---|
| 401 | do 140 j=2,nz |
|---|
| 402 | iemax(j)=1+max0(ismax(j-1),ismax(j)) |
|---|
| 403 | 140 continue |
|---|
| 404 | iemax(nz1)=1+ismax(nz) |
|---|
| 405 | c set er and ez to zero |
|---|
| 406 | do 150 i=1,im2 |
|---|
| 407 | er(i)=0.0 |
|---|
| 408 | ez(i)=0.0 |
|---|
| 409 | 150 continue |
|---|
| 410 | c sum up fields |
|---|
| 411 | do 200 js=1,nz |
|---|
| 412 | c js1=js+1 |
|---|
| 413 | ism=ismax(js) |
|---|
| 414 | if(ism.eq.0) go to 200 |
|---|
| 415 | 160 continue |
|---|
| 416 | do 190 is=1,ism |
|---|
| 417 | l=(js-1)*nr+is |
|---|
| 418 | a1=aa(l) |
|---|
| 419 | if (a1.eq.0.) go to 190 |
|---|
| 420 | c l=(is-1)*im3 need to modify this because ers and ezs are |
|---|
| 421 | c loaded in different order |
|---|
| 422 | l=(is-1)*nr1 |
|---|
| 423 | c***** do 170 je=1,js |
|---|
| 424 | do 170 je=1,js-1 |
|---|
| 425 | c k1=l+(js-je)*nr1 need to modify this because ers and |
|---|
| 426 | c ezs are loaded in different order |
|---|
| 427 | k1=l+(js-je)*im4 |
|---|
| 428 | m1=(je-1)*nr1 |
|---|
| 429 | iem=iemax(je) |
|---|
| 430 | if(iem.le.1 .or. k1.ge.maxdim-nr1)goto 170 |
|---|
| 431 | do 168 ie=1,iem |
|---|
| 432 | c n=m1+ie |
|---|
| 433 | c k=k1+ie |
|---|
| 434 | er(m1+ie)=er(m1+ie)+a1*ers(k1+ie) |
|---|
| 435 | ez(m1+ie)=ez(m1+ie)-a1*ezs(k1+ie) |
|---|
| 436 | 168 continue |
|---|
| 437 | 170 continue |
|---|
| 438 | c***** do 180 je=js1,nz1 |
|---|
| 439 | do 180 je=js,nz |
|---|
| 440 | c k1=l+(je-js1)*nr1 need to modify this because ers and |
|---|
| 441 | c ezs are loaded in different order |
|---|
| 442 | c***** k1=l+(je-js1)*im4 |
|---|
| 443 | k1=l+(je-js)*im4 |
|---|
| 444 | m1=(je-1)*nr1 |
|---|
| 445 | iem=iemax(je) |
|---|
| 446 | if(iem.le.1 .or. k1.ge.maxdim-nr1) go to 180 |
|---|
| 447 | do 178 ie=1,iem |
|---|
| 448 | c n=m1+ie |
|---|
| 449 | c k=k1+ie |
|---|
| 450 | er(m1+ie)=er(m1+ie)+a1*ers(k1+ie) |
|---|
| 451 | ez(m1+ie)=ez(m1+ie)+a1*ezs(k1+ie) |
|---|
| 452 | 178 continue |
|---|
| 453 | 180 continue |
|---|
| 454 | 190 continue |
|---|
| 455 | 200 continue |
|---|
| 456 | c***** |
|---|
| 457 | do 300 i=1,nr1 |
|---|
| 458 | ez(im3+i)=ez(i) |
|---|
| 459 | er(im3+1)=er(i) |
|---|
| 460 | 300 continue |
|---|
| 461 | c***** |
|---|
| 462 | c---evaluate and apply space charge impulse. |
|---|
| 463 | 201 continue |
|---|
| 464 | do 240 np=n1,ngood |
|---|
| 465 | bgfstar=cord(6,np) |
|---|
| 466 | cbug if(cord(5,np).le.-dist*bt) go to 235 |
|---|
| 467 | cbug if(cord(5,np).lt.0.and.(opt.eq.4. .or. opt.eq.5))go to 235 |
|---|
| 468 | cbug if(cord(5,np).le.dist*bt)then |
|---|
| 469 | cbug if(cord(5,np).gt.0)then |
|---|
| 470 | cbug c3=.5*cord(5,np)/(bt*gt)+.5*dist/gt |
|---|
| 471 | cbug else |
|---|
| 472 | cbug c3=.5*(dist*bt+cord(5,np))**2/(bt*gt*(dist*bt)) |
|---|
| 473 | cbug endif |
|---|
| 474 | cbug else |
|---|
| 475 | cbug c3=dist/gt |
|---|
| 476 | cbug endif |
|---|
| 477 | c goto 235 no space charge kicks until particle is emitted from the cathode |
|---|
| 478 | if(cord(5,np).le.0) go to 235 |
|---|
| 479 | c turn-off screening at dist*bt, skip screening for metallic cathodes |
|---|
| 480 | if(cord(5,np).le.dist*bt.and.(opt.ne.4. .or. opt.ne.5))then |
|---|
| 481 | if(cord(5,np).gt.0)then |
|---|
| 482 | c3=.5*cord(5,np)/(bt*gt)+.5*dist/gt |
|---|
| 483 | else |
|---|
| 484 | c3=.5*(dist*bt+cord(5,np))**2/(bt*gt*(dist*bt)) |
|---|
| 485 | endif |
|---|
| 486 | else |
|---|
| 487 | c3=dist/gt |
|---|
| 488 | endif |
|---|
| 489 | zstar=gt*(cord(5,np)-zc) |
|---|
| 490 | xstar=cord(1,np) |
|---|
| 491 | ystar=cord(3,np) |
|---|
| 492 | r=sqrt(repsq*(xstar-xbar)**2+epsq*(ystar-ybar)**2) |
|---|
| 493 | if(np.eq.1) go to 202 |
|---|
| 494 | if(r.lt.0.000001) r=.000001 |
|---|
| 495 | xor=(xstar-xbar)*xfac/r |
|---|
| 496 | yor=(ystar-ybar)*yfac/r |
|---|
| 497 | go to 203 |
|---|
| 498 | 202 continue |
|---|
| 499 | xor=0. |
|---|
| 500 | yor=0. |
|---|
| 501 | 203 continue |
|---|
| 502 | z=zstar |
|---|
| 503 | 205 continue |
|---|
| 504 | if(abs(z).le.hl) go to 210 |
|---|
| 505 | if (nip.le.0) go to 235 |
|---|
| 506 | if (pl.ne.dzmax) go to 235 |
|---|
| 507 | z=z-sign(pl,z) |
|---|
| 508 | go to 205 |
|---|
| 509 | c interpolate impulse within mesh. |
|---|
| 510 | c 210 if(r.gt.rmax) go to 220 |
|---|
| 511 | 210 continue |
|---|
| 512 | if(amod(opt,2.).eq.0..and.r.le.rmax)then |
|---|
| 513 | rb=(r*r/dr) |
|---|
| 514 | i=1.+rb |
|---|
| 515 | a=(r-rm(i))/(rm(i+1)-rm(i)) |
|---|
| 516 | elseif(r.le.rmax)then |
|---|
| 517 | rb=r/dr |
|---|
| 518 | i=1.0+rb |
|---|
| 519 | a=rb-float(i-1) |
|---|
| 520 | endif |
|---|
| 521 | if(r.gt.rmax.and.rwall.ne.0)then |
|---|
| 522 | a=(r-rmax)/(rwall-rmax) |
|---|
| 523 | i=nr1 |
|---|
| 524 | c---if rwall is zero arbitraily asume an effective rwall of 2*rmax to |
|---|
| 525 | c calculate a and i |
|---|
| 526 | elseif(r.gt.rmax)then |
|---|
| 527 | a=(r-rmax)/rmax |
|---|
| 528 | i=nr1 |
|---|
| 529 | endif |
|---|
| 530 | b=1.0-a |
|---|
| 531 | c****** zb=(z+hl)/dz |
|---|
| 532 | zb=(z+hl-.5*dz)/dz |
|---|
| 533 | if(zb.lt.0.)zb=zb+float(nz) |
|---|
| 534 | c***** |
|---|
| 535 | j=1.0+zb |
|---|
| 536 | c=zb-float(j-1) |
|---|
| 537 | d=1.0-c |
|---|
| 538 | l=i+(j-1)*nr1 |
|---|
| 539 | m=l+nr1 |
|---|
| 540 | if(r.gt.rmax.and.r.lt.rwall)then |
|---|
| 541 | c the following two lines assume all the charge is within rmax. |
|---|
| 542 | cbgr=c3*(d*er(l)+c*er(m))*rmax/r |
|---|
| 543 | cbgz=c3*(d*b*ez(l)+c*b*ez(m)) |
|---|
| 544 | elseif(r.lt.rwall .or. rwall.eq.0.)then |
|---|
| 545 | cbgr=c3*(d*(a*er(l+1)+b*er(l))+c*(a*er(m+1)+b*er(m))) |
|---|
| 546 | cbgz=c3*(d*(a*ez(l+1)+b*ez(l))+c*(a*ez(m+1)+b*ez(m))) |
|---|
| 547 | endif |
|---|
| 548 | c---The following makes a correction to the kinetic energy for the |
|---|
| 549 | c--- space charge potential energy. |
|---|
| 550 | if(cord(5,np).lt.dist*bt.and.cord(5,np).gt.0..and.rwall.ne.0. |
|---|
| 551 | * .and.cord(6,np).gt.0..and.rm(iemax(j)).gt.0.)then |
|---|
| 552 | iem=min0(iemax(j),iemax(j+1)) |
|---|
| 553 | cbgz=cbgz-(log(rwall)-log(rm(iem)))*(d*er(iem+(j-1)*nr1)+ |
|---|
| 554 | * c*er(iem+j*nr1)) |
|---|
| 555 | cbgz=cbgz-.5*(d*(a*er(l+1)+b*er(l))+c*(a*er(m+1)+b*er(m)) |
|---|
| 556 | * +d*er(l+1)+c*er(m+1))*(rm(i+1)-r) |
|---|
| 557 | do 215 k=i+2,iem |
|---|
| 558 | cbgz=cbgz-.5*(d*(er(k+(j-1)*nr1)+er(k-1+(j-1)*nr1)) |
|---|
| 559 | * +c*(er(k+j*nr1)+er(k-1+j*nr1)))*(rm(k)-rm(k-1)) |
|---|
| 560 | 215 continue |
|---|
| 561 | endif |
|---|
| 562 | c estimate impulse based on point charge at xbar,ybar,zbar. |
|---|
| 563 | c220 z=(zstar-zbar)*(1.+bt*cord(6,np)/gam(np)) |
|---|
| 564 | c---the following lines are for a point charge at beam centroid. |
|---|
| 565 | c220 z=(zstar-zbar) |
|---|
| 566 | c d=sqrt(z**2+r**2) |
|---|
| 567 | c rod3=r/d**3 |
|---|
| 568 | c zod3=z/d**3 |
|---|
| 569 | c if(nip.eq.0) go to 228 |
|---|
| 570 | c include neighboring bunches. |
|---|
| 571 | c do 226 i=1,nip |
|---|
| 572 | c xi=i |
|---|
| 573 | c do 226 j=1,2 |
|---|
| 574 | c s=z+xi*pl |
|---|
| 575 | c d=sqrt(s**2+r**2) |
|---|
| 576 | c rod3=rod3+r/d**3 |
|---|
| 577 | c zod3=zod3+s/d**3 |
|---|
| 578 | c 226 xi=-xi |
|---|
| 579 | c evaluate impulse. |
|---|
| 580 | c 228 cbgr=eng*c1*c3*rod3/(4.*pi) |
|---|
| 581 | c cbgz=eng*c1*c4*zod3/(4.*pi) |
|---|
| 582 | c---following two lines are for line charge with density qol(j) |
|---|
| 583 | c 220 cbgz=0. |
|---|
| 584 | c cbgr=c1*c3*qol(j)/(2.*pi*r**2) |
|---|
| 585 | c apply impulse |
|---|
| 586 | 230 continue |
|---|
| 587 | dbgx=cord(2,np)+cbgr*xor |
|---|
| 588 | cord(2,np)=dbgx |
|---|
| 589 | dbgy=cord(4,np)+cbgr*yor |
|---|
| 590 | cord(4,np)=dbgy |
|---|
| 591 | bgfstar=cord(6,np)+cbgz |
|---|
| 592 | 231 continue |
|---|
| 593 | gfstar=sqrt(1.d0+dbgx*dbgx+dbgy*dbgy+bgfstar*bgfstar) |
|---|
| 594 | cord(6,np)=gt*(bgfstar+bt*gfstar) |
|---|
| 595 | gam(np)=gt*(gfstar+bt*bgfstar) |
|---|
| 596 | go to 240 |
|---|
| 597 | 235 continue |
|---|
| 598 | dbgx=cord(2,np) |
|---|
| 599 | dbgy=cord(4,np) |
|---|
| 600 | bgfstar=cord(6,np) |
|---|
| 601 | go to 231 |
|---|
| 602 | ctest du bug |
|---|
| 603 | cbug write(89,*) cbgr*xor*gam(np)/dist, |
|---|
| 604 | cbug * cbgr*yor*gam(np)/dist, |
|---|
| 605 | cbug * cbgz/dist, |
|---|
| 606 | cbug * cord(1,np),cord(3,np),cord(5,np) |
|---|
| 607 | c |
|---|
| 608 | 240 continue |
|---|
| 609 | return |
|---|
| 610 | c |
|---|
| 611 | c very simple inflexible 3d space charge calculator |
|---|
| 612 | c uses point-to-point interaction with finite radius |
|---|
| 613 | c as suggested by G. Boicourt |
|---|
| 614 | c relevant scheff parameters are: |
|---|
| 615 | c sce(1) = beam current, amps |
|---|
| 616 | c sce(3) = impact parameter, cm (and flag for point-to-point) |
|---|
| 617 | c |
|---|
| 618 | 400 continue |
|---|
| 619 | if (dist.eq.0.) then |
|---|
| 620 | write(nnout,*) ' Space charge subroutine number ',iprog |
|---|
| 621 | write(nnout,*) ' Point to point interaction with finite radius' |
|---|
| 622 | write(nnout,*) ' Impact parameter =',sce(3),' cm' |
|---|
| 623 | return |
|---|
| 624 | endif |
|---|
| 625 | call rtrans(gt,bt) |
|---|
| 626 | a3 = sce(3)**2 |
|---|
| 627 | points = npoints |
|---|
| 628 | q=sce(1)/(freq*1.e6*points) |
|---|
| 629 | con=8.984e5*q*dist/(gt*erest) |
|---|
| 630 | do 410 i=1,ngood |
|---|
| 631 | dxp(i) = 0. |
|---|
| 632 | dyp(i) = 0. |
|---|
| 633 | dzp(i) = 0. |
|---|
| 634 | 410 continue |
|---|
| 635 | do 450 np=1,ngood |
|---|
| 636 | c |
|---|
| 637 | c only do calculation once for each pair of points |
|---|
| 638 | c |
|---|
| 639 | cdir$ IVDEP |
|---|
| 640 | do 430 mp=np+1,ngood |
|---|
| 641 | x = cord(1,np) - cord(1,mp) |
|---|
| 642 | y = cord(3,np) - cord(3,mp) |
|---|
| 643 | z =(cord(5,np) - cord(5,mp))*gt |
|---|
| 644 | d = x*x + y*y + z*z + a3 |
|---|
| 645 | d3i= 1./(d*sqrt(d)) |
|---|
| 646 | dxp(np) = dxp(np) + x * d3i |
|---|
| 647 | dyp(np) = dyp(np) + y * d3i |
|---|
| 648 | dzp(np) = dzp(np) + z * d3i |
|---|
| 649 | dxp(mp) = dxp(mp) - x * d3i |
|---|
| 650 | dyp(mp) = dyp(mp) - y * d3i |
|---|
| 651 | dzp(mp) = dzp(mp) - z * d3i |
|---|
| 652 | 430 continue |
|---|
| 653 | c |
|---|
| 654 | c increment xprime, yprime and energy |
|---|
| 655 | c |
|---|
| 656 | rr=amax1(0.,cord(5,np))/(cord(5,np)+1.e-20) |
|---|
| 657 | cord(2,np) = cord(2,np) + con * dxp(np) * rr |
|---|
| 658 | cord(4,np) = cord(4,np) + con * dyp(np) * rr |
|---|
| 659 | bgfstar=gt*(cord(6,np)-bt*gam(np))+con*dzp(np)*rr |
|---|
| 660 | gfstar=sqrt(1.+cord(2,np)**2+cord(4,np)**2+bgfstar**2) |
|---|
| 661 | cord(6,np)=gt*(bgfstar+bt*gfstar) |
|---|
| 662 | gam(np)=gt*(gfstar+bt*bgfstar) |
|---|
| 663 | 450 continue |
|---|
| 664 | return |
|---|
| 665 | end |
|---|
| 666 | c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* |
|---|
![(please configure the [header_logo] section in trac.ini)](/PSPA/chrome/site/pspa.gif)
