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

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

parmela pspa initial

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