source: PSPA/parmelaPSPA/trunk/input.f @ 18

Last change on this file since 18 was 18, checked in by lemeur, 12 years ago

corr. bug dans distrir.f

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