| 1 | subroutine out6dp(iz0,iplace) |
|---|
| 2 | c print emittances |
|---|
| 3 | c iz0 = 0 -- include all good particles |
|---|
| 4 | c iz0 = 1 -- include only good particles with cord(7,i) > 0 |
|---|
| 5 | c |
|---|
| 6 | c optcon(1) = 6 -- to invoke OUT6 via the OUTPUT card |
|---|
| 7 | c optcon(2) ne 0 -- use canonical momentum rather than |
|---|
| 8 | c ordinary momentum in emittance calculation |
|---|
| 9 | c optcon(3) ne 0 -- printer plot of x vs gamma*betax |
|---|
| 10 | c < 0 -- also print projection plots |
|---|
| 11 | c optcon(4) ne 0 -- plot y vs gamma*betay |
|---|
| 12 | c optcon(5) ne 0 -- plot z vs r |
|---|
| 13 | c optcon(6) ne 0 -- plot z vs gamma*beta |
|---|
| 14 | c optcon(7) ne 0 -- plot z vs del(gamma*beta)/<gamma*beta> |
|---|
| 15 | c |
|---|
| 16 | c OUT6 is called with IZ0 = 0 from INPUT if INTYPE = 10, If want |
|---|
| 17 | c plots at this time, the OUTPUT card must precede the INPUT card |
|---|
| 18 | c |
|---|
| 19 | c |
|---|
| 20 | c----------------------------------------------------------------------- |
|---|
| 21 | implicit double precision (b,g,r,t,x,y,z) |
|---|
| 22 | c warning variables en simple precision dans les commons |
|---|
| 23 | real zloc ! systeme |
|---|
| 24 | real gam ! coords |
|---|
| 25 | real beami,brhof,radian,twopi,z0,zlimit ! const |
|---|
| 26 | save |
|---|
| 27 | c |
|---|
| 28 | include 'param_sz.h' |
|---|
| 29 | include 'cfldscom.h' |
|---|
| 30 | include 'constcom.h' |
|---|
| 31 | include 'coordcom.h' |
|---|
| 32 | include 'flagcom.h' |
|---|
| 33 | include 'ncordscom.h' |
|---|
| 34 | include 'outcom.h' |
|---|
| 35 | include 'syscom.h' |
|---|
| 36 | include 'tstepcom.h' |
|---|
| 37 | include 'ucom.h' |
|---|
| 38 | include 'upawcom.h' |
|---|
| 39 | include 'var_char.h' |
|---|
| 40 | c |
|---|
| 41 | double precision adp |
|---|
| 42 | common/ddat/ctime,ddat |
|---|
| 43 | cbm common/hist2/xmax,bxmax,ymax,bymax,zminc,zmaxc,rmax,gbmin, |
|---|
| 44 | cbm * gbmax,gbav,iiz0,ne |
|---|
| 45 | cbm common/hist3/gbxx(imaa),gbyy(imaa) |
|---|
| 46 | dimension gbxx(imaa),gbyy(imaa) |
|---|
| 47 | common/out6/nbphase,nbremesh |
|---|
| 48 | common/outpar/noutp,outpar(0:lmx,10),moutp(10) |
|---|
| 49 | c-------------------------------------------------------------------------- |
|---|
| 50 | c* |
|---|
| 51 | c need double precision for emittance calc |
|---|
| 52 | c |
|---|
| 53 | data phold/0./,neold/0/,cp/1./,sp/0./ |
|---|
| 54 | c hfac = Compton wavelength times critical magnetic field |
|---|
| 55 | data hfac/1703./ |
|---|
| 56 | c |
|---|
| 57 | iiz0=iz0 |
|---|
| 58 | c copy parameters in case OUT6 not called from OUTPUT |
|---|
| 59 | if(iplace.eq.0) then |
|---|
| 60 | navv=nnout |
|---|
| 61 | else |
|---|
| 62 | navv=nav |
|---|
| 63 | endif |
|---|
| 64 | do 5 i = 1,lmx |
|---|
| 65 | optcon(i) = outpar(i,7) |
|---|
| 66 | 5 continue |
|---|
| 67 | xav = 0. |
|---|
| 68 | xsqav = 0. |
|---|
| 69 | bxav = 0. |
|---|
| 70 | bxsqav = 0. |
|---|
| 71 | xbxav = 0. |
|---|
| 72 | xbxsav = 0. |
|---|
| 73 | yav = 0. |
|---|
| 74 | ysqav = 0. |
|---|
| 75 | byav = 0. |
|---|
| 76 | bysqav = 0. |
|---|
| 77 | ybyav = 0. |
|---|
| 78 | ybysav = 0. |
|---|
| 79 | zav = 0. |
|---|
| 80 | zsqav = 0. |
|---|
| 81 | bzav = 0. |
|---|
| 82 | bzsqav = 0. |
|---|
| 83 | zbzav = 0. |
|---|
| 84 | zbzsav = 0. |
|---|
| 85 | gav = 0. |
|---|
| 86 | gsqav = 0. |
|---|
| 87 | bav = 0. |
|---|
| 88 | bsqav = 0. |
|---|
| 89 | gbav = 0. |
|---|
| 90 | gbsqav = 0. |
|---|
| 91 | zmaxc = 0. |
|---|
| 92 | zminc = 99999. |
|---|
| 93 | adp = 0. |
|---|
| 94 | xmax = 0. |
|---|
| 95 | bxmax = 0. |
|---|
| 96 | ymax = 0. |
|---|
| 97 | bymax = 0. |
|---|
| 98 | rmax = 0. |
|---|
| 99 | gmax = 0. |
|---|
| 100 | gmin = 99999. |
|---|
| 101 | gbmax = 0. |
|---|
| 102 | gbmin = 99999. |
|---|
| 103 | if (ngood.eq.0) go to 110 |
|---|
| 104 | c###################################################################### |
|---|
| 105 | do 100 i = 1,ngood |
|---|
| 106 | ne = cord(7,i) |
|---|
| 107 | c skip if particle is still upstream of first element |
|---|
| 108 | if (ne.eq.0.and.iz0.eq.1) go to 100 |
|---|
| 109 | adp = adp + 1. |
|---|
| 110 | x = cord(1,i) |
|---|
| 111 | if (abs(x).gt.xmax) xmax = abs(x) |
|---|
| 112 | xsq = x**2 |
|---|
| 113 | y = cord(3,i) |
|---|
| 114 | if (abs(y).gt.ymax) ymax = abs(y) |
|---|
| 115 | ysq = y**2 |
|---|
| 116 | r = dsqrt(xsq + ysq) |
|---|
| 117 | if (r.gt.rmax) rmax = r |
|---|
| 118 | z = cord(5,i) |
|---|
| 119 | if (z.gt.zmaxc) zmaxc = z |
|---|
| 120 | if (z.lt.zminc) zminc = z |
|---|
| 121 | ax = 0. |
|---|
| 122 | ay = 0. |
|---|
| 123 | if (optcon(2).eq.0.) go to 20 ! ================================ |
|---|
| 124 | c calculate (transverse) vector potential for canonical momentum |
|---|
| 125 | c actually calc eA/mc**2 |
|---|
| 126 | c skip if not yet in first element |
|---|
| 127 | if (ne.eq.0) go to 20 ! ========================================= |
|---|
| 128 | nt = ntype(ne) |
|---|
| 129 | c skip if not in an rf cell (tho other elements can have nonzero |
|---|
| 130 | c vector potential...) |
|---|
| 131 | if (nt.ne.7) go to 20 |
|---|
| 132 | if (ne.ne.neold) then |
|---|
| 133 | neold = ne |
|---|
| 134 | e0 = el(5,ne) |
|---|
| 135 | s0 = el(9,ne) |
|---|
| 136 | c0 = el(10,ne) |
|---|
| 137 | efac = e0*wavel*(freq/el(4,ne))/twopi/erest |
|---|
| 138 | nc=el(6,ne) |
|---|
| 139 | c spzoff = spzloc(ne) - hcll(nc) |
|---|
| 140 | spzoff = zloc(ne) - hcll(nc) |
|---|
| 141 | if (el(8,ne).gt.0.) spzoff = spzoff + hcll(nc) |
|---|
| 142 | hz = el(11,ne)/hfac/2. |
|---|
| 143 | endif |
|---|
| 144 | ph=wt*radian*el(4,ne)/freq |
|---|
| 145 | if(ph.eq.phold) go to 10 |
|---|
| 146 | sp=sin(ph) |
|---|
| 147 | cp=cos(ph) |
|---|
| 148 | c s = sp*c0+cp*s0 |
|---|
| 149 | c = cp*c0-sp*s0 |
|---|
| 150 | phold=ph |
|---|
| 151 | write(navv,12) e0,s0,c0,wt,c,efac,hz |
|---|
| 152 | 12 format (' e0,s0,c0,wt,c,efac,hz =',3f8.4,f8.2,3f8.4) |
|---|
| 153 | 10 continue |
|---|
| 154 | if (r.gt.0.) then |
|---|
| 155 | spzc = z - spzoff |
|---|
| 156 | call cfield(r,spzc,nc,ez,er,spbt) |
|---|
| 157 | c rf field |
|---|
| 158 | avec = efac*er*c |
|---|
| 159 | ax = avec*x/r |
|---|
| 160 | ay = avec*y/r |
|---|
| 161 | c static mag field |
|---|
| 162 | ax = ax - y*hz |
|---|
| 163 | ay = ay + x*hz |
|---|
| 164 | endif |
|---|
| 165 | 20 continue ! ===================================================== |
|---|
| 166 | c write(navv,11) r,spzc,ez,er,spbt,x,y,ax,cord(2,i),ay,cord(4,i) |
|---|
| 167 | 11 format (' r,zc,ez,er,bt =',5f8.4/ |
|---|
| 168 | . ' x,y,ax,bx,ay,by =',6f8.4) |
|---|
| 169 | c |
|---|
| 170 | c bx = (canonical momentum)/mc, if 2nd parameter on OUTPUT card ne 0 |
|---|
| 171 | bx = cord(2,i)+ ax |
|---|
| 172 | gbxx(i) = bx |
|---|
| 173 | if (abs(bx).gt.bxmax) bxmax = abs(bx) |
|---|
| 174 | bxsq = bx**2 |
|---|
| 175 | xbx = x*bx |
|---|
| 176 | xbxsq = xbx**2 |
|---|
| 177 | xav = xav + x |
|---|
| 178 | xsqav = xsqav + xsq |
|---|
| 179 | bxav = bxav + bx |
|---|
| 180 | bxsqav = bxsqav + bxsq |
|---|
| 181 | xbxav = xbxav + xbx |
|---|
| 182 | xbxsav = xbxsav + xbxsq |
|---|
| 183 | c |
|---|
| 184 | by = cord(4,i) + ay |
|---|
| 185 | gbyy(i) = by |
|---|
| 186 | if (abs(by).gt.bymax) bymax = abs(by) |
|---|
| 187 | bysq = by**2 |
|---|
| 188 | yby = y*by |
|---|
| 189 | ybysq = yby**2 |
|---|
| 190 | yav = yav + y |
|---|
| 191 | ysqav = ysqav + ysq |
|---|
| 192 | byav = byav + by |
|---|
| 193 | bysqav = bysqav + bysq |
|---|
| 194 | ybyav = ybyav + yby |
|---|
| 195 | ybysav = ybysav + ybysq |
|---|
| 196 | c |
|---|
| 197 | zsq = z**2 |
|---|
| 198 | bz = cord(6,i) |
|---|
| 199 | bzsq = bz**2 |
|---|
| 200 | zbz = z*bz |
|---|
| 201 | zbzsq = zbz**2 |
|---|
| 202 | zav = zav + z |
|---|
| 203 | zsqav = zsqav + zsq |
|---|
| 204 | bzav = bzav + bz |
|---|
| 205 | bzsqav = bzsqav + bzsq |
|---|
| 206 | zbzav = zbzav + zbz |
|---|
| 207 | zbzsav = zbzsav + zbzsq |
|---|
| 208 | c |
|---|
| 209 | c roundoff trouble extracting beta from gam(i) when gamma near 1 |
|---|
| 210 | gb = sqrt(cord(2,i)**2 + cord(4,i)**2 + cord(6,i)**2) |
|---|
| 211 | gbsq = gb**2 |
|---|
| 212 | gsq = 1. + gbsq |
|---|
| 213 | g = dsqrt(gsq) |
|---|
| 214 | if (g.gt.gmax) gmax = g |
|---|
| 215 | if (g.lt.gmin) gmin = g |
|---|
| 216 | b = gb/g |
|---|
| 217 | c 20/03/2001 |
|---|
| 218 | if(b.gt.1.) then |
|---|
| 219 | b = 1. |
|---|
| 220 | endif |
|---|
| 221 | bsq = b**2 |
|---|
| 222 | if (gb.gt.gbmax) gbmax = gb |
|---|
| 223 | if (gb.lt.gbmin) gbmin = gb |
|---|
| 224 | gav = gav + g |
|---|
| 225 | gsqav = gsqav + gsq |
|---|
| 226 | bav = bav + b |
|---|
| 227 | bsqav = bsqav + bsq |
|---|
| 228 | gbav = gbav + gb |
|---|
| 229 | gbsqav = gbsqav + gbsq |
|---|
| 230 | 100 continue |
|---|
| 231 | c##################################################################### |
|---|
| 232 | 110 continue |
|---|
| 233 | c 1000. rad to mrad |
|---|
| 234 | c 1.e6=(1.e3)**2 |
|---|
| 235 | c --------------------------------- x -------------------------------- |
|---|
| 236 | if (adp.gt.0.) then |
|---|
| 237 | xav = xav/adp |
|---|
| 238 | xsqav = xsqav/adp |
|---|
| 239 | bxav = bxav/adp*1000. |
|---|
| 240 | bxsqav = bxsqav/adp*1.e6 |
|---|
| 241 | xbxav = xbxav/adp*1000. |
|---|
| 242 | xbxsav = xbxsav/adp*1.e6 |
|---|
| 243 | endif |
|---|
| 244 | xrms = 0. |
|---|
| 245 | bxrms = 0. |
|---|
| 246 | xbrms = 0. |
|---|
| 247 | xemit = 0. |
|---|
| 248 | xsqav = xsqav - xav**2 |
|---|
| 249 | if (xsqav.gt.0.) xrms = dsqrt(xsqav) |
|---|
| 250 | bxsqav = bxsqav - bxav**2 |
|---|
| 251 | if (bxsqav.gt.0.) bxrms = dsqrt(bxsqav) |
|---|
| 252 | temp = xbxsav - xbxav**2 |
|---|
| 253 | c if(temp.gt.0) xbrms = sqrt(temp) |
|---|
| 254 | xbxav = xbxav - xav*bxav |
|---|
| 255 | temp = xsqav*bxsqav - xbxav**2 |
|---|
| 256 | if (temp.gt.0.) xemit = dsqrt(temp) |
|---|
| 257 | c --------------------------------- y -------------------------------- |
|---|
| 258 | if (adp.gt.0.) then |
|---|
| 259 | yav = yav/adp |
|---|
| 260 | ysqav = ysqav/adp |
|---|
| 261 | byav = byav/adp*1000. |
|---|
| 262 | bysqav = bysqav/adp*1.e6 |
|---|
| 263 | ybyav = ybyav/adp*1000. |
|---|
| 264 | ybysav = ybysav/adp*1.e6 |
|---|
| 265 | endif |
|---|
| 266 | yrms = 0. |
|---|
| 267 | byrms = 0. |
|---|
| 268 | ybrms = 0. |
|---|
| 269 | yemit = 0. |
|---|
| 270 | ysqav = ysqav - yav**2 |
|---|
| 271 | if (ysqav.gt.0.) yrms = dsqrt(ysqav) |
|---|
| 272 | bysqav = bysqav - byav**2 |
|---|
| 273 | if (bysqav.gt.0.) byrms = dsqrt(bysqav) |
|---|
| 274 | temp = ybysav - ybyav**2 |
|---|
| 275 | c if(temp.gt.0) ybrms = sqrt(temp) |
|---|
| 276 | ybyav = ybyav - yav*byav |
|---|
| 277 | temp = ysqav*bysqav - ybyav**2 |
|---|
| 278 | if (temp.gt.0.) yemit = dsqrt(temp) |
|---|
| 279 | c rrms = 0.5*(xrms + yrms) |
|---|
| 280 | c --------------------------------- z -------------------------------- |
|---|
| 281 | if (adp.gt.0.) then |
|---|
| 282 | zav = zav/adp |
|---|
| 283 | zsqav = zsqav/adp |
|---|
| 284 | bzav = bzav/adp |
|---|
| 285 | bzsqav = bzsqav/adp |
|---|
| 286 | zbzav = zbzav/adp |
|---|
| 287 | zbzsav = zbzsav/adp |
|---|
| 288 | endif |
|---|
| 289 | zrms = 0. |
|---|
| 290 | zemit = 0. |
|---|
| 291 | zsqav = zsqav - zav**2 |
|---|
| 292 | if (zsqav.gt.0.) zrms = dsqrt(zsqav) |
|---|
| 293 | bzsqav = bzsqav - bzav**2 |
|---|
| 294 | temp = zbzsav - zbzav**2 |
|---|
| 295 | zbzav = zbzav - zav*bzav |
|---|
| 296 | temp = zsqav*bzsqav - zbzav**2 |
|---|
| 297 | c 10. cm to mm |
|---|
| 298 | if (temp.gt.0.) zemit = 10.*dsqrt(temp) |
|---|
| 299 | c --------------------------------- bg -------------------------------- |
|---|
| 300 | if (adp.gt.0.) then |
|---|
| 301 | gav = gav/adp |
|---|
| 302 | gsqav = gsqav/adp |
|---|
| 303 | bav = bav/adp |
|---|
| 304 | bsqav = bsqav/adp |
|---|
| 305 | gbav = gbav/adp |
|---|
| 306 | gbsqav = gbsqav/adp |
|---|
| 307 | endif |
|---|
| 308 | c write(navv,199) gav,gsqav,bav,bsqav,gbav,gbsqav |
|---|
| 309 | 199 format (' g & b aves ',1pg11.3,5g11.3) |
|---|
| 310 | grms = 0. |
|---|
| 311 | brms = 0. |
|---|
| 312 | gbrms = 0. |
|---|
| 313 | temp = gsqav - gav**2 |
|---|
| 314 | if (temp.gt.0.) grms = dsqrt(temp) |
|---|
| 315 | temp = bsqav - bav**2 |
|---|
| 316 | if (temp.gt.0.) brms = dsqrt(temp) |
|---|
| 317 | temp = gbsqav - gbav**2 |
|---|
| 318 | if (temp.gt.0.) gbrms = dsqrt(temp) |
|---|
| 319 | c --------------------------------- write -------------------------------- |
|---|
| 320 | iadp=adp |
|---|
| 321 | write(navv,*)'---------' |
|---|
| 322 | write(navv,190) wt,adp,cord(5,1) |
|---|
| 323 | 190 format (' phase =',f11.2,' degrees, ',f6.0,' active particles', |
|---|
| 324 | 1' zrp = ',f10.4,' cm') |
|---|
| 325 | c convert emittances to mm-mrad (cm-mrad to mm-mrad) |
|---|
| 326 | xemit = xemit*10. |
|---|
| 327 | yemit = yemit*10. |
|---|
| 328 | c zemit=zemit*(360.*1000*erest)/wavel ! [keV-deg] |
|---|
| 329 | c for now, average x and y quantities |
|---|
| 330 | c xrms = 0.5*(xrms + yrms) |
|---|
| 331 | c bxrms = 0.5*(bxrms + byrms) |
|---|
| 332 | c xbxav = 0.5*(xbxav + ybyav) |
|---|
| 333 | c xemit = 0.5*(xemit + yemit) |
|---|
| 334 | write(navv,200) xrms,bxrms,xbxav,xemit |
|---|
| 335 | c changer les f8.4 en f8.2 pour les emit f7.3 en f7.2 pour les gb(x/y)rms |
|---|
| 336 | 200 format (' xrms = ',f7.4,' gbxrms = ',f7.2,' xgbxav = ',f8.2, |
|---|
| 337 | . ' xemit = ',f8.2,' mm-mrad') |
|---|
| 338 | xyemit = 0.5*(xemit + yemit) |
|---|
| 339 | write(navv,210) yrms,byrms,ybyav,yemit,xyemit |
|---|
| 340 | 210 format (' yrms = ',f7.4,' gbyrms = ',f7.2,' ygbyav = ',f8.2, |
|---|
| 341 | . ' yemit = ',f8.2,' <',f8.2,'>') |
|---|
| 342 | write(navv,220) zminc,zmaxc,zemit |
|---|
| 343 | 220 format (' zmin = ',f8.2,' zmax = ',f8.2,' zemit = ',1pg11.3) |
|---|
| 344 | write(navv,230) zav,zrms,gav,grms |
|---|
| 345 | 230 format (' zav = ',g12.5,' +- ',g12.5,' gav = ',g12.5,' +- ',g12.5) |
|---|
| 346 | write(navv,231) bav,brms,gbav,gbrms |
|---|
| 347 | 231 format (' bav = ',g12.5,' +- ',g12.5,' gbav= ',g12.5,' +- ',g12.5) |
|---|
| 348 | write(nsemit)wt,adp,xrms,bxrms,xbxav,xemit,yrms,byrms,ybyav,yemit, |
|---|
| 349 | *zav,zrms,zminc,zmaxc,zemit |
|---|
| 350 | c |
|---|
| 351 | c histogram setup |
|---|
| 352 | c |
|---|
| 353 | do 9996 i=1,ngood |
|---|
| 354 | spxx=cord(1,i) |
|---|
| 355 | if ((ne.eq.0).and.(i.eq.1)) then |
|---|
| 356 | spxxp=0. |
|---|
| 357 | else |
|---|
| 358 | spxxp=cord(2,i)/cord(6,i) |
|---|
| 359 | endif |
|---|
| 360 | spbgx=gbxx(i) |
|---|
| 361 | spyy=cord(3,i) |
|---|
| 362 | if ((ne.eq.0).and.(i.eq.1)) then |
|---|
| 363 | spyyp=0. |
|---|
| 364 | else |
|---|
| 365 | spyyp=cord(4,i)/cord(6,i) |
|---|
| 366 | endif |
|---|
| 367 | spbgy=gbyy(i) |
|---|
| 368 | spzz=cord(5,i) |
|---|
| 369 | spbgz=cord(6,i) |
|---|
| 370 | phi=wt |
|---|
| 371 | wz=(gam(i)-1.)*erest |
|---|
| 372 | c unit ndes2 at constant phase wt |
|---|
| 373 | cbm write(ndes2)spxx,spxxp,spbgx,spyy,spyyp,spbgy,spzz, |
|---|
| 374 | cbm 1 spbgz,phi,wz,ne,i,ngood,i,phizero(i) |
|---|
| 375 | write(ndes2)spxx,spxxp,spbgx,spyy,spyyp,spbgy,spzz, |
|---|
| 376 | 1 spbgz,phi,wz,ne,i,iadp,i,phizero(i),ksi1(i),ksi2(i),ksi3(i) |
|---|
| 377 | if(ne.eq.0) then |
|---|
| 378 | spzz=0. |
|---|
| 379 | phi=phizero(i) |
|---|
| 380 | c unit ndes1 end of element z=cte |
|---|
| 381 | cbm write(ndes1)spxx,spxxp,spbgx,spyy,spyyp,spbgy,spzz, |
|---|
| 382 | cbm 1 spbgz,phi,wz,ne,i,ngood,i,phizero(i) |
|---|
| 383 | write(ndes1,*)spxx,spxxp,spbgx,spyy,spyyp,spbgy,spzz, |
|---|
| 384 | 1 spbgz,phi,wz,ne,i,iadp,i,phizero(i),ksi1(i),ksi2(i),ksi3(i) |
|---|
| 385 | endif |
|---|
| 386 | 9996 continue |
|---|
| 387 | c skip if no hists desired |
|---|
| 388 | ihist=optcon(3)+optcon(4)+optcon(5)+optcon(6)+optcon(7) |
|---|
| 389 | if (ihist.eq.0) return |
|---|
| 390 | cbm 18/07/2k call histogram2 |
|---|
| 391 | return |
|---|
| 392 | end |
|---|
| 393 | c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* |
|---|
![(please configure the [header_logo] section in trac.ini)](/PSPA/chrome/site/pspa.gif)
