subroutine input(lonnom, work_path) c-------last modification date 03/27/95 c generalized input routine c type 0, read np particles from file 'parin' and store c in cor c type 1, np particles are generated randomly in three c phase plane ellipses c type 2, np particles are generated randomly in a six c dimensional ellipse c type 3, np particles are generated uniformily on the c perimiter of one of the phase plane ellipses c type 4, particle number np is positioned at six specified c cordinates c type 5, np particles are generated randomly in a four c dimensional transverse hyperspace with random c phase and energy spread within an ellipse. c type 6, np particles are generated randomly in a four c dimensional transverse hyperspace with uniform c phase and random energy spread. c type 7, np particles are generated in a kapchinskij- c vladimerskij distribution transversely with c random phase and energy spread within an ellipse. c type -n,random number generator reset to beginning c type 8, np particles are gererated randomly with uniform c distribution in real space (x,y,z), then xp, yp, c and zp are chosen from within each phase-plane c ellipse. .z,zp are then converted to phi,w. c type 9, gausian distribution in x, y and z. xp, yp, and energy c type 4 spread are zero. parameters are ntype,number of c particles,sigmar,rmax,sigmaz(deg),zmax(deg). c this used to be the following. c rectangular array in one of the three phase planes. c ntype,kind,hl,hr,nh,vb,vt,nv c type 10 np particles are generated at a photocathode, c see comments after line 10000 c type 11 same as type 10 uses Hammersley sequence c see comments after line 11000 c type 12 np particles are generated from output data from c EGUN data. parameters are ntype, np, number of rays, c phase spread, egun mesh size in cm. The file egun.dat c contains the data. The format of the data is: c ray#,r,rdot,zdot,tdot,r(cathode),i/r(cathode)/2pi. c type 14 same as type 4 with cord(7,i)=0. for i=1,np c type 15 hot cathode microwave electron gun c type 19 rectangular array c type 100*n, as type n, but xp and yp are adjusted to c force the angular momentum to be zero. c valid for types 1, 2, 5, and 6. c vv(1)=type c vv(2)=no of particles or particle no c vv(3-8)=transverse ellipse parameters c vv(9+10)=phase and energy spread (half) c vv(11+16)=displacements of input beam c phase spread or delta phase is entered in degrees c c------------------------------------------------------------------------ c------------------------------------------------------------------------ save c include 'param_sz.h' include 'param_paw.h' include 'constcom.h' include 'coordcom.h' include 'misccom.h' include 'ncordscom.h' include 'pcordcom.h' include 'upawcom.h' include 'ucom.h' c character*256 nomfic character*256 work_path character*256 nomf integer lonnom c parameter (igun=31) common/cf/fmax common/cout7/corout7(6,imaa),bgs,nelal,nbufl,iout7 common/intype/intype ! mcd common/jones/ajones,zjones,zcath(imaa) common/out6/nbphase,nbremesh common/pawc/h(nxpaw) common/part/nparticle(imaa) common/phase/uphmin,uphmax,uph(imaa),sigdeg,ipas common/ux0/ux0min,ux0max,ux0(imaa) common/vx0/vy0min,vy0max,vy0(imaa) common/dist/distp(imaa),xr(imaa),yr(imaa) common/wtstep/wtstep,nstep common/hb2/rmax common/hom/xhom(imaa),yhom(imaa) dimension xquiet(imaa),yquiet(imaa),anglex(imaa) dimension x1(imaa),y1(imaa) dimension p1(imaa),p2(imaa),p3(imaa),p4(imaa) dimension iwork1(imaa) dimension regun(igun),rdot(igun),zdot(igun),tdot(igun), * rcath(igun),ioverr(igun),charge(igun) double precision dppi,dpepsout real f(41) external fsch c-------------------------------------------------------------------------- c* data rnumb/0.0/ data f/0.,.0398,.0793,.1179,.1554,.1915,.2258,.258,.2881,.3159, * .3413,.3643,.3849,.4032,.4192,.4332,.4452,.4554,.4641,.4713, * .4772,.4821,.4861,.4893,.4918,.4938,.4953,.4965,.4974,.4981, * .4987,.4990,.4993,.4995,.4997,.49975,.4998,.49985,.4999,.49995, * .5/ print *, ' PARMELA ENTREE INPUT ' do 1 i=1,imaa ! mcd zcath(i)=0. ! mcd 1 continue ! mcd r1=1.0 r2=1.0 r3=1.0 r4=1.0 r5=1.0 r6=1.0 phis=0. es=w0 c fetch first random number c if (rnumb.eq.0.0) dum=sandsv (rnumb) gg=es/erest bgs=sqrt(gg*(2.+gg)) bs=sqrt(gg*(2.+gg))/(1.+gg) ntype=abs(vv(1)) intype=ntype if(ip.eq.1 .or. ip.eq.2)write(nnout,2) ntype 2 format (' input type ',t40,i8) c c---- type 0 read np particle c if(ntype.eq.0)then nomf = nomfic(work_path(1:lonnom),'parin.input0') print *, ' parmela file : fichier parin : ', nomf open(unit=3,file=nomf,status='old',err=10) npoints=vv(2)+1 cor(5,1)=0. cor(6,1)=0. tcharge=0. do 11 j=2,npoints read(3,*,err=12,end=12) cor(1,j),cor(2,j),cor(3,j), * cor(4,j),cor(5,j),cor(6,j),weight(j) weight(j)=abs(weight(j)) tcharge=tcharge+weight(j) npart=j cor(5,1)=cor(5,1)+cor(5,j) cor(6,1)=cor(6,1)+cor(6,j) write(nnout,*) j,cor(1,j),cor(2,j),cor(3,j),cor(4,j), * cor(5,j),cor(6,j),weight(j) ! 02/92 11 continue 12 close(3) cor(5,1)=cor(5,1)/(npart-1) cor(6,1)=cor(6,1)/(npart-1) cor(1,1)=0. cor(2,1)=0. cor(3,1)=0. cor(4,1)=0. weight(1)=0. npoints=npart ngood=npoints write(nnout,*)' total charge= ',tcharge,' coul' return 10 continue call appendparm stop ' Abnormal stop parin not found' endif c if (ntype.gt.100) ntype=ntype/100 c np1=npoints+1 c if(ntype.eq.4 .or. ntype.eq.14)go to 400 c c if(ntype.eq.9)go to 900 c if(ntype.eq.19)go to 19000 c npoints=vv(2)+npoints if(npoints.gt.imaa)npoints=imaa c if(ntype.eq.3)go to 300 c if(ntype.eq.10)go to 10000 c if(ntype.eq.11)go to 11000 c if(ntype.eq.12)go to 12000 c if(ntype.eq.15)go to 15000 c if(ntype.eq.5)then rtest=1.2/float(npoints)**.2 rtest2=rtest**2 itest=0 endif ax=vv(3) bx=vv(4) ex=vv(5) ay=vv(6) by=vv(7) ey=vv(8) if(ntype.ne.8)go to 15 c c---- type 8. get ellipse parameters for z, zp. c az=vv(9) bz=vv(10) ez=vv(11) n1=12 gz=(1.+az**2)/bz zmx=sqrt(ez/gz) zpmx=ez/zmx dz=-az/gz go to 16 15 continue n1=11 dphi=vv(9)*radian de=vv(10) c calculate ellipse parameters 16 continue gx=(1.0+ax**2)/bx gy=(1.0+ay**2)/by xmx=sqrt(ex/gx) xpmx=sqrt(ex*gx) ymx=sqrt(ey/gy) dy=-ay/gy ypmx=sqrt(ey*gy) dx=-ax/gx if (vv(1).lt.0.0) call sandst (rnumb) c c types 1,2,5,6,and 7 c xnp=vv(2)-1. if(xnp.lt.1.)xnp=1. area=pi/vv(2) 198 continue dr5=sqrt((1.-sqrt(1-area**2))/2.) itest=0 istart=2 199 continue do 60 np=np1,npoints if(ntype.ne.8)go to 19 200 continue r1=2.*ranf()-1. r3=2.*ranf()-1. r5=2.*ranf()-1. rr=r1**2+r3**2+r5**2 if (rr.gt.1.)go to 200 201 continue r2=2.*ranf()-1. r4=2.*ranf()-1. r6=2.*ranf()-1. if(r2**2+r4**2+r6**2.gt.1.)go to 201 sq=sqrt(1.-rr) r2=r2*sq r4=r4*sq r6=r6*sq go to 50 19 continue if(ntype.eq.7)go to 45 20 continue r1=2.0*ranf()-1.0 r2=2.0*ranf()-1.0 if ((r1**2+r2**2).gt.1.0) go to 20 30 continue r3=2.0*ranf()-1.0 r4=2.0*ranf()-1.0 if ((r3**2+r4**2).gt.1.0) go to 30 if (ntype.lt.5) go to 40 if (r1**2+r2**2+r3**2+r4**2.gt.1) go to 20 40 continue if (ntype.eq.6) r5=2.0*float(np-np1)/xnp-1.0 if (ntype.eq.5) r6=sin(abs(acos(r5)))*(ranf()*2.-1) c--- this code is to try to make the distribution in real space more c--- uniform for ntype 5. if(ntype.eq.5)then itest=itest+1 if(itest.gt.10)then write(ndiag,*)'rtest too big' rtest=.9*rtest rtest2=rtest**2 r5=1.0 write(nnout,*)'rtest now=',rtest go to 198 endif do 53 i=istart,np-1 dxx=cord(1,i)-r1 dpx=cord(2,i)-r1 dyy=cord(3,i)-r3 dpy=cord(4,i)-r4 dzz=cord(5,i)-r5 if(dxx**2+dpx**2+dyy**2+dpy**2+dzz**2.lt.rtest2)go to 20 if(dzz.gt.rtest)istart=i+1 53 continue itest=0 cord(1,np)=r1 cord(2,np)=r2 cord(3,np)=r3 cord(4,np)=r4 cord(5,np)=r5 endif if (ntype.eq.6) r6=2.0*ranf()-1.0 if (ntype.ge.6) go to 50 if (ntype.ne.2) go to 50 if (r1**2+r2**2+r3**2+r4**2+r5**2+r6**2.gt.1.0) go to 20 go to 50 c c type 7 k-v distribution c 45 continue a=2.0*ranf()-1.0 a=acos(a)/2.0 b=twopi*(ranf()-0.5) c=twopi*(ranf()-0.5) r1=cos(a)*cos(b) r3=cos(a)*sin(b) r2=sin(a)*cos(c) r4=sin(a)*sin(c) c makes phase and energy distribution into an ellipse c 46 r5=2.*ranf()-1. c r6=2.*ranf()-1. c if((r5*r5+r6*r6).gt.1.) go to 46 46 continue r6=sin(abs(acos(r5)))*(ranf()*2.-1.) 50 continue if (vv(1).le.10.) go to 52 if (ntype.gt.6) go to 52 if (ntype.eq.3 .or. ntype.eq.4 .or. ntype.eq.14) go to 52 c---force zero transverse angular momentum. rr=sqrt(r1**2+r3**2) rp=sqrt(r2**2+r4**2)*sign(1.,r2) r2=rp*r1/rr r4=rp*r3/rr 52 continue cor(1,np)=r1*xmx+dx*r2*xpmx cor(2,np)=r2*xpmx cor(3,np)=r3*ymx+dy*r4*ypmx cor(4,np)=r4*ypmx if(ntype.ne.8)go to 55 zp=r6*zpmx z=r5*zmx+dz*zp cor(5,np)=-z*twopi/(bs*wavel)+phis cor(6,np)=2.*es*zp+es go to 60 55 continue cor(5,np)=r5*dphi+phis cor(6,np)=r6*de+es if(ntype.ne.5.and.ntype.ne.7)go to 60 r5=r5-dr5 if(r5.le.-1)go to 65 sac=sin(acos(r5)) dr5=.5*area*sac if((r5-dr5).le.-1.)go to 60 dr5=area/(sac+sin(acos(r5-dr5))) 60 continue go to 67 65 continue npoints=np 67 continue if (nn.ge.n1) go to 99998 go to 99999 c c type 3 c 300 continue do 301 n=np1,npoints do 302 m=1,4 cor(m,n)=0.0 302 continue cor(5,n)=phis cor(6,n)=es 301 continue npp=vv(2) j=vv(3) j=2*j gg=(1.+vv(4)**2)/vv(5) dd=-vv(4)/gg xyzmx=sqrt(vv(6)/gg) xyzpmx=vv(6)/xyzmx da=twopi/npp kk=np1-1 do 303 np=1,npp ang=float(np-1)*da xyzp=xyzpmx*sin(ang) xyz=xyzmx*cos(ang)+xyzp*dd if(j.eq.6)xyz=xyz*radian kk=kk+1 cor(j-1,kk)=cor(j-1,kk)+xyz cor(j,kk)=xyzp+cor(j,kk) 303 continue n1=7 if (nn.gt.n1) go to 99998 go to 99999 c c type 4 and 14 c 400 continue c input ntype x xp y yp phi w ksi1 ksi2 ksi3 n4 c the first 7 parameters must be on all cards input c n4 is the total number of input cards c the first card must have 11 parameters c corrige le 10/12/2008 (nn.eq.8)-->(nn.eq.11) if(nn.eq.11) then maxc=vv(11) ncarte=1 else ncarte=ncarte+1 endif npoints=np1 if(npoints.gt.imaa)npoints=imaa do 401 n=1,4 cor(n,npoints)=vv(n+1) 401 continue cor(5,npoints)=vv(6)*radian+phis cor(6,npoints)=vv(7)+es phizero(npoints)=vv(6) go to 99999 c c---- type 9. gausian array. c 900 continue np1=npoints+1 npoints=npoints+vv(2) if(npoints.gt.imaa)npoints=imaa sigmar=vv(3) rmax=vv(4) sigmaz=vv(5) zmax=vv(6) ss=-1. rnorm=1.41421356/sigmaz zzmax=rnorm*zmax if(zzmax.gt.4.)zzmax=3.999999999 iz=zzmax*10. a=zzmax*10.-iz fmax=f(iz+1)*(1.-a)+f(iz+2)*a do 901 i=np1,npoints if(sigmar.gt.10.*rmax)then r1=sqrt(ranf())*rmax else 902 continue r1=sigmar*sqrt(-log(ranf())) if(r1.gt.rmax)go to 902 endif theata=twopi*ranf() cor(1,i)=sin(theata)*r1 cor(3,i)=cos(theata)*r1 if(sigmaz.gt.10.*zmax)then cor(5,i)=2.*radian*zmax*(float(i-np1)/float(npoints-np1)-.5) else rand=fmax*float(i-np1)/float(npoints-np1) do 903 j=2,41 if(rand.le.f(j))go to 904 903 continue 904 continue r2=zmax/zzmax*(float(j-2)+(rand-f(j-1))/(f(j)-f(j-1)))*.1 ss=-ss cor(5,i)=r2*radian*ss endif cor(2,i)=0. cor(4,i)=0. cor(6,i)=es 901 continue go to 99999 c c type 19. rectangular array. c 19000 kind=vv(2) hl=vv(3) hr=vv(4) if(kind.eq.3)hl=hl*radian if(kind.eq.3)hr=hr*radian nh=vv(5) vb=vv(6) vt=vv(7) nv=vv(8) k1=2*kind-1 k2=k1+1 npoints=np1+nh*nv-1 if (npoints.gt.imaa) npoints=imaa do 19001 np=np1,npoints cor(5,np)=phis cor(6,np)=es do 19001 n=1,4 cor(n,np)=0. 19001 continue hz=hl if (k1.eq.5) hz=hz+phis vz=vb if (k2.eq.6) vz=vz+es dh=0. if (nh.gt.1) dh=(hr-hl)/(nh-1) dv=0. if (nv.gt.1) dv=(vt-vb)/(nv-1) np=np1-1 hh=hz-dh do 19002 i=1,nh hh=hh+dh v=vz-dv do 19002 j=1,nv np=np+1 if (np.gt.imaa) go to 19003 v=v+dv cor(k1,np)=hh cor(k2,np)=v 19002 continue 19003 continue n1=9 if(nn.lt.n1)go to 99999 c c displace input beam c 99998 continue n2=n1+4 if(nn.ge.n2)vv(n2)=vv(n2)*radian do 99997 n=n1,nn k=n-n1+1 do 99997 np=np1,npoints cor(k,np)=cor(k,np)+vv(n) 99997 continue go to 99999 c---convert from x,xp,y,yp,phi,w to x,bgx,y,bgy,z,bgz 99999 continue ngood=npoints nbufl=npoints c ii=2,npoints reference particle is treated in main program if (ntype.ne.4 .or. ntype.eq.14) then do 99990 ii=2,npoints phizero(ii)=cor(5,ii) do 99991 iii=1,6 corout7(iii,ii)=cor(iii,ii) 99991 continue 99990 continue nelal=0 iout7=iout7+1 call outlaldp iout7=iout7+1 endif g=w0/erest zcon=-wavel/twopi do 99992 np=np1,npoints phi=cor(5,np)/radian phizero(np)=phi wz=cor(6,np) g=wz/erest bgz=sqrt(g*(2.+g)) dz=zcon*cor(5,np)*bgz/(1.+g) bgx=cor(2,np)*bgz bgy=cor(4,np)*bgz cor(1,np)=cor(1,np)+cor(2,np)*(z0+dz) cor(2,np)=bgx cor(3,np)=cor(3,np)+cor(4,np)*(z0+dz) cor(4,np)=bgy cor(5,np)=z0+dz cor(6,np)=sqrt(bgz**2-bgx**2-bgy**2) nparticle(np)=np 99992 continue c for emittance summary, copy parameters into cord & gam c note that officially this is done at line 210 of the main program do 99993 i = np1,npoints cord(7,i)=0. do 99994 j = 1,6 cord(j,i) = cor(j,i) 99994 continue if(ntype.eq.14) cord(7,i)=0.0 gam(i) = sqrt(1. + cord(2,i)**2 + cord(4,i)**2 + cord(6,i)**2) 99993 continue write(nnout,*) ' ncarte ',ncarte,' maxc ',maxc if((ntype.eq.4 .or. ntype.eq.14).and.(ncarte.lt.maxc)) return if(ip.ne.0) then 8337 continue cbm call ou6dp normal call out6dp(0,0) endif return ! return input 1 a 9 c c---type 10 = particles generated at a photocathode c 10000 continue c last modification 03/08/90 c vv(2) = no. of particles c vv(3) = sigma of laser pulse duration in picoseconds c vv(4) = half width of laser pulse. i.e. generate a gaussian c time profile with sigma = vv(3), but truncate to +-vv(4) c vv(5) = sigma of laser spot in cm - assumes gaussian radial profile c vv(6) = max radius of photocathode. i.e. truncate radial dist at vv(6) c vv(7) = mean K.E. of electrons leaving the photocathode, in MeV c This should be the same as paramter W0 on the RUN card. c vv(8) = sigma of energy-spread distribution at the photocathode (MeV). c vv(9) = rf step size in deg. This must be the same as parmater DWT c on the START card c vv(10) = type of random (new) 1 by rannor c 2 x,y by rannor phi by didtribution c vv(11) = spota -- Angle of laser beam with the normal at the cathode c (degrees), if non-zero you need vv(12) c vv(12) = spotd -- Distance from focal point to cathode (cm) c if no zero you need vv(11) c vv(13) = number of chenals for histo (new) c vv(13) nonzero you need vv(11) and vv(13) egal or not zero c vv(14) = mu -- the cathode shape parameter if Mike Jones. Must have c 0 < mu < .5. If mu is nonzero, the field shape of the cell c with element no. 1 will be taken from Jones' analytic form, c overriding the fourier coefficients, if any c vv(15) = length of cell with cathode as front surface (cm). Must be c the same a L on the first CELL card c c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ if (nn.lt.10) then write(nnout,10001)nn 10001 format (' Too few parameters for input type 10 - nn = ',i2, 1 ' parameters only 10 are needed') call appendparm stop ' Abnormal stop input 10 a ' endif if (nn.lt.11.and.ip.eq.2) then write(nnout,10002)ip,nn 10002 format (' Too few parameters for input type 10 - ip = ', 1 i1,' nn = ',i2,' parameters only',/,' 11 are needed ') call appendparm stop ' Abnormal stop input 10 b ' endif tsig = vv(3) tmax = vv(4) if (tsig.le.0.) then write(nnout,10003) 10003 format (' sigma of laser pulse must be > 0') call appendparm stop ' Abnormal stop input 10 c ' endif tfac = tmax/tsig rsig = vv(5) rmax = vv(6) rmaxsq = rmax**2 ecath = vv(7) decath = vv(8) c convert laser pulse sigma to degrees sigdeg = tsig*freq*360./1.e6 c generate particles during +- tmax in time c set the starting clock back by an integer number of time steps, c corresponding to at least tmax degmax = tmax*freq*360./1.e6 dwt = vv(9) nstep = int(degmax/dwt) if (degmax.ne.dwt*float(nstep)) nstep = nstep + 1 wtstep = nstep*dwt iran=vv(10) if (ip.eq.1 .or. ip.eq.2) . write(nnout,10004) vv(2),tsig,sigdeg,tmax,rsig,rmax,ecath, . decath,dwt,wtstep 10004 format ( . ' Generating ',f10.0,' particles from a photocathode'/ . ' sigma of laser pulse =',t40,1pg12.3,' psec'/ . ' sigma of laser pulse =',t40,g12.3,' degrees'/ . ' tmax of laser pulse =',t40,g12.3,' psec'/ . ' sigma of laser spot =',t40,g12.3,' cm'/ . ' rmax of photocathode =',t40,g12.3,' cm'/ . ' mean electron energy at cathode =',t40,g12.3,' MeV'/ . ' sigma of energy spread at cathode =',t40,g12.3,' MeV'/ . ' integration step size =',t40,g12.3,' degrees'/ . ' clock offset for earliest particle =',t40,g12.3,' degrees') if (vv(10).eq.1) write(nnout,10019) iran 10019 format (' generated by rannor only ran = ',t45,i2) if (vv(10).eq.2) write(nnout,10020) iran 10020 format (' generated by distriphase-rannor ran = ',t45,i2) if (vv(10).gt.2) then write(nnout,*) ' Input 10 ran.gt.2 ran = ',iran call appendparm stop ' Abnormal stop input 10 d ' endif c---- laser beam angle if(nn.ge.11) then spota=vv(11)*radian spotd=vv(12) if (ip.eq.1 .or. ip.eq.2) then write(nnout,10025) spota/radian,spotd 10025 format( . ' angle of injection of laser light =',t40,g12.3,' degrees'/ . ' distance focal point - cathode =',t40,g12.3,' cm') endif endif c---- if (nn.ge.14) then ajones = vv(14) if (ajones.lt.0.0 .or. ajones.gt.0.499) then write(nnout,10005) ajones 10005 format (' cathode shape parameter = ',1pg12.3, . ' is out of range') call appendparm stop ' Abnormal stop input 10 e ' endif zjones = 0. if (ip.eq.1 .or. ip.eq.2 .or. ip.eq.4) . write(nnout,10006) ajones 10006 format (' cathode shape parameter =',t40,1pg12.3) endif if (nn.ge.15) then zjones = vv(15) if (ip.eq.1 .or. ip.eq.2 .or. ip.eq.4) . write(nnout,10007) zjones 10007 format (' gun cavity length =',t40,1pg12.3,' cm') endif c----- c c offset z of reference particle so it reaches z = 0 at center of pulse bgr = cor(6,1) bzr = bgr/sqrt(1. + bgr**2) cor(5,1) = -wtstep*wavel*bzr/360. c----- uphmin=0. uphmax=0. ux0min=0. ux0max=0. vy0min=0. vy0max=0. c----------------------------- ran = 1 ---- if(vv(10).eq.1) then do 10008 i = np1,npoints c choose time (phase) of generation, relative to center of the pulse 10009 continue call rannor(u,v) if (abs(u).gt.tfac) go to 10009 phase = u*sigdeg c choose intial energy at cathode 10010 continue call rannor(u,v) e = ecath + u*decath if (e.le.0.) go to 10010 betcat = sqrt(2.*e/erest) gamma = 1./sqrt(1. - betcat**2) c choose initial momentum, isotropic, but betaz > 0 cose = sandom(d) sine = sqrt(1. - cose**2) phi = twopi*sandom(d) bx = betcat*sine*cos(phi) by = betcat*sine*sin(phi) bz = betcat*cose cor(2,i) = gamma*bx cor(4,i) = gamma*by cor(6,i) = gamma*bz c choose (x,y,z) 10011 continue call rannor(u,v) x0 = u*rsig y0 = v*rsig rsq = x0**2 + y0**2 if (rsq.gt.rmaxsq) go to 10011 c---------------------------------- if((spota.eq.0.).or.(spotd.eq.0)) go to 10021 anglex10 = atan(x0/spotd) delta = spotd*(1-(cos(spota)/cos(spota+(anglex10)))) delta = delta*360./wavel x0= spotd*sin(anglex10)/cos(anglex10+spota) phase=phase+delta 10021 continue c--------------------------------------- uph(i)=phase if(uph(i).gt.uphmax)uphmax=uph(i) if(uph(i).lt.uphmin)uphmin=uph(i) c--------------------------------------- ux0(i)=x0 if(ux0(i).gt.ux0max)ux0max=ux0(i) if(ux0(i).lt.ux0min)ux0min=ux0(i) vy0(i)=y0 if(vy0(i).gt.vy0max)vy0max=vy0(i) if(vy0(i).lt.vy0min)vy0min=vy0(i) c----- z0 = 0. if (ajones.gt.0.) then c calculate the shape of the cathode surface z0 = cathod(rsq) zcath(i) = z0 endif c----- c shift initial (x,y,z) so that the particle will emerge from c the photocathode at (x0,y0,z0) at time = phase toff = (wtstep + phase)*wavel/360. cor(1,i) = x0 - toff*bx cor(3,i) = y0 - toff*by cor(5,i) = z0 - toff*bz 10008 continue else ! vv(10).ne.1 c---- call distriphase(degmax,sigdeg,npoints-1) call distrir(rmax,rsig,npoints-1) c---- do 10012 i = np1,npoints c choose time (phase) of generation, relative to center of the pulse c---- u=distp(i-1) c---- phase=u c---- c choose intial energy at cathode 10013 continue call rannor(u,v) c---- e = ecath + u*decath if (e.le.0.) go to 10013 betcat = sqrt(2.*e/erest) gamma = 1./sqrt(1. - betcat**2) c choose initial momentum, isotropic, but betaz > 0 cose = sandom(d) sine = sqrt(1. - cose**2) phi = twopi*sandom(d) bx = betcat*sine*cos(phi) by = betcat*sine*sin(phi) bz = betcat*cose cor(2,i) = gamma*bx cor(4,i) = gamma*by cor(6,i) = gamma*bz c choose (x,y,z) c---- c---- if(vv(10).eq.2) then 10014 continue call rannor(u,v) x0=u*rsig y0=v*rsig rsq=x0**2+y0**2 if(rsq.gt.rmaxsq) go to 10014 else ! vv(10) .ne.2 x0 = xr(i-1) y0 = yr(i-1) endif ! vv(10) .eq. 2 c---------------------------------- if((spota.eq.0.).or.(spotd.eq.0)) go to 10022 anglex10 = atan(x0/spotd) delta = spotd*(1-(cos(spota)/cos(spota+(anglex10)))) delta = delta*360./wavel x0= spotd*sin(anglex10)/cos(anglex10+spota) phase=phase+delta 10022 continue c--------------------------------------- uph(i)=phase if(uph(i).gt.uphmax)uphmax=uph(i) if(uph(i).lt.uphmin)uphmin=uph(i) c---- ux0(i)=x0 if(ux0(i).gt.ux0max)ux0max=ux0(i) if(ux0(i).lt.ux0min)ux0min=ux0(i) vy0(i)=y0 if(vy0(i).gt.vy0max)vy0max=vy0(i) if(vy0(i).lt.vy0min)vy0min=vy0(i) c---- z0 = 0. if (ajones.gt.0.) then rsq = x0**2 + y0**2 c calculate the shape of the cathode surface z0 = cathod(rsq) zcath(i) = z0 endif c shift initial (x,y,z) so that the particle will emerge from c the photocathode at (x0,y0,z0) at time = phase toff = (wtstep + phase)*wavel/360. cor(1,i) = x0 - toff*bx cor(3,i) = y0 - toff*by cor(5,i) = z0 - toff*bz 10012 continue c---------- endif ! vv(1).eq.1 c---------- if(ip.eq.2) then open(unit=ninput,file='input10',status='unknown') write(ninput,*) npoints,sigdeg,rmax do 10015 i=1,npoints write(ninput,*) uph(i),ux0(i),vy0(i) 10015 continue endif cbm 18/07/2k call histogram1 c---------- c for emittance summary, copy parameters into cord & gam c note that officially this is done at line 210 of the main program ngood = npoints do 10016 i = 1,npoints do 10017 j = 1,6 cord(j,i) = cor(j,i) 10017 continue phizero(i)=uph(i) cord(7,i) = 0. gam(i) = sqrt(1. + cord(2,i)**2 + cord(4,i)**2 + cord(6,i)**2) 10016 continue write(nnout,10018) 10018 format (/' emittance summary for generated particles:') iout7=1 call outlaldp iout7=iout7+1 cbm call out6dp pour type 10 call out6dp(0,0) return c c--- type 11 c c---type 11 = same as input 10 with smooth distributions c uses hammersley sequence c 11000 continue c vv(2) = no. of particles c vv(3) = sigma of laser pulse duration in picoseconds c vv(4) = half width of laser pulse. i.e. generate a gaussian c time profile with sigma = vv(3), but truncate to +-vv(4 c NOT EFFECTIVE in this version. The gaussian is c truncated at roughly +-2.5 sigma. c vv(5) = sigma of laser spot in cm - assumes gaussian radial c profile NOT EFFECTIVE. c vv(6) = max radius of photocathode. i.e. truncate radial dist a c vv(7) = mean K.E. of electrons leaving the photocathode, in MeV c This should be the same as paramter W0 on the RUN card. c vv(8) = sigma of energy-spread distribution at the photocathode c vv(9) = rf step size in deg. This must be the same as parmater c on the START card c--------------- c The following parameters define the way the particles are c generated in each phase space and the correlation between c them. The program transform the decimal number which reper c the particle (j=1,2,....,N) into its equivalent in another c base, inverses each of the digit and reconverts the number c into decimal. This leads for example in base 3: (1/3,2/3,1/9, c 1/3+1/9,....). The base must be a prime number. c If the base is 0 the program just produces a random set using c the usual random fortran generator: RANF(). c vv(10) = bit reversal base for x distribution c if vv(10)<0 don't use bit reversal but produce a c uniform set in x and y. c vv(11)= bit reversal base for y distribution c vv(12)= bit reversal base for phi distribution c if vv(12)<0 call distriphase as in INPUT 10 c vv(13)= bit reversal base for W0 distribution c Not in use in this version. Same as INPUT 10 c vv(14) = spota -- Angle of laser beam with the normal at the cathode c (degrees), if non-zero you need vv(15) c vv(15) = spotd -- Distance from focal point to cathode (cm) c if no zero you need vv(14) c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ tsig = vv(3) tmax = vv(4) rsig = vv(5) rmax = vv(6) rmaxsq = rmax**2 ecath = vv(7) decath = vv(8) nbasex=vv(10) nbasey=vv(11) nbasep=vv(12) nbasew=vv(13) if(nbasep.ge.0) tmax=2.5 tfac=tmax*tsig c convert laser pulse sigma to degrees sigdeg = tsig*freq*360./1.e6 c generate particles during +- tmax in time c set the starting clock back by an integer number of time steps, c corresponding to at least tmax degmax = tfac*freq*360./1.e6 dwt = vv(9) nstep = int(degmax/dwt) if (degmax.ne.dwt*float(nstep)) nstep = nstep + 1 wtstep = nstep*dwt c c offset z of reference particle so it reaches z = 0 at center of cg it will reach z0 not 0 bgr = cor(6,1) bzr = bgr/sqrt(1. + bgr**2) cg cor(5,1) = -wtstep*wavel*bzr/360. cor(5,1) = z0-wtstep*wavel*bzr/360. c----- if(nbasex.lt.0)then call distrirhom(rmax,npoints-1,flag,nhom) npoints=nhom+1 else call loduni(nbasex,npoints-1,p1) call loduni(nbasey,npoints-1,p2) j=1 do 11002 i=np1,npoints x1(i)=(2*p1(i)-1)*rmax y1(i)=(2*p2(i)-1)*rmax r=sqrt(x1(i)**2+y1(i)**2) if(r.le.rmax)then xquiet(j)=x1(i) yquiet(j)=y1(i) j=j+1 endif 11002 continue npoints=j+1 endif c---- laser beam angle if(nn.ge.14) then spota=vv(14)*radian spotd=vv(15) do 11020 j=1,npoints anglex(j) = atan(xquiet(j)/spotd) xquiet(j)= spotd*sin(anglex(j))/cos(anglex(j)+spota) 11020 continue endif c---- uphmin=0. uphmax=0. if(nbasep.lt.0)then call distriphase(degmax,sigdeg,npoints-1) else call lodgau(nbasep,npoints-1,p1,p2,p3,p4,iwork1) endif do 11003 i=np1,npoints c choose phase if(nbasep.lt.0)then phase=distp(i-1) else phase=p1(i-1)*sigdeg endif c---- if((spota.eq.0.).or.(spotd.eq.0)) go to 11023 delta = spotd*(1-(cos(spota)/cos(spota+(anglex(i))))) delta = delta*360./wavel phase=phase+delta 11023 continue uph(i)=phase phizero(i)=uph(i) if(uph(i).gt.uphmax)uphmax=uph(i) if(uph(i).lt.uphmin)uphmin=uph(i) c choose intial energy at cathode 11004 continue call rannor(u,v) e = ecath + u*decath if (e.le.0.) go to 11004 betcat = sqrt(2.*e/erest) cbm warning if 1.-betcat**2 < 0 !!!!! gamma = 1./sqrt(1. - betcat**2) c choose initial momentum, isotropic, but betaz > 0 cose = sandom(d) sine = sqrt(1. - cose**2) phi = twopi*sandom(d) bx = betcat*sine*cos(phi) by = betcat*sine*sin(phi) bz = betcat*cose cor(2,i) = gamma*bx cor(4,i) = gamma*by cor(6,i) = gamma*bz c c choose (x,y,z) c if(nbasex.lt.0)then x0=xhom(i-1) y0=yhom(i-1) else x0=xquiet(i-1) y0=yquiet(i-1) endif c----------- ux0(i)=x0 if(ux0(i).gt.ux0max)ux0max=ux0(i) if(ux0(i).lt.ux0min)ux0min=ux0(i) vy0(i)=y0 if(vy0(i).gt.vy0max)vy0max=vy0(i) if(vy0(i).lt.vy0min)vy0min=vy0(i) c----------- c c shift initial (x,y,z) so that the particle will emerge from c the photocathode at (x0,y0,z0) at time = phase toff = (wtstep + phase)*wavel/360. cor(1,i) = x0 - toff*bx cor(3,i) = y0 - toff*by cor(5,i) = z0 - toff*bz 11003 continue c---------- c for emittance summary, copy parameters into cord & gam c note that officially this is done at line 210 of the main program ngood = npoints do 11016 i = 1,npoints do 11017 j = 1,6 cord(j,i) = cor(j,i) 11017 continue phizero(i)=uph(i) cord(7,i) = 0. gam(i) = sqrt(1. + cord(2,i)**2 + cord(4,i)**2 + cord(6,i)**2) 11016 continue c---- c---- paw en test cbm 18/07/2k call histogram1 if (ip.eq.1 .or. ip.eq.2) . write(nnout,11001) vv(2),ngood,tsig,sigdeg,tfac,rmax, . ecath,decath,dwt,wtstep,nbasex,nbasey,nbasep,nbasew 11001 format ( . ' Generating ',f5.0,' particles from a photocathode',i5, . ' really active particles '/ . ' Quiet start with Hammersley sequence'/ . ' sigma of laser pulse =',t40,1pg12.3,' psec'/ . ' sigma of laser pulse =',t40,g12.3,' degrees'/ . ' tmax of laser pulse = +- 2.5 sigma = ',t40,g12.3,' psec'/ . ' rmax of photocathode =',t40,g12.3,' cm'/ . ' mean electron energy at cathode =',t40,g12.3,' MeV'/ . ' sigma of energy spread at cathode =',t40,g12.3,' MeV'/ . ' integration step size =',t40,g12.3,' degrees'/ . ' clock offset for earliest particle =',t40,g12.3,' degrees'/ . ' bit reversal base for x =',i3 / . ' bit reversal base for y =',i3 / . ' bit reversal base for phi =',i3 / . ' bit reversal base for W0 =',i3 ) if(nn.ge.14) then write(nnout,11025) spota/radian,spotd 11025 format( . ' angle of injection of laser light =',t40,g12.3,' degrees'/ . ' distance focal point - cathode =',t40,g12.3,' cm') endif write(nnout,10018) 11018 format (/' emittance summary for generated particles:') iout7=1 call outlaldp iout7=iout7+1 cbm call out6dp pour type 11 call out6dp(0,0) return ! input 11 c c--- type 12 EGUN c 12000 continue open(unit=3,file='egun.dat',status='old',err=12001) go to 12002 12001 continue call appendparm stop ' Abnormal stop file egun.dat not found ' 12002 continue nray=vv(3) 12003 continue read(3,*)i,regun(i),rdot(i),zdot(i),tdot(i),rcath(i),ioverr(i) if(i.lt.nray)go to 12003 close(3) do 12004 i=1,nray charge(i)=ioverr(i)*rcath(i)*twopi if(i.ge.2)charge(i)=charge(i)+charge(i-1) 12004 continue do 12005 i=1,nray charge(i)=charge(i)/charge(nray) 12005 continue r2=-1.-2./(npoints-2) do 12006 np=2,npoints r1=ranf() r2=r2+2./(npoints-2) r3=ranf()*twopi do 12007 j=1,nray if(r1.le.charge(j))go to 12008 12007 continue c use the jth ray for this particle. determine bgx,bgy,bgz from c rdot,zdot,tdot 12008 continue beta=sqrt(rdot(j)**2+zdot(j)**2+tdot(j)**2) gamma=1./sqrt(1-beta**2) bgz=gamma*zdot(j) cor(6,np)=bgz c determine a z position cor(5,np)=z0+wavel/twopi*vv(4)*r2*radian*zdot(j) c determin x and y x=cos(r3)*regun(j)*vv(5) y=sin(r3)*regun(j)*vv(5) bgx=(rdot(j)*cos(r3)+tdot(j)*sin(r3))*gamma bgy=(rdot(j)*sin(r3)-tdot(j)*cos(r3))*gamma cor(1,np)=x+bgx/bgz*cor(5,np) cor(2,np)=bgx cor(3,np)=y+bgy/bgz*cor(5,np) cor(4,np)=bgy 12006 continue return c c--- type 15 hot cathode microwave electron gun c 15000 continue c megin for hot-cathode microwave electron gun (meg) c H. Liu Aug. 29, 1989 IHEP, Academia Sinica c modified again on oct 16 1989 c vv(2) = no. of particles c vv(3) = E0 ( in MV/m av axial e field, same as in cathode cell card) c vv(4)= initial phase of reference particle (degres) c vv(5) = SIGMAr of cathode (cm) c vv(6) = radius of a thermionic cathode emitter (cm) c vv(7) = temperature (Kelvin) c vv(8) = SIGMAe for energy distribution (MeV) c vv(9) = dwt in deg. c vv(10)= work function (eV) c vv(11)= storing flag if(nn.lt.10) then write(nnout,15001) 15001 format('too few parameters for input type 15 and z0.gt.0.') call appendparm stop ' Abnormal stop input 15 ' endif c e0 = vv(3) rpphase=vv(4) rsig=vv(5) rmax=vv(6) rmaxsq=rmax**2 tempk= 0.001*vv(7) ! (in kilo-Kelvin) decath=vv(8) dwt=vv(9) wf=vv(10) ecath=w0 c ----------------------------- phiran=1. rran=1. if(nn.gt.11) then phiran=vv(12) numphi=vv(13) rran=vv(14) endif write(nnout,*) 'phi random =',phiran,' r random =',rran c ----------------------------- cforj=0.44016*sqrt(e0)/tempk fmax=exp(cforj) c -------------------------------------- if(phiran.eq.0.) then phidiv=180./numphi fsum=0.0 do 15002 i=1,numphi phiarg=(i-0.5)*phidiv fphi=fmax**(sqrt(sin(phiarg*radian))) fsum=fsum+fphi 15002 continue ithphi=1 staphi=0. do 15003 i=1,numphi phinex=(i-0.5)*phidiv fphi=fmax**(sqrt(sin(phinex*radian))) nphidiv=npoints*fphi/fsum subphi=phidiv/nphidiv do 15004 j=1,nphidiv ijk=ithphi+j if(ijk.gt.npoints) go to 15006 phinew=staphi+(j-1)*subphi c print*,ijk,phinew cor(5,ijk)=phinew 15004 continue ithphi=ijk staphi=phinew+subphi 15003 continue c -------- npwant=npoints-ijk if(npwant.gt.0) then do 15007 i=1,npwant phinew=(i-0.5)*180./npwant cor(5,ijk+i)=phinew 15007 continue endif 15006 continue endif c -------------- c The number of background data is calculated c according to iphi*numr*(numr+1)/2. So when c iphi=4 and numr=32, there are 2112 pairs of c data produced. Usually npoints is less than c 2000. If the last circle of particle c distribution is required to be exactly c close, set npoints equal to 2*numr*(numr+1). if(rran.eq.0.) then numra=32 iphi=4 isum0=1 do 15008 i=1,numra sqri=float(i) numphi=iphi*i deltaphi=360./numphi do 15009 j=1,numphi isum=isum0+j phij=(j-1.)*deltaphi phirad=phij*radian cor(1,isum)=sqri*cos(phirad) cor(3,isum)=sqri*sin(phirad) if(isum.eq.npoints) go to 15010 15009 continue isum0=isum 15008 continue 15010 continue factor=rmax/sqri do 15011 i=2,npoints cor(1,i)=factor*cor(1,i) cor(3,i)=factor*cor(3,i) 15011 continue endif c -------------- c -------------------------------------- acj0=120.4*tempk**2*1.e+6*exp(-wf/(8.62*1.e-5*tempk*1000.)) acjmax=acj0*fmax c aci=acj*pi*rmaxsq c ---------------------------------- dppi=pi qint=GPINDP(0.d0,dppi,1.d-6,dpepsout,fsch,1) QCHARGE=0.5*rmaxsq*acj0*qint*1000./freq c -------------------------------------- nstep=int(rpphase/dwt) if(rpphase.ne.dwt*float(nstep)) nstep=nstep+1 wtstep=nstep*dwt c esc=0. escrf=esc/e0 if(ip.eq.1 .or. ip.eq.2) write(nnout,15012) (vv(i),i=1,11), . ecath,acj0,acjmax,fmax,qcharge,qcharge*1.e+10/1.6,nstep,wtstep 15012 format( . ' intype = ',t40,f12.3,' '/ . ' generating ',f6.0,' particles from a thermionic cathode'/ . ' average axial e field =',t40,g12.3,' MV/m'/ . ' initial phase of r.p. =',t40,g12.3,'degrees'/ . ' SIGMAr for emitter =',t40,g12.3,' cm'/ . ' radius of emitter =',t40,g12.3,'cm'/ . ' operating temperature =',t40,g12.3, 'kilo-Kelvin'/ . ' SIGMAe of energy distribution =',t40,g12.3,'MeV'/ . ' rf step size dwt =',t40,g12.3,' degrees'/ . ' work function of emitter =',t40,g12.3,'eV'/ . ' storing flag =',t40,g12.3,' '/ . ' mean kinetic energy W0 =',t40,g12.3,' MeV'/ . ' current density J0 =',t40,g12.3,' A/cm**2'/ . ' current density Jmax =',t40,g12.3, ' A/cm**2'/ . ' Fmax (Jmax/J0) =',t40,g12.3,' '/ . ' total charge amount in the bunch =',t40,g12.3,' nC'/ . ' number of electrons in the bunch =',t40,g12.3, ' '/ . ' nstep= ',t40,i6,' '/ . ' wtstep =',t40,g12.3,' degrees'/) c zcath(1)=z0 c bgr=cor(6,1) bzr=bgr/sqrt(1.+bgr**2) cor(5,1)=zcath(1)-wtstep*wavel*bzr/360. c ----------------------- c-------------------------------------------------------------------------- uphmin=0. uphmax=0. ux0min=0. ux0max=0. vy0min=0. vy0max=0. c-------------------------------------------------------------------------- do 15013 i=np1,npoints if(phiran.ne.0.) then c choose initial phases 15014 ph1=pi*ranf() se=sin(ph1)+escrf if (se.le.0.) go to 15014 fph1=fmax**sqrt(se) ph2=rn32() fph2=fmax*ph2 if(fph2.gt.fph1) go to 15014 phase=ph1/radian else phase=cor(5,i) endif c---- uph(i)=phase if(uph(i).gt.uphmax)uphmax=uph(i) if(uph(i).lt.uphmin)uphmin=uph(i) c---- c choose initial energy at cathode 15015 call norran(u) e=ecath+u*decath if(e.le.0.) go to 15015 betcat=sqrt(2.*e/erest) gamma=1./sqrt(1.-betcat**2) c choose initial momentum,isotropic, but betaz>0 sand1=ranf() theta=(sand1-0.5)*pi cose=cos(theta) sine=sqrt(1.-cose**2) phi=twopi*sandom(d) bx=betcat*sine*cos(phi) by=betcat*sine*sin(phi) bz=betcat*cose c ------------------------------------ cor(2,i)=gamma*bx cor(4,i)=gamma*by cor(6,i)=gamma*bz c choose (x,y,z) if(rran.ne.0.) then 15016 call rannor(u,v) x0=u*rsig y0=v*rsig rsq=x0**2+y0**2 c if(rsq.gt.rmaxsq .or. rsq.lt.rinsq) go to 15016 if(rsq.gt.rmaxsq) go to 15016 else x0=cor(1,i) y0=cor(3,i) endif c---- ux0(i)=x0 if(ux0(i).gt.ux0max)ux0max=ux0(i) if(ux0(i).lt.ux0min)ux0min=ux0(i) vy0(i)=y0 if(vy0(i).gt.vy0max)vy0max=vy0(i) if(vy0(i).lt.vy0min)vy0min=vy0(i) c---- c c shift x,y,z zcath(i)=zcath(1) cliu toff=(wtstep+phase)*wavel/360. toff=phase*wavel/360. cor(1,i)=x0-toff*bx cor(3,i)=y0-toff*by cor(5,i)=zcath(i)-toff*bz c -------------------- phizero(i)=phase c ----------------------------------- 15013 continue c-------------------------------------------------------------------------- if(ip.eq.2) then open(unit=ninput,file='input15',status='unknown') write(ninput,*) npoints,sigdeg,rmax do 15017 i=1,npoints write(ninput,*) uph(i),ux0(i),vy0(i) 15017 continue endif c--- paw en test cbm 18/07/2k call histogram1 c-------------------------------------------------------------------------- c ---------------------------- c for emittance summary, copy parameters into cord & gam c note that officially this is done at line 210 of the main program ngood = npoints c --------------- do 15018 i = 1,npoints do 15019 j = 1,6 cord(j,i) = cor(j,i) 15019 continue cord(7,i) = 0. gam(i) = sqrt(1. + cord(2,i)**2 + cord(4,i)**2 + cord(6,i)**2) 15018 continue if (ip.ne.0) then write(nav,*) ' ' write (nav,15020) 15020 format (/' emittance summary for generated particles:') cbm call out6dp pour type 15 call out6dp(0,0) endif return ! return input 15 end c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*