1 | subroutine input(lonnom, work_path) |
---|
2 | c-------last modification date 03/27/95 |
---|
3 | c generalized input routine |
---|
4 | c type 0, read np particles from file 'parin' and store |
---|
5 | c in cor |
---|
6 | c type 1, np particles are generated randomly in three |
---|
7 | c phase plane ellipses |
---|
8 | c type 2, np particles are generated randomly in a six |
---|
9 | c dimensional ellipse |
---|
10 | c type 3, np particles are generated uniformily on the |
---|
11 | c perimiter of one of the phase plane ellipses |
---|
12 | c type 4, particle number np is positioned at six specified |
---|
13 | c cordinates |
---|
14 | c type 5, np particles are generated randomly in a four |
---|
15 | c dimensional transverse hyperspace with random |
---|
16 | c phase and energy spread within an ellipse. |
---|
17 | c type 6, np particles are generated randomly in a four |
---|
18 | c dimensional transverse hyperspace with uniform |
---|
19 | c phase and random energy spread. |
---|
20 | c type 7, np particles are generated in a kapchinskij- |
---|
21 | c vladimerskij distribution transversely with |
---|
22 | c random phase and energy spread within an ellipse. |
---|
23 | c type -n,random number generator reset to beginning |
---|
24 | c type 8, np particles are gererated randomly with uniform |
---|
25 | c distribution in real space (x,y,z), then xp, yp, |
---|
26 | c and zp are chosen from within each phase-plane |
---|
27 | c ellipse. .z,zp are then converted to phi,w. |
---|
28 | c type 9, gausian distribution in x, y and z. xp, yp, and energy |
---|
29 | c type 4 spread are zero. parameters are ntype,number of |
---|
30 | c particles,sigmar,rmax,sigmaz(deg),zmax(deg). |
---|
31 | c this used to be the following. |
---|
32 | c rectangular array in one of the three phase planes. |
---|
33 | c ntype,kind,hl,hr,nh,vb,vt,nv |
---|
34 | c type 10 np particles are generated at a photocathode, |
---|
35 | c see comments after line 10000 |
---|
36 | c type 11 same as type 10 uses Hammersley sequence |
---|
37 | c see comments after line 11000 |
---|
38 | c type 12 np particles are generated from output data from |
---|
39 | c EGUN data. parameters are ntype, np, number of rays, |
---|
40 | c phase spread, egun mesh size in cm. The file egun.dat |
---|
41 | c contains the data. The format of the data is: |
---|
42 | c ray#,r,rdot,zdot,tdot,r(cathode),i/r(cathode)/2pi. |
---|
43 | c type 14 same as type 4 with cord(7,i)=0. for i=1,np |
---|
44 | c type 15 hot cathode microwave electron gun |
---|
45 | c type 19 rectangular array |
---|
46 | c type 100*n, as type n, but xp and yp are adjusted to |
---|
47 | c force the angular momentum to be zero. |
---|
48 | c valid for types 1, 2, 5, and 6. |
---|
49 | c vv(1)=type |
---|
50 | c vv(2)=no of particles or particle no |
---|
51 | c vv(3-8)=transverse ellipse parameters |
---|
52 | c vv(9+10)=phase and energy spread (half) |
---|
53 | c vv(11+16)=displacements of input beam |
---|
54 | c phase spread or delta phase is entered in degrees |
---|
55 | c |
---|
56 | c------------------------------------------------------------------------ |
---|
57 | c------------------------------------------------------------------------ |
---|
58 | save |
---|
59 | c |
---|
60 | include 'param_sz.h' |
---|
61 | include 'param_paw.h' |
---|
62 | include 'constcom.h' |
---|
63 | include 'coordcom.h' |
---|
64 | include 'misccom.h' |
---|
65 | include 'ncordscom.h' |
---|
66 | include 'pcordcom.h' |
---|
67 | include 'upawcom.h' |
---|
68 | include 'ucom.h' |
---|
69 | c |
---|
70 | character*256 nomfic |
---|
71 | character*256 work_path |
---|
72 | character*256 nomf |
---|
73 | integer lonnom |
---|
74 | c |
---|
75 | parameter (igun=31) |
---|
76 | common/cf/fmax |
---|
77 | common/cout7/corout7(6,imaa),bgs,nelal,nbufl,iout7 |
---|
78 | common/intype/intype ! mcd |
---|
79 | common/jones/ajones,zjones,zcath(imaa) |
---|
80 | common/out6/nbphase,nbremesh |
---|
81 | common/pawc/h(nxpaw) |
---|
82 | common/part/nparticle(imaa) |
---|
83 | common/phase/uphmin,uphmax,uph(imaa),sigdeg,ipas |
---|
84 | common/ux0/ux0min,ux0max,ux0(imaa) |
---|
85 | common/vx0/vy0min,vy0max,vy0(imaa) |
---|
86 | common/dist/distp(imaa),xr(imaa),yr(imaa) |
---|
87 | common/wtstep/wtstep,nstep |
---|
88 | common/hb2/rmax |
---|
89 | common/hom/xhom(imaa),yhom(imaa) |
---|
90 | dimension xquiet(imaa),yquiet(imaa),anglex(imaa) |
---|
91 | dimension x1(imaa),y1(imaa) |
---|
92 | dimension p1(imaa),p2(imaa),p3(imaa),p4(imaa) |
---|
93 | dimension iwork1(imaa) |
---|
94 | dimension regun(igun),rdot(igun),zdot(igun),tdot(igun), |
---|
95 | * rcath(igun),ioverr(igun),charge(igun) |
---|
96 | double precision dppi,dpepsout |
---|
97 | real f(41) |
---|
98 | external fsch |
---|
99 | c-------------------------------------------------------------------------- |
---|
100 | c* |
---|
101 | data rnumb/0.0/ |
---|
102 | data f/0.,.0398,.0793,.1179,.1554,.1915,.2258,.258,.2881,.3159, |
---|
103 | * .3413,.3643,.3849,.4032,.4192,.4332,.4452,.4554,.4641,.4713, |
---|
104 | * .4772,.4821,.4861,.4893,.4918,.4938,.4953,.4965,.4974,.4981, |
---|
105 | * .4987,.4990,.4993,.4995,.4997,.49975,.4998,.49985,.4999,.49995, |
---|
106 | * .5/ |
---|
107 | print *, ' PARMELA ENTREE INPUT ' |
---|
108 | do 1 i=1,imaa ! mcd |
---|
109 | zcath(i)=0. ! mcd |
---|
110 | 1 continue ! mcd |
---|
111 | r1=1.0 |
---|
112 | r2=1.0 |
---|
113 | r3=1.0 |
---|
114 | r4=1.0 |
---|
115 | r5=1.0 |
---|
116 | r6=1.0 |
---|
117 | phis=0. |
---|
118 | es=w0 |
---|
119 | c fetch first random number |
---|
120 | c if (rnumb.eq.0.0) dum=sandsv (rnumb) |
---|
121 | gg=es/erest |
---|
122 | bgs=sqrt(gg*(2.+gg)) |
---|
123 | bs=sqrt(gg*(2.+gg))/(1.+gg) |
---|
124 | ntype=abs(vv(1)) |
---|
125 | intype=ntype |
---|
126 | if(ip.eq.1 .or. ip.eq.2)write(nnout,2) ntype |
---|
127 | 2 format (' input type ',t40,i8) |
---|
128 | c |
---|
129 | c---- type 0 read np particle |
---|
130 | c |
---|
131 | if(ntype.eq.0)then |
---|
132 | nomf = nomfic(work_path(1:lonnom),'parin.input0') |
---|
133 | print *, ' parmela file : fichier parin : ', nomf |
---|
134 | |
---|
135 | open(unit=3,file=nomf,status='old',err=10) |
---|
136 | npoints=vv(2)+1 |
---|
137 | cor(5,1)=0. |
---|
138 | cor(6,1)=0. |
---|
139 | tcharge=0. |
---|
140 | do 11 j=2,npoints |
---|
141 | read(3,*,err=12,end=12) cor(1,j),cor(2,j),cor(3,j), |
---|
142 | * cor(4,j),cor(5,j),cor(6,j),weight(j) |
---|
143 | weight(j)=abs(weight(j)) |
---|
144 | tcharge=tcharge+weight(j) |
---|
145 | npart=j |
---|
146 | cor(5,1)=cor(5,1)+cor(5,j) |
---|
147 | cor(6,1)=cor(6,1)+cor(6,j) |
---|
148 | write(nnout,*) j,cor(1,j),cor(2,j),cor(3,j),cor(4,j), |
---|
149 | * cor(5,j),cor(6,j),weight(j) ! 02/92 |
---|
150 | 11 continue |
---|
151 | 12 close(3) |
---|
152 | cor(5,1)=cor(5,1)/(npart-1) |
---|
153 | cor(6,1)=cor(6,1)/(npart-1) |
---|
154 | cor(1,1)=0. |
---|
155 | cor(2,1)=0. |
---|
156 | cor(3,1)=0. |
---|
157 | cor(4,1)=0. |
---|
158 | weight(1)=0. |
---|
159 | npoints=npart |
---|
160 | ngood=npoints |
---|
161 | write(nnout,*)' total charge= ',tcharge,' coul' |
---|
162 | return |
---|
163 | 10 continue |
---|
164 | call appendparm |
---|
165 | stop ' Abnormal stop parin not found' |
---|
166 | endif |
---|
167 | c |
---|
168 | if (ntype.gt.100) ntype=ntype/100 |
---|
169 | c |
---|
170 | np1=npoints+1 |
---|
171 | c |
---|
172 | if(ntype.eq.4 .or. ntype.eq.14)go to 400 |
---|
173 | c |
---|
174 | c |
---|
175 | if(ntype.eq.9)go to 900 |
---|
176 | c |
---|
177 | if(ntype.eq.19)go to 19000 |
---|
178 | c |
---|
179 | npoints=vv(2)+npoints |
---|
180 | if(npoints.gt.imaa)npoints=imaa |
---|
181 | c |
---|
182 | if(ntype.eq.3)go to 300 |
---|
183 | c |
---|
184 | if(ntype.eq.10)go to 10000 |
---|
185 | c |
---|
186 | if(ntype.eq.11)go to 11000 |
---|
187 | c |
---|
188 | if(ntype.eq.12)go to 12000 |
---|
189 | c |
---|
190 | if(ntype.eq.15)go to 15000 |
---|
191 | c |
---|
192 | if(ntype.eq.5)then |
---|
193 | rtest=1.2/float(npoints)**.2 |
---|
194 | rtest2=rtest**2 |
---|
195 | itest=0 |
---|
196 | endif |
---|
197 | ax=vv(3) |
---|
198 | bx=vv(4) |
---|
199 | ex=vv(5) |
---|
200 | ay=vv(6) |
---|
201 | by=vv(7) |
---|
202 | ey=vv(8) |
---|
203 | if(ntype.ne.8)go to 15 |
---|
204 | c |
---|
205 | c---- type 8. get ellipse parameters for z, zp. |
---|
206 | c |
---|
207 | az=vv(9) |
---|
208 | bz=vv(10) |
---|
209 | ez=vv(11) |
---|
210 | n1=12 |
---|
211 | gz=(1.+az**2)/bz |
---|
212 | zmx=sqrt(ez/gz) |
---|
213 | zpmx=ez/zmx |
---|
214 | dz=-az/gz |
---|
215 | go to 16 |
---|
216 | 15 continue |
---|
217 | n1=11 |
---|
218 | dphi=vv(9)*radian |
---|
219 | de=vv(10) |
---|
220 | c calculate ellipse parameters |
---|
221 | 16 continue |
---|
222 | gx=(1.0+ax**2)/bx |
---|
223 | gy=(1.0+ay**2)/by |
---|
224 | xmx=sqrt(ex/gx) |
---|
225 | xpmx=sqrt(ex*gx) |
---|
226 | ymx=sqrt(ey/gy) |
---|
227 | dy=-ay/gy |
---|
228 | ypmx=sqrt(ey*gy) |
---|
229 | dx=-ax/gx |
---|
230 | if (vv(1).lt.0.0) call sandst (rnumb) |
---|
231 | c |
---|
232 | c types 1,2,5,6,and 7 |
---|
233 | c |
---|
234 | xnp=vv(2)-1. |
---|
235 | if(xnp.lt.1.)xnp=1. |
---|
236 | area=pi/vv(2) |
---|
237 | 198 continue |
---|
238 | dr5=sqrt((1.-sqrt(1-area**2))/2.) |
---|
239 | itest=0 |
---|
240 | istart=2 |
---|
241 | 199 continue |
---|
242 | do 60 np=np1,npoints |
---|
243 | if(ntype.ne.8)go to 19 |
---|
244 | 200 continue |
---|
245 | r1=2.*ranf()-1. |
---|
246 | r3=2.*ranf()-1. |
---|
247 | r5=2.*ranf()-1. |
---|
248 | rr=r1**2+r3**2+r5**2 |
---|
249 | if (rr.gt.1.)go to 200 |
---|
250 | 201 continue |
---|
251 | r2=2.*ranf()-1. |
---|
252 | r4=2.*ranf()-1. |
---|
253 | r6=2.*ranf()-1. |
---|
254 | if(r2**2+r4**2+r6**2.gt.1.)go to 201 |
---|
255 | sq=sqrt(1.-rr) |
---|
256 | r2=r2*sq |
---|
257 | r4=r4*sq |
---|
258 | r6=r6*sq |
---|
259 | go to 50 |
---|
260 | 19 continue |
---|
261 | if(ntype.eq.7)go to 45 |
---|
262 | 20 continue |
---|
263 | r1=2.0*ranf()-1.0 |
---|
264 | r2=2.0*ranf()-1.0 |
---|
265 | if ((r1**2+r2**2).gt.1.0) go to 20 |
---|
266 | 30 continue |
---|
267 | r3=2.0*ranf()-1.0 |
---|
268 | r4=2.0*ranf()-1.0 |
---|
269 | if ((r3**2+r4**2).gt.1.0) go to 30 |
---|
270 | if (ntype.lt.5) go to 40 |
---|
271 | if (r1**2+r2**2+r3**2+r4**2.gt.1) go to 20 |
---|
272 | 40 continue |
---|
273 | if (ntype.eq.6) r5=2.0*float(np-np1)/xnp-1.0 |
---|
274 | if (ntype.eq.5) r6=sin(abs(acos(r5)))*(ranf()*2.-1) |
---|
275 | c--- this code is to try to make the distribution in real space more |
---|
276 | c--- uniform for ntype 5. |
---|
277 | if(ntype.eq.5)then |
---|
278 | itest=itest+1 |
---|
279 | if(itest.gt.10)then |
---|
280 | write(ndiag,*)'rtest too big' |
---|
281 | rtest=.9*rtest |
---|
282 | rtest2=rtest**2 |
---|
283 | r5=1.0 |
---|
284 | write(nnout,*)'rtest now=',rtest |
---|
285 | go to 198 |
---|
286 | endif |
---|
287 | do 53 i=istart,np-1 |
---|
288 | dxx=cord(1,i)-r1 |
---|
289 | dpx=cord(2,i)-r1 |
---|
290 | dyy=cord(3,i)-r3 |
---|
291 | dpy=cord(4,i)-r4 |
---|
292 | dzz=cord(5,i)-r5 |
---|
293 | if(dxx**2+dpx**2+dyy**2+dpy**2+dzz**2.lt.rtest2)go to 20 |
---|
294 | if(dzz.gt.rtest)istart=i+1 |
---|
295 | 53 continue |
---|
296 | itest=0 |
---|
297 | cord(1,np)=r1 |
---|
298 | cord(2,np)=r2 |
---|
299 | cord(3,np)=r3 |
---|
300 | cord(4,np)=r4 |
---|
301 | cord(5,np)=r5 |
---|
302 | endif |
---|
303 | if (ntype.eq.6) r6=2.0*ranf()-1.0 |
---|
304 | if (ntype.ge.6) go to 50 |
---|
305 | if (ntype.ne.2) go to 50 |
---|
306 | if (r1**2+r2**2+r3**2+r4**2+r5**2+r6**2.gt.1.0) go to 20 |
---|
307 | go to 50 |
---|
308 | c |
---|
309 | c type 7 k-v distribution |
---|
310 | c |
---|
311 | 45 continue |
---|
312 | a=2.0*ranf()-1.0 |
---|
313 | a=acos(a)/2.0 |
---|
314 | b=twopi*(ranf()-0.5) |
---|
315 | c=twopi*(ranf()-0.5) |
---|
316 | r1=cos(a)*cos(b) |
---|
317 | r3=cos(a)*sin(b) |
---|
318 | r2=sin(a)*cos(c) |
---|
319 | r4=sin(a)*sin(c) |
---|
320 | c makes phase and energy distribution into an ellipse |
---|
321 | c 46 r5=2.*ranf()-1. |
---|
322 | c r6=2.*ranf()-1. |
---|
323 | c if((r5*r5+r6*r6).gt.1.) go to 46 |
---|
324 | 46 continue |
---|
325 | r6=sin(abs(acos(r5)))*(ranf()*2.-1.) |
---|
326 | 50 continue |
---|
327 | if (vv(1).le.10.) go to 52 |
---|
328 | if (ntype.gt.6) go to 52 |
---|
329 | if (ntype.eq.3 .or. ntype.eq.4 .or. ntype.eq.14) go to 52 |
---|
330 | c---force zero transverse angular momentum. |
---|
331 | rr=sqrt(r1**2+r3**2) |
---|
332 | rp=sqrt(r2**2+r4**2)*sign(1.,r2) |
---|
333 | r2=rp*r1/rr |
---|
334 | r4=rp*r3/rr |
---|
335 | 52 continue |
---|
336 | cor(1,np)=r1*xmx+dx*r2*xpmx |
---|
337 | cor(2,np)=r2*xpmx |
---|
338 | cor(3,np)=r3*ymx+dy*r4*ypmx |
---|
339 | cor(4,np)=r4*ypmx |
---|
340 | if(ntype.ne.8)go to 55 |
---|
341 | zp=r6*zpmx |
---|
342 | z=r5*zmx+dz*zp |
---|
343 | cor(5,np)=-z*twopi/(bs*wavel)+phis |
---|
344 | cor(6,np)=2.*es*zp+es |
---|
345 | go to 60 |
---|
346 | 55 continue |
---|
347 | cor(5,np)=r5*dphi+phis |
---|
348 | cor(6,np)=r6*de+es |
---|
349 | if(ntype.ne.5.and.ntype.ne.7)go to 60 |
---|
350 | r5=r5-dr5 |
---|
351 | if(r5.le.-1)go to 65 |
---|
352 | sac=sin(acos(r5)) |
---|
353 | dr5=.5*area*sac |
---|
354 | if((r5-dr5).le.-1.)go to 60 |
---|
355 | dr5=area/(sac+sin(acos(r5-dr5))) |
---|
356 | 60 continue |
---|
357 | go to 67 |
---|
358 | 65 continue |
---|
359 | npoints=np |
---|
360 | 67 continue |
---|
361 | if (nn.ge.n1) go to 99998 |
---|
362 | go to 99999 |
---|
363 | c |
---|
364 | c type 3 |
---|
365 | c |
---|
366 | 300 continue |
---|
367 | do 301 n=np1,npoints |
---|
368 | do 302 m=1,4 |
---|
369 | cor(m,n)=0.0 |
---|
370 | 302 continue |
---|
371 | cor(5,n)=phis |
---|
372 | cor(6,n)=es |
---|
373 | 301 continue |
---|
374 | npp=vv(2) |
---|
375 | j=vv(3) |
---|
376 | j=2*j |
---|
377 | gg=(1.+vv(4)**2)/vv(5) |
---|
378 | dd=-vv(4)/gg |
---|
379 | xyzmx=sqrt(vv(6)/gg) |
---|
380 | xyzpmx=vv(6)/xyzmx |
---|
381 | da=twopi/npp |
---|
382 | kk=np1-1 |
---|
383 | do 303 np=1,npp |
---|
384 | ang=float(np-1)*da |
---|
385 | xyzp=xyzpmx*sin(ang) |
---|
386 | xyz=xyzmx*cos(ang)+xyzp*dd |
---|
387 | if(j.eq.6)xyz=xyz*radian |
---|
388 | kk=kk+1 |
---|
389 | cor(j-1,kk)=cor(j-1,kk)+xyz |
---|
390 | cor(j,kk)=xyzp+cor(j,kk) |
---|
391 | 303 continue |
---|
392 | n1=7 |
---|
393 | if (nn.gt.n1) go to 99998 |
---|
394 | go to 99999 |
---|
395 | c |
---|
396 | c type 4 and 14 |
---|
397 | c |
---|
398 | 400 continue |
---|
399 | c input ntype x xp y yp phi w ksi1 ksi2 ksi3 n4 |
---|
400 | c the first 7 parameters must be on all cards input |
---|
401 | c n4 is the total number of input cards |
---|
402 | c the first card must have 11 parameters |
---|
403 | c corrige le 10/12/2008 (nn.eq.8)-->(nn.eq.11) |
---|
404 | if(nn.eq.11) then |
---|
405 | maxc=vv(11) |
---|
406 | ncarte=1 |
---|
407 | else |
---|
408 | ncarte=ncarte+1 |
---|
409 | endif |
---|
410 | npoints=np1 |
---|
411 | if(npoints.gt.imaa)npoints=imaa |
---|
412 | do 401 n=1,4 |
---|
413 | cor(n,npoints)=vv(n+1) |
---|
414 | 401 continue |
---|
415 | cor(5,npoints)=vv(6)*radian+phis |
---|
416 | cor(6,npoints)=vv(7)+es |
---|
417 | |
---|
418 | phizero(npoints)=vv(6) |
---|
419 | go to 99999 |
---|
420 | c |
---|
421 | c---- type 9. gausian array. |
---|
422 | c |
---|
423 | 900 continue |
---|
424 | np1=npoints+1 |
---|
425 | npoints=npoints+vv(2) |
---|
426 | if(npoints.gt.imaa)npoints=imaa |
---|
427 | sigmar=vv(3) |
---|
428 | rmax=vv(4) |
---|
429 | sigmaz=vv(5) |
---|
430 | zmax=vv(6) |
---|
431 | ss=-1. |
---|
432 | rnorm=1.41421356/sigmaz |
---|
433 | zzmax=rnorm*zmax |
---|
434 | if(zzmax.gt.4.)zzmax=3.999999999 |
---|
435 | iz=zzmax*10. |
---|
436 | a=zzmax*10.-iz |
---|
437 | fmax=f(iz+1)*(1.-a)+f(iz+2)*a |
---|
438 | do 901 i=np1,npoints |
---|
439 | if(sigmar.gt.10.*rmax)then |
---|
440 | r1=sqrt(ranf())*rmax |
---|
441 | else |
---|
442 | 902 continue |
---|
443 | r1=sigmar*sqrt(-log(ranf())) |
---|
444 | if(r1.gt.rmax)go to 902 |
---|
445 | endif |
---|
446 | theata=twopi*ranf() |
---|
447 | cor(1,i)=sin(theata)*r1 |
---|
448 | cor(3,i)=cos(theata)*r1 |
---|
449 | if(sigmaz.gt.10.*zmax)then |
---|
450 | cor(5,i)=2.*radian*zmax*(float(i-np1)/float(npoints-np1)-.5) |
---|
451 | else |
---|
452 | rand=fmax*float(i-np1)/float(npoints-np1) |
---|
453 | do 903 j=2,41 |
---|
454 | if(rand.le.f(j))go to 904 |
---|
455 | 903 continue |
---|
456 | 904 continue |
---|
457 | r2=zmax/zzmax*(float(j-2)+(rand-f(j-1))/(f(j)-f(j-1)))*.1 |
---|
458 | ss=-ss |
---|
459 | cor(5,i)=r2*radian*ss |
---|
460 | endif |
---|
461 | cor(2,i)=0. |
---|
462 | cor(4,i)=0. |
---|
463 | cor(6,i)=es |
---|
464 | 901 continue |
---|
465 | go to 99999 |
---|
466 | c |
---|
467 | c type 19. rectangular array. |
---|
468 | c |
---|
469 | 19000 kind=vv(2) |
---|
470 | hl=vv(3) |
---|
471 | hr=vv(4) |
---|
472 | if(kind.eq.3)hl=hl*radian |
---|
473 | if(kind.eq.3)hr=hr*radian |
---|
474 | nh=vv(5) |
---|
475 | vb=vv(6) |
---|
476 | vt=vv(7) |
---|
477 | nv=vv(8) |
---|
478 | k1=2*kind-1 |
---|
479 | k2=k1+1 |
---|
480 | npoints=np1+nh*nv-1 |
---|
481 | if (npoints.gt.imaa) npoints=imaa |
---|
482 | do 19001 np=np1,npoints |
---|
483 | cor(5,np)=phis |
---|
484 | cor(6,np)=es |
---|
485 | do 19001 n=1,4 |
---|
486 | cor(n,np)=0. |
---|
487 | 19001 continue |
---|
488 | hz=hl |
---|
489 | if (k1.eq.5) hz=hz+phis |
---|
490 | vz=vb |
---|
491 | if (k2.eq.6) vz=vz+es |
---|
492 | dh=0. |
---|
493 | if (nh.gt.1) dh=(hr-hl)/(nh-1) |
---|
494 | dv=0. |
---|
495 | if (nv.gt.1) dv=(vt-vb)/(nv-1) |
---|
496 | np=np1-1 |
---|
497 | hh=hz-dh |
---|
498 | do 19002 i=1,nh |
---|
499 | hh=hh+dh |
---|
500 | v=vz-dv |
---|
501 | do 19002 j=1,nv |
---|
502 | np=np+1 |
---|
503 | if (np.gt.imaa) go to 19003 |
---|
504 | v=v+dv |
---|
505 | cor(k1,np)=hh |
---|
506 | cor(k2,np)=v |
---|
507 | 19002 continue |
---|
508 | 19003 continue |
---|
509 | n1=9 |
---|
510 | if(nn.lt.n1)go to 99999 |
---|
511 | c |
---|
512 | c displace input beam |
---|
513 | c |
---|
514 | 99998 continue |
---|
515 | n2=n1+4 |
---|
516 | if(nn.ge.n2)vv(n2)=vv(n2)*radian |
---|
517 | do 99997 n=n1,nn |
---|
518 | k=n-n1+1 |
---|
519 | do 99997 np=np1,npoints |
---|
520 | cor(k,np)=cor(k,np)+vv(n) |
---|
521 | 99997 continue |
---|
522 | go to 99999 |
---|
523 | c---convert from x,xp,y,yp,phi,w to x,bgx,y,bgy,z,bgz |
---|
524 | 99999 continue |
---|
525 | ngood=npoints |
---|
526 | nbufl=npoints |
---|
527 | c ii=2,npoints reference particle is treated in main program |
---|
528 | if (ntype.ne.4 .or. ntype.eq.14) then |
---|
529 | do 99990 ii=2,npoints |
---|
530 | phizero(ii)=cor(5,ii) |
---|
531 | do 99991 iii=1,6 |
---|
532 | corout7(iii,ii)=cor(iii,ii) |
---|
533 | 99991 continue |
---|
534 | 99990 continue |
---|
535 | nelal=0 |
---|
536 | iout7=iout7+1 |
---|
537 | call outlaldp |
---|
538 | iout7=iout7+1 |
---|
539 | endif |
---|
540 | g=w0/erest |
---|
541 | zcon=-wavel/twopi |
---|
542 | do 99992 np=np1,npoints |
---|
543 | phi=cor(5,np)/radian |
---|
544 | phizero(np)=phi |
---|
545 | wz=cor(6,np) |
---|
546 | g=wz/erest |
---|
547 | bgz=sqrt(g*(2.+g)) |
---|
548 | dz=zcon*cor(5,np)*bgz/(1.+g) |
---|
549 | bgx=cor(2,np)*bgz |
---|
550 | bgy=cor(4,np)*bgz |
---|
551 | cor(1,np)=cor(1,np)+cor(2,np)*(z0+dz) |
---|
552 | cor(2,np)=bgx |
---|
553 | cor(3,np)=cor(3,np)+cor(4,np)*(z0+dz) |
---|
554 | cor(4,np)=bgy |
---|
555 | cor(5,np)=z0+dz |
---|
556 | cor(6,np)=sqrt(bgz**2-bgx**2-bgy**2) |
---|
557 | nparticle(np)=np |
---|
558 | 99992 continue |
---|
559 | c for emittance summary, copy parameters into cord & gam |
---|
560 | c note that officially this is done at line 210 of the main program |
---|
561 | do 99993 i = np1,npoints |
---|
562 | cord(7,i)=0. |
---|
563 | do 99994 j = 1,6 |
---|
564 | cord(j,i) = cor(j,i) |
---|
565 | 99994 continue |
---|
566 | if(ntype.eq.14) cord(7,i)=0.0 |
---|
567 | gam(i) = sqrt(1. + cord(2,i)**2 + cord(4,i)**2 + cord(6,i)**2) |
---|
568 | 99993 continue |
---|
569 | write(nnout,*) ' ncarte ',ncarte,' maxc ',maxc |
---|
570 | if((ntype.eq.4 .or. ntype.eq.14).and.(ncarte.lt.maxc)) return |
---|
571 | if(ip.ne.0) then |
---|
572 | 8337 continue |
---|
573 | cbm call ou6dp normal |
---|
574 | call out6dp(0,0) |
---|
575 | endif |
---|
576 | return ! return input 1 a 9 |
---|
577 | c |
---|
578 | c---type 10 = particles generated at a photocathode |
---|
579 | c |
---|
580 | 10000 continue |
---|
581 | c last modification 03/08/90 |
---|
582 | c vv(2) = no. of particles |
---|
583 | c vv(3) = sigma of laser pulse duration in picoseconds |
---|
584 | c vv(4) = half width of laser pulse. i.e. generate a gaussian |
---|
585 | c time profile with sigma = vv(3), but truncate to +-vv(4) |
---|
586 | c vv(5) = sigma of laser spot in cm - assumes gaussian radial profile |
---|
587 | c vv(6) = max radius of photocathode. i.e. truncate radial dist at vv(6) |
---|
588 | c vv(7) = mean K.E. of electrons leaving the photocathode, in MeV |
---|
589 | c This should be the same as paramter W0 on the RUN card. |
---|
590 | c vv(8) = sigma of energy-spread distribution at the photocathode (MeV). |
---|
591 | c vv(9) = rf step size in deg. This must be the same as parmater DWT |
---|
592 | c on the START card |
---|
593 | c vv(10) = type of random (new) 1 by rannor |
---|
594 | c 2 x,y by rannor phi by didtribution |
---|
595 | c vv(11) = spota -- Angle of laser beam with the normal at the cathode |
---|
596 | c (degrees), if non-zero you need vv(12) |
---|
597 | c vv(12) = spotd -- Distance from focal point to cathode (cm) |
---|
598 | c if no zero you need vv(11) |
---|
599 | c vv(13) = number of chenals for histo (new) |
---|
600 | c vv(13) nonzero you need vv(11) and vv(13) egal or not zero |
---|
601 | c vv(14) = mu -- the cathode shape parameter if Mike Jones. Must have |
---|
602 | c 0 < mu < .5. If mu is nonzero, the field shape of the cell |
---|
603 | c with element no. 1 will be taken from Jones' analytic form, |
---|
604 | c overriding the fourier coefficients, if any |
---|
605 | c vv(15) = length of cell with cathode as front surface (cm). Must be |
---|
606 | c the same a L on the first CELL card |
---|
607 | c |
---|
608 | c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
609 | |
---|
610 | if (nn.lt.10) then |
---|
611 | write(nnout,10001)nn |
---|
612 | 10001 format (' Too few parameters for input type 10 - nn = ',i2, |
---|
613 | 1 ' parameters only 10 are needed') |
---|
614 | call appendparm |
---|
615 | stop ' Abnormal stop input 10 a ' |
---|
616 | endif |
---|
617 | if (nn.lt.11.and.ip.eq.2) then |
---|
618 | write(nnout,10002)ip,nn |
---|
619 | 10002 format (' Too few parameters for input type 10 - ip = ', |
---|
620 | 1 i1,' nn = ',i2,' parameters only',/,' 11 are needed ') |
---|
621 | call appendparm |
---|
622 | stop ' Abnormal stop input 10 b ' |
---|
623 | endif |
---|
624 | tsig = vv(3) |
---|
625 | tmax = vv(4) |
---|
626 | if (tsig.le.0.) then |
---|
627 | write(nnout,10003) |
---|
628 | 10003 format (' sigma of laser pulse must be > 0') |
---|
629 | call appendparm |
---|
630 | stop ' Abnormal stop input 10 c ' |
---|
631 | endif |
---|
632 | tfac = tmax/tsig |
---|
633 | rsig = vv(5) |
---|
634 | rmax = vv(6) |
---|
635 | rmaxsq = rmax**2 |
---|
636 | ecath = vv(7) |
---|
637 | decath = vv(8) |
---|
638 | c convert laser pulse sigma to degrees |
---|
639 | sigdeg = tsig*freq*360./1.e6 |
---|
640 | c generate particles during +- tmax in time |
---|
641 | c set the starting clock back by an integer number of time steps, |
---|
642 | c corresponding to at least tmax |
---|
643 | degmax = tmax*freq*360./1.e6 |
---|
644 | dwt = vv(9) |
---|
645 | nstep = int(degmax/dwt) |
---|
646 | if (degmax.ne.dwt*float(nstep)) nstep = nstep + 1 |
---|
647 | wtstep = nstep*dwt |
---|
648 | iran=vv(10) |
---|
649 | if (ip.eq.1 .or. ip.eq.2) |
---|
650 | . write(nnout,10004) vv(2),tsig,sigdeg,tmax,rsig,rmax,ecath, |
---|
651 | . decath,dwt,wtstep |
---|
652 | 10004 format ( |
---|
653 | . ' Generating ',f10.0,' particles from a photocathode'/ |
---|
654 | . ' sigma of laser pulse =',t40,1pg12.3,' psec'/ |
---|
655 | . ' sigma of laser pulse =',t40,g12.3,' degrees'/ |
---|
656 | . ' tmax of laser pulse =',t40,g12.3,' psec'/ |
---|
657 | . ' sigma of laser spot =',t40,g12.3,' cm'/ |
---|
658 | . ' rmax of photocathode =',t40,g12.3,' cm'/ |
---|
659 | . ' mean electron energy at cathode =',t40,g12.3,' MeV'/ |
---|
660 | . ' sigma of energy spread at cathode =',t40,g12.3,' MeV'/ |
---|
661 | . ' integration step size =',t40,g12.3,' degrees'/ |
---|
662 | . ' clock offset for earliest particle =',t40,g12.3,' degrees') |
---|
663 | if (vv(10).eq.1) write(nnout,10019) iran |
---|
664 | 10019 format (' generated by rannor only ran = ',t45,i2) |
---|
665 | if (vv(10).eq.2) write(nnout,10020) iran |
---|
666 | 10020 format (' generated by distriphase-rannor ran = ',t45,i2) |
---|
667 | if (vv(10).gt.2) then |
---|
668 | write(nnout,*) ' Input 10 ran.gt.2 ran = ',iran |
---|
669 | call appendparm |
---|
670 | stop ' Abnormal stop input 10 d ' |
---|
671 | endif |
---|
672 | c---- laser beam angle |
---|
673 | if(nn.ge.11) then |
---|
674 | spota=vv(11)*radian |
---|
675 | spotd=vv(12) |
---|
676 | if (ip.eq.1 .or. ip.eq.2) then |
---|
677 | write(nnout,10025) spota/radian,spotd |
---|
678 | 10025 format( |
---|
679 | . ' angle of injection of laser light =',t40,g12.3,' degrees'/ |
---|
680 | . ' distance focal point - cathode =',t40,g12.3,' cm') |
---|
681 | endif |
---|
682 | endif |
---|
683 | c---- |
---|
684 | if (nn.ge.14) then |
---|
685 | ajones = vv(14) |
---|
686 | if (ajones.lt.0.0 .or. ajones.gt.0.499) then |
---|
687 | write(nnout,10005) ajones |
---|
688 | 10005 format (' cathode shape parameter = ',1pg12.3, |
---|
689 | . ' is out of range') |
---|
690 | call appendparm |
---|
691 | stop ' Abnormal stop input 10 e ' |
---|
692 | endif |
---|
693 | zjones = 0. |
---|
694 | if (ip.eq.1 .or. ip.eq.2 .or. ip.eq.4) |
---|
695 | . write(nnout,10006) ajones |
---|
696 | 10006 format (' cathode shape parameter =',t40,1pg12.3) |
---|
697 | endif |
---|
698 | if (nn.ge.15) then |
---|
699 | zjones = vv(15) |
---|
700 | if (ip.eq.1 .or. ip.eq.2 .or. ip.eq.4) |
---|
701 | . write(nnout,10007) zjones |
---|
702 | 10007 format (' gun cavity length =',t40,1pg12.3,' cm') |
---|
703 | endif |
---|
704 | |
---|
705 | |
---|
706 | c----- |
---|
707 | c |
---|
708 | c offset z of reference particle so it reaches z = 0 at center of pulse |
---|
709 | bgr = cor(6,1) |
---|
710 | bzr = bgr/sqrt(1. + bgr**2) |
---|
711 | cor(5,1) = -wtstep*wavel*bzr/360. |
---|
712 | c----- |
---|
713 | uphmin=0. |
---|
714 | uphmax=0. |
---|
715 | ux0min=0. |
---|
716 | ux0max=0. |
---|
717 | vy0min=0. |
---|
718 | vy0max=0. |
---|
719 | c----------------------------- ran = 1 ---- |
---|
720 | if(vv(10).eq.1) then |
---|
721 | do 10008 i = np1,npoints |
---|
722 | c choose time (phase) of generation, relative to center of the pulse |
---|
723 | 10009 continue |
---|
724 | call rannor(u,v) |
---|
725 | if (abs(u).gt.tfac) go to 10009 |
---|
726 | phase = u*sigdeg |
---|
727 | c choose intial energy at cathode |
---|
728 | 10010 continue |
---|
729 | call rannor(u,v) |
---|
730 | e = ecath + u*decath |
---|
731 | if (e.le.0.) go to 10010 |
---|
732 | betcat = sqrt(2.*e/erest) |
---|
733 | gamma = 1./sqrt(1. - betcat**2) |
---|
734 | c choose initial momentum, isotropic, but betaz > 0 |
---|
735 | cose = sandom(d) |
---|
736 | sine = sqrt(1. - cose**2) |
---|
737 | phi = twopi*sandom(d) |
---|
738 | bx = betcat*sine*cos(phi) |
---|
739 | by = betcat*sine*sin(phi) |
---|
740 | bz = betcat*cose |
---|
741 | cor(2,i) = gamma*bx |
---|
742 | cor(4,i) = gamma*by |
---|
743 | cor(6,i) = gamma*bz |
---|
744 | c choose (x,y,z) |
---|
745 | 10011 continue |
---|
746 | call rannor(u,v) |
---|
747 | x0 = u*rsig |
---|
748 | y0 = v*rsig |
---|
749 | rsq = x0**2 + y0**2 |
---|
750 | if (rsq.gt.rmaxsq) go to 10011 |
---|
751 | c---------------------------------- |
---|
752 | if((spota.eq.0.).or.(spotd.eq.0)) go to 10021 |
---|
753 | anglex10 = atan(x0/spotd) |
---|
754 | delta = spotd*(1-(cos(spota)/cos(spota+(anglex10)))) |
---|
755 | delta = delta*360./wavel |
---|
756 | x0= spotd*sin(anglex10)/cos(anglex10+spota) |
---|
757 | phase=phase+delta |
---|
758 | 10021 continue |
---|
759 | c--------------------------------------- |
---|
760 | uph(i)=phase |
---|
761 | if(uph(i).gt.uphmax)uphmax=uph(i) |
---|
762 | if(uph(i).lt.uphmin)uphmin=uph(i) |
---|
763 | c--------------------------------------- |
---|
764 | ux0(i)=x0 |
---|
765 | if(ux0(i).gt.ux0max)ux0max=ux0(i) |
---|
766 | if(ux0(i).lt.ux0min)ux0min=ux0(i) |
---|
767 | vy0(i)=y0 |
---|
768 | if(vy0(i).gt.vy0max)vy0max=vy0(i) |
---|
769 | if(vy0(i).lt.vy0min)vy0min=vy0(i) |
---|
770 | c----- |
---|
771 | z0 = 0. |
---|
772 | if (ajones.gt.0.) then |
---|
773 | c calculate the shape of the cathode surface |
---|
774 | z0 = cathod(rsq) |
---|
775 | zcath(i) = z0 |
---|
776 | endif |
---|
777 | c----- |
---|
778 | c shift initial (x,y,z) so that the particle will emerge from |
---|
779 | c the photocathode at (x0,y0,z0) at time = phase |
---|
780 | toff = (wtstep + phase)*wavel/360. |
---|
781 | cor(1,i) = x0 - toff*bx |
---|
782 | cor(3,i) = y0 - toff*by |
---|
783 | cor(5,i) = z0 - toff*bz |
---|
784 | 10008 continue |
---|
785 | else ! vv(10).ne.1 |
---|
786 | c---- |
---|
787 | call distriphase(degmax,sigdeg,npoints-1) |
---|
788 | call distrir(rmax,rsig,npoints-1) |
---|
789 | c---- |
---|
790 | do 10012 i = np1,npoints |
---|
791 | c choose time (phase) of generation, relative to center of the pulse |
---|
792 | c---- |
---|
793 | u=distp(i-1) |
---|
794 | c---- |
---|
795 | phase=u |
---|
796 | c---- |
---|
797 | c choose intial energy at cathode |
---|
798 | 10013 continue |
---|
799 | call rannor(u,v) |
---|
800 | c---- |
---|
801 | e = ecath + u*decath |
---|
802 | if (e.le.0.) go to 10013 |
---|
803 | betcat = sqrt(2.*e/erest) |
---|
804 | gamma = 1./sqrt(1. - betcat**2) |
---|
805 | c choose initial momentum, isotropic, but betaz > 0 |
---|
806 | cose = sandom(d) |
---|
807 | sine = sqrt(1. - cose**2) |
---|
808 | phi = twopi*sandom(d) |
---|
809 | bx = betcat*sine*cos(phi) |
---|
810 | by = betcat*sine*sin(phi) |
---|
811 | bz = betcat*cose |
---|
812 | cor(2,i) = gamma*bx |
---|
813 | cor(4,i) = gamma*by |
---|
814 | cor(6,i) = gamma*bz |
---|
815 | c choose (x,y,z) |
---|
816 | c---- |
---|
817 | c---- |
---|
818 | if(vv(10).eq.2) then |
---|
819 | 10014 continue |
---|
820 | call rannor(u,v) |
---|
821 | x0=u*rsig |
---|
822 | y0=v*rsig |
---|
823 | rsq=x0**2+y0**2 |
---|
824 | if(rsq.gt.rmaxsq) go to 10014 |
---|
825 | else ! vv(10) .ne.2 |
---|
826 | x0 = xr(i-1) |
---|
827 | y0 = yr(i-1) |
---|
828 | endif ! vv(10) .eq. 2 |
---|
829 | c---------------------------------- |
---|
830 | if((spota.eq.0.).or.(spotd.eq.0)) go to 10022 |
---|
831 | anglex10 = atan(x0/spotd) |
---|
832 | delta = spotd*(1-(cos(spota)/cos(spota+(anglex10)))) |
---|
833 | delta = delta*360./wavel |
---|
834 | x0= spotd*sin(anglex10)/cos(anglex10+spota) |
---|
835 | phase=phase+delta |
---|
836 | 10022 continue |
---|
837 | c--------------------------------------- |
---|
838 | uph(i)=phase |
---|
839 | if(uph(i).gt.uphmax)uphmax=uph(i) |
---|
840 | if(uph(i).lt.uphmin)uphmin=uph(i) |
---|
841 | c---- |
---|
842 | ux0(i)=x0 |
---|
843 | if(ux0(i).gt.ux0max)ux0max=ux0(i) |
---|
844 | if(ux0(i).lt.ux0min)ux0min=ux0(i) |
---|
845 | vy0(i)=y0 |
---|
846 | if(vy0(i).gt.vy0max)vy0max=vy0(i) |
---|
847 | if(vy0(i).lt.vy0min)vy0min=vy0(i) |
---|
848 | c---- |
---|
849 | z0 = 0. |
---|
850 | if (ajones.gt.0.) then |
---|
851 | rsq = x0**2 + y0**2 |
---|
852 | c calculate the shape of the cathode surface |
---|
853 | z0 = cathod(rsq) |
---|
854 | zcath(i) = z0 |
---|
855 | endif |
---|
856 | c shift initial (x,y,z) so that the particle will emerge from |
---|
857 | c the photocathode at (x0,y0,z0) at time = phase |
---|
858 | toff = (wtstep + phase)*wavel/360. |
---|
859 | cor(1,i) = x0 - toff*bx |
---|
860 | cor(3,i) = y0 - toff*by |
---|
861 | cor(5,i) = z0 - toff*bz |
---|
862 | 10012 continue |
---|
863 | |
---|
864 | c---------- |
---|
865 | endif ! vv(1).eq.1 |
---|
866 | c---------- |
---|
867 | |
---|
868 | if(ip.eq.2) then |
---|
869 | open(unit=ninput,file='input10',status='unknown') |
---|
870 | write(ninput,*) npoints,sigdeg,rmax |
---|
871 | do 10015 i=1,npoints |
---|
872 | write(ninput,*) uph(i),ux0(i),vy0(i) |
---|
873 | 10015 continue |
---|
874 | endif |
---|
875 | cbm 18/07/2k call histogram1 |
---|
876 | c---------- |
---|
877 | c for emittance summary, copy parameters into cord & gam |
---|
878 | c note that officially this is done at line 210 of the main program |
---|
879 | ngood = npoints |
---|
880 | do 10016 i = 1,npoints |
---|
881 | do 10017 j = 1,6 |
---|
882 | cord(j,i) = cor(j,i) |
---|
883 | 10017 continue |
---|
884 | phizero(i)=uph(i) |
---|
885 | cord(7,i) = 0. |
---|
886 | gam(i) = sqrt(1. + cord(2,i)**2 + cord(4,i)**2 + cord(6,i)**2) |
---|
887 | 10016 continue |
---|
888 | write(nnout,10018) |
---|
889 | 10018 format (/' emittance summary for generated particles:') |
---|
890 | iout7=1 |
---|
891 | call outlaldp |
---|
892 | iout7=iout7+1 |
---|
893 | cbm call out6dp pour type 10 |
---|
894 | call out6dp(0,0) |
---|
895 | return |
---|
896 | c |
---|
897 | c--- type 11 |
---|
898 | c |
---|
899 | c---type 11 = same as input 10 with smooth distributions |
---|
900 | c uses hammersley sequence |
---|
901 | c |
---|
902 | 11000 continue |
---|
903 | c vv(2) = no. of particles |
---|
904 | c vv(3) = sigma of laser pulse duration in picoseconds |
---|
905 | c vv(4) = half width of laser pulse. i.e. generate a gaussian |
---|
906 | c time profile with sigma = vv(3), but truncate to +-vv(4 |
---|
907 | c NOT EFFECTIVE in this version. The gaussian is |
---|
908 | c truncated at roughly +-2.5 sigma. |
---|
909 | c vv(5) = sigma of laser spot in cm - assumes gaussian radial |
---|
910 | c profile NOT EFFECTIVE. |
---|
911 | c vv(6) = max radius of photocathode. i.e. truncate radial dist a |
---|
912 | c vv(7) = mean K.E. of electrons leaving the photocathode, in MeV |
---|
913 | c This should be the same as paramter W0 on the RUN card. |
---|
914 | c vv(8) = sigma of energy-spread distribution at the photocathode |
---|
915 | c vv(9) = rf step size in deg. This must be the same as parmater |
---|
916 | c on the START card |
---|
917 | c--------------- |
---|
918 | c The following parameters define the way the particles are |
---|
919 | c generated in each phase space and the correlation between |
---|
920 | c them. The program transform the decimal number which reper |
---|
921 | c the particle (j=1,2,....,N) into its equivalent in another |
---|
922 | c base, inverses each of the digit and reconverts the number |
---|
923 | c into decimal. This leads for example in base 3: (1/3,2/3,1/9, |
---|
924 | c 1/3+1/9,....). The base must be a prime number. |
---|
925 | c If the base is 0 the program just produces a random set using |
---|
926 | c the usual random fortran generator: RANF(). |
---|
927 | c vv(10) = bit reversal base for x distribution |
---|
928 | c if vv(10)<0 don't use bit reversal but produce a |
---|
929 | c uniform set in x and y. |
---|
930 | c vv(11)= bit reversal base for y distribution |
---|
931 | c vv(12)= bit reversal base for phi distribution |
---|
932 | c if vv(12)<0 call distriphase as in INPUT 10 |
---|
933 | c vv(13)= bit reversal base for W0 distribution |
---|
934 | c Not in use in this version. Same as INPUT 10 |
---|
935 | c vv(14) = spota -- Angle of laser beam with the normal at the cathode |
---|
936 | c (degrees), if non-zero you need vv(15) |
---|
937 | c vv(15) = spotd -- Distance from focal point to cathode (cm) |
---|
938 | c if no zero you need vv(14) |
---|
939 | c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
---|
940 | tsig = vv(3) |
---|
941 | tmax = vv(4) |
---|
942 | rsig = vv(5) |
---|
943 | rmax = vv(6) |
---|
944 | rmaxsq = rmax**2 |
---|
945 | ecath = vv(7) |
---|
946 | decath = vv(8) |
---|
947 | nbasex=vv(10) |
---|
948 | nbasey=vv(11) |
---|
949 | nbasep=vv(12) |
---|
950 | nbasew=vv(13) |
---|
951 | if(nbasep.ge.0) tmax=2.5 |
---|
952 | tfac=tmax*tsig |
---|
953 | c convert laser pulse sigma to degrees |
---|
954 | sigdeg = tsig*freq*360./1.e6 |
---|
955 | c generate particles during +- tmax in time |
---|
956 | c set the starting clock back by an integer number of time steps, |
---|
957 | c corresponding to at least tmax |
---|
958 | degmax = tfac*freq*360./1.e6 |
---|
959 | dwt = vv(9) |
---|
960 | nstep = int(degmax/dwt) |
---|
961 | if (degmax.ne.dwt*float(nstep)) nstep = nstep + 1 |
---|
962 | wtstep = nstep*dwt |
---|
963 | c |
---|
964 | c offset z of reference particle so it reaches z = 0 at center of |
---|
965 | cg it will reach z0 not 0 |
---|
966 | bgr = cor(6,1) |
---|
967 | bzr = bgr/sqrt(1. + bgr**2) |
---|
968 | cg cor(5,1) = -wtstep*wavel*bzr/360. |
---|
969 | cor(5,1) = z0-wtstep*wavel*bzr/360. |
---|
970 | c----- |
---|
971 | if(nbasex.lt.0)then |
---|
972 | call distrirhom(rmax,npoints-1,flag,nhom) |
---|
973 | npoints=nhom+1 |
---|
974 | else |
---|
975 | call loduni(nbasex,npoints-1,p1) |
---|
976 | call loduni(nbasey,npoints-1,p2) |
---|
977 | j=1 |
---|
978 | do 11002 i=np1,npoints |
---|
979 | x1(i)=(2*p1(i)-1)*rmax |
---|
980 | y1(i)=(2*p2(i)-1)*rmax |
---|
981 | r=sqrt(x1(i)**2+y1(i)**2) |
---|
982 | if(r.le.rmax)then |
---|
983 | xquiet(j)=x1(i) |
---|
984 | yquiet(j)=y1(i) |
---|
985 | j=j+1 |
---|
986 | endif |
---|
987 | 11002 continue |
---|
988 | npoints=j+1 |
---|
989 | endif |
---|
990 | c---- laser beam angle |
---|
991 | if(nn.ge.14) then |
---|
992 | spota=vv(14)*radian |
---|
993 | spotd=vv(15) |
---|
994 | do 11020 j=1,npoints |
---|
995 | anglex(j) = atan(xquiet(j)/spotd) |
---|
996 | xquiet(j)= spotd*sin(anglex(j))/cos(anglex(j)+spota) |
---|
997 | 11020 continue |
---|
998 | endif |
---|
999 | c---- |
---|
1000 | uphmin=0. |
---|
1001 | uphmax=0. |
---|
1002 | if(nbasep.lt.0)then |
---|
1003 | call distriphase(degmax,sigdeg,npoints-1) |
---|
1004 | else |
---|
1005 | call lodgau(nbasep,npoints-1,p1,p2,p3,p4,iwork1) |
---|
1006 | endif |
---|
1007 | do 11003 i=np1,npoints |
---|
1008 | c choose phase |
---|
1009 | if(nbasep.lt.0)then |
---|
1010 | phase=distp(i-1) |
---|
1011 | else |
---|
1012 | phase=p1(i-1)*sigdeg |
---|
1013 | endif |
---|
1014 | c---- |
---|
1015 | if((spota.eq.0.).or.(spotd.eq.0)) go to 11023 |
---|
1016 | delta = spotd*(1-(cos(spota)/cos(spota+(anglex(i))))) |
---|
1017 | delta = delta*360./wavel |
---|
1018 | phase=phase+delta |
---|
1019 | 11023 continue |
---|
1020 | uph(i)=phase |
---|
1021 | phizero(i)=uph(i) |
---|
1022 | if(uph(i).gt.uphmax)uphmax=uph(i) |
---|
1023 | if(uph(i).lt.uphmin)uphmin=uph(i) |
---|
1024 | c choose intial energy at cathode |
---|
1025 | 11004 continue |
---|
1026 | call rannor(u,v) |
---|
1027 | e = ecath + u*decath |
---|
1028 | if (e.le.0.) go to 11004 |
---|
1029 | betcat = sqrt(2.*e/erest) |
---|
1030 | cbm warning if 1.-betcat**2 < 0 !!!!! |
---|
1031 | gamma = 1./sqrt(1. - betcat**2) |
---|
1032 | c choose initial momentum, isotropic, but betaz > 0 |
---|
1033 | cose = sandom(d) |
---|
1034 | sine = sqrt(1. - cose**2) |
---|
1035 | phi = twopi*sandom(d) |
---|
1036 | bx = betcat*sine*cos(phi) |
---|
1037 | by = betcat*sine*sin(phi) |
---|
1038 | bz = betcat*cose |
---|
1039 | cor(2,i) = gamma*bx |
---|
1040 | cor(4,i) = gamma*by |
---|
1041 | cor(6,i) = gamma*bz |
---|
1042 | c |
---|
1043 | c choose (x,y,z) |
---|
1044 | c |
---|
1045 | if(nbasex.lt.0)then |
---|
1046 | x0=xhom(i-1) |
---|
1047 | y0=yhom(i-1) |
---|
1048 | else |
---|
1049 | x0=xquiet(i-1) |
---|
1050 | y0=yquiet(i-1) |
---|
1051 | endif |
---|
1052 | c----------- |
---|
1053 | ux0(i)=x0 |
---|
1054 | if(ux0(i).gt.ux0max)ux0max=ux0(i) |
---|
1055 | if(ux0(i).lt.ux0min)ux0min=ux0(i) |
---|
1056 | vy0(i)=y0 |
---|
1057 | if(vy0(i).gt.vy0max)vy0max=vy0(i) |
---|
1058 | if(vy0(i).lt.vy0min)vy0min=vy0(i) |
---|
1059 | c----------- |
---|
1060 | c |
---|
1061 | c shift initial (x,y,z) so that the particle will emerge from |
---|
1062 | c the photocathode at (x0,y0,z0) at time = phase |
---|
1063 | toff = (wtstep + phase)*wavel/360. |
---|
1064 | cor(1,i) = x0 - toff*bx |
---|
1065 | cor(3,i) = y0 - toff*by |
---|
1066 | cor(5,i) = z0 - toff*bz |
---|
1067 | 11003 continue |
---|
1068 | c---------- |
---|
1069 | c for emittance summary, copy parameters into cord & gam |
---|
1070 | c note that officially this is done at line 210 of the main program |
---|
1071 | ngood = npoints |
---|
1072 | do 11016 i = 1,npoints |
---|
1073 | do 11017 j = 1,6 |
---|
1074 | cord(j,i) = cor(j,i) |
---|
1075 | 11017 continue |
---|
1076 | phizero(i)=uph(i) |
---|
1077 | cord(7,i) = 0. |
---|
1078 | gam(i) = sqrt(1. + cord(2,i)**2 + cord(4,i)**2 + cord(6,i)**2) |
---|
1079 | 11016 continue |
---|
1080 | c---- |
---|
1081 | c---- paw en test |
---|
1082 | cbm 18/07/2k call histogram1 |
---|
1083 | if (ip.eq.1 .or. ip.eq.2) |
---|
1084 | . write(nnout,11001) vv(2),ngood,tsig,sigdeg,tfac,rmax, |
---|
1085 | . ecath,decath,dwt,wtstep,nbasex,nbasey,nbasep,nbasew |
---|
1086 | 11001 format ( |
---|
1087 | . ' Generating ',f5.0,' particles from a photocathode',i5, |
---|
1088 | . ' really active particles '/ |
---|
1089 | . ' Quiet start with Hammersley sequence'/ |
---|
1090 | . ' sigma of laser pulse =',t40,1pg12.3,' psec'/ |
---|
1091 | . ' sigma of laser pulse =',t40,g12.3,' degrees'/ |
---|
1092 | . ' tmax of laser pulse = +- 2.5 sigma = ',t40,g12.3,' psec'/ |
---|
1093 | . ' rmax of photocathode =',t40,g12.3,' cm'/ |
---|
1094 | . ' mean electron energy at cathode =',t40,g12.3,' MeV'/ |
---|
1095 | . ' sigma of energy spread at cathode =',t40,g12.3,' MeV'/ |
---|
1096 | . ' integration step size =',t40,g12.3,' degrees'/ |
---|
1097 | . ' clock offset for earliest particle =',t40,g12.3,' degrees'/ |
---|
1098 | . ' bit reversal base for x =',i3 / |
---|
1099 | . ' bit reversal base for y =',i3 / |
---|
1100 | . ' bit reversal base for phi =',i3 / |
---|
1101 | . ' bit reversal base for W0 =',i3 ) |
---|
1102 | if(nn.ge.14) then |
---|
1103 | write(nnout,11025) spota/radian,spotd |
---|
1104 | 11025 format( |
---|
1105 | . ' angle of injection of laser light =',t40,g12.3,' degrees'/ |
---|
1106 | . ' distance focal point - cathode =',t40,g12.3,' cm') |
---|
1107 | endif |
---|
1108 | write(nnout,10018) |
---|
1109 | 11018 format (/' emittance summary for generated particles:') |
---|
1110 | iout7=1 |
---|
1111 | call outlaldp |
---|
1112 | iout7=iout7+1 |
---|
1113 | cbm call out6dp pour type 11 |
---|
1114 | call out6dp(0,0) |
---|
1115 | return ! input 11 |
---|
1116 | c |
---|
1117 | c--- type 12 EGUN |
---|
1118 | c |
---|
1119 | 12000 continue |
---|
1120 | open(unit=3,file='egun.dat',status='old',err=12001) |
---|
1121 | go to 12002 |
---|
1122 | 12001 continue |
---|
1123 | call appendparm |
---|
1124 | stop ' Abnormal stop file egun.dat not found ' |
---|
1125 | 12002 continue |
---|
1126 | nray=vv(3) |
---|
1127 | 12003 continue |
---|
1128 | read(3,*)i,regun(i),rdot(i),zdot(i),tdot(i),rcath(i),ioverr(i) |
---|
1129 | if(i.lt.nray)go to 12003 |
---|
1130 | close(3) |
---|
1131 | do 12004 i=1,nray |
---|
1132 | charge(i)=ioverr(i)*rcath(i)*twopi |
---|
1133 | if(i.ge.2)charge(i)=charge(i)+charge(i-1) |
---|
1134 | 12004 continue |
---|
1135 | do 12005 i=1,nray |
---|
1136 | charge(i)=charge(i)/charge(nray) |
---|
1137 | 12005 continue |
---|
1138 | r2=-1.-2./(npoints-2) |
---|
1139 | do 12006 np=2,npoints |
---|
1140 | r1=ranf() |
---|
1141 | r2=r2+2./(npoints-2) |
---|
1142 | r3=ranf()*twopi |
---|
1143 | do 12007 j=1,nray |
---|
1144 | if(r1.le.charge(j))go to 12008 |
---|
1145 | 12007 continue |
---|
1146 | c use the jth ray for this particle. determine bgx,bgy,bgz from |
---|
1147 | c rdot,zdot,tdot |
---|
1148 | 12008 continue |
---|
1149 | beta=sqrt(rdot(j)**2+zdot(j)**2+tdot(j)**2) |
---|
1150 | gamma=1./sqrt(1-beta**2) |
---|
1151 | bgz=gamma*zdot(j) |
---|
1152 | cor(6,np)=bgz |
---|
1153 | c determine a z position |
---|
1154 | cor(5,np)=z0+wavel/twopi*vv(4)*r2*radian*zdot(j) |
---|
1155 | c determin x and y |
---|
1156 | x=cos(r3)*regun(j)*vv(5) |
---|
1157 | y=sin(r3)*regun(j)*vv(5) |
---|
1158 | bgx=(rdot(j)*cos(r3)+tdot(j)*sin(r3))*gamma |
---|
1159 | bgy=(rdot(j)*sin(r3)-tdot(j)*cos(r3))*gamma |
---|
1160 | cor(1,np)=x+bgx/bgz*cor(5,np) |
---|
1161 | cor(2,np)=bgx |
---|
1162 | cor(3,np)=y+bgy/bgz*cor(5,np) |
---|
1163 | cor(4,np)=bgy |
---|
1164 | 12006 continue |
---|
1165 | return |
---|
1166 | c |
---|
1167 | c--- type 15 hot cathode microwave electron gun |
---|
1168 | c |
---|
1169 | 15000 continue |
---|
1170 | c megin for hot-cathode microwave electron gun (meg) |
---|
1171 | c H. Liu Aug. 29, 1989 IHEP, Academia Sinica |
---|
1172 | c modified again on oct 16 1989 |
---|
1173 | c vv(2) = no. of particles |
---|
1174 | c vv(3) = E0 ( in MV/m av axial e field, same as in cathode cell card) |
---|
1175 | c vv(4)= initial phase of reference particle (degres) |
---|
1176 | c vv(5) = SIGMAr of cathode (cm) |
---|
1177 | c vv(6) = radius of a thermionic cathode emitter (cm) |
---|
1178 | c vv(7) = temperature (Kelvin) |
---|
1179 | c vv(8) = SIGMAe for energy distribution (MeV) |
---|
1180 | c vv(9) = dwt in deg. |
---|
1181 | c vv(10)= work function (eV) |
---|
1182 | c vv(11)= storing flag |
---|
1183 | if(nn.lt.10) then |
---|
1184 | write(nnout,15001) |
---|
1185 | 15001 format('too few parameters for input type 15 and z0.gt.0.') |
---|
1186 | call appendparm |
---|
1187 | stop ' Abnormal stop input 15 ' |
---|
1188 | endif |
---|
1189 | c |
---|
1190 | e0 = vv(3) |
---|
1191 | rpphase=vv(4) |
---|
1192 | rsig=vv(5) |
---|
1193 | rmax=vv(6) |
---|
1194 | rmaxsq=rmax**2 |
---|
1195 | tempk= 0.001*vv(7) ! (in kilo-Kelvin) |
---|
1196 | decath=vv(8) |
---|
1197 | dwt=vv(9) |
---|
1198 | wf=vv(10) |
---|
1199 | ecath=w0 |
---|
1200 | c ----------------------------- |
---|
1201 | phiran=1. |
---|
1202 | rran=1. |
---|
1203 | if(nn.gt.11) then |
---|
1204 | phiran=vv(12) |
---|
1205 | numphi=vv(13) |
---|
1206 | rran=vv(14) |
---|
1207 | endif |
---|
1208 | write(nnout,*) 'phi random =',phiran,' r random =',rran |
---|
1209 | c ----------------------------- |
---|
1210 | cforj=0.44016*sqrt(e0)/tempk |
---|
1211 | fmax=exp(cforj) |
---|
1212 | c -------------------------------------- |
---|
1213 | if(phiran.eq.0.) then |
---|
1214 | phidiv=180./numphi |
---|
1215 | fsum=0.0 |
---|
1216 | do 15002 i=1,numphi |
---|
1217 | phiarg=(i-0.5)*phidiv |
---|
1218 | fphi=fmax**(sqrt(sin(phiarg*radian))) |
---|
1219 | fsum=fsum+fphi |
---|
1220 | 15002 continue |
---|
1221 | ithphi=1 |
---|
1222 | staphi=0. |
---|
1223 | do 15003 i=1,numphi |
---|
1224 | phinex=(i-0.5)*phidiv |
---|
1225 | fphi=fmax**(sqrt(sin(phinex*radian))) |
---|
1226 | nphidiv=npoints*fphi/fsum |
---|
1227 | subphi=phidiv/nphidiv |
---|
1228 | do 15004 j=1,nphidiv |
---|
1229 | ijk=ithphi+j |
---|
1230 | if(ijk.gt.npoints) go to 15006 |
---|
1231 | phinew=staphi+(j-1)*subphi |
---|
1232 | c print*,ijk,phinew |
---|
1233 | cor(5,ijk)=phinew |
---|
1234 | 15004 continue |
---|
1235 | ithphi=ijk |
---|
1236 | staphi=phinew+subphi |
---|
1237 | 15003 continue |
---|
1238 | c -------- |
---|
1239 | npwant=npoints-ijk |
---|
1240 | if(npwant.gt.0) then |
---|
1241 | do 15007 i=1,npwant |
---|
1242 | phinew=(i-0.5)*180./npwant |
---|
1243 | cor(5,ijk+i)=phinew |
---|
1244 | 15007 continue |
---|
1245 | endif |
---|
1246 | 15006 continue |
---|
1247 | endif |
---|
1248 | c -------------- |
---|
1249 | c The number of background data is calculated |
---|
1250 | c according to iphi*numr*(numr+1)/2. So when |
---|
1251 | c iphi=4 and numr=32, there are 2112 pairs of |
---|
1252 | c data produced. Usually npoints is less than |
---|
1253 | c 2000. If the last circle of particle |
---|
1254 | c distribution is required to be exactly |
---|
1255 | c close, set npoints equal to 2*numr*(numr+1). |
---|
1256 | if(rran.eq.0.) then |
---|
1257 | numra=32 |
---|
1258 | iphi=4 |
---|
1259 | isum0=1 |
---|
1260 | do 15008 i=1,numra |
---|
1261 | sqri=float(i) |
---|
1262 | numphi=iphi*i |
---|
1263 | deltaphi=360./numphi |
---|
1264 | do 15009 j=1,numphi |
---|
1265 | isum=isum0+j |
---|
1266 | phij=(j-1.)*deltaphi |
---|
1267 | phirad=phij*radian |
---|
1268 | cor(1,isum)=sqri*cos(phirad) |
---|
1269 | cor(3,isum)=sqri*sin(phirad) |
---|
1270 | if(isum.eq.npoints) go to 15010 |
---|
1271 | 15009 continue |
---|
1272 | isum0=isum |
---|
1273 | 15008 continue |
---|
1274 | 15010 continue |
---|
1275 | factor=rmax/sqri |
---|
1276 | do 15011 i=2,npoints |
---|
1277 | cor(1,i)=factor*cor(1,i) |
---|
1278 | cor(3,i)=factor*cor(3,i) |
---|
1279 | 15011 continue |
---|
1280 | endif |
---|
1281 | c -------------- |
---|
1282 | c -------------------------------------- |
---|
1283 | acj0=120.4*tempk**2*1.e+6*exp(-wf/(8.62*1.e-5*tempk*1000.)) |
---|
1284 | acjmax=acj0*fmax |
---|
1285 | c aci=acj*pi*rmaxsq |
---|
1286 | c ---------------------------------- |
---|
1287 | dppi=pi |
---|
1288 | qint=GPINDP(0.d0,dppi,1.d-6,dpepsout,fsch,1) |
---|
1289 | QCHARGE=0.5*rmaxsq*acj0*qint*1000./freq |
---|
1290 | c -------------------------------------- |
---|
1291 | nstep=int(rpphase/dwt) |
---|
1292 | if(rpphase.ne.dwt*float(nstep)) nstep=nstep+1 |
---|
1293 | wtstep=nstep*dwt |
---|
1294 | c |
---|
1295 | esc=0. |
---|
1296 | escrf=esc/e0 |
---|
1297 | if(ip.eq.1 .or. ip.eq.2) write(nnout,15012) (vv(i),i=1,11), |
---|
1298 | . ecath,acj0,acjmax,fmax,qcharge,qcharge*1.e+10/1.6,nstep,wtstep |
---|
1299 | 15012 format( |
---|
1300 | . ' intype = ',t40,f12.3,' '/ |
---|
1301 | . ' generating ',f6.0,' particles from a thermionic cathode'/ |
---|
1302 | . ' average axial e field =',t40,g12.3,' MV/m'/ |
---|
1303 | . ' initial phase of r.p. =',t40,g12.3,'degrees'/ |
---|
1304 | . ' SIGMAr for emitter =',t40,g12.3,' cm'/ |
---|
1305 | . ' radius of emitter =',t40,g12.3,'cm'/ |
---|
1306 | . ' operating temperature =',t40,g12.3, 'kilo-Kelvin'/ |
---|
1307 | . ' SIGMAe of energy distribution =',t40,g12.3,'MeV'/ |
---|
1308 | . ' rf step size dwt =',t40,g12.3,' degrees'/ |
---|
1309 | . ' work function of emitter =',t40,g12.3,'eV'/ |
---|
1310 | . ' storing flag =',t40,g12.3,' '/ |
---|
1311 | . ' mean kinetic energy W0 =',t40,g12.3,' MeV'/ |
---|
1312 | . ' current density J0 =',t40,g12.3,' A/cm**2'/ |
---|
1313 | . ' current density Jmax =',t40,g12.3, ' A/cm**2'/ |
---|
1314 | . ' Fmax (Jmax/J0) =',t40,g12.3,' '/ |
---|
1315 | . ' total charge amount in the bunch =',t40,g12.3,' nC'/ |
---|
1316 | . ' number of electrons in the bunch =',t40,g12.3, ' '/ |
---|
1317 | . ' nstep= ',t40,i6,' '/ |
---|
1318 | . ' wtstep =',t40,g12.3,' degrees'/) |
---|
1319 | c |
---|
1320 | zcath(1)=z0 |
---|
1321 | c |
---|
1322 | bgr=cor(6,1) |
---|
1323 | bzr=bgr/sqrt(1.+bgr**2) |
---|
1324 | cor(5,1)=zcath(1)-wtstep*wavel*bzr/360. |
---|
1325 | c ----------------------- |
---|
1326 | c-------------------------------------------------------------------------- |
---|
1327 | uphmin=0. |
---|
1328 | uphmax=0. |
---|
1329 | ux0min=0. |
---|
1330 | ux0max=0. |
---|
1331 | vy0min=0. |
---|
1332 | vy0max=0. |
---|
1333 | c-------------------------------------------------------------------------- |
---|
1334 | do 15013 i=np1,npoints |
---|
1335 | if(phiran.ne.0.) then |
---|
1336 | c choose initial phases |
---|
1337 | 15014 ph1=pi*ranf() |
---|
1338 | se=sin(ph1)+escrf |
---|
1339 | if (se.le.0.) go to 15014 |
---|
1340 | fph1=fmax**sqrt(se) |
---|
1341 | ph2=rn32() |
---|
1342 | fph2=fmax*ph2 |
---|
1343 | if(fph2.gt.fph1) go to 15014 |
---|
1344 | phase=ph1/radian |
---|
1345 | else |
---|
1346 | phase=cor(5,i) |
---|
1347 | endif |
---|
1348 | c---- |
---|
1349 | uph(i)=phase |
---|
1350 | if(uph(i).gt.uphmax)uphmax=uph(i) |
---|
1351 | if(uph(i).lt.uphmin)uphmin=uph(i) |
---|
1352 | c---- |
---|
1353 | c choose initial energy at cathode |
---|
1354 | 15015 call norran(u) |
---|
1355 | e=ecath+u*decath |
---|
1356 | if(e.le.0.) go to 15015 |
---|
1357 | betcat=sqrt(2.*e/erest) |
---|
1358 | gamma=1./sqrt(1.-betcat**2) |
---|
1359 | c choose initial momentum,isotropic, but betaz>0 |
---|
1360 | sand1=ranf() |
---|
1361 | theta=(sand1-0.5)*pi |
---|
1362 | cose=cos(theta) |
---|
1363 | sine=sqrt(1.-cose**2) |
---|
1364 | phi=twopi*sandom(d) |
---|
1365 | bx=betcat*sine*cos(phi) |
---|
1366 | by=betcat*sine*sin(phi) |
---|
1367 | bz=betcat*cose |
---|
1368 | c ------------------------------------ |
---|
1369 | cor(2,i)=gamma*bx |
---|
1370 | cor(4,i)=gamma*by |
---|
1371 | cor(6,i)=gamma*bz |
---|
1372 | c choose (x,y,z) |
---|
1373 | if(rran.ne.0.) then |
---|
1374 | 15016 call rannor(u,v) |
---|
1375 | x0=u*rsig |
---|
1376 | y0=v*rsig |
---|
1377 | rsq=x0**2+y0**2 |
---|
1378 | c if(rsq.gt.rmaxsq .or. rsq.lt.rinsq) go to 15016 |
---|
1379 | if(rsq.gt.rmaxsq) go to 15016 |
---|
1380 | else |
---|
1381 | x0=cor(1,i) |
---|
1382 | y0=cor(3,i) |
---|
1383 | endif |
---|
1384 | c---- |
---|
1385 | ux0(i)=x0 |
---|
1386 | if(ux0(i).gt.ux0max)ux0max=ux0(i) |
---|
1387 | if(ux0(i).lt.ux0min)ux0min=ux0(i) |
---|
1388 | vy0(i)=y0 |
---|
1389 | if(vy0(i).gt.vy0max)vy0max=vy0(i) |
---|
1390 | if(vy0(i).lt.vy0min)vy0min=vy0(i) |
---|
1391 | c---- |
---|
1392 | c |
---|
1393 | c shift x,y,z |
---|
1394 | zcath(i)=zcath(1) |
---|
1395 | cliu toff=(wtstep+phase)*wavel/360. |
---|
1396 | toff=phase*wavel/360. |
---|
1397 | cor(1,i)=x0-toff*bx |
---|
1398 | cor(3,i)=y0-toff*by |
---|
1399 | cor(5,i)=zcath(i)-toff*bz |
---|
1400 | c -------------------- |
---|
1401 | phizero(i)=phase |
---|
1402 | c ----------------------------------- |
---|
1403 | 15013 continue |
---|
1404 | c-------------------------------------------------------------------------- |
---|
1405 | if(ip.eq.2) then |
---|
1406 | open(unit=ninput,file='input15',status='unknown') |
---|
1407 | write(ninput,*) npoints,sigdeg,rmax |
---|
1408 | do 15017 i=1,npoints |
---|
1409 | write(ninput,*) uph(i),ux0(i),vy0(i) |
---|
1410 | 15017 continue |
---|
1411 | endif |
---|
1412 | c--- paw en test |
---|
1413 | cbm 18/07/2k call histogram1 |
---|
1414 | c-------------------------------------------------------------------------- |
---|
1415 | c ---------------------------- |
---|
1416 | c for emittance summary, copy parameters into cord & gam |
---|
1417 | c note that officially this is done at line 210 of the main program |
---|
1418 | ngood = npoints |
---|
1419 | c --------------- |
---|
1420 | do 15018 i = 1,npoints |
---|
1421 | do 15019 j = 1,6 |
---|
1422 | cord(j,i) = cor(j,i) |
---|
1423 | 15019 continue |
---|
1424 | cord(7,i) = 0. |
---|
1425 | gam(i) = sqrt(1. + cord(2,i)**2 + cord(4,i)**2 + cord(6,i)**2) |
---|
1426 | 15018 continue |
---|
1427 | if (ip.ne.0) then |
---|
1428 | write(nav,*) ' ' |
---|
1429 | write (nav,15020) |
---|
1430 | 15020 format (/' emittance summary for generated particles:') |
---|
1431 | cbm call out6dp pour type 15 |
---|
1432 | call out6dp(0,0) |
---|
1433 | endif |
---|
1434 | return ! return input 15 |
---|
1435 | end |
---|
1436 | c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* |
---|