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

Last change on this file since 315 was 315, checked in by lemeur, 11 years ago

fichier parin.input0

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