program ICOOL c c Monte carlo simulation of muon ionization cooling c character*4 namepath ! namepath = ' ' call ICOOLS(namepath) ! end c ********************************************************* subroutine ICOOLS(namepath) c c Monte carlo simulation of muon ionization cooling c implicit double precision (a-h,o-z) include 'icommon.inc' character*(*) namepath c version = '3.29f' ! pathname = namepath ! move I/O path specifier into common ! pi = 4.d0*ATAN2(1.d0,1.d0) twopi = 2.d0 * pi piov2 = pi / 2.d0 permfp = 4.d-7 * pi dbleps = 1d0 100 dbleps = dbleps*0.5d0 if (1d0+dbleps.gt.1d0) goto 100 ! call SETUP if( iflag .ne. 0 ) go to 900 ! if( lelms ) call ELMS_INIT c call SIMULATE c call FINISH c 900 continue end c ********************************************************* BLOCK DATA ! ! Initialize common block data in F77-compatible manner ! implicit double precision (a-h,o-z) include 'icommon.inc' ! ! backgrn2 data bftag/'NONE'/ data ztotalbkg/0.d0/ ! ! beamdat data nbeamtyp/1/,bmalt/.false./ ! ! cav4 data nzom4/20/,nrom4/20/ ! const data cc/2.99792458d8/,ee/0.299792458d0/, + avagn/6.02214199d23/,relcl/2.817940285d-13/ data j01/2.4048255576957727686216318793264d0/ data mass/0.51099892d-3, 0.105658367d0, 0.13957018d0, ! GeV/c^2 & 0.493677d0, 0.93827201d0, & 1.8756128d0, 2.809413d0, 6.535366/ ! d,He3,Li7 data chrg/6*1d0,2d0,3d0/ ! ! control data beep/.false./,bgen/.true./,bunchcut/1.d6/, & dectrk/.false./,diagref/.false./,epsf/0.05d0/,epsreq/1d-3/, & epsstep/1d-6/,forcerp/.true./, & fsav/.false./,fsavset/.false./, & ffcr/.false./,f9dp/4/,goodtrack/.true./,magconf/0/,mapdef/0/, & neighbor/.false./,neutrino/0/,nnudk/1/, & nprnt/-1/,npskip/0/,nsections/1/,ntuple/.false./,nuthmin/0./, & nuthmax/3.141593/, & output1/.false./,phantom/.false./,phasemodel/1/,prlevel/1/, & prnmax/300/,pzmintrk/0.001/, & rfdiag/0/,rfphase/0/,rnseed/-1/, & rtuple/.false./,rtuplen/5/,run_env/.false./, & scalestep/1./,spin/.false./, & spinmatter/0/, & spintrk/0/, & stepmin/1d-5/,stepmax/1.d0/,steprk/.true./, & summary/.true./,termout/.true./, & timelim/1e9/, & varstep/.true./ ! ! covars data ncovar/0/ ! ! dendat data dirdenp/-1/ ! ! dendatc data matdens/' '/,matdenp/' '/ ! ! dvars data labvar/' X',' Y',' Z',' Px',' Py',' Pz',' ct',' Pm', & ' E',' KE', ' X''',' Y''',' TH',' R','PHI', & ' Pr','Pph',' Lz',' L2',' S', & ' Pt',' ',' ',' ',' ',' ',' ',' A2', & ' R2','HEL',' Bx',' By',' Bz',' Ex',' Ey',' Ez', & ' Sx',' Sy',' Sz','phs'/ ! ! emits data nemit/0/,emevt/100*0/,emevcut/100*0/ data emitx/100*0./,emity/100*0./,emitz/100*0./,emit6/100*0./ data sigmacut/.true./,pxycorr/.false./,edetail/.false./ data sig_cut/4./,discorr/.false./,ipzcor/0/ ! ! hist data nhist/0/,hauto/.true./,hcprn/0/ data ihcon/20*0/,ihund/20*0/,ihovf/20*0/,hval/1000*0/ ! ! intdat data ldedx/.true./,lscatter/.true./,lbrem/.false./, & ldray/.true./,trainsc/.false./, + linteract/.false./,ldecay/.true./,lspace/.false./, & lstrag/.true./,lelms/.false./,lsamcs/.false./ data delev/2/,straglev/4/,scatlev/6/,spacelev/1/, + bremlev/1/,intlev/1/,declev/1/,draylev/1/ data facfms/1./,facmms/1./,dcute/0.003d0/,dcutm/0.003d0/, & parbunsc/4e12/,pdelev5/0.200/,gfactsc/.false./, & frfbunsc/201./,wanga/90.1/,wangb/3.35/,wangc/1.22/, & wangd/4.66/,wangpmx/1.5/,wangfmx/13.706/, & fastdecay/.false./,elmscor/0/ ! ! intdatc ! data elmsdir/' '/,elmsroot/' '/ ! ! matdata data nmat/19/,matid/'LH','GH','LHE','GHE','LI','BE','B','C', & 'AL','TI','FE','CU','W','HG','PB','LIH','AM','CH2','SS'/ ! data order: Z rho __ Lr Xo X1 Sa Sc Sm ! : nel { Zel(i), Ael(i), IPel(i), fel(i) ,i=1,nel} data matdat/ & 1.,0.071,0.,866., 0.4759,1.9215,0.1348,-3.2632,5.6249, & 1, 1.,1.008,21.8,1., 8*0., ! Liq H & 1.,8.38d-5,0.,7.31d5, 0.,0.,0.,0.,0., & 1, 1.,1.008,19.2,1., 8*0., ! Gas H & 2.,0.125,0.,756., 0.48,1.93,0.134,-3.3,5.7, & 1, 2.,4.003,42.3,1., 8*0., ! Liq He & 2.,1.66d-4,0.,5.67d5, 0.,0.,0.,0.,0., & 1, 2.,4.003,41.8,1., 8*0., ! Gas He & 3.,0.534,0.,155., 0.1304,1.6397,0.9514,-3.1221,2.4993, & 1, 3.,6.94,40.0,1., 8*0., ! Li & 4.,1.848,0.,35.3, 0.0592,1.6922,0.8039,-2.7847,2.4339, & 1, 4.,9.01,63.7,1., 8*0., ! Be & 5.,2.370,0.,22.2, 0.0305,1.9688,0.5622,-2.8477,2.4512, & 1, 5.,10.811,76.0,1., 8*0., ! B & 6.,2.265,0.,18.8, -0.0178,2.3415,0.2614,-2.8680,2.8697, & 1, 6.,12.011,78.0,1., 8*0., ! C & 13.,2.699,0.,8.9, 0.1708,3.0127,0.0802,-4.2395,3.6345, & 1, 13.,26.982,166.,1., 8*0., ! Al & 22.,4.54,0.,3.56, 0.0957,3.0386,0.1566,-4.4450,3.0302, & 1, 22.,47.88,233.,1., 8*0., ! Ti & 26.,7.87,0.,1.76, -0.0012,3.1531,0.1468,-4.2911,2.9632, & 1, 26.,55.845,286.,1., 8*0., ! Fe & 29.,8.96,0.,1.43, -0.0254,3.2792,0.1434,-4.4190,2.9044, & 1, 29.,63.546,322.,1., 8*0., ! Cu & 74.,19.3,0.,0.35, 0.2167,3.4960,0.1551,-5.406,2.845, & 1, 74.,183.85,727.,1., 8*0., ! W & 80.,13.546,0.,0.475, 0.2756,3.7275,0.1101,-5.9605,3.0519, & 1, 80.,200.59,800.,1., 8*0., ! Hg & 82.,11.35,0.,0.561, 0.3776,3.8073,0.0936,-6.2018,3.1608, & 1, 82.,207.21,823.,1., 8*0., ! Pb & 4.,0.780,0.,102., -0.0988,1.4515,0.9057,-2.3580,2.5849, & 2, 3.,6.94,40.0,1., 1.,1.008,21.8,1., 4*0., ! LiH & 4.,2.1,0.,19.2, 5*0., 2, 4.,9.01,63.7,0.64, & 13.,26.982,166.,0.36, 4*0., ! AlBeMet & 8.,0.94,57.4,47.9, 0.1370,2.5177,0.1211,-3.0016,3.4292, & 2, 6.,12.001,78.,1., 1.,1.008,21.8,2., 4*0., ! CH2 (polyethylene) & 26.,8.00,0.,1.76, 5*0., 3, 26.,55.85,275.,0.72, & 24.,52.00,260.,0.19, 28.,58.69,310.,0.09/ ! 304 SS ! ! polar data polx/100*0./,poly/100*0./,polz/100*0./,helic/100*0./ ! ! refpart data xpref/3*0.d0/,ppref/3*0.d0/,tpref/0.d0/,trefmean/0.d0/ data t2refmean/0d0/ data iptypref/0/,ipflgref/0/,phmodref/0/ ! ! regndat data nfparm/15/,ngparm/10/ ! ! rfdat data phase/0./ ! ! rhist data nrhist/0/,rhauto/.true./,rhprin/.false./ ! ! scat data iscon/20*0/,isxund/20*0/,isyund/20*0/,isxovf/20*0/ data isyovf/20*0/,sval/23000*0/,nscat/0/ data sauto/.true./ ! ! simul data krad/1/,nfail/0/ data htrk/0.d0/,hptrk/0.d0/,gtrk/0.d0/,gptrk/0.d0/,tortion/0.d0/ data dsffag/0.d0/ data isffile/-1/ data firstsreg/.true./ ! ! taper data ncelltap/-1/,jctap/0/ ! ! zhist data nzhist/0/,izhknt/220*0/,zhauto/.true./,zhprin/.false./ ! end c ******************************************************** subroutine BACKGET(x,t,e,b) ! ! Get background field ! implicit double precision (a-h,o-z) include 'icommon.inc' ! save jx,jy,jz dimension x(3),e(3),b(3) data eps/1.d-6/ ! ! --------------------------------------------------------- do i=1,3 e(i) = 0. b(i) = 0. end do ! if( x(3).lt.slogr .or. x(3).gt.slogr+ztotalbkg+eps ) go to 900 ! if( bftag .eq. 'RFUS' ) then call BACKRFGET(x,t,e,b) ! else x4 = x(1) y4 = x(2) z4 = x(3) call HUNT(xgr,nxgr,x4,jx) ! returns lower grid corner call HUNT(ygr,nygr,y4,jy) call HUNT(sgr,nsgr,z4,jz) ! if( jx .eq. 0 ) jx = 1 ! handle case when particle if( jy .eq. 0 ) jy = 1 ! is exactly on bottom edge if( jz .eq. 0 ) jz = 1 ! of the grid ! if( interbkg .le. 1 ) then ! linear u = (x4-xgr(jx)) / dxgr v = (y4-ygr(jy)) / dygr w = (z4-sgr(jz)) / (dsgr*(1.+hfldbkg*x4)) ! call TRILINEAR(u,v,w,jx,jy,jz,bxgr,mxxg,mxyg,mxsg,b(1)) call TRILINEAR(u,v,w,jx,jy,jz,bygr,mxxg,mxyg,mxsg,b(2)) call TRILINEAR(u,v,w,jx,jy,jz,bsgr,mxxg,mxyg,mxsg,b(3)) ! if( bentbkg ) then htrk = htrkgr(jz) + w*(htrkgr(jz+1)-htrkgr(jz)) gtrk = gtrkgr(jz) + w*(gtrkgr(jz+1)-gtrkgr(jz)) else htrk = 0.d0 gtrk = 0.d0 end if ! else if( interbkg.eq.2 .or. interbkg.eq.3 ) then ! if( bentbkg ) then kpt = interbkg + 1 k = MIN( MAX(jz-(kpt-1)/2,1), nsgr+1-kpt ) call POLINT(sgr(k),htrkgr(k),kpt,z4,htrk,dy) call POLINT(sgr(k),hptrkbkg(k),kpt,z4,hptrk,dy) ordtrk = ordtrkbkg end if ! ! Find offset into grid arrays kx = MIN( MAX(jx-interbkg/2,1), nxgr-interbkg ) ky = MIN( MAX(jy-interbkg/2,1), nygr-interbkg ) kz = MIN( MAX(jz-interbkg/2,1), nsgr-interbkg ) ! call TRIPOLY(kx,ky,kz,xgr,ygr,sgr,bxgr,mxxg,mxyg,mxsg, & interbkg,x4,y4,z4,htrk,b(1),df) call TRIPOLY(kx,ky,kz,xgr,ygr,sgr,bygr,mxxg,mxyg,mxsg, & interbkg,x4,y4,z4,htrk,b(2),df) call TRIPOLY(kx,ky,kz,xgr,ygr,sgr,bsgr,mxxg,mxyg,mxsg, & interbkg,x4,y4,z4,htrk,b(3),df) ! end if ! end of test on interbkg end if ! end of test on bftag ! --------------------------------------------------------- ! 900 continue return end c ********************************************************* subroutine BACKRFGET(x,t,e,b) ! ! Get background RF field ! implicit double precision (a-h,o-z) include 'icommon.inc' ! real*8 x(3),t,e(3),b(3) ! dum = t ! defeat annoying compiler warnings dum = x(1) dum = e(1) dum = b(1) ! return end c ********************************************************* subroutine BACKSET(iarg) ! ! Generate background field grid ! implicit double precision (a-h,o-z) include 'icommon.inc' ! dimension x(3),e(3),b(3),xloc(3) character*10 fname ! if( iarg .eq. 0 ) then ! initialize ! do ix=1,nxgr xgr(ix) = (ix-1) * dxgr + xlogr end do do iy=1,nygr ygr(iy) = (iy-1) * dygr + ylogr end do do iz=1,nsgr sgr(iz) = (iz-1) * dsgr end do ! do ix=1,nxgr do iy=1,nygr do iz=1,nsgr bxgr(ix,iy,iz) = 0. ! clear background array bygr(ix,iy,iz) = 0. bsgr(ix,iy,iz) = 0. end do end do end do end if ! ! --------------------------------------------------------- if( iarg .eq. 1 ) then ! add background field if( bftag .eq. 'STUS' ) then ! read in static user field mfile = bfparm(2) bscfac = bfparm(5) write(fname,40) mfile 40 format('for0',i2,'.dat') open(unit=mfile, & file=pathname(1:LASTNB(pathname)) // fname, & status='unknown',iostat=ioc) if( ioc .ne. 0 ) then write(2,*)'Couldnt open BACKGROUND user static field' iflag = -1 go to 900 end if ! read(mfile,*) !skip title record read(mfile,*) mx,my,mz,hfldbkg if( mx.ne.nxgr .or. my.ne.nygr .or. mz.ne.nsgr ) then write(2,*)'Dimension mismatch BACKGROUND user static field' iflag = -1 go to 900 end if ! ! if( hfldbkg .ne. 0.d0 ) ncurved = 1 ! read(mfile,*) (xgr(j),j=1,nxgr) read(mfile,*) (ygr(j),j=1,nygr) read(mfile,*) (sgr(j),j=1,nsgr) ! nrecords = mx*my*mz ! do i=1,nrecords read(mfile,*)ix,iy,iz,bx,by,bz bxgr(ix,iy,iz) = bxgr(ix,iy,iz) + bx * bscfac bygr(ix,iy,iz) = bygr(ix,iy,iz) + by * bscfac bsgr(ix,iy,iz) = bsgr(ix,iy,iz) + bz * bscfac end do ! close(mfile) ! ......................................................... else if( bftag .eq. 'RFUS' ) then ! read in RF user field mfile = bfparm(1) freq = bfparm(2) * 1.e6 phas0 = bfparm(3) * pi / 180. write(fname,50) mfile 50 format('for0',i2,'.dat') open(unit=mfile, & file=pathname(1:LASTNB(pathname)) // fname, & status='unknown',iostat=ioc) if( ioc .ne. 0 ) then write(2,*)'Couldnt open BACKGROUND user rf field' iflag = -1 go to 900 end if ! read(mfile,*) !skip title record read(mfile,*) mx,my,mz if( mx.ne.nxgr .or. my.ne.nygr .or. mz.ne.nsgr ) then write(2,*)'Dimension mismatch BACKGROUND user rf field' iflag = -1 go to 900 end if ! nrecords = mx*my*mz ! do i=1,nrecords read(mfile,*)ix,iy,iz,bx,by,bz,ex,ey,ez bxgr(ix,iy,iz) = bxgr(ix,iy,iz) + bx bygr(ix,iy,iz) = bygr(ix,iy,iz) + by bsgr(ix,iy,iz) = bsgr(ix,iy,iz) + bz exgr(ix,iy,iz) = exgr(ix,iy,iz) + ex eygr(ix,iy,iz) = eygr(ix,iy,iz) + ey esgr(ix,iy,iz) = esgr(ix,iy,iz) + ez end do ! close(mfile) ! ......................................................... else ! use predefined field types ! do 250 ix=1,nxgr x(1) = xgr(ix) ! do 230 iy=1,nygr x(2) = ygr(iy) r = SQRT(x(1)**2 + x(2)**2) if( r .gt. rmaxbkg ) go to 230 ! do 210 iz=1,nsgr ! The grid zbkg is defined in relative coordinates x(3) = sgr(iz) if( x(3).lt.zminbkg .or. x(3).gt.zmaxbkg ) go to 210 ! ! Translate into magnet local coordinate system xloc(1) = x(1) xloc(2) = x(2) xloc(3) = x(3) - zoffbkg ! call FIELDFNC(bftag,bfparm,xloc,t,e,b) ! bxgr(ix,iy,iz) = bxgr(ix,iy,iz) + b(1) bygr(ix,iy,iz) = bygr(ix,iy,iz) + b(2) bsgr(ix,iy,iz) = bsgr(ix,iy,iz) + b(3) ! 210 continue ! 230 continue ! 250 continue end if ! end of test on bftag end if ! end of test on iarg=1 ! 900 continue return end ! ********************************************************* subroutine BEG_OF_REGION(iscell) ! ! Handle tasks required at the beginning of a region ! implicit double precision (a-h,o-z) include 'icommon.inc' integer k,rc,EDGE_ACCEL_13 double precision r,bs,pr,pt logical iscell if (.not.iscell) then r = sqrt(xp(1)*xp(1)+xp(2)*xp(2)) k = 0 100 k = k + 1 if (k.le.nrreg .and. (r.lt.rlow(k) .or. r.gt.rhigh(k))) goto 100 end if if (iscell .or. k.le.nrreg) then if ( .not.iscell .and. ftag(k).eq.'ACCE' ) then if (jpar.eq.0) call ACCEL_PHASE(fparm(1,k)) if (iflag.eq.-1) return if ( fparm(1,k).eq.13 .and. $ (fparm(5,k).eq.0 .or. fparm(5,k).eq.2) ) then rc = EDGE_ACCEL_13(xp,pp,tp,pcharge/ee,pmass, $ twopi*fparm(2,k)*1d6,phase,1d-3*fparm(3,k)) if (rc.ne.0) iflag = -86 if (rtuple) call OUTPUT end if ! else if (.not.iscell.and.ftag(k).eq.'SOL'.and.fparm(1,k).eq.8) $ then pr = pp(1)*pp(1)+pp(2)*pp(2)+pp(3)*pp(3) bs = fparm(4,k)*0.25d0*pcharge*pcharge/sqrt(pr) pp(1) = pp(1) + bs*xp(1) pp(2) = pp(2) + bs*xp(2) tp = tp - 0.5d0*bs*r*r*sqrt(pr+pmass*pmass)/(cc*pr) if (fparm(3,k).eq.0.or.fparm(3,k).eq.2) then bs = 0.5d0*pcharge*fparm(2,k) pp(1) = pp(1) + bs*xp(2) pp(2) = pp(2) - bs*xp(1) end if pt = pr-pp(1)*pp(1)-pp(2)*pp(2) if (pt.lt.0d0) then iflag = -86 else pp(3) = sqrt(pt) end if if (rtuple) call OUTPUT ! else if(.not.iscell .and. nbffag.eq.1d0) then by = fparm(2,k) r0 = fparm(5,k) th2 = slen / r0 dsffag = 2d0 * xp(1) * SIN(0.5d0*th2) call VMAG(pp,pm) dpx = pcharge * by * dsffag pp(1) = pp(1) + dpx pp(3) = SQRT( pm*pm - pp(1)**2 - pp(2)**2 ) end if end if end ! ********************************************************* subroutine BIPOLYINT(jx,jy,xg,yg,fg,mp,np,kk,x,y,f,df) ! ! Polynomial interpolation on 2-dimensional grid ! ! cf. Num. Rec., 2nd ed, p.103,118; uses POLINT ! ! kk: polynomial order {1,2,3} ! implicit double precision (a-h,o-z) dimension fg(mp,np),xg(mp),yg(np),fxtmp(4),fytmp(4) ! kx = MIN( MAX(jx-kk/2,1), mp-kk ) ky = MIN( MAX(jy-kk/2,1), np-kk ) ! do j=1,kk+1 ! loop over rows ! do k=1,kk+1 ! save row in temporary array fytmp(k) = fg(kx+j-1,ky+k-1) end do ! call POLINT(yg(ky),fytmp,kk+1,y,fxtmp(j),df) end do ! call POLINT(xg(kx),fxtmp,kk+1,x,f,df) ! return end c ********************************************************* logical function BLANK(ch) ! ! Return true if this character is a space or tab ! character*1 ch,spaceChar, tabChar parameter (spaceChar = ' ') ! ! Ascii character code for horizontal tab tabChar = char(9) ! if( ch.eq.spaceChar .or. ch.eq.tabChar ) then blank = .true. else blank = .false. end if ! end c ********************************************************* subroutine BORIS_STEP(y,neqn,x) c ! Propagate particle one step (hdid) thru field ! Recommend new step size (hnextf) based on field variation c implicit double precision (a-h,o-z) include 'icommon.inc' c real*8 myu,mycrossE,mycrossB1,mycrossB2,myfront,mydz2,mylambda real*8 myefld dimension y(neqn),myu(3),v(3),myefld(3),delspin(3) ! ! X = s ! Y = {x,y,t,x',y',Ps; Sx,Sy,Sz} ! ps =y(6) px = ps * y(4) py = ps * y(5) pm = DSQRT( px**2 + py**2 + ps**2 ) gamma = DSQRT(1.d0 + pm*pm/pmass/pmass) h = hdid h2 = 0.5d0 * hdid mydz2 = 0.5d0 * pcharge * h ! q*dz/2 c c Peter Stoltz added this section, 06/20/01 c Advance momenta/energy using Boris push (Cary's notes) c First store px,py,U/c in the array v: c myu(1) = px myu(2) = py myu(3) = gamma*pmass c c First, do a half-push of the positions with the c old velocity c y(1) = y(1) + 0.5d0*h*y(4) y(2) = y(2) + 0.5d0*h*y(5) c This uses the above def. of velocity c Don't forget the units of ps are not MKS right now y(3) = y(3) + 0.5d0*h*myu(3)/(ps*cc) x = x + 0.5d0*h c c Need fields, so first construct the vector x,y,z xp(1) = y(1) xp(2) = y(2) xp(3) = x c Then pass to FIELD call FIELD( xp(1),y(3) ) do i=1,3 ! convert to E/c myefld(i) = efld(i)/cc end do c c half-step the part that direct-pushes p_z c so, use 0.5*h*q=mydz2 myu(1) = myu(1) - mydz2*bfld(2) myu(2) = myu(2) + mydz2*bfld(1) myu(3) = myu(3) + mydz2*myefld(3) c Recalculate p_s (will then be constant in the next step) ps2 = myu(3)*myu(3)-(myu(1)*myu(1)+myu(2)*myu(2)+pmass*pmass) if (ps2 .gt. 0.d0) then ps = DSQRT(ps2) else ! particle reflected before reaches end of step y(6) = -y(6) go to 900 end if c qfact = mydz2/ps ! 1/2 * dz * charge / pz mylambda=qfact*qfact*(bfld(3)*bfld(3)- + (myefld(1)*myefld(1)+myefld(2)*myefld(2))) c c Bump using the matrix that doesn't change p_z mycrossE=myefld(1)*myefld(2) mycrossB1=myefld(1)*bfld(3) mycrossB2=myefld(2)*bfld(3) myfront=2.d0*qfact/(1.d0+mylambda) v(1)=myfront*( + myu(1)*qfact*(myefld(1)*myefld(1) - bfld(3)*bfld(3))+ + myu(2)*(bfld(3) + qfact*mycrossE) + myu(3)*(myefld(1) + + qfact*mycrossB2)) v(2)=myfront*(myu(1)*(qfact*mycrossE - bfld(3)) + + myu(2)*qfact*(myefld(2)*myefld(2) - bfld(3)*bfld(3))+ + myu(3)*(myefld(2) - qfact*mycrossB1)) v(3)=myfront*(myu(1)*(myefld(1) - qfact*mycrossB2) + + myu(2)*(myefld(2)+qfact*mycrossB1)+ + myu(3)*qfact*(myefld(1)*myefld(1)+myefld(2)*myefld(2))) c c Bump the second half of the step that direct-pushes p_z c recall that the v's above are a delta-v myu(1) = myu(1) + v(1) - mydz2*bfld(2) myu(2) = myu(2) + v(2) + mydz2*bfld(1) myu(3) = myu(3) + v(3) + mydz2*myefld(3) ! c Recalculate p_z ps2 = myu(3)*myu(3)-(myu(1)*myu(1)+myu(2)*myu(2)+pmass*pmass) if (ps2 .gt. 0.d0) then ps = DSQRT(ps2) else ! particle reflected before reaches end of step y(6) = -y(6) go to 900 end if c c Put correct t into y array y(3) = y(3) + 0.5d0*h*myu(3)/(ps*cc) c c Put dx/ds and dy/ds into the y-array y(4) = myu(1)/ps y(5) = myu(2)/ps c Bump positions using new velocities y(1) = y(1)+0.5d0*h*y(4) y(2) = y(2)+0.5d0*h*y(5) c c Put new p_s into y array y(6) = ps c if (neqn .eq. 9) then bmt2 = axbmt(1)*axbmt(1)+axbmt(2)*axbmt(2)+ & axbmt(3)*axbmt(3) spinlambda=h2*h2*bmt2 spincross1=axbmt(2)*axbmt(3) spincross2=axbmt(3)*axbmt(1) spincross3=axbmt(1)*axbmt(2) spinfront=h/(1.d0+spinlambda) delspin(1) = spinfront*( & -y(7)*h2*(axbmt(2)*axbmt(2) + axbmt(3)*axbmt(3)) & +y(8)*(axbmt(3) + h2*spincross3) & +y(9)*(-axbmt(2) + h2*spincross2)) delspin(2) = spinfront*( & -y(8)*h2*(axbmt(3)*axbmt(3) + axbmt(1)*axbmt(1)) & +y(9)*(axbmt(1) + h2*spincross1) & +y(7)*(-axbmt(3) + h2*spincross3)) delspin(3) = spinfront*( & -y(9)*h2*(axbmt(1)*axbmt(1) + axbmt(2)*axbmt(2)) & +y(7)*(axbmt(2) + h2*spincross2) & +y(8)*(-axbmt(1) + h2*spincross1)) y(7)=y(7) + delspin(1) y(8)=y(8) + delspin(2) y(9)=y(9) + delspin(3) end if ! c Put correct z into x x = x + 0.5d0*h c 900 continue return end c ********************************************************* subroutine CHKTIME(tlim) c c Check elapsed time from start of execution c time calculated in minutes c implicit double precision (a-h,o-z) include 'icommon.inc' ! integer*2 iyr,imon,iday,ihr,imin,isec c call ICOOLDATE(iyr,imon,iday) call ICOOLTIME(ihr,imin,isec) t0 = ihr0*60 + imin0 t = ihr*60 + imin if(iday .ne. iday0) t = t + (iday-iday0)*1440 dt = t - t0 if(dt .gt. tlim) then iflag = -1 write(2,*) 'Exceeded time limit' end if c return end c ********************************************************* function CURVED(chf,kk) ! ! Determine if this field uses a curved coordinate system ! ! kk = 1:region 2:cell 3:background ! implicit double precision (a-h,o-z) include 'icommon.inc' ! logical curved character(4) chf ! curved = .false. ! if( chf.eq.'BROD' .or. chf.eq.'BSOL' .or. chf.eq.'DIP' ) then curved = .true. go to 900 end if ! if( kk .eq. 1 ) then ! region if(chf.eq.'HDIP' .and. fparm(3,krad).eq.1 ) curved = .true. if(chf.eq.'STUS' .and. fparm(3,krad).eq.1 ) curved = .true. if(chf.eq.'G43D' .and. fparm(3,krad).eq.1 ) curved = .true. if(chf.eq.'HELI' .and. fparm(1,krad).eq.6 ) curved = .true. if(chf.eq.'EFLD' .and. fparm(12,krad).ne.0d0 ) curved = .true. ! else if( kk .eq. 2 ) then ! cell if(chf.eq.'HDIP' .and. cfparm(3).eq.1 ) curved = .true. if(chf.eq.'STUS' .and. cfparm(3).eq.1 ) curved = .true. if(chf.eq.'G43D' .and. cfparm(3).eq.1 ) curved = .true. if(chf.eq.'HELI' .and. cfparm(1).eq.6 ) curved = .true. if(chf.eq.'EFLD' .and. cfparm(12).ne.0d0 ) curved = .true. ! else if( kk .eq. 3 ) then ! background if(chf.eq.'HDIP' .and. bfparm(3).eq.1 ) curved = .true. if(chf.eq.'STUS' .and. bfparm(3).eq.1 ) curved = .true. if(chf.eq.'G43D' .and. bfparm(3).eq.1 ) curved = .true. if(chf.eq.'HELI' .and. bfparm(1).eq.6 ) curved = .true. if(chf.eq.'EFLD' .and. bfparm(12).ne.0d0 ) curved = .true. ! end if ! 900 continue end c ********************************************************* SUBROUTINE DCROSS(A,B,C) c c Vector cross product ! implicit double precision (a-h,o-z) dimension A(3),B(3),C(3) c C(1)=A(2)*B(3)-A(3)*B(2) C(2)=A(3)*B(1)-A(1)*B(3) C(3)=A(1)*B(2)-A(2)*B(1) c RETURN END c ********************************************************* subroutine DERIV(s,y,dyds) c c Compute values of derivatives dYdS for the variables Y at s c implicit double precision (a-h,o-z) include 'icommon.inc' c dimension y(15),dyds(15),xscr(3) ! ! Use space stepping ! Y = {x, y, t, Px/Ps, Py/Ps, Ps; Sx, Sy, Sz} ! xscr(1) = y(1) xscr(2) = y(2) xscr(3) = s call FIELD( xscr,y(3) ) ! ! Watch out for crazy stepping results if( ABS(y(1)).gt.100. .or. ABS(y(2)).gt.100. .or. & ABS(y(4)*y(6)).gt.1000. .or. ABS(y(5)*y(6)).gt.1000. ) then iflag = -76 go to 900 end if ! arg = y(4)**2+y(5)**2+(pmass/y(6))**2+1.0d0 if( arg .lt. 0.d0 ) then iflag = -12 go to 900 end if eonps = SQRT(arg) ! 1 / betaZ ! if( ncurved .eq. 0 ) then ! straight section dyds(1) = y(4) dyds(2) = y(5) dyds(3) = eonps/cc ! Load xscr with dPi/dt / Vz xscr(1) = pcharge*(eonps*efld(1)/cc + y(5)*bfld(3) - bfld(2)) xscr(2) = pcharge*(eonps*efld(2)/cc - y(4)*bfld(3) + bfld(1)) xscr(3) = pcharge* & (eonps*efld(3)/cc + y(4)*bfld(2) - y(5)*bfld(1)) dyds(4) = (xscr(1)-y(4)*xscr(3))/y(6) dyds(5) = (xscr(2)-y(5)*xscr(3))/y(6) dyds(6) = xscr(3) ! else ! curved section ! c Version 020821a c c Report all suspected bugs to J. Scott Berg c or c c $Date: 2002/08/21 19:22:21 $ c $Revision: 1.6 $ c $Source: /home/jsberg/cvsroot/icool/integ/deriv1.f,v $ ! hs = 1.0d0 + htrk*y(1) + gtrk*y(2) dyds(1) = hs*y(4) + tortion*y(2) dyds(2) = hs*y(5) - tortion*y(1) dyds(3) = hs*eonps/cc hs = hs*pcharge xscr(1) = (htrk + tortion*y(5))*y(6) + $ hs*(eonps*efld(1)/cc + y(5)*bfld(3) - bfld(2)) xscr(2) = (gtrk - tortion*y(4))*y(6) + $ hs*(eonps*efld(2)/cc + bfld(1) - y(4)*bfld(3)) xscr(3) = hs*(y(4)*efld(1) + y(5)*efld(2) + efld(3))/cc dyds(6) = eonps*xscr(3) - y(4)*xscr(1) - y(5)*xscr(2) dyds(4) = (xscr(1)-y(4)*dyds(6))/y(6) dyds(5) = (xscr(2)-y(5)*dyds(6))/y(6) ! end if ! end of test on curved ! ! ......................................................... if( spintrk.gt.0 .and. ABS(iptyp).eq.2 ) then call DCROSS( y(7),axbmt,dyds(7) ) end if ! ! ENVELOPE EQUATIONS ADDITION if (run_env .and. (ievt.eq.1)) then call ENV_DERIV(y,dyds) end if ! 900 continue return end c ********************************************************* function DERIVTAB(xa,f,npts,jt,x,m) ! ! Compute derivative numerically from tabulated function f ! ! The maximum value of JT must be NPTS-1 ! implicit double precision (a-h,o-z) dimension xa(npts),f(npts) ! dx = xa(jt+1)-xa(jt) ! if( m .eq. 1 ) then ! first derivative if( jt .eq. 1 )then djt = (f(2) - f(1)) / dx djt1 = (f(jt+2) - f(jt)) / (2. * dx) else if( jt .eq. npts-1 ) then djt = (f(npts) - f(npts-2)) / (2.*dx) djt1 = (f(npts) - f(npts-1)) / dx else djt = (f(jt+1) - f(jt-1)) / (2. * dx) djt1 = (f(jt+2) - f(jt)) / (2. * dx) end if ! else ! second derivative if( jt .eq. 1) then djt1 = (f(jt+2) - 2.*f(jt+1) + f(jt)) / (dx * dx) djt = djt1 else if( jt .eq. npts-1 ) then djt = (f(jt+1) - 2.*f(jt) + f(jt-1)) / (dx * dx) djt1 = djt else djt = (f(jt+1) - 2.*f(jt) + f(jt-1)) / (dx * dx) djt1 = (f(jt+2) - 2.*f(jt+1) + f(jt)) / (dx * dx) end if end if ! derivtab = djt + (djt1-djt) * (x-xa(jt)) / dx ! end c ********************************************************* SUBROUTINE DIPROT(X,PX,Y,PY,T,E,L,Q,M,B,TH,RC) C TRACK FROM ONE PLANE TO ANOTHER IN A UNIFORM DIPOLE FIELD C X,PX,Y,PY,T,E : (I/O) PHASE SPACE COORDINATES C L : (O) ARC LENGTH TRAVERSED C Q : (I) PARTICLE CHARGE C M : (I) PARTICLE MASS C B : (I) FIELD C TH : (I) ANGLE BETWEEN PLANES C RC : (O) RETURN CODE C COPYRIGHT (C) 2010 J. SCOTT BERG IMPLICIT NONE INCLUDE 'icommon.inc' DOUBLE PRECISION X,PX,Y,PY,T,E,L,Q,M,B,TH INTEGER RC DOUBLE PRECISION P,PH,PS,PHI0,DPHI,H,DDPHI,DPHI0,DDPHI0,SINTH,PS1 DOUBLE PRECISION DL INTRINSIC ABS,ASIN,SIN,SQRT P = E*E-M*M PH = P-PY*PY PS = PH-PX*PX IF (PS.LE.0D0) THEN RC = -1 RETURN ENDIF P = SQRT(P) PH = SQRT(PH) PS = SQRT(PS) H = Q*B/PH PHI0 = ASIN(PX/PH) SINTH = SIN(TH) IF (ABS(SIN(PHI0+TH)-H*X*SINTH).GT.1D0) THEN RC = -2 RETURN ENDIF DPHI = PHI0+TH-ASIN(SIN(PHI0+TH)-H*X*SINTH) C NEWTON ITERATE TO GET DPHI. THIS IS NECESSARY TO AVOID ROUNDOFF ERRORS C FOR SMALL FIELDS DDPHI = ATAN(1D0) 1 DDPHI0 = ABS(DDPHI) DPHI0 = DPHI IF (COS(PHI0+TH-DPHI).EQ.0D0) THEN RC = -3 RETURN ENDIF DDPHI = (H*X*SINTH-2*SIN(0.5D0*DPHI)*COS(PHI0+TH-0.5D0*DPHI))/ $ COS(PHI0+TH-DPHI) DPHI = DPHI+DDPHI IF (DPHI.NE.DPHI0.AND.DDPHI.NE.0D0.AND.ABS(DDPHI).LT.DDPHI0) $ GOTO 1 IF (DPHI.NE.DPHI0.AND.DDPHI.NE.0D0) DPHI = DPHI0 C FIND ARC LENGTH IN BEND. FANCY STUFF TO HANDLE CASE WHERE THE BEND FIELD C IS NEARLY ZERO. PS1 = PS*COS(TH-0.5D0*DPHI)-PX*SIN(TH-0.5D0*DPHI) IF (PS1.EQ.0D0) THEN RC = -4 RETURN ENDIF X = X/PS1 DL = 0.5D0*Q*B*X*SINTH IF (ABS(DL).GT.1D0) THEN RC = -5 RETURN ENDIF IF (ABS(DL).LT.SQRT(1.2D1*DBLEPS)) THEN DL = 1D0 ELSE DL = ASIN(DL)/DL ENDIF DL = DL*X*SINTH X = X*(PS*COS(0.5*DPHI)+PX*SIN(0.5*DPHI)) PX = PX*COS(TH-DPHI)+PS*SIN(TH-DPHI) Y = Y + PY*DL T = T + E*DL/CC L = L + DL*P RC = 0 RETURN END c ********************************************************* SUBROUTINE DIPEDG(X,PX,Y,PY,T,E,Q,M,DB,RC) C COPYRIGHT (C) 2010 J. SCOTT BERG IMPLICIT NONE DOUBLE PRECISION X,PX,Y,PY,T,E,DB,Q,M INTEGER RC DOUBLE PRECISION ZI(6),ZM(6),V0(6),V1(6),ERR(6),EE INTEGER I LOGICAL DONE EXTERNAL DIPVEC ZI(1) = X ZI(2) = PX ZI(3) = Y ZI(4) = PY ZI(5) = T ZI(6) = E CALL DIPVEC(ZI(1),ZI(2),ZI(3),ZI(4),ZI(5),ZI(6),Q,M,DB,V0,RC) IF (RC.NE.0) THEN RC = RC - 10 RETURN ENDIF DO 10000 I=1,6 ZM(I) = ZI(I) + V0(I) 10000 CONTINUE CALL DIPVEC(ZM(1),ZM(2),ZM(3),ZM(4),ZM(5),ZM(6),Q,M,DB,V1,RC) IF (RC.NE.0) THEN RC = RC - 20 RETURN ENDIF DO 10010 I=1,6 ZM(I) = ZI(I) + 0.25D0*(V0(I)+V1(I)) ERR(I) = 0.25D0*ABS(V1(I)-V0(I)) 10010 CONTINUE 10020 CALL DIPVEC(ZM(1),ZM(2),ZM(3),ZM(4),ZM(5),ZM(6),Q,M,DB,V0,RC) IF (RC.NE.0) THEN RC = RC - 30 RETURN ENDIF DONE = .TRUE. DO 10030 I=1,6 EE = ABS(ZM(I)-ZI(I)-0.5D0*V0(I)) DONE = DONE.AND.(EE.EQ.0D0.OR.EE.GE.ERR(I)) IF (EE.NE.0D0.AND.EE.LT.ERR(I)) ERR(I) = EE ZM(I) = ZI(I) + 0.5D0*V0(I) 10030 CONTINUE IF (.NOT.DONE) GOTO 10020 X = X + V0(1) PX = PX + V0(2) Y = Y + V0(3) PY = PY + V0(4) T = T + V0(5) E = E + V0(6) RETURN END c ********************************************************* SUBROUTINE DIPVEC(X,PX,Y,PY,T,E,Q,M,DB,V,RC) C COPYRIGHT (C) 2010 J. SCOTT BERG IMPLICIT NONE INCLUDE 'icommon.inc' DOUBLE PRECISION X,PX,Y,PY,T,E,M,Q,DB,V(6) INTEGER RC DOUBLE PRECISION PH2,PS2,C0,C1 ! x = x ! suppress compiler warnings for unused variables t = t ! PH2 = E*E-M*M-PY*PY PS2 = PH2-PX*PX IF (PS2.LT.0D0) THEN RC = -1 RETURN ENDIF C0 = Y*Q*DB/SQRT(PS2) C1 = 0.5D0*Y*C0/PS2 V(1) = C1*PH2 V(2) = 0D0 V(3) = C1*PX*PY V(4) = -C0*PX V(5) = C1*PX*E/CC V(6) = 0D0 RC = 0 RETURN END c ********************************************************* subroutine DISP_COORD(d,fi) ! ! Displace particle coordinates by d in the direction fi ! implicit double precision (a-h,o-z) include 'icommon.inc' ! dimension x(2) ! ! d=d*ran1(idum) ! tild angle ! fi=fi*ran1(idum) ! direction of rotation respect to axis x in plane (x-y) cosfi = COS(fi) sinfi = SIN(fi) do i=1,2 x(i) = xp(i) end do ! xp(1)=x(1)-d*cosfi xp(2)=x(2)-d*sinfi ! return end c ********************************************************* real*8 function DOT3(a,b) ! ! Double precision vector dot product ! implicit double precision (a-h,o-z) dimension a(3),b(3) ! dot3 = 0.d0 ! do i=1,3 dot3 = dot3 + a(i)*b(i) end do ! return end c ********************************************************* subroutine END_OF_REGION(iscell) ! ! Handle tasks required at the end of a region ! implicit double precision (a-h,o-z) include 'icommon.inc' integer k,rc,EDGE_ACCEL_13 DOUBLE PRECISION E,PS2,TH,L double precision r,bs,pr,pt logical iscell ! if (.not.iscell) then r = sqrt(xp(1)*xp(1)+xp(2)*xp(2)) k = 0 100 k = k + 1 if (k.le.nrreg .and. (r.lt.rlow(k) .or. r.gt.rhigh(k))) goto 100 end if ! if (iscell .or. k.le.nrreg) then if( (iscell .and. cftag.eq.'DIP' .and. cfparm(1).eq.3 .and. $ cfparm(2).lt.0. .and. cfparm(6).eq.1 ) .or. $ (.not.iscell .and. ftag(k).eq.'DIP' .and. & fparm(1,k).eq.3 .and. fparm(2,k).lt.0. .and. & fparm(6,k).eq.1) ) then ! radial sector FFAG reverse bend field xp(1) = -xp(1) ! restore x pp(1) = -pp(1) if (rtuple) call OUTPUT else if ( .not.iscell .and. ftag(k).eq.'ACCE' .and. $ fparm(1,k).eq.13 .and. $ (fparm(5,k).eq.0 .or. fparm(5,k).eq.1) ) then rc = EDGE_ACCEL_13(xp,pp,tp,pcharge/ee,pmass, $ twopi*fparm(2,k)*1d6,phase,-1d-3*fparm(3,k)) if (rc.ne.0) iflag = -86 if (rtuple) call OUTPUT ! else if (.not.iscell.and.ftag(k).eq.'SOL'.and.fparm(1,k).eq.8) $ then pr = pp(1)*pp(1)+pp(2)*pp(2)+pp(3)*pp(3) if (fparm(3,k).eq.0.or.fparm(3,k).eq.1) then bs = 0.5d0*pcharge*fparm(2,k) pp(1) = pp(1) - bs*xp(2) pp(2) = pp(2) + bs*xp(1) end if bs = fparm(5,k)*0.25d0*pcharge*pcharge/sqrt(pr) pp(1) = pp(1) + bs*xp(1) pp(2) = pp(2) + bs*xp(2) tp = tp - 0.5d0*bs*r*r*sqrt(pr+pmass*pmass)/(cc*pr) pt = pr-pp(1)*pp(1)-pp(2)*pp(2) if (pt.lt.0d0) then iflag = -86 else pp(3) = sqrt(pt) end if if (rtuple) call OUTPUT ! ELSE IF (.NOT.ISCELL.AND.FTAG(K).EQ.'DIP'.AND.FPARM(1,K).EQ.5) $ THEN E = SQRT(PP(1)**2+PP(2)**2+PP(3)**2+PMASS**2) TH = TWOPI*FPARM(5,K)/3.6D2 C BASIC PHYSICS: UNBEND TO THE REAL EDGE, DO THE VERTICAL EDGE FOCUSING, C THEN DRIFT TO THE ENTRANCE OF A SECTOR MAGNET CALL DIPROT(XP(1),PP(1),XP(2),PP(2),TP,E,L, $ PCHARGE,PMASS,FPARM(2,K),-TH,RC) IF (RC.LT.0) THEN IFLAG = -86 RETURN ENDIF IF (FPARM(6,K).EQ.0.OR.FPARM(6,K).EQ.1) THEN CALL DIPEDG(XP(1),PP(1),XP(2),PP(2),TP,E, $ PCHARGE,PMASS,-FPARM(2,K),RC) IF (RC.LT.0) THEN IFLAG = -86 RETURN ENDIF ENDIF CALL DIPROT(XP(1),PP(1),XP(2),PP(2),TP,E,L, $ PCHARGE,PMASS,0D0,TH,RC) IF (RC.LT.0) THEN IFLAG = -86 RETURN ENDIF PS2 = E**2-PMASS**2-PP(1)**2-PP(2)**2 IF (PS2.LE.0D0) THEN IFLAG = -86 ELSE PP(3) = SQRT(PS2) ENDIF IF (RTUPLE) CALL OUTPUT ! ELSE IF (.NOT.ISCELL.AND.FTAG(K).EQ.'DIP'.AND.FPARM(1,K).EQ.5) $ THEN E = SQRT(PP(1)**2+PP(2)**2+PP(3)**2+PMASS**2) TH = TWOPI*FPARM(4,K)/3.6D2 C BASIC PHYSICS: DRIFT TO THE REAL EDGE, DO THE VERTICAL EDGE FOCUSING, C THEN TRACK IN A UNIFORM FIELD BACK TO THE ENTRANCE OF A SECTOR MAGNET CALL DIPROT(XP(1),PP(1),XP(2),PP(2),TP,E,L, $ PCHARGE,PMASS,0D0,TH,RC) IF (RC.LT.0) THEN IFLAG = -86 RETURN ENDIF IF (FPARM(6,K).EQ.0.OR.FPARM(6,K).EQ.2) THEN CALL DIPEDG(XP(1),PP(1),XP(2),PP(2),TP,E, $ PCHARGE,PMASS,FPARM(2,K),RC) IF (RC.LT.0) THEN IFLAG = -86 RETURN ENDIF ENDIF CALL DIPROT(XP(1),PP(1),XP(2),PP(2),TP,E,L, $ PCHARGE,PMASS,FPARM(2,K),-TH,RC) IF (RC.LT.0) THEN IFLAG = -86 RETURN ENDIF PS2 = E**2-PMASS**2-PP(1)**2-PP(2)**2 IF (PS2.LE.0D0) THEN IFLAG = -86 ELSE PP(3) = SQRT(PS2) ENDIF IF (RTUPLE) CALL OUTPUT ! else if(.not.iscell .and. nbffag.eq.1d0) then by = fparm(2,k) r0 = fparm(5,k) th2 = slen / r0 dsffag = 2d0 * xp(1) * SIN(0.5d0*th2) call VMAG(pp,pm) dpx = pcharge * by * dsffag pp(1) = pp(1) + dpx pp(3) = SQRT( pm*pm - pp(1)**2 - pp(2)**2 ) end if end if ! end c ********************************************************* subroutine ENV_DERIV(y,dydx) c c envelope dynamic equations c implicit double precision (a-h,o-z) include 'icommon.inc' c dimension y(15),dydx(15),ctau(5) real*8 msdiff,kappa data ctau/0., 658.654, 7.8045, 3.713, 0./ !c ps = y(6) bg = ps / pmass gamma = DSQRT(1+bg**2) bet = bg / gamma vs = bet * cc if (mtag(krad) .ne. 'VAC') then c Determine mean value of energy loss c call ENV_ELOSS(ps,dedxm) dps = 0.001 ps2 = ps + dps call ENV_ELOSS(ps2,dedx2) deslope = (dedx2-dedxm)/(dps*bet) ! ! multiple scatter term, pD/2mc ! uses Gaussian fit; tail contrib unimportant? msdiff = (100. /(2* ps * pmass * mparm(4,krad))) & * (0.0136 * cc / vs)**2 else deslope = 0 dedxm = 0 msdiff = 0 end if ! ! focussing strength kappa kappa = pcharge * bfld(3) / (2 * ps) ! ! rest of the derivatives t_emit = y(7) t_beta = y(8) t_alpha = y(9) t_angm = y(10) t_N = y(11) t_gauss = y(12) t_eL = y(13) c replace dPz/dz with Ez + mean energy loss terms dydx(6) = pcharge * efld(3) / vs - cc * dedxm / vs dydx(7) = t_beta * msdiff - t_emit * cc * dedxm / (vs * ps) dydx(8) = - 2 * t_alpha + t_beta * pcharge * efld(3) / (vs * ps) & - (t_beta**2 / t_emit) * msdiff dydx(9) = -(1+t_alpha**2 + (t_beta * kappa - t_angm)**2)/t_beta & + 2 * kappa * (t_beta * kappa - t_angm) & - (t_alpha * t_beta / t_emit) * msdiff dydx(10) = t_beta * kappa * cc * dedxm / (vs * ps) & - (t_angm * t_beta / t_emit) * msdiff if (env_decay) then dydx(11) = - t_N / (ctau(ABS(iptyp)) * bg) else dydx(11) = 0 end if dydx(12) = (6. * t_beta / t_emit) * msdiff dydx(13) = - t_eL * deslope c return end c ********************************************************* subroutine ENV_ELOSS(pm,dedxm) c c calculates avg energy loss dedxm -- same as DEDX c implicit double precision (a-h,o-z) include 'icommon.inc' c real*8 mel,ip_mat ! data bbcon/5.099102e-32/ ! [GeV m^2] ! mel = mass(1) e = SQRT(pm**2 + pmass**2) gam = e / pmass bet = pm / pmass / gam bg = bet * gam bet2 = bet * bet zm = mparm(1,krad) am = aelcomp(1,krad) dm = mparm(2,krad) densel = zm * densat ! [# cm^-3] c1 = bbcon * densel * 1.e6 ! [GeV/m] c3 = 0.5 * c1 ! [GeV/m] ip_mat = ipelcomp(1,krad) x0 = mparm(5,krad) x1 = mparm(6,krad) sa = mparm(7,krad) sc = mparm(8,krad) sm = mparm(9,krad) ! c Determine mean value of energy loss c dens = 0. if( delev .gt. 1 ) then ! density effect correction if( bg .le. 0. ) then iflag = -25 go to 900 end if xs = LOG10( bg ) if( xs.ge.x0 .and. xs.le.x1 )then dens = 4.6052 * xs + sa * (x1 - xs)**sm + sc else if( xs .gt. x1 )then dens = 4.6052 * xs + sc end if end if ! tmax = 2.*mel*bg*bg/(1.+2.*gam*mel/pmass+(mel/pmass)**2) ! Bethe-Bloch mean loss if( delev.eq.1 .or. delev.eq.2 ) then arg = 2.*mel*bg**2*tmax*1e18 / ip_mat**2 if( arg .le. 0. ) then iflag = -25 go to 900 end if dedxm = c1/bet**2 * ( 0.5*LOG(arg) - bet**2 - dens/2. ) c demean = dedxm * hstp end if ! ! Restricted Bethe-Bloch loss ! cf. Review of Particle Properties, 1998, eq. 23.2,6 if( delev .eq. 3 ) then if( ABS(iptyp) .eq. 1 ) then tcut = dcute else tcut = dcutm end if tupper = MIN(tcut,tmax) ! arg = 2.*mel*bg**2*tupper*1e18 / ip_mat**2 if( arg .le. 0. ) then iflag = -25 go to 900 end if dedxm = c1 / bet**2 * ( 0.5*LOG(arg) - & 0.5*(1.+tupper/tmax)*bet**2 - dens/2. ) c demean = dedxm * hstp end if c 900 continue return end c ********************************************************* subroutine ENV_MAKE_BEAM c c reset ievt=1 particle, on axis c implicit double precision (a-h,o-z) include 'icommon.inc' c ievt = 1 ipnum = 1 ipflg = 0 evtwt = 1. xp(1) = 0. xp(2) = 0. xp(3) = 0. pp(1) = 0. pp(2) = 0. pp(3) = p1bt(3,1) iptyp = bmtype(1) pmass = mass( ABS(iptyp) ) pcharge = SIGN( ee, DBLE(iptyp) ) * chrg( ABS(iptyp) ) bg = pp(3) / pmass gamma = DSQRT(1+bg**2) bet = bg / gamma vs = bet * cc tp = x1bt(3,1) / vs ! initialize envelope parameters env_alpha = 0. env_angm = 0. env_beta = p1bt(3,1) * x2bt(1,1) / p2bt(1,1) env_emit = x2bt(1,1) * p2bt(1,1) / pmass env_pz = p1bt(3,1) ! env_N = 1.0 env_gauss = 100. env_eL = x2bt(3,1) * p2bt(3,1) / pmass do i=1,nbcorr(1) ic = corrtyp(i,1) if( ic .eq. 1 ) then ! angular momentum Bc = corr1(i,1) env_angm = 0.5 * pcharge * Bc * x2bt(1,1) / p2bt(1,1) end if end do c return end c ********************************************************* subroutine ENV_SCRAPE c implicit double precision (a-h,o-z) include 'icommon.inc' c f_old = 1 - (1 + env_gauss/2.)*EXP(-env_gauss/2.) g_old = 1 - (1 + env_gauss/2. + env_gauss**2/8.) & * EXP(-env_gauss/2.) amp_sigma = env_emit * f_old / g_old amp_max = amp_sigma * env_gauss r_max = rhigh(nrreg) amp_cut = (env_pz/pmass) * r_max**2 / env_beta ! make amplitude cut larger to allow for "elliptical" orbits amp_cut = amp_cut * (SQRT(3.) + amp_cut/(8.*amp_sigma)) & / (1. + amp_cut/(8.*amp_sigma)) if (amp_cut .lt. amp_max) then env_gauss = env_gauss * amp_cut / amp_max f_new = 1 - (1 + env_gauss/2.)*EXP(-env_gauss/2.) g_new = 1 - (1 + env_gauss/2. + env_gauss**2/8.) & * EXP(-env_gauss/2.) env_emit = amp_sigma / (f_new/g_new) env_N = env_N * f_new/f_old end if c return end c ********************************************************* subroutine ENV_SETUP(title,ioc) c c predetermined options for envelope code c implicit double precision (a-h,o-z) include 'icommon.inc' c integer ioc character*80 title c ! beam moments output file header open(unit=11, & file=pathname(1:LASTNB(pathname)) // 'for011.dat', & status='unknown',iostat=ioc) write(11,'(a1,a79)') '#',title write(11,200) 200 format(t1,'regn',t10,'z',t24,'Bz',t38,'t',t52,'pz',t66,'eN', & t80,'beta',t94,'alpha',t108,'ang mom',t122,'N',t136,'gauss', & t150,'eL') c bgen = .true. if( npart .gt. 0 ) then env_n = npart else env_n = 1 end if npart=1 spintrk=0 steprk = .true. c return end c ********************************************************* subroutine FIELD(xyz,tt) c c Return EM field at location xyz and time tt c implicit double precision (a-h,o-z) include 'icommon.inc' c dimension xyz(3),e(3),b(3) ! do i=1,3 ! initialize; also zero field for drift space efld(i) = 0.d0 bfld(i) = 0.d0 end do ! call FIELD_CELL(xyz,tt,e,b) do i=1,3 efld(i) = efld(i) + e(i) bfld(i) = bfld(i) + b(i) end do ! call FIELD_BKG(xyz,tt,e,b) do i=1,3 efld(i) = efld(i) + e(i) bfld(i) = bfld(i) + b(i) end do ! call FIELD_REG(xyz,tt,e,b) do i=1,3 efld(i) = efld(i) + e(i) bfld(i) = bfld(i) + b(i) end do ! end c ********************************************************* subroutine FIELDFNC(tag,fpar,x,t,e,b) ! ! Select appropriate field function ! ! The position x must be given in terms of magnet local system ! implicit double precision (a-h,o-z) dimension x(3),e(3),b(3),fpar(15) character*4 tag ! if( tag .eq. 'ACCE' ) then ! accelerator call ACCEL(fpar,x,t,e,b) ! else if( tag .eq. 'BACK' ) then ! previously defined background call BACKGET(x,t,e,b) ! else if( tag .eq. 'BLOC' ) then ! solenoidal current block call BLOCK(fpar,x,e,b) ! else if( tag .eq. 'BROD' ) then ! bent axial current rod call BROD(fpar,x,e,b) ! else if( tag .eq. 'BSOL' ) then ! bent solenoid call BENTSOL(fpar,x,e,b) ! else if( tag .eq. 'COIL' ) then ! circular coil call COIL(fpar,x,e,b) ! else if( tag .eq. 'DIP' ) then ! vertical sector dipole call DIPOLE(fpar,x,e,b) ! else if( tag .eq. 'EFLD' ) then ! electric field region call EFIELD(fpar,x,e,b) ! else if( tag .eq. 'FOFO' ) then ! solenoid FOFO lattice call FOFO(fpar,x,e,b) ! else if ( tag .eq. 'HDIP' ) then ! horizontal sector dipole call DIPOLE_H(fpar,x,e,b) ! else if ( tag .eq. 'HELI' ) then ! helix call HELIX(fpar,x,e,b) ! else if ( tag .eq. 'HORN' ) then ! horn call HORN(fpar,x,e,b) ! else if( tag .eq. 'KICK' ) then ! kicker call KICKER(fpar,x,t,e,b) ! else if ( tag .eq. 'NONE' ) then ! no field do i=1,3 e(i) = 0.d0 b(i) = 0.d0 end do ! else if ( tag .eq. 'QUAD' ) then ! quad call QUADRUPOLE(fpar,x,e,b) c else if( tag .eq. 'ROD' ) then ! axial current rod call ROD(fpar,x,e,b) ! else if ( tag .eq. 'SEX' ) then ! sextupole call SEXTUPOLE(fpar,x,e,b) c else if( tag .eq. 'SHEE' ) then ! solenoidal sheet call SHEET(fpar,x,e,b) c else if( tag .eq. 'SOL' ) then ! solenoid call SOLENOID(fpar,x,e,b) c else if( tag .eq. 'SQUA' ) then ! skew quadrupole call QUAD_SKEW(fpar,x,e,b) c else if( tag .eq. 'STUS' ) then ! static 3D user field call STUS(fpar,x,e,b) c else if( tag .eq. 'WIG' ) then ! tapered solenoid call WIGGLER(fpar,x,e,b) c end if ! end of test for field tags ! return end ! ********************************************************* subroutine FIELD_BKG(xyz,tt,e,b) c c Return background EM field at location xyz and time tt c implicit double precision (a-h,o-z) include 'icommon.inc' c dimension xyz(3),x(3),e(3),b(3) ! do i=1,3 e(i) = 0.d0 b(i) = 0.d0 x(i) = xyz(i) end do ! if( bftag .ne. 'NONE' ) then kbbkg = 1 kbcell = 0 x(3) = xyz(3) - sbg0 ! transform to magnet coordinate system x(3) = MAX( x(3),0.d0 ) ! call BACKGET(x,tt,e,b) ! end if ! end ! ********************************************************* subroutine FIELD_CELL(xyz,tt,e,b) c c Return cell EM field at location xyz and time tt c implicit double precision (a-h,o-z) include 'icommon.inc' c dimension xyz(3),x(3),e(3),b(3) ! do i=1,3 e(i) = 0.d0 b(i) = 0.d0 x(i) = xyz(i) end do ! if( cftag .ne. 'NONE' ) then kbcell = 1 ! inform field routines that this is a cell field kbbkg = 0 x(3) = xyz(3) - sc0 ! transform to magnet coordinate system x(3) = MAX( x(3),0.d0 ) ! call FIELDFNC(cftag,cfparm,x,tt,e,b) ! do i=1,3 ! check symmetry of cell fields b(i) = cellslp * b(i) end do ! if( cftag.eq.'DIP' .or. cftag.eq.'BSOL' .or. & cftag.eq.'WIG' ) then htrk = cellslp * htrk hptrk = cellslp * hptrk end if ! end if ! end ! ********************************************************* subroutine FIELD_REG(xyz,tt,e,b) c c Return region EM field at location xyz and time tt c implicit double precision (a-h,o-z) include 'icommon.inc' c dimension xyz(3),x(3),e(3),b(3),e1(3),b1(3),fpar(15) character*4 tag ! do i=1,3 e(i) = 0.d0 b(i) = 0.d0 x(i) = xyz(i) end do ! if( ft .ne. 'NONE' ) then kbcell = 0 kbbkg = 0 x(3) = xyz(3) - sr0 ! transform to magnet coordinate system x(3) = MAX( x(3),0.d0 ) ! call FIELDFNC(ft,fparm(1,krad),x,tt,e,b) ! end if ! ! Neighbor region fields if( neighbor ) then do j=1,nregnb if( j7rec .eq. i7regnb(j) ) then if( j .gt. 1 ) then ! get previous region tag = tagnb(j-1) do i=1,15 fpar(i) = fparnb(i,j-1) end do x(3) = xyz(3) - sr0nb(j-1) ! transform to magnet system call FIELDFNC(tag,fpar,x,tt,e1,b1) do i=1,3 e(i) = e(i) + e1(i) b(i) = b(i) + b1(i) end do end if ! if( j .lt. nregnb ) then ! get following region tag = tagnb(j+1) do i=1,15 fpar(i) = fparnb(i,j+1) end do x(3) = xyz(3) - sr0nb(j+1) ! transform to magnet system call FIELDFNC(tag,fpar,x,tt,e1,b1) do i=1,3 e(i) = e(i) + e1(i) b(i) = b(i) + b1(i) end do end if end if end do end if ! end c ********************************************************* subroutine FINISH c c Take care of end of simulation tasks c implicit double precision (a-h,o-z) include 'icommon.inc' c integer*2 iyr,imon,iday,ihr,imin,isec ! if( ffcr ) call FORMFEED ! call OUT_EMIT ! emittances if( spin ) call OUT_POL ! polarization call OUT_STATS ! statistics call OUT_HIST ! histograms call OUT_SCAT ! scatterplots call OUT_ZHIST ! z-histories call OUT_RHIST ! r-histoties c ! write(2,50) ndeltas,nflips !50 format('NDELTAS,NFLIPS: ',2i10) ! write(2,*)' ICOOL finished' call ICOOLDATE(iyr,imon,iday) call ICOOLTIME(ihr,imin,isec) t0 = ihr0*3600 + imin0*60 + isec0 t = ihr*3600 + imin*60 + isec if(iday .ne. iday0) t = t + (iday-iday0)*86400 dt = t - t0 nh = dt / 3600. nm = (dt - 3600.*nh) / 60. ns = dt - 3600.*nh - 60.*nm if( termout ) write(*,85)nh,nm,ns 85 format(' Total elapsed time = ',i4,' hours,',i4,' minutes,', + i4,' seconds') write(2,85)nh,nm,ns c close(1) close(2) close(3) close(4) close(7) close(8) close(9) ! ! ENVELOPE EQUATIONS ADDITION if (run_env) close(11) ! if( rfdiag .gt. 19 ) close(rfdiag) if( rfphase .gt. 19 ) close(rfphase) c if( beep ) then call SOUND(400,200) call SOUND(4,200) call SOUND(400,200) call SOUND(4,200) call SOUND(400,200) call SOUND(4,200) call SOUND(200,600) end if c return end ! ********************************************************* subroutine FPSIH(x,f,df) ! ! Define transcendental equation and derivative for determining ! phase of the head of the alpha bucket ! implicit double precision (a-h,o-z) include 'icommon.inc' ! sps = SIN(psis) cps = COS(psis) ! f = COS(x) + x*sps - (pi-psis)*sps + cps df = sps - SIN(x) ! return end c ********************************************************* subroutine FSTEP(y,neqn,x) c ! Propagate particle one step (hdid) thru field ! Recommend new step size (hnextf) based on field variation c implicit double precision (a-h,o-z) include 'icommon.inc' c parameter (grow = -0.20d0, shrink = -0.25d0, fcor = 1.d0/15.d0, + safety = 0.9d0, errcon = 1.89d-4 , nmax = 15) ! dimension y(neqn),dydx(nmax),yscal(nmax),yt(nmax) dimension ysav(nmax),dysav(nmax) ! h = hdid call DERIV(x,y,dydx) if( iflag .ne. 0 ) go to 900 ! xsav = x ! save initial values ! do i=1,neqn ysav(i) = y(i) dysav(i) = dydx(i) yscal(i) = ABS( y(i) ) + ABS( hdid*dydx(i) ) if( yscal(i) .eq. 0.d0 ) yscal(i) = 1.d0 end do ! if( varstep ) then c Take two half-steps .................................. hh = 0.5d0 * h call RKSTEP(ysav,dysav,neqn,xsav,hh,yt) if( iflag .ne. 0 ) go to 900 x = xsav + hh call DERIV(x,yt,dydx) if( iflag .ne. 0 ) go to 900 call RKSTEP(yt,dydx,neqn,x,hh,y) if( iflag .ne. 0 ) go to 900 x = xsav + h c c Take one large step .................................. call RKSTEP(ysav,dysav,neqn,xsav,h,yt) if( iflag .ne. 0 ) go to 900 c ......................................................... errmax = 0.d0 c do i=1,neqn yt(i) = y(i) - yt(i) ! yt holds error errmax = MAX( errmax, ABS( yt(i)/yscal(i) ) ) end do c errmax = errmax / epsreq ! scale relative to required tolerance ! if( errmax .gt. 1.d0 ) then ! try smaller steps hold = h h = safety * h * (errmax**shrink) if( h .lt. 0.1d0*hold ) then h = 0.1d0 * hold ! no smaller than factor of 10 end if hnextf = h ! else ! this step OK; try bigger one next time ! if( errmax .gt. errcon ) then hnextf = safety * h * (errmax**grow) else hnextf = 5.d0 * h ! no bigger than factor of 5 end if end if c do i=1,neqn ! correct 5th order truncation error y(i) = y(i) + yt(i) * fcor end do ! else ! fixed step size ....................... ! call RKSTEP(ysav,dysav,neqn,xsav,h,y) x = xsav + h ! end if ! 900 continue end c ********************************************************* function GAUSS(xm,sd) c c Return a random value from a gaussian distribution c with mean xm and standard deviation sd c uses Numerical Recipe function GASDEV c c get random value from gaussian distribution with mean 0 c and standard deviation 1 ! implicit double precision (a-h,o-z) include 'icommon.inc' ! gauss = GASDEV(irnarg) c c renormalize to requested mean and standard deviation gauss = sd*gauss + xm c return end c ********************************************************* subroutine GO_REGION c c Propagate particle thru specific s-region c implicit double precision (a-h,o-z) include 'icommon.inc' c ! logical first dimension y(15),v(3),bet(3),exv(3),v1(3),bl(3),bt(3),pol1(3) c data a/1.165923d-3/ ! muon magnetic moment anomaly ! data a/0./ ! test muon anomaly <<<<<<<<<<<<<<<<<<< ! c The units used are: t[s] x[m] v[m/s] E[GV/m] B[T] c p[Gev/c] m[GeV/c^2] E[GeV] c call VMAG(pp,pm) gamma = SQRT( 1.d0 + pm*pm/pmass/pmass ) vz = pp(3) / pmass / gamma * cc if( vz .le. 0.d0 ) then iflag = -32 go to 900 end if ! call BEG_OF_REGION(.false.) if (iflag.ne.0) go to 900 ! ! It defeats the step size adjustment if you put ! this inside the fine step loop below! ! ft = ftag(1) hregn = slen / 10.d0 hnextf = zstep hnextde = zstep hnextms = zstep jknt = 0 ! first = .false. ! if( varstep ) then ! adaptive step size maxsteps = MAX( 100.*slen/zstep,1. ) sacc = xp(3) ! first = .true. ! else ! fixed step size if( scalestep*zstep .le. slen ) then maxsteps = NINT( slen / (scalestep*zstep) ) else maxsteps = 1 end if fixstep = slen / DBLE(maxsteps) z0 = xpint(3) if( z0 .gt. sr0 ) then nskip = MIN(NINT(z0/fixstep),maxsteps) else nskip = 0 end if end if ! ihwin = 1 sft = cftag krold = -1 ! flag when we change radial regions ! ! --------------------------------------------------------- ! do 500 iz=1,maxsteps ! loop over axial steps ! do 450 ir=1,nrreg ! loop over radial regions ! if( mgeom(1) .eq. 'HWIN' ) then ! hemi window patch ft = ftag(ihwin) mt = 'xxx' ! IN_HEMIWIN will set this later mgt = mgeom(1) krad = 1 go to 150 end if ! rp = SQRT( xp(1)**2 + xp(2)**2 ) if( rp .lt. rlow(ir) .or. rp .gt. rhigh(ir) ) go to 450 ! if( ir .ne. krold ) then ! new radial region hregn = MIN( slen/10.d0,(rhigh(ir)-rlow(ir))/10.d0 ) krad = ir mt = mtag(ir) ft = ftag(ir) mgt = mgeom(ir) ! if( mt .ne. 'VAC') call LOAD_MATERIAL ! krold = ir end if ! 150 continue ! enter here for patches that skip normal start ! ! Check each step for density profile adjustments if( mt.ne.'VAC' .and. dirdenp.ne.-1 ) then call LOAD_MATERIAL end if ! ! Save current particle info (beginning of step) ......... do i=1,3 xp0(i) = xp(i) pp0(i) = pp(i) pol0(i) = pol(i) end do ! tp0 = tp call VMAG( pp,pm0 ) ! save initial magnitude ! if( spintrk.gt.0 .and. ABS(iptyp).eq.2 ) then ! save current BMT direction gamma = SQRT(1.d0+pm0*pm0/pmass/pmass) do i=1,3 v(i) = pp(i) / pmass / gamma * cc end do call VMAG(v,vm) do i=1,3 v1(i) = v(i) / vm ! unit velocity vector bet(i) = v(i) / cc end do ! call FIELD( xp(1),tp ) call DCROSS( efld,bet,exv ) term = DOT3( bfld,v1 ) do i=1,3 bl(i) = term * v1(i) bt(i) = bfld(i) - bl(i) end do ! term0 = DOT3(bfld,v1) do i=1,3 term = (a*gamma+1.d0)*bfld(i) term2 = -term0 * (gamma-1.d0)*a*v1(i) term3 = gamma*( a+1.d0/(gamma+1.d0) )*exv(i)/cc axbmt(i) = pcharge/pp(3) * (term+term2+term3) end do end if ! ......................................................... 170 continue ! enter here if step needs to be redone ! jknt = jknt + 1 if( jknt .gt. maxsteps ) then iflag = -33 go to 900 end if ! ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! Determine next step size ! call STEPSIZE(pm0) ! if( iflag .lt. 0 ) go to 900 ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! if( phantom ) then ! move parallel to axis for field grid do i=1,3 xp(i) = xp0(i) pp(i) = pp0(i) end do xp(3) = xp0(3) + hdid call VMAG(pp,pm) gamma = SQRT( 1.d0 + pm*pm/pmass/pmass ) vz = pp(3) / pmass / gamma * cc tp = tp0 + hdid/vz call FIELD( xp(1),tp ) go to 350 end if ! ! Move particle thru field (if any) if( jfieldon.eq.0 .and. .not.neighbor) then ! drift xp(1) = xp0(1) + pp0(1)/pp0(3)*hdid xp(2) = xp0(2) + pp0(2)/pp0(3)*hdid xp(3) = xp0(3) + hdid ! do i=1,3 pp(i) = pp0(i) bfld(i) = 0.d0 efld(i) = 0.d0 end do ! call VMAG(pp,pm) gamma = SQRT( 1.d0 + pm*pm/pmass/pmass ) vz = pp(3) / pmass / gamma * cc tp = tp0 + hdid/vz ! if( spintrk.gt.0 .and. ABS(iptyp).eq.2 ) then do i=1,3 pol(i) = pol0(i) end do end if ! .... ..................................................... else ! swim thru field ! ! Compute initial values and derivatives for FSTEP x = xp0(3) ! s y(1) = xp0(1) ! x y(2) = xp0(2) ! y y(3) = tp0 ! t y(4) = pp0(1) / pp0(3) ! x' y(5) = pp0(2) / pp0(3) ! y' y(6) = pp0(3) ! Ps ! if( spintrk.gt.0 .and. ABS(iptyp).eq.2 ) then do i=1,3 y(i+6) = pol0(i) end do neqn = 9 else neqn = 6 end if ! ! ENVELOPE EQUATIONS ADDITION if (run_env .and. (ievt.eq.1)) then neqn = 13 y(6) = env_pz y(7) = env_emit y(8) = env_beta y(9) = env_alpha y(10) = env_angm y(11) = env_N y(12) = env_gauss y(13) = env_eL end if ! ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c Take a step thru appropriate field c if( steprk ) then call FSTEP(y,neqn,x) else call BORIS_STEP(y,neqn,x) end if c if( iflag .lt. 0 ) go to 900 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! ! Update particle variables with integration results ! xp(1) = y(1) xp(2) = y(2) xp(3) = x pp(1) = y(6) * y(4) pp(2) = y(6) * y(5) pp(3) = y(6) tp = y(3) ! ! ENVELOPE EQUATIONS ADDITION if (run_env .and. (ievt.eq.1)) then env_pz = y(6) env_emit = y(7) env_beta = y(8) env_alpha = y(9) env_angm = y(10) env_N = y(11) env_gauss = y(12) env_eL = y(13) end if ! ! Dipole vertical focusing from rotated pole faces (NOT end of region) if( edgekick ) then call VMAG(pp,pm) pp(2) = pp(2) + dpyedge arg = pm**2 - pp(1)**2 - pp(2)**2 if( arg .gt. 0. ) then pp(3) = SQRT(arg) else iflag = -11 go to 900 end if edgekick = .false. end if ! ! Watch out for crazy stepping results if( ABS(xp(1)).gt.100. .or. ABS(xp(2)).gt.100. .or. & ABS(pp(1)).gt.1000. .or. ABS(pp(2)).gt.1000. ) then iflag = -76 go to 900 end if ! if( spintrk.gt.0 .and. ABS(iptyp).eq.2 ) then do i=1,3 pol(i) = y(i+6) end do ! Require component parallel to BMT direction be conserved call VMAG(axbmt,bmt) ! convert axbmt to unit vector if( bmt .ne. 0. ) then do i=1,3 axbmt(i) = axbmt(i) / bmt end do else iflag = -62 go to 900 end if ! term = DOT3(pol0,axbmt) if( term .ge. 1.d0 ) then ! spin along bmt axis do i=1,3 pol(i) = axbmt(i) end do go to 300 else if( term .le. -1.d0 ) then ! spin against bmt axis do i=1,3 pol(i) = -axbmt(i) end do go to 300 end if ! General case do i=1,3 pol1(i) = pol(i) - term*axbmt(i) end do ! call VMAG(pol1,r1) if( r1 .gt. 0.d0 ) then term2 = 1. - term**2 if( term2 .gt. 0.d0 ) then do i=1,3 pol(i) = term*axbmt(i) + SQRT(term2)/r1*pol1(i) end do else iflag = -62 go to 900 end if else ! r1 = 0 => pol=pol0 do i=1,3 pol(i) = pol0(i) end do end if ! end of test on r1 ! end if ! end of test on spintrk > 0 and muon 300 continue ! ! Ensure conservation of momentum call VMAG( pp,pm ) if( ft.ne.'ACCE' .and. ft.ne.'KICK' ) then do i=1,3 pp(i) = pp(i) * pm0 / pm end do end if ! end if ! end of test on jfieldon=0 ! 350 continue ! call FIELD( xp(1),tp ) ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! ! Check for material-independent processes ! ! Use average p from the field step for these processes pavg = 0.5d0 * (pm0 + pm) ! if ( ldecay ) call DECAY(pavg) ! particle decay ! if( iflag .lt. 0 ) go to 900 ! ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! c Simulate material-dependent processes c if (run_env .and. (ievt.eq.1)) then call ENV_SCRAPE ! else if ( mt .ne. 'VAC' ) then ! call MSTEP(pavg) ! if( iflag .eq. 1 ) go to 170 ! redo this step if( iflag .lt. 0 ) go to 900 end if end if ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! sarc = sarc + sdid ! if( varstep ) then sacc = sacc + hdid if( sacc .gt. sr0+slen ) go to 550 ! safety valve ! if( sacc .lt. xpint(3) ) then ! go to 500 ! else ! if( first ) then ! do j=1,3 ! xp(j) = xpint(j) ! pp(j) = ppint(j) ! end do ! first = .false. ! end if ! end if ! else ! fixed steps if( nskip .gt. 0 ) then nskip = nskip - 1 if( nskip .eq. 0 ) then ! reset position and momentum do j=1,3 ! keep common time for all particles xp(j) = xpint(j) pp(j) = ppint(j) end do end if go to 500 end if end if ! ! Report particle info if( jpar .le. nprnt ) then if( jknt .le. prnmax ) then if( prlevel .gt. 1 ) then write(2,425)jknt,tp,hdid,hnextf,hnextde,hnextms,xp,pp 425 format('KNT,TP,HDID,HF,HE,HS:',i5,5e10.3,/, & ' XP,PP: ',6e12.4) if( spintrk.gt.0 .and. jpar.gt.0 ) then if( ABS(iptyp) .eq. 2 ) then hlab = HELICITY( pol ) write(2,426)pol,hlab 426 format(' Spin,Hlab:',3f7.3,5x,f8.4) else write(2,427) pol 427 format(' Spin:',3f7.3) end if end if end if c if( prlevel .gt. 2 ) then write(2,430)efld,bfld 430 format(' E,B: ',6e12.4) end if c if( prlevel.gt.3 .and. .not.(xp(1).eq.0..and.xp(2).eq.0.)) & then phi = ATAN2( xp(2),xp(1) ) rp = sqrt( xp(1)**2 + xp(2)**2 ) sph = sin(phi) cph = cos(phi) pr = pp(1)*cph + pp(2)*sph pph = pp(2)*cph - pp(1)*sph br = bfld(1)*cph + bfld(2)*sph bph = bfld(2)*cph - bfld(1)*sph write(2,435) rp,phi*180./pi,pr,pph,br,bph 435 format(' R,PHI,Pr,Pp,Br,Bp:',6e10.3) end if c end if end if ! ! Test if Pz is below minimum allowed for tracking if( pp(3) .le. pzmintrk ) then iflag = -43 go to 900 end if ! ......................................................... go to 480 ! step was successful; go on c 450 continue ! end of loop over radial regions c ......................................................... ! ! Particle radius did not match any radial region iflag = -23 go to 900 c 480 continue ! ! Load diagnostic information if( nzhist.gt.0 .and. jpar.ge.0 ) call FILL_ZHIST ! if( rtuple .and. MOD(jknt,rtuplen).eq.0 ) call OUTPUT ! if( ftag(krad) .eq. 'ACCE' ) then model = fparm(1,krad) kfile = fparm(7,krad) if( model.eq.8 .and. kfile.gt.19 ) then dt = tp - t0_il if( dt .ge. dt_il ) then write(kfile,490) jsrg,jpar,xp(3),tp,efld(3) 490 format(2i5,3e13.5) t0_il = tp end if end if end if ! --------------------------------------------------------- if( varstep ) then ! check if we are at end of region sac = xp(3) - sr0 ! accumulated length along s for this region ds = slen - sac ! remaining distance if( ds .le. epsstep ) go to 550 ! stop tracking end if ! 500 continue ! End of loop over axial steps ! --------------------------------------------------------- ! 550 continue call END_OF_REGION(.false.) ! 900 continue if( jpar .le. nprnt ) then write(2,'(i5,a,i5)')jknt,' steps in region ',jsrg end if ! return end c ********************************************************* subroutine HRDEND(x,y,t,px,py,p,q,m,n,bx,by) c c Compute phase space evolution due to end fields c c x : horizontal position (m) c y : vertical position (m) c t : time (s) c px : Horizontal transverse momentum (GeV/c) c py : Vertical transverse momentum (GeV/c) c p : Total momentum (GeV/c) c q : Charge, in units of e c m : Mass (GeV/c^2) c n : Number of multipoles (1=dipole, 2=quad, etc.) c bx : Array of d^{k-1}B_x/dx^{k-1}, k=1..n (k=1: dipole; k=2: quad) c by : Array of d^{k-1}B_y/dx^{k-1}, k=1..n (T/m^{n-1}) c c Version 030429a c c Report all suspected bugs to J. Scott Berg c or c c $Date: 2003/04/28 19:41:28 $ c $Revision: 1.7 $ c $Source: /home/jsberg/cvsroot/icool/ends/ends.f,v $ c implicit none double precision c parameter (c=2.99792458d8) integer n integer i double precision x,y,t,px,py,p,q,m,bx(n),by(n) double precision u(2,2),tr,ti,tr2,ti2,tt,r2,xf,yf r2 = x*x+y*y ! tr = 0.25d0*q*c/p tr = 0.25d-9*q*c/p ! use p in GeV/c ti = 0d0 xf = tr*by(1)*r2 yf = -tr*bx(1)*r2 u(1,1) = tr*(by(1)*x + bx(1)*y) u(1,2) = -tr*(bx(1)*x - by(1)*y) if (n.gt.1) then u(1,1) = u(1,1) + tr*by(2)*r2 u(1,2) = u(1,2) - tr*bx(2)*r2 endif ti = y*tr tr = x*tr if (n.gt.1) then tr2 = 0.5d0*tr*r2 ti2 = 0.5d0*ti*r2 xf = xf + by(2)*tr2 - bx(2)*ti2 yf = yf - by(2)*ti2 - bx(2)*tr2 if (n.gt.2) then u(1,1) = u(1,1) + by(3)*tr2 - bx(3)*ti2 u(1,2) = u(1,2) - by(3)*ti2 - bx(3)*tr2 endif endif u(2,1) = by(1)*ti + bx(1)*tr do 100 i=2,n-2 tt = (x*tr-y*ti)/i ti = (x*ti+y*tr)/i tr = tt tr2 = tr*r2/(i+1) ti2 = ti*r2/(i+1) c tr + i ti = (q/4p)(x+iy)^i/i! c tr2 + i ti2 = (q/4p)(x+iy)^i (x^2+y^2)/(i+1)! xf = xf + by(i+1)*tr2 - bx(i+1)*ti2 - by(i-1)*tr + bx(i-1)*ti yf = yf - by(i+1)*ti2 - bx(i+1)*tr2 - by(i-1)*ti - bx(i-1)*tr u(1,1) = u(1,1) + by(i+2)*tr2 - bx(i+2)*ti2 u(1,2) = u(1,2) - by(i+2)*ti2 - bx(i+2)*tr2 u(2,1) = u(2,1) + by(i)*ti + bx(i)*tr 100 continue if (n.gt.2) then tt = (x*tr-y*ti)/(n-1) ti = (x*ti+y*tr)/(n-1) tr = tt tr2 = tr*r2/n ti2 = ti*r2/n xf = xf + by(n)*tr2 - bx(n)*ti2 - by(n-2)*tr + bx(n-2)*ti yf = yf - by(n)*ti2 - bx(n)*tr2 - by(n-2)*ti - bx(n-2)*tr u(2,1) = u(2,1) + by(n-1)*ti + bx(n-1)*tr endif if (n.gt.1) then tt = (x*tr-y*ti)/n ti = (x*ti+y*tr)/n tr = tt xf = xf - by(n-1)*tr + bx(n-1)*ti yf = yf - by(n-1)*ti - bx(n-1)*tr u(2,1) = u(2,1) + by(n)*ti + bx(n)*tr endif tt = (x*tr-y*ti)/(n+1) ti = (x*ti+y*tr)/(n+1) tr = tt xf = xf - by(n)*tr + bx(n)*ti yf = yf - by(n)*ti - bx(n)*tr tt = 2d0*u(2,1) u(2,1) = u(1,2) - tt u(1,2) = u(1,2) + tt u(2,2) = 1d0 - u(1,1) u(1,1) = 1d0 + u(1,1) tt = 1d0/(u(1,1)*u(2,2)-u(1,2)*u(2,1)) tr = (px*u(2,2) - py*u(2,1))*tt py = (py*u(1,1) - px*u(1,2))*tt px = tr t = t + (px*xf+py*yf)*sqrt(p*p+m*m)/(p*p*c) x = x + xf y = y + yf end c ********************************************************* INTEGER FUNCTION LASTNB(STRING) CHARACTER*(*) STRING INTEGER I I = LEN(STRING) 10 IF (STRING(I:I).NE.' ') GO TO 20 I = I - 1 if (I.EQ.0) GO TO 20 GO TO 10 20 LASTNB = I RETURN END c ********************************************************* subroutine LOAD_MATERIAL ! ! Load parameters for this material ! implicit double precision (a-h,o-z) include 'icommon.inc' ! ! Load material parameters do i=1,nmat if( mt .eq. matid(i) ) then if( mt .eq. matdens ) then ! DENS command case mparm(2,krad) = matdat(2,i) * densfac mparm(4,krad) = matdat(4,i) / densfac ! else if( mt .eq. matdenp ) then ! DENP command case if( dirdenp .eq. 1 ) then x = xp(1) fact = adenp + bdenp*x + cdenp*x*x + ddenp*x*x*x mparm(2,krad) = matdat(2,i) * fact mparm(4,krad) = matdat(4,i) / fact else if( dirdenp .eq. 2 ) then x = xp(2) fact = adenp + bdenp*x + cdenp*x*x + ddenp*x*x*x mparm(2,krad) = matdat(2,i) * fact mparm(4,krad) = matdat(4,i) / fact else if( dirdenp .eq. 3 ) then x = SQRT( xp(1)**2 + xp(2)**2 ) fact = adenp + bdenp*x + cdenp*x*x + ddenp*x*x*x mparm(2,krad) = matdat(2,i) * fact mparm(4,krad) = matdat(4,i) / fact end if ! else ! normal case mparm(2,krad) = matdat(2,i) mparm(4,krad) = matdat(4,i) end if ! mparm(1,krad) = matdat(1,i) ! mparm(3,krad) = matdat(3,i) do j=5,9 mparm(j,krad) = matdat(j,i) end do ! ! Watch out for user setting density to 0 if( mparm(2,krad) .eq. 0d0 ) mt = 'VAC' ! nelcomp(krad) = matdat(10,i) do j=1,nelcomp(krad) zelcomp(j,krad) = matdat(7+4*j,i) aelcomp(j,krad) = matdat(8+4*j,i) ipelcomp(j,krad) = matdat(9+4*j,i) felcomp(j,krad) = matdat(10+4*j,i) end do go to 155 ! end if end do 155 continue ! ! Load some material dependent quantities for Moliere scattering ! Ref. Geant3 manual, p. 234 ! FELCOMP is the number of atoms of this type in a compound zsmcs = 0.d0 zemcs = 0.d0 ! if( scatlev.eq.6 .and. ABS(iptyp).ne.1 ) then ! Fano, heavy do i=1,nelcomp(krad) zsmcs = zsmcs + felcomp(i,krad) * & zelcomp(i,krad)**2 zemcs = zemcs + felcomp(i,krad) * & zelcomp(i,krad)**2 * & LOG( zelcomp(i,krad)**(-0.6667) ) end do else do i=1,nelcomp(krad) zsmcs = zsmcs + felcomp(i,krad) * & zelcomp(i,krad) * (zelcomp(i,krad)+1.) zemcs = zemcs + felcomp(i,krad) * & zelcomp(i,krad) * (zelcomp(i,krad)+1.) * & LOG( zelcomp(i,krad)**(-0.6667) ) end do end if ! wtmol = 0d0 do i=1,nelcomp(krad) ! molecular weight wtmol = wtmol + felcomp(i,krad)*aelcomp(i,krad) end do ! Geant3 manual, p. 233, eq. 4 chiccmcs = 0.156915d0 * facfms * zsmcs * mparm(2,krad) / wtmol ! densat = avagn * mparm(2,krad) / wtmol ! Ref. PDG p. 270, eq. 27.5 drnorm = twopi * relcl**2 * mparm(1,krad) & * (ABS(pcharge)/ee)**2 ! [cm^2] ! end ! ********************************************************* subroutine LOREN4(beta,x,t) ! ! Lorentz transformation of 4-vector along z axis ! implicit double precision (a-h,o-z) dimension x(3) ! gamma = 1d0 / SQRT(1d0 - beta**2) zp = x(3) tp = t ! t = gamma * (tp - beta*zp) x(3) = gamma * (zp - beta*tp) ! end ! ********************************************************* subroutine LORENEBG(bet,gamma,e,b) ! ! General transform of rest-frame electromagnetic field to LAB frame ! implicit real*8 (a-h,o-z) include 'icommon.inc' dimension e(3),b(3),bet(3),erf(3),brf(3),crbb(3),crbe(3) ! do i=1,3 erf(i) = e(i) brf(i) = b(i) end do f = gamma**2 / (gamma+1d0) dbb = DOT3(bet,brf) dbe = DOT3(bet,erf) call DCROSS(bet,brf,crbb) call DCROSS(bet,erf,crbe) ! do i=1,3 e(i) = gamma*(erf(i)+crbb(i)*cc) - f*dbe*bet(i) b(i) = gamma*(brf(i)-crbe(i)/cc) - f*dbb*bet(i) end do ! end c ********************************************************* subroutine MAKE_BEAM c c Generate incident beam particles c implicit double precision (a-h,o-z) include 'icommon.inc' save nbad data nbad/0/ c ievt = jpar ipnum = 1 ipflg = 0 evtwt = 1.d0 c c Assign particle mass randomly according to specified beam fractions rn = RAN1(irnarg) xlim = 0.d0 c do i=1,nbeamtyp xlim = xlim + fracbt(i) if( rn .le. xlim ) then iptyp = bmtype(i) kk = i go to 100 end if end do c 100 continue if( bmalt ) then if( MOD(ievt,2) .eq. 0 ) iptyp = -iptyp end if pmass = mass( ABS(iptyp) ) pcharge = SIGN( ee, DBLE(iptyp) ) * chrg( ABS(iptyp) ) c -------------------------------------------------------- 120 continue if( bdistyp(kk) .eq. 1 ) then ! gaussian do i=1,3 xp(i) = GAUSS( x1bt(i,kk),x2bt(i,kk) ) end do c do i=1,3 pp(i) = GAUSS( p1bt(i,kk),p2bt(i,kk) ) end do ! ........................................................ else if( bdistyp(kk) .eq. 2 ) then ! uniform circle ! Weight the radius to prevent enhancement near center (YF) r = SQRT( RANV(x1bt(1,kk)**2,x2bt(1,kk)**2) ) phi = RANV(x1bt(2,kk),x2bt(2,kk)) xp(3) = RANV(x1bt(3,kk),x2bt(3,kk)) pr = SQRT( RANV(p1bt(1,kk)**2,p2bt(1,kk)**2) ) pphi = RANV(p1bt(2,kk),p2bt(2,kk)) pp(3) = RANV(p1bt(3,kk),p2bt(3,kk)) ! ! Get independent phi angle for the momentum (YF) phip = RANV(x1bt(2,kk),x2bt(2,kk)) cph = COS(phip) sph = SIN(phip) ! xp(1) = r * COS(phi) xp(2) = r * SIN(phi) pp(1) = pr*cph - pphi*sph pp(2) = pr*sph + pphi*cph end if c --------------------------------------------------------- ! ! Beam correlations do i=1,nbcorr(kk) ic = corrtyp(i,kk) if( ic .eq. 1 ) then ! angular momentum call VMAG(pp,pm) c = corr1(i,kk) pp(1) = pp(1) + 0.5d0 * pcharge * xp(2) * c pp(2) = pp(2) - 0.5d0 * pcharge * xp(1) * c arg = pm**2 - pp(1)**2 - pp(2)**2 if( arg .gt. 0.d0 ) then pp(3) = SQRT( arg ) else nbad = nbad + 1 if( nbad .lt. 100 ) then go to 120 ! try again else write(2,*) 'Trouble setting angular momentum correlation' iflag = -13 go to 900 end if end if ! else if( ic .eq. 2 ) then ! Palmer amplitude r2 = xp(1)**2 + xp(2)**2 th2 = (pp(1)/pp(3))**2 + (pp(2)/pp(3))**2 c = corr1(i,kk) bp = corr2(i,kk) dpz = c * ( r2/(bp*bp) + th2 ) pp(3) = pp(3) + dpz ! else if( ic.gt.2 .and. ic.lt.6 ) then ! rf bucket call RF_BUCKET(i,kk) ! else if( ic.gt.5 .and. ic.lt.8 ) then ! Twiss parameters call TWISS(i,kk) ! else if( ic .eq. 9 ) then ! Compensate solenoid path length with extra transverse momentum b0 = corr1(i,kk) if( b0 .gt. 0.d0 ) then arg = pp(3)**2*(1.d0-b0*b0)-(pmass*b0)**2 if( arg .gt. 0.d0 ) then pt = SQRT( arg ) / b0 phi = corr2(i,kk) * pi / 180.d0 if( phi .eq. 0.d0 ) phi = RANV(0.d0,twopi) pp(1) = pt * COS(phi) pp(2) = pt * SIN(phi) end if end if ! else if( ic .eq. 10 ) then ! Balbekov energy-amplitude eref = corr1(i,kk) babs = corr2(i,kk) sige = corr3(i,kk) r2 = xp(1)**2 + xp(2)**2 pt2 = pp(1)**2 + pp(2)**2 ab = (pt2 + r2*(ee*babs/2.d0)**2) / pmass**2 de = GAUSS( 0.d0,sige ) e = eref*SQRT(1.d0+ab) + de pm2 = e**2 - pmass**2 pp(3) = SQRT( pm2 - pt2 ) ! else if( ic .eq. 11 ) then ! dispersion d = corr1(i,kk) pref = corr2(i,kk) idir = NINT(corr3(i,kk)) del = (pp(3)-pref) / pref fac = d * del ! if( idir .eq. 1 ) then ! 1 => x xp(1) = xp(1) + fac else if( idir .eq. 2) then ! 2 => y xp(2) = xp(2) + fac else if( idir .eq. 3) then ! 3 => x' call VMAG(pp,pm0) xprm = pp(1)/pp(3) + fac pp(1) = xprm * pp(3) arg = pm0**2 - pp(1)**2 - pp(2)**2 if( arg .gt. 0. ) then pp(3) = SQRT(arg) else nbad = nbad + 1 if( nbad .lt. 100 ) then go to 120 ! try again else write(2,*) 'Trouble setting dispersion correlation' iflag = -16 go to 900 end if end if else if( idir .eq. 4) then ! 4 => y' call VMAG(pp,pm0) yprm = pp(2)/pp(3) + fac pp(2) = yprm * pp(3) arg = pm0**2 - pp(1)**2 - pp(2)**2 if( arg .gt. 0. ) then pp(3) = SQRT(arg) else nbad = nbad + 1 if( nbad .lt. 100 ) then go to 120 ! try again else write(2,*) 'Trouble setting dispersion correlation' iflag = -16 go to 900 end if end if end if ! end of test on idir ! end if ! end of test on corrtyp end do ! end of loop on nbcorr ! --------------------------------------------------------- ! Convert initial bunch length into time spread call VMAG(pp,pm) gamma = SQRT(1.d0 + pm*pm/pmass/pmass) vz = pp(3) / pmass / gamma * cc tp = xp(3) / vz if( bmalt ) then if( MOD(ievt,2) .eq. 0 ) tp = tp + dtalt end if xp(3) = 0.d0 ! do i=1,3 pol(i) = 0.d0 end do POL(3) = 1.d0 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! 900 continue return end c ********************************************************* subroutine MSTEP(pm) ! ! Simulate processes in a step thru material medium ! implicit double precision (a-h,o-z) include 'icommon.inc' ! logical in_wedge,in_pwedge,in_aspwedge,in_hemiwin,in_asrwedge logical in_ring,in_absorb,in_nia character*4 blnk data blnk/' '/ ! pm0 = pm ! original magnitude of momentum hstp = sdid ! hstp can be changed in wedges ! in_absorb = .false. if( mgt .eq. 'CBLOCK' ) then ! normal absorber in_absorb = .true. ! else if( mgt .eq. 'WEDGE' ) then mt = mtag(krad) call LOAD_MATERIAL in_absorb = IN_WEDGE(hstp) if( .not. in_absorb ) then ! test for external material if( mtag2(krad).ne.blnk .and. mtag2(krad).ne.'VAC' ) then mt = mtag2(krad) call LOAD_MATERIAL in_absorb = .true. end if end if ! else if( mgt .eq. 'PWEDGE' ) then in_absorb = IN_PWEDGE(hstp) else if( mgt .eq. 'HWIN' ) then in_absorb = IN_HEMIWIN() else if( mgt .eq. 'ASPW' ) then in_absorb = IN_ASPWEDGE(hstp) else if( mgt .eq. 'RING' ) then in_absorb = IN_RING(hstp) else if( mgt .eq. 'ASRW' ) then in_absorb = IN_ASRWEDGE(hstp) else if( mgt .eq. 'NIA' ) then in_absorb = IN_NIA(hstp) end if ! if( iflag .ne. 0 ) go to 900 if( hstp .eq. 0. ) go to 900 ! --------------------------------------------------------- if( in_absorb ) then if( lelms .and. mt.eq.'LH' ) then ! use ELMS call ELMS(hstp,pm) if( iflag .ne. 0 ) go to 900 ! else if( lsamcs ) then ! use SAMCS call SAMCS_DRIVER(hstp,pm) if( iflag .ne. 0 ) go to 900 ! else if( ldray ) then ! use delta rays ! Continous processes for energy transfers below DCUTx delev = 2 ! straglev=5 needs mean dE/dx straglev = 5 call DEDX(hstp,pm) if( iflag .ne. 0 ) go to 900 ! pavg = 0.5 * (pm + pm0) if( lscatter ) call SCATTER(hstp,pavg) ! Multiple scattering if( iflag .ne. 0 ) go to 900 ! Discrete processes for energy transfers above DCUTx call DELTA_RAY(hstp,pavg) if( iflag .ne. 0 ) go to 900 ! ......................................................... else ! uncorrelated electromagnetic processes if( ldedx ) call DEDX(hstp,pm) ! ionization energy loss if( iflag .ne. 0 ) go to 900 ! ! pm is returned as the magnitude of the momentum after energy loss pavg = 0.5 * (pm + pm0) ! if( lscatter ) call SCATTER(hstp,pavg) ! Multiple scattering if( iflag .ne. 0 ) go to 900 ! end if ! pavg = 0.5 * (pm + pm0) if ( linteract ) call INTERACTION(hstp,pavg) ! nuclear interaction if( iflag .ne. 0 ) go to 900 ! end if ! --------------------------------------------------------- ! 900 continue end ! ********************************************************* subroutine NEW_PARTICLE(ip,ifl) ! ! Set up data for new particle ! implicit double precision (a-h,o-z) include 'icommon.inc' ! ifl = 0 ! --------------------------------------------------------- ! if( ip .eq. 0 ) then ! initialize REFERENCE particle ievt = 0 ipnum = 0 iptyp = 2 ! default to a muon refpar = 2 ipflg = 0 do i=1,3 xp(i) = 0.d0 pp(i) = 0.d0 end do tp = 0.d0 pmass = mass(2) pcharge = ee evtwt = 1.d0 ! if( .not. bgen )then ! take data from FOR003.DAT beam file tp = tbref xp(3) = zbref pp(3) = pbref refpar = typref ! pmass = mass( ABS(typref) ) ! pcharge = SIGN( ee,DBLE(typref) ) * chrg(ABS(typref)) end if ! ......................................................... else if( ip .eq. -1 ) then ! initialize 2nd REFERENCE particle ievt = -1 ipnum = 0 iptyp = 2 ! default to a muon ref2par = 2 ipflg = 0 do i=1,3 xp(i) = 0.d0 pp(i) = 0.d0 end do tp = 0.d0 pmass = mass(2) pcharge = ee evtwt = 1.d0 ! if( .not. bgen )then ! take data from FOR003.DAT beam file tp = tb2ref xp(3) = zb2ref pp(3) = pb2ref ref2par = typref ! pmass = mass( ABS(typref) ) ! pcharge = SIGN( ee,DBLE(typref) ) * chrg(ABS(typref)) end if ! ! ......................................................... else ! initialize NORMAL particles ! if( bgen ) then ! generate new particle if (run_env) then call ENV_MAKE_BEAM else call MAKE_BEAM if( iflag .ne. 0 ) then ifl = -1 go to 900 end if end if ! else ! read in an old one read(3,*,iostat=ioc)ievt,ipnum,iptyp,ipflg,tp,evtwt,xp,pp, & pol ! if( ioc .ne. 0 ) then ifl = -1 go to 900 end if ! if( ipflg .ne. 0 ) then if( goodtrack ) then ! skip bad tracks ifl = -1 go to 900 else ! use bad tracks ipflg = 0 end if end if ! pmass = mass( ABS(iptyp) ) pcharge = SIGN( ee,DBLE(iptyp) ) * chrg(ABS(iptyp)) if( declev .eq. 5 ) pol(1) = 0. end if ! end of test on BGEN ! end if ! end of test on ip=0 ! --------------------------------------------------------- ! if( ip .le. nprnt ) then write(2,105)ip,iptyp,xp,pp,tp 105 format(' INITIAL:',2i4,/,7e11.4) ! if( spin .and. ip.gt.0 ) then write(2,110)pol 110 format(' SPIN:',3f8.4) end if ! end of test on spin end if ! end of test on ip.le.nprnt ! do i=1,3 ! save initial particle properties xpint(i) = xp(i) ppint(i) = pp(i) polinit(i) = pol(i) end do tpint = tp ! xp(3) = 0d0 ! 900 continue return end c ********************************************************* subroutine OUTPUT ! ! Write particle information at this location to postprocessor file ! ! The header for this write statement is in subr SETUP ! implicit double precision (a-h,o-z) include 'icommon.inc' character*80 fmt ! ipn = MIN(ipnum,100) ! Boris step has staggered fields/particles, so get local field if( (.not.steprk) .and. (jsrg.ne.1) ) call FIELD(xp,tp) ! if( f9dp .eq. 4 ) then fmt = '(i7,2i4,i5,i6,1x,1p,18e12.4)' else if( f9dp .eq. 6 ) then fmt = '(i7,2i4,i5,i6,1x,1p,18e14.6)' else if( f9dp .eq. 8 ) then fmt = '(i7,2i4,i5,i6,1x,1p,18e16.8)' else if( f9dp .eq. 10 ) then fmt = '(i7,2i4,i5,i6,1x,1p,18e18.10)' else if( f9dp .eq. 12 ) then fmt = '(i7,2i4,i5,i6,1x,1p,18e20.12)' else if( f9dp .eq. 14 ) then fmt = '(i7,2i4,i5,i6,1x,1p,18e22.14)' else if( f9dp .eq. 16 ) then fmt = '(i7,2i4,i5,i6,1x,1p,18e24.16)' else if( f9dp .eq. 17 ) then fmt = '(i7, 2i4, i5, i6, 1x, 1p, 18e25.16E3)' else fmt = '(i7,2i4,i5,i6,1x,1p,18e12.4)' ! default end if ! if( jsrg .eq. 1 ) then write(9,fmt)ievt,ipn,iptyp,ipflg,jsrg,tpint,xpint,ppint, & bfld,evtwt,efld,sarc,polinit else write(9,fmt)ievt,ipn,iptyp,ipflg,jsrg,tp,xp,pp,bfld,evtwt, & efld,sarc,pol end if ! if (run_env .and. (ievt.eq.1) ) then write(11,200)jsrg,xp(3),bfld(3),tp,env_pz,env_emit, & env_beta,env_alpha,env_angm,env_N,env_gauss,env_eL 200 format(i5,11e14.5) end if ! return end ! ********************************************************* subroutine PHMODEL(iskip) ! ! Set reference particle properties ! implicit double precision (a-h,o-z) include 'icommon.inc' ! iskip = 0 ! C Reference particle phase model 2 C Set the phase of the cavity such that if the phase were C zero, the reference particle would not be accelerated C Use the actual cavity fields and frequency if ( phmodref.eq.2 .and. ftag(1).eq.'ACCE' .and. & (fparm(1,1).eq.1d0 .or. fparm(1,1).eq.2d0) .or. $ fparm(1,1).eq.13d0) then C To implement other accel models, modify findzc and/or C this call to change range of integration and/or C the longituinally-dependent position trefmean = tp - FINDZC(twopi*fparm(2,1)*1d6,fparm(3,1)*1d-3, & slen,zstep,pp(3),pmass) end if ! ! --------------------------------------------------------- if( phmodref .eq. 3 ) then ! assume constant ref particle velocity ! do this for every region for plotting purposes u = pp(3) / pmass betam = u / SQRT(1.d0+u*u) tnow = tp dt = slen / (betam*cc) if (varstep) then tp = tnow + dt xp(3) = xp(3) + slen sarc = sarc + slen else if( scalestep*zstep .le. slen ) then maxsteps = NINT( slen / (scalestep*zstep) ) else maxsteps = 1 end if fixstep = slen / DBLE(maxsteps) if (rtuple.and. $ (ftag(1).eq.'ACCE'.and.fparm(1,1).eq.13.and. $ (fparm(5,1).eq.0.or.fparm(5,1).eq.2).or. $ ftag(1).eq.'SOL'.and.fparm(1,1).eq.8)) then CALL FIELD(xp,tp) CALL OUTPUT end if do 310 iz=1,maxsteps xp(3) = xp(3) + fixstep sarc = sarc + fixstep tp = tp + fixstep/(betam*cc) if( rtuple .and. MOD(iz-1,rtuplen).eq.0 ) then call FIELD(xp,tp) call OUTPUT end if 310 continue if (rtuple.and. $ (ftag(1).eq.'ACCE'.and.fparm(1,1).eq.13.and. $ (fparm(5,1).eq.0.or.fparm(5,1).eq.1).or. $ ftag(1).eq.'SOL'.and.fparm(1,1).eq.8.or. $ ftag(1).eq.'DIP'.and.fparm(1,1).eq.3.and. $ fparm(2,1).lt.0d0.and.fparm(6,1).eq.1)) then CALL FIELD(xp,tp) CALL OUTPUT end if end if if( jpar .eq. 0 ) then trefmean = tnow + 0.5d0*dt tpref = tp ! needed for accel, model=6 logic else if( jpar .eq. -1 ) then t2refmean = tnow + 0.5d0*dt tp2ref = tp end if ft = ftag(1) pzrefmean = pp(3) call ACCEL_PHASE(fparm) ! n.b. call FIELD after ACCEL-PHASE to get ACCEL(11) to work call FIELD(xp,tp) ! get field value for an OUTPUT iskip = 1 ! skip tracking end if ! end of test on phmoderef = 3 ! ! --------------------------------------------------------- if( phmodref.eq.4 .and. ftag(1).eq.'ACCE' ) then ! n.b. leave ftag in this test or ref par doesn't track thru absorbers ! C Logic here is as follows: with a variable step size, just compute C the step all at once. For a fixed step size, we must mach what C is happening for non-reference particles. If the gradient is C zero or really small, the time of flight is computed based on C the average of the inverse of the velocity. If the gradient is C non zero, the time-of-flight is is just dp/g, where dp is the C change in momentum and g is the gradient C pz0 = pp(3) e0 = SQRT( pz0**2 + pmass**2 ) tnow = tp if (jpar .eq. 0) then a1 = 1d-3*v3ref ! [GeV/m] else if (jpar .eq. -1) then a1 = 1d-3*v3ref2 end if ! if (varstep) then xp(1) = xp(1) + pp(1)/pz0*slen xp(2) = xp(2) + pp(2)/pz0*slen xp(3) = xp(3) + slen sarc = sarc + slen e1 = e0 + a1*slen pp(3) = sqrt(e1**2-pmass**2) tp = tnow + 0.5d0*(e0/pz0+e1/pp(3))*slen/cc else if( scalestep*zstep .le. slen ) then maxsteps = NINT( slen / (scalestep*zstep) ) else maxsteps = 1 end if fixstep = slen / DBLE(maxsteps) if (rtuple.and.fparm(1,1).eq.13.and. $ (fparm(5,1).eq.0.or.fparm(5,1).eq.2)) then CALL FIELD(xp,tp) CALL OUTPUT end if ! do 330 iz=1,maxsteps xp(1) = xp(1) + pp(1)/pz0*fixstep xp(2) = xp(2) + pp(2)/pz0*fixstep xp(3) = xp(3) + fixstep sarc = sarc + fixstep ! e1 = e0 + gr*fixstep e1 = e0 + a1*fixstep pz1 = sqrt(e1**2-pmass**2) ! if (0.5d0*(gr*fixstep*pmass*pp(3)**(-2))**2.lt. ! $ dbleps) then tp = tp + 0.5d0*(e0/pp(3)+e1/pz1)*fixstep/cc ! else ! tp = tp + (pz1-pp(3))/(gr*cc) ! end if e0 = e1 pp(3) = pz1 if( rtuple .and. MOD(iz-1,rtuplen).eq.0 ) then call FIELD(xp,tp) call OUTPUT end if 330 continue ! if (rtuple.and.fparm(1,1).eq.13.and. $ (fparm(5,1).eq.0.or.fparm(5,1).eq.1)) then CALL FIELD(xp,tp) CALL OUTPUT end if end if ! end of test on varstep pzrefmean = 0.5d0*(pz0+pp(3)) if (jpar .eq. 0) then trefmean = 0.5d0*(tnow+tp) tpref = tp else if (jpar .eq. -1) then t2refmean = 0.5d0*(tnow+tp) tp2ref = tp end if ft = ftag(1) if( forcerp ) then ! force ref particle to axis call VMAG(pp,pm) pp(1) = 0.0d0 pp(2) = 0.0d0 pp(3) = pm xp(1) = 0.0d0 xp(2) = 0.0d0 end if call ACCEL_PHASE(fparm) ! n.b. call FIELD after ACCEL-PHASE to get ACCEL(11) to work call FIELD(xp,tp) ! get field value for an OUTPUT iskip = 1 ! skip tracking end if ! end of test on phmoderef = 4 ! --------------------------------------------------------- if( phmodref.eq.5 .and. ftag(1).eq.'ACCE' ) then ! n.b. leave ftag in this test or ref par doesn't track thru absorbers ! pz0 = pp(3) e0 = SQRT( pz0**2 + pmass**2 ) tnow = tp zend = xp(3) + slen ! actual position at end of cavity zacc = zend - z0pmref ! distance of end from last reset if (jpar .eq. 0) then a0 = v1ref ! [GeV] a1 = 1d-3*v2ref ! [GeV/m] a2 = 1d-3*v3ref ! [GeV/m^2] else if (jpar .eq. -1) then a0 = v1ref2 ! [GeV] a1 = 1d-3*v2ref2 ! [GeV/m] a2 = 1d-3*v3ref2 ! [GeV/m^2] end if ! if (varstep) then xp(1) = xp(1) + pp(1)/pz0*slen xp(2) = xp(2) + pp(2)/pz0*slen xp(3) = zend sarc = sarc + slen e1 = a0 + a1*zacc + a2*zacc**2 pp(3) = sqrt(e1**2-pmass**2) tp = tnow + 0.5d0*(e0/pz0+e1/pp(3))*slen/cc else if( scalestep*zstep .le. slen ) then maxsteps = NINT( slen / (scalestep*zstep) ) else maxsteps = 1 end if fixstep = slen / DBLE(maxsteps) if (rtuple.and.fparm(1,1).eq.13.and. $ (fparm(5,1).eq.0.or.fparm(5,1).eq.2)) then CALL FIELD(xp,tp) CALL OUTPUT end if ! do 350 iz=1,maxsteps xp(1) = xp(1) + pp(1)/pz0*fixstep xp(2) = xp(2) + pp(2)/pz0*fixstep xp(3) = xp(3) + fixstep sarc = sarc + fixstep ! e1 = e0 + gr*fixstep z1 = xp(3) - z0pmref e1 = a0 + a1*z1 + a2*z1*z1 pz1 = SQRT(e1**2-pmass**2) ! if (0.5d0*(gr*fixstep*pmass*pp(3)**(-2))**2.lt. ! $ dbleps) then tp = tp + 0.5d0*(e0/pp(3)+e1/pz1)*fixstep/cc ! else ! tp = tp + (pz1-pp(3))/(gr*cc) ! end if e0 = e1 pp(3) = pz1 if( rtuple .and. MOD(iz-1,rtuplen).eq.0 ) then call FIELD(xp,tp) call OUTPUT end if 350 continue ! if (rtuple.and.fparm(1,1).eq.13.and. $ (fparm(5,1).eq.0.or.fparm(5,1).eq.1)) then CALL FIELD(xp,tp) CALL OUTPUT end if end if ! end of test on varstep pzrefmean = 0.5d0*(pz0+pp(3)) if (jpar .eq. 0) then trefmean = 0.5d0*(tnow+tp) tpref = tp else if (jpar .eq. -1) then t2refmean = 0.5d0*(tnow+tp) tp2ref = tp end if ft = ftag(1) if( forcerp ) then ! force ref particle to axis call VMAG(pp,pm) pp(1) = 0.0d0 pp(2) = 0.0d0 pp(3) = pm xp(1) = 0.0d0 xp(2) = 0.0d0 end if call ACCEL_PHASE(fparm) ! n.b. call FIELD after ACCEL-PHASE to get ACCEL(11) to work call FIELD(xp,tp) ! get field value for an OUTPUT iskip = 1 ! skip tracking end if ! end of test on phmoderef = 5 ! --------------------------------------------------------- if( phmodref .eq. 6) then ! pz0 = pp(3) e0 = SQRT( pz0**2 + pmass**2 ) tnow = tp zend = xp(3) + slen ! actual position at end of cavity zacc = zend - z0pmref ! distance of end from last reset if (jpar .eq. 0) then a0 = v1ref ! [GeV] a1 = 1d-3*v2ref ! [GeV/m] a2 = 1d-3*v3ref ! [GeV/m^2] else if (jpar .eq. -1) then a0 = v1ref2 ! [GeV] a1 = 1d-3*v2ref2 ! [GeV/m] a2 = 1d-3*v3ref2 ! [GeV/m^2] end if ! if (varstep) then xp(1) = xp(1) + pp(1)/pz0*slen xp(2) = xp(2) + pp(2)/pz0*slen xp(3) = zend sarc = sarc + slen e1 = a0 + a1*zacc + a2*zacc**2 pp(3) = sqrt(e1**2-pmass**2) tp = tnow + 0.5d0*(e0/pz0+e1/pp(3))*slen/cc else if( scalestep*zstep .le. slen ) then maxsteps = NINT( slen / (scalestep*zstep) ) else maxsteps = 1 end if fixstep = slen / DBLE(maxsteps) if (rtuple.and.fparm(1,1).eq.13.and. $ (fparm(5,1).eq.0.or.fparm(5,1).eq.2)) then CALL FIELD(xp,tp) CALL OUTPUT end if ! do 370 iz=1,maxsteps xp(1) = xp(1) + pp(1)/pz0*fixstep xp(2) = xp(2) + pp(2)/pz0*fixstep xp(3) = xp(3) + fixstep sarc = sarc + fixstep ! e1 = e0 + gr*fixstep z1 = xp(3) - z0pmref e1 = a0 + a1*z1 + a2*z1*z1 pz1 = SQRT(e1**2-pmass**2) ! if (0.5d0*(gr*fixstep*pmass*pp(3)**(-2))**2.lt. ! $ dbleps) then tp = tp + 0.5d0*(e0/pp(3)+e1/pz1)*fixstep/cc ! else ! tp = tp + (pz1-pp(3))/(gr*cc) ! end if e0 = e1 pp(3) = pz1 if( rtuple .and. MOD(iz-1,rtuplen).eq.0 ) then call FIELD(xp,tp) call OUTPUT end if 370 continue ! if (rtuple.and.fparm(1,1).eq.13.and. $ (fparm(5,1).eq.0.or.fparm(5,1).eq.1)) then CALL FIELD(xp,tp) CALL OUTPUT end if end if ! end of test on varstep pzrefmean = 0.5d0*(pz0+pp(3)) if (jpar .eq. 0) then trefmean = 0.5d0*(tnow+tp) tpref = tp else if (jpar .eq. -1) then t2refmean = 0.5d0*(tnow+tp) tp2ref = tp end if ft = ftag(1) if( forcerp ) then ! force ref particle to axis call VMAG(pp,pm) pp(1) = 0.0d0 pp(2) = 0.0d0 pp(3) = pm xp(1) = 0.0d0 xp(2) = 0.0d0 end if call ACCEL_PHASE(fparm) ! n.b. call FIELD after ACCEL-PHASE to get ACCEL(11) to work call FIELD(xp,tp) ! get field value for an OUTPUT iskip = 1 ! skip tracking end if ! end of test on phmoderef = 6 ! end ! ********************************************************* subroutine PKICK(gamma,length,e,b) ! ! Kick particle momentum array PP due to EM fields ! implicit double precision (a-h,o-z) include 'icommon.inc' real*8 length ! dimension e(3),b(3),bet(3),dp(3) ! do i=1,3 bet(i) = pp(i) / pmass / gamma end do ! term = length / bet(3) dp(1) = pcharge * ( e(1)/cc + bet(2)*b(3) - bet(3)*b(2) ) * term dp(2) = pcharge * ( e(2)/cc - bet(1)*b(3) + bet(3)*b(1) ) * term dp(3) = pcharge * ( e(3)/cc + bet(1)*b(2) - bet(2)*b(1) ) * term ! do i=1,3 pp(i) = pp(i) + dp(i) end do ! return end c ********************************************************* subroutine PSEUDOREGION(i7,ip) ! ! Special purpose "region" commands ! implicit double precision (a-h,o-z) include 'icommon.inc' ! save jx,jy dimension bxar(3),byar(3) character*10 fname character*4 grtype,mat ! idum = i7 ! ........................................................ if( jpsr .eq. -1 ) then ! rotation plane ang = val1ps iaxis = ival1ps iclass = ival2ps if( iclass.eq.0 .or. (iclass.eq.1 .and. ip.le.0) .or. & (iclass.eq.2 .and. ip.gt.0) ) then ! call ROTATE_COORD(ang,iaxis) ! end if if( ip .le. nprnt ) then write(2,75) ang*180./pi,iaxis,xp,pp 75 format(' ROT=',f7.2,i3,6e11.4) end if ! ! ........................................................ ! Save output plane information ! else if( jpsr.eq.-2 .and. .not.(ntuple.or.rtuple) ) then ! ! loutput = .true. ! ! ........................................................ else if( jpsr.eq.-3 .and. ip.gt.0 ) then ! aperture plane iapertype = ival1ps aperlim1 = val1ps aperlim2 = val2ps ! if( iapertype .eq. 1 ) then ! elliptical cut if( ABS(xp(1)) .ge. aperlim1 ) then iflag = -28 go to 900 end if yell = aperlim2 * SQRT(1. - (xp(1)/aperlim1)**2) if( ABS(xp(2)) .ge. yell ) then iflag = -28 go to 900 end if ! else if( iapertype .eq. 2 ) then ! rectangular cut idir = val3ps if( idir .eq. 1 ) then if( xp(1).le.aperlim1 .or. xp(1).ge.aperlim2 ) then iflag = -28 go to 900 end if else if( xp(2).le.aperlim1 .or. xp(2).ge.aperlim2 ) then iflag = -28 go to 900 end if end if ! end of test on IDIR ! else if( iapertype .eq. 3 ) then ! normal quad cut parab = val3ps r = SQRT( xp(1)**2 + xp(2)**2 ) if( r .lt. aperlim1 ) go to 900 if( r .gt. aperlim2 ) then iflag = -28 go to 900 end if ! Rotate particle to c.s. along pole with origin at vertex th = pi / 4.d0 cth = COS(th) sth = SIN(th) xq = cth*xp(1) + sth*xp(2) yq = -sth*xp(1) + cth*xp(2) arg = 4.d0*parab*(xq-aperlim1) if( arg. ge. 0.d0 ) then ypole = SQRT(arg) if( ABS(yq) .le. ypole ) then iflag = -28 go to 900 end if end if th = 3.d0*pi / 4.d0 cth = COS(th) sth = SIN(th) xq = cth*xp(1) + sth*xp(2) yq = -sth*xp(1) + cth*xp(2) arg = 4.d0*parab*(xq-aperlim1) if( arg. ge. 0.d0 ) then ypole = SQRT(arg) if( ABS(yq) .le. ypole ) then iflag = -28 go to 900 end if end if th = 5.d0*pi / 4.d0 cth = COS(th) sth = SIN(th) xq = cth*xp(1) + sth*xp(2) yq = -sth*xp(1) + cth*xp(2) arg = 4.d0*parab*(xq-aperlim1) if( arg. ge. 0.d0 ) then ypole = SQRT(arg) if( ABS(yq) .le. ypole ) then iflag = -28 go to 900 end if end if th = 7.d0*pi / 4.d0 cth = COS(th) sth = SIN(th) xq = cth*xp(1) + sth*xp(2) yq = -sth*xp(1) + cth*xp(2) arg = 4.d0*parab*(xq-aperlim1) if( arg. ge. 0.d0 ) then ypole = SQRT(arg) if( ABS(yq) .le. ypole ) then iflag = -28 go to 900 end if end if ! else if( iapertype .eq. 4 ) then ! skew quad cut parab = val3ps r = SQRT( xp(1)**2 + xp(2)**2 ) if( r .lt. aperlim1 ) go to 900 if( r .gt. aperlim2 ) then iflag = -28 go to 900 end if arg = 4.d0*parab*(xp(1)-aperlim1) if( arg. ge. 0.d0 ) then ypole = SQRT(arg) if( ABS(xp(2)) .le. ypole ) then iflag = -28 go to 900 end if end if ! Rotate particle to c.s. along pole with origin at vertex th = piov2 cth = COS(th) sth = SIN(th) xq = cth*xp(1) + sth*xp(2) yq = -sth*xp(1) + cth*xp(2) arg = 4.d0*parab*(xq-aperlim1) if( arg. ge. 0.d0 ) then ypole = SQRT(arg) if( ABS(yq) .le. ypole ) then iflag = -28 go to 900 end if end if th = pi cth = COS(th) sth = SIN(th) xq = cth*xp(1) + sth*xp(2) yq = -sth*xp(1) + cth*xp(2) arg = 4.d0*parab*(xq-aperlim1) if( arg. ge. 0.d0 ) then ypole = SQRT(arg) if( ABS(yq) .le. ypole ) then iflag = -28 go to 900 end if end if th = 3.d0*piov2 cth = COS(th) sth = SIN(th) xq = cth*xp(1) + sth*xp(2) yq = -sth*xp(1) + cth*xp(2) arg = 4.d0*parab*(xq-aperlim1) if( arg. ge. 0.d0 ) then ypole = SQRT(arg) if( ABS(yq) .le. ypole ) then iflag = -28 go to 900 end if end if ! end if ! ! ........................................................ ! Background field needed else if( jpsr.eq.-4 ) then if(termout) write(*,*)' Resetting background field grid' call BACKSET(0) ! reset background field map sbg0 = saccum ! ! ........................................................ ! Background field region else if( jpsr.eq.-5 ) then if(termout) write(*,*)' Computing background field: ',bftag call BACKSET(1) ! add field to background if( iflag .ne. 0 ) then write(2,*)' Error in computing background field:',bftag go to 900 end if ! ! ........................................................ ! End background field region else if( jpsr.eq.-6 ) then mfile = ival1ps if( mfile.gt.19 .and. mfile.lt.100 ) then ! Save field grid to file write(fname,120) mfile 120 format('for0',i2,'.dat') open(unit=mfile, & file=pathname(1:LASTNB(pathname)) // fname, & status='unknown',iostat=ioc) if( ioc .ne. 0 ) then iflag = -50 go to 900 end if write(mfile,'(a)') 'Background field map' write(mfile,'(3i8,2x,e12.4)') nxgr,nygr,nsgr,hfldbkg do i=1,nxgr write(mfile,'(5e12.4)') xgr(i) end do do j=1,nygr write(mfile,'(5e12.4)') ygr(j) end do do k=1,nsgr write(mfile,'(5e12.4)') sgr(k) end do do k=1,nsgr do j=1,nygr do i=1,nxgr write(mfile,'(3i4,3e12.4)') i,j,k,bxgr(i,j,k), & bygr(i,j,k),bsgr(i,j,k) end do end do end do close(mfile) end if ! if( bentbkg ) then ordtrkbkg = 1 call HUNT(xgr,nxgr,0.d0,jx) if( jx .eq. 0 ) jx = 1 call HUNT(ygr,nygr,0.d0,jy) if( jy .eq. 0 ) jy = 1 do i=1,nsgr htrkgr(i) = pcharge / prefbkg * bygr(jx,jy,i) gtrkgr(i) = pcharge / prefbkg * bxgr(jx,jy,i) if( i .eq. 1 ) then dby = (bygr(jx,jy,2) - bygr(jx,jy,1)) / dsgr else if( i.gt.1 .and. i.lt.nsgr ) then dby = (bygr(jx,jy,i+1) - bygr(jx,jy,i-1)) / (2.*dsgr) else if( i .eq. nsgr ) then dby = (bygr(jx,jy,nsgr) - bygr(jx,jy,nsgr-1)) / dsgr end if hptrkbkg(i) = pcharge / prefbkg * dby if( hptrkbkg(i) .ne. 0. ) ordtrkbkg = 2. end do end if ! ........................................................ ! User transport matrix else if( jpsr.eq.-7 ) then if( jpar .gt. 0 ) then call TRANSPORT end if ! ........................................................ ! Initialize reference particle else if( jpsr.eq.-8 .and. ip.eq.0 ) then ! call READ_EVT(0) ! initial data set by SECT command if( forcerp ) then do i=1,2 xp(i) = 0.0d0 pp(i) = 0.0d0 end do end if if( v1ref.ne.0d0 .and. (phmodref.eq.3.or.phmodref.eq.4) ) & pp(3) = v1ref if( v2ref.ne.0d0 .and. (phmodref.eq.3.or.phmodref.eq.4) ) & tp = v2ref if( phmodref.eq.5 .or. phmodref.eq.6 ) then z0pmref = xp(3) ! reset starting position end if iptyp = refpar ipflg = 0 pmass = mass( ABS(refpar) ) pcharge = SIGN( ee, DBLE(iptyp) ) * chrg( ABS(refpar) ) ppint(3) = pp(3) tpint = tp call WRITE_EVT(0) ! sets xpref,ppref,tpref ! if( ip .le. nprnt ) then write(2,125)ip,iptyp,xp,pp,tp 125 format(' REF PARTICLE REDEFINED:',2i4,/,7e11.4) end if ! ! ........................................................ ! The DUMMY command has jsec = -9 ! ........................................................ else if( jpsr .eq. -10 ) then ! tilted field psi = GAUSS(0.d0,val1ps) if( val2ps .lt. 0. ) then fi = RANV(0.d0,360.d0) else fi = val2ps end if psi = psi*pi/180.d0 fi = fi*pi/180.d0 ! call TILT_COORD(psi,fi) ! if( ip .le. nprnt ) then write(2,175) psi*180./pi,fi*180./pi,xp,pp 175 format(' TILT=',2f7.2,2x,6e11.4) end if ! ........................................................ else if( jpsr .eq. -11 ) then ! displaced field d = GAUSS(0.d0,val1ps) if( val2ps .lt. 0. ) then fi = RANV(0.d0,360.d0) else fi = val2ps end if fi = fi*pi/180.d0 ! call DISP_COORD(d,fi) ! if( ip .le. nprnt ) then write(2,185) d,fi*180./pi,xp,pp 185 format(' DISP=',2f7.2,2x,6e11.4) end if ! ........................................................ else if( jpsr.eq.-12 .and. ip.eq.-1 ) then ! 2nd reference command ! call READ_EVT(-1) ! initial data set by SECT command if( forcerp ) then do i=1,2 xp(i) = 0.0d0 pp(i) = 0.0d0 end do end if if( v1ref2.ne.0d0 .and. (phmodref.eq.3.or.phmodref.eq.4) ) & pp(3) = v1ref2 if( v2ref2.ne.0d0 .and. (phmodref.eq.3.or.phmodref.eq.4) ) & tp = v2ref2 iptyp = ref2par ipflg = 0 pmass = mass( ABS(iptyp) ) pcharge = SIGN( ee, DBLE(iptyp) ) * chrg( ABS(iptyp) ) ppint(3) = pp(3) call WRITE_EVT(-1) ! if( ip .le. nprnt ) then write(2,190)ip,iptyp,xp,pp,tp 190 format(' 2nd REF PARTICLE REDEFINED:',2i4,/,7e11.4) end if ! ! --------------------------------------------------------- else if( jpsr.eq.-13 .and. ip.gt.0 ) then ! cut on variable ivar = ival1ps ivtyp = ival2ps val = val1ps ! if( ivar .le. 3 ) then if( ivtyp .eq. 1 ) then if( xp(ivar) .lt. val ) then iflag = -15 go to 900 end if else if( xp(ivar) .gt. val ) then iflag = -15 go to 900 end if end if ! else if( ivar.ge.4 .and. ivar.le.6 ) then if( ivtyp .eq. 1 ) then if( pp(ivar-3) .lt. val ) then iflag = -15 go to 900 end if else if( pp(ivar-3) .gt. val ) then iflag = -15 go to 900 end if end if ! else if( ivar .eq. 7 ) then if( ivtyp .eq. 1 ) then if( tp .lt. val ) then iflag = -15 go to 900 end if else if( tp .gt. val ) then iflag = -15 go to 900 end if end if ! else if( ivar .eq. 8 ) then call VMAG(pp,pm0) if( ivtyp .eq. 1 ) then if( pm0 .lt. val ) then iflag = -15 go to 900 end if else if( pm0 .gt. val ) then iflag = -15 go to 900 end if end if ! else if( ivar .eq. 9 ) then call VMAG(pp,pm0) e = SQRT(pm0**2 + pmass**2) if( ivtyp .eq. 1 ) then if( e .lt. val ) then iflag = -15 go to 900 end if else if( e .gt. val ) then iflag = -15 go to 900 end if end if ! else if( ivar .eq. 10 ) then call VMAG(pp,pm0) ek = SQRT(pm0**2 + pmass**2) - pmass if( ivtyp .eq. 1 ) then if( ek .lt. val ) then iflag = -15 go to 900 end if else if( ek .gt. val ) then iflag = -15 go to 900 end if end if ! else if( ivar .eq. 11 ) then call VMAG(pp,pm) xprm = pp(1)/pm if( ivtyp .eq. 1 ) then if( xprm .lt. val ) then iflag = -15 go to 900 end if else if( xprm .gt. val ) then iflag = -15 go to 900 end if end if ! else if( ivar .eq. 12 ) then call VMAG(pp,pm) yprm = pp(2)/pm if( ivtyp .eq. 1 ) then if( yprm .lt. val ) then iflag = -15 go to 900 end if else if( yprm .gt. val ) then iflag = -15 go to 900 end if end if ! end if ! end of test on ivar ! ! --------------------------------------------------------- else if( jpsr .eq. -14 ) then ! fringe field kick and other edge effects grtype = matps model = ival1ps p1 = val1ps p2 = val2ps p3 = val3ps p4 = val4ps ! call VMAG(pp,pm) ! ........................................................ if( grtype .eq. 'SOL' ) then bs = p1 rp = SQRT( xp(1)**2 + xp(2)**2 ) dpf = 0.5d0 * ee * bs * rp * SIGN(1.d0,DBLE(iptyp)) if( rp .gt. 0.d0 ) then phi = ATAN2( xp(2),xp(1) ) else phi = 0.d0 end if dpx = -dpf * SIN(phi) dpy = dpf * COS(phi) pp(1) = pp(1) + dpx pp(2) = pp(2) + dpy ! ! ........................................................ else if( grtype .eq. 'DIP' ) then x = xp(1) y = xp(2) px = pp(1) py = pp(2) q = SIGN(1.d0,pcharge) n = 1 bxar(1) = 0.d0 byar(1) = p1 ! [T] ! call HRDEND(x,y,tp,px,py,pm,q,pmass,n,bxar,byar) ! xp(1) = x xp(2) = y pp(1) = px pp(2) = py ! ! ........................................................ else if( grtype .eq. 'HDIP' ) then x = xp(1) y = xp(2) px = pp(1) py = pp(2) q = SIGN(1.d0,pcharge) n = 1 bxar(1) = p1 ! [T] byar(1) = 0.d0 ! call HRDEND(x,y,tp,px,py,pm,q,pmass,n,bxar,byar) ! xp(1) = x xp(2) = y pp(1) = px pp(2) = py ! ! ........................................................ else if( grtype .eq. 'QUAD' ) then x = xp(1) y = xp(2) px = pp(1) py = pp(2) q = SIGN(1.d0,pcharge) n = 2 bxar(1) = 0.d0 bxar(2) = 0.d0 byar(1) = 0.d0 byar(2) = p1 ! [T/m] ! call HRDEND(x,y,tp,px,py,pm,q,pmass,n,bxar,byar) ! xp(1) = x xp(2) = y pp(1) = px pp(2) = py ! ! ........................................................ else if( grtype .eq. 'SQUA' ) then x = xp(1) y = xp(2) px = pp(1) py = pp(2) q = SIGN(1.d0,pcharge) n = 2 bxar(1) = 0.d0 bxar(2) = p1 ! [T/m] byar(1) = 0.d0 byar(2) = 0.d0 ! call HRDEND(x,y,tp,px,py,pm,q,pmass,n,bxar,byar) ! xp(1) = x xp(2) = y pp(1) = px pp(2) = py ! ! ........................................................ else if( grtype .eq. 'SEX' ) then x = xp(1) y = xp(2) px = pp(1) py = pp(2) q = SIGN(1.d0,pcharge) n = 3 bxar(1) = 0.d0 bxar(2) = 0.d0 bxar(3) = 0.d0 byar(1) = 0.d0 byar(2) = 0.d0 byar(3) = p1 ! [T/m^2] ! call HRDEND(x,y,tp,px,py,pm,q,pmass,n,bxar,byar) ! xp(1) = x xp(2) = y pp(1) = px pp(2) = py ! ! ........................................................ else if( grtype .eq. 'BSOL' ) then if( p3 .eq. 0.d0 ) then ! entrance face solenoid bs = p1 rp = SQRT( xp(1)**2 + xp(2)**2 ) dpf = 0.50d0 * ee * bs * rp * SIGN(1.d0,DBLE(iptyp)) if( rp .gt. 0.d0 ) then phi = ATAN2( xp(2),xp(1) ) else phi = 0.d0 end if dpx = -dpf * SIN(phi) dpy = dpf * COS(phi) pp(1) = pp(1) + dpx pp(2) = pp(2) + dpy end if ! x = xp(1) ! dipole y = xp(2) px = pp(1) py = pp(2) q = SIGN(1.d0,pcharge) n = 1 bxar(1) = 0.d0 byar(1) = p2 ! [T] ! call HRDEND(x,y,tp,px,py,pm,q,pmass,n,bxar,byar) ! xp(1) = x xp(2) = y pp(1) = px pp(2) = py ! if( p3 .eq. 1.d0 ) then ! exit face solenoid bs = p1 rp = SQRT( xp(1)**2 + xp(2)**2 ) dpf = 0.50d0 * ee * bs * rp * SIGN(1.d0,DBLE(iptyp)) if( rp .gt. 0.d0 ) then phi = ATAN2( xp(2),xp(1) ) else phi = 0.d0 end if dpx = -dpf * SIN(phi) dpy = dpf * COS(phi) pp(1) = pp(1) + dpx pp(2) = pp(2) + dpy end if ! ! ........................................................ else if( grtype .eq. 'FACE' ) then ! p1: pole face rotation angle; p2: radius of curvature ! p3: flag to turn on momentum correction ! p4: flag to turn on horizontal focusing pref = ppref(3) if( pref .ne. 0.d0 ) then del = (pm - pref) / pref else del = 0.d0 end if fac = TAN(p1*pi/180.d0) / p2 if( p3 .ne. 0.d0 ) fac = fac / (1.d0 + del) if( p4 .ne. 0.d0 ) then xprm = fac*xp(1) + pp(1)/pp(3) pp(1) = xprm * pp(3) end if yprm = -fac*xp(2) + pp(2)/pp(3) pp(2) = yprm * pp(3) ! else if( grtype .eq. 'DIP3' ) then ! edge effects for hard-edge radial-sector negative-bend ! p1: angle between reference orbit and pole face [deg] ! p2: hard-edge dipole field strength [T] ! p3: entrance (1) or exit (2) flag phi = p1 * pi / 180.d0 cphi = COS(phi) sphi = SIN(phi) dsffag = 2.d0 * sphi / slen ! arclength correction factor x0 = xp(1) ! if( p3 .eq. 1.d0 ) then ! entrance xp(1) = x0 * cphi ! coordinate system rotation pp(1) = pp(1) * cphi pp(1) = pp(1) + pcharge*p2*x0*sphi ! focus correction else ! exit pp(1) = pp(1) + pcharge*p2*x0*sphi ! focus correction xp(1) = x0 * cphi ! coordinate system rotation pp(1) = pp(1) * cphi end if ! end if ! end of test on edge type ! ........................................................ arg = pm**2 - pp(1)**2 - pp(2)**2 if( arg .gt. 0.d0 ) then pp(3) = SQRT(arg) else iflag = -86 go to 900 end if ! ! --------------------------------------------------------- else if( jpsr.eq.-15 ) then !.and. ! grid definition ! & ( jpar.eq.-1 .or. ! & ( jpar.eq.0 .and. ref2par.eq.0 ) .or. ! & ( jpar.eq.1 .and. refpar.eq.0) ) ) then jgrid = ival1ps grtype = matps call GRID(jgrid,grtype,fparps) ! ! --------------------------------------------------------- else if( jpsr .eq. -16 ) then ! reset particle times if( ip .gt. 0 ) then tp = tpref end if ! ! --------------------------------------------------------- else if( jpsr .eq. -17 ) then ! change variable ivar = ival1ps val = val1ps iclass = ival2ps if( iclass.eq.0 .or. (iclass.eq.1 .and. ip.le.0) .or. & (iclass.eq.2 .and. ip.gt.0) ) then ! if( ivar .le. 3 ) then xp(ivar) = xp(ivar) + val ! else if( ivar.ge.4 .and. ivar.le.6 ) then pp(ivar-3) = pp(ivar-3) + val ! else if( ivar .eq. 7 ) then tp = tp + val ! else if( ivar .eq. 8 ) then call VMAG(pp,pm0) pm = pm0 + val do i=1,3 pp(i) = pp(i) * pm / pm0 end do ! else if( ivar .eq. 9 ) then call VMAG(pp,pm0) e = SQRT(pm0**2 + pmass**2) + val pm = SQRT(e*e - pmass**2) do i=1,3 pp(i) = pp(i) * pm / pm0 end do ! else if( ivar .eq. 10 ) then call VMAG(pp,pm0) ek = SQRT(pm0**2 + pmass**2) - pmass + val e = ek + pmass pm = SQRT(e*e - pmass**2) do i=1,3 pp(i) = pp(i) * pm / pm0 end do ! else if( ivar .eq. 11 ) then call VMAG(pp,pm) xprm = val + pp(1)/pm pp(1) = xprm * pm arg = pm**2 - pp(1)**2 - pp(2)**2 if( arg .gt. 0.d0 ) then pp(3) = SQRT(arg) else iflag = -86 go to 900 end if ! else if( ivar .eq. 12 ) then call VMAG(pp,pm) yprm = val + pp(2)/pm pp(2) = yprm * pm arg = pm**2 - pp(1)**2 - pp(2)**2 if( arg .gt. 0.d0 ) then pp(3) = SQRT(arg) else iflag = -86 go to 900 end if ! ! else if( ivar .eq. 99 ) then ! patch for HK scaling dipole work ! n.b. ROTATE does this transformation, but also changes x and z ! cv = cos(val*pi/180.) ! sv = sin(val*pi/180.) ! px = cv*pp(1) + sv*pp(3) ! pz = -sv*pp(1) + cv*pp(3) ! pp(1) = px ! if( pz .gt. 0. ) then ! pp(3) = pz ! else ! iflag = -86 ! go to 900 ! end if ! end if ! end of test on ivar end if ! end of test on iclass ! --------------------------------------------------------- ! else if( jpsr .eq. -18 ) then ! change material density ! mat = matps ! fact = val1ps ! do i=1,nmat ! if( mat .eq. matid(i) ) then ! matdat(2,i) = matdat(2,i) * fact ! density ! if( fact .ne. 0. ) then ! matdat(4,i) = matdat(4,i) / fact ! radiation length ! end if ! end if ! end do ! --------------------------------------------------------- else if( jpsr .eq. -19 ) then ! random field kick if( ip .le. 0 ) go to 900 call VMAG(pp,pm) dpx = 0.d0 dpy = 0.d0 mat = rkmat(jrkval) fact = rkval(jrkval) azm = rkphi(jrkval) rkaz = rkazm(jrkval) ! if( mat .eq. 'SOL' ) then ! length = rkaz alp = 0.5d0*ee*fact / pp(3) c = COS(alp*rkaz) s = SIN(alp*rkaz) xpr0 = pp(1) / pp(3) ypr0 = pp(2) / pp(3) x0 = xp(1) y0 = xp(2) xp(1) = c*c*x0 + s*c/alp*xpr0 + s*c*y0 + s*s/alp*ypr0 xpr = -alp*s*c*x0 + c*c*xpr0 - alp*s*s*y0 + s*c*ypr0 xp(2) = -s*c*x0 - s*s/alp*xpr0 + c*c*y0 + s*c/alp*ypr0 ypr = alp*s*s*x0 - s*c*xpr0 - alp*s*c*y0 + c*c*ypr0 dpx = xpr*pp(3) - pp(1) dpy = ypr*pp(3) - pp(2) ! else ! not solenoid if( azm .gt. -0.5d0 ) then ! use fixed azimuth for kick direction if( mat .eq. 'DIP' ) then dpx = -pcharge * fact else if( mat .eq. 'HDIP' ) then dpy = pcharge * fact else if( mat .eq. 'QUAD' ) then dpx = -pcharge * xp(1) * fact dpy = pcharge * xp(2) * fact else if( mat .eq. 'SQUA' ) then dpx = pcharge * xp(2) * fact dpy = pcharge * xp(1) * fact else if( mat .eq. 'SEX' ) then dpx = pcharge * (xp(2)**2-xp(1)**2) * fact dpy = 2.d0*pcharge * xp(1)*xp(2) * fact end if ! end of test of field type ! else ! use random azimuthal direction if( mat.eq.'DIP' .or. mat.eq.'HDIP' ) then dpx = -pcharge * fact * SIN(rkaz) dpy = pcharge * fact * COS(rkaz) else if( mat.eq.'QUAD' .or. mat.eq.'SQUA' ) then dpx = -pcharge * fact & * (xp(1)*SIN(2.d0*rkaz)-xp(2)*COS(2.d0*rkaz)) dpy = pcharge * fact & * (xp(2)*SIN(2.d0*rkaz)+xp(1)*COS(2.d0*rkaz)) else if( mat .eq. 'SEX' ) then dpx = -pcharge * fact & * ( (xp(1)**2-xp(2)**2)*SIN(3.d0*rkaz) & -2.*xp(1)*xp(2)*COS(3.d0*rkaz) ) dpy = pcharge * fact & * ( 2.*xp(1)*xp(2)*SIN(3.d0*rkaz) & +(xp(1)**2-xp(2)**2)*COS(3.d0*rkaz) ) ! end if end if ! end of test on RKAZ > -0.5 end if ! end of test on mat=SOL ! pp(1) = pp(1) + dpx pp(2) = pp(2) + dpy arg = pm**2 - pp(1)**2 - pp(2)**2 if( arg .gt. 0.d0 ) then pp(3) = SQRT(arg) else iflag = -14 go to 900 end if ! end if ! end of test on jsec ! ! --------------------------------------------------------- 900 continue return end c ********************************************************* function RANV(xlo,xhi) ! ! Return uniform distribution of variable over specified range ! implicit double precision (a-h,o-z) include 'icommon.inc' ! rn = RAN1(irnarg) RANV = (xhi-xlo)*rn + xlo ! return end ! ********************************************************* subroutine READ_EVT(ip) ! ! Read particle data from common arrays or unit(8) ! implicit double precision (a-h,o-z) include 'icommon.inc' ! if( ip .le. mxpar ) then if( ip .eq. 0 ) then ! reference particle do i=1,3 xp(i) = xpref(i) pp(i) = ppref(i) pol(i) = 0.d0 xpint(i) = 0d0 ppint(i) = 0d0 polinit(i) = 0d0 end do ppint(3) = ppref(3) tp = tpref sarc = sarcref ievt = 0 iptyp = iptypref ipflg = ipflgref evtwt = 1.d0 if( iptyp .ne. 0. ) then pmass = mass( ABS(iptypref) ) pcharge = SIGN( ee, DBLE(iptyp) ) * chrg( ABS(iptypref) ) else pmass = 0.d0 end if ! ......................................................... else if( ip .eq. -1 ) then ! 2nd reference particle do i=1,3 xp(i) = xp2ref(i) pp(i) = pp2ref(i) pol(i) = 0.d0 xpint(i) = 0d0 ppint(i) = 0d0 polinit(i) = 0d0 end do ppint(3) = ppref(3) tp = tp2ref sarc = sarc2ref ievt = -1 iptyp = iptyp2ref ipflg = ipflg2ref evtwt = 1.d0 if( iptyp .ne. 0. ) then pmass = mass( ABS(iptyp) ) pcharge = SIGN( ee, DBLE(iptyp) ) * chrg( ABS(iptyp) ) else pmass = 0.d0 end if ! ......................................................... else ! normal particles do i=1,3 xp(i) = xpar(i,ip) pp(i) = ppar(i,ip) pol(i) = polpar(i,ip) xpint(i) = xpintpar(i,ip) ppint(i) = ppintpar(i,ip) polinit(i) = polinitpar(i,ip) end do ! tp = tpar(ip) tpint = tpintpar(ip) sarc = sarcpar(ip) bfld(3) = bzfldpar(ip) ievt = ievtpar(ip) ipnum = ipnumpar(ip) iptyp = iptypar(ip) ipflg = ipflgpar(ip) evtwt = evtwtpar(ip) pmass = mass( ABS(iptyp) ) pcharge = SIGN( ee, DBLE(iptyp) ) * chrg( ABS(iptyp) ) end if ! end of test on ip = 0 ! else ! ip > mxpar read(8,rec=ip-mxpar) xp,pp,tp,pol,sarc,evtwt, & ievt,ipnum,iptyp,ipflg,bfld(3), & xpint,ppint,tpint,polinit pmass = mass( ABS(iptyp) ) pcharge = SIGN( ee, DBLE(iptyp) ) * chrg( ABS(iptyp) ) end if ! return end c ********************************************************* subroutine REG_SUMMARY ! ! Print summary of regions ! Fill tables for RKICK, NEIGHBOR ! ! nrlines: total # of region lines ! nplines: number of pre-looping lines ! implicit double precision (a-h,o-z) include 'icommon.inc' ! dimension fpar(15) logical*1 cell1,rep1,repincell,bkgon character*80 probtitle character*4 cmd,ctag,grtype,mat data rkscmn/1./,rkscsd/1./ data rkickt/0./,rkickl/0./,rphi/0./ ! if( summary ) then write(7,358) 358 format(t1,20('*'),t30,'REGION SUMMARY',t58,20('*')) write(7,360) 360 format(t3, 'IZ', t7, 'IS', t11, 'IC', t14, 'CTAG', t19, 'BTAG', + t26, 'SLEN', t39, 'ZLO', t52, 'ZHI', t62, 'IR', t65, 'FTAG', + t70, 'MTAG', t75, 'MGEOM') end if ! i7 = 0 cell1 = .false. ! set true first time thru cell rep1 = .false. ! set true first time thru repetition repincell = .false. ! set true if there is a repeat in a cell nclines = 0 nreplines = 0 jcel = 0 kcel = 0 cftag = 'NONE' ft = 'NONE' bftag = 'NONE' bkgon = .false. ipstart = 1 ngoodp = 0 nrkval = 0 jrkval = 0 nregnb = 0 shift = 0. ! if( .not.bgen ) then ! determine shift from beam restart read(3,'(a80)',iostat=ioc) probtitle if( ioc .ne. 0 ) then write(2,*)'Error reading particle info from FOR003.DAT' iflag = -1 go to 900 end if read(3,*) zbref,pbref,tbref,zb2ref,pb2ref,tb2ref,typref if( ioc .ne. 0 ) then write(2,*) 'Cant read beam reference particle FOR003.DAT' iflag = -1 go to 900 end if shift = zbref rewind(3) end if zlo = shift zlostart = shift ! -------------------------------------------------------- ! do 450 is=1,nsections jsec = is ss0 = zlo ! FIX THIS IF WE USE IT ! ! Rewind input to top of region data in the input file if( is .eq. 1 ) then nlines = nrlines else ! we have already read the list before if( is .eq. 2 ) nreg = i7 - i7b ! save # of regions in a section if( nplines .gt. 0 ) then nlines = nrlines - nplines else write(2,*) 'Error in REG_SUMMARY (missing BEGS command?)' if(termout) & write(*,*) 'Error in REG_SUMMARY (missing BEGS command?)' iflag = -1 go to 900 end if end if ! do i=1,nlines backspace (unit=1) end do c ......................................................... 400 continue ! loop over region commands c read(1,'(a4)',iostat=ioc)cmd ! ........................................................ if( cmd .eq. 'CELL' )then read(1,*,iostat=ioc) ncells read(1,*) cellflip read(1,'(a4)',iostat=ioc) cftag read(1,*,iostat=ioc)(cfparm(j),j=1,nfparm) jcel = jcel + 1 ! jcel counts new CELL boundaries kcel = 1 ! kcel counts individual cells inside CELL boundary ! Thus if kcel .ne. 0, then we are currently inside a CELL boundary sc0 = zlo ! initial position of this CELL cellslp = 1. ! cell field symmetry cell1 = .true. if( summary ) write(7,'(t20,a,i5)')' CELLS = ',ncells ! ! ........................................................ else if( cmd .eq. 'REPE' )then read(1,*,iostat=ioc) nrep rep1 = .true. if( cell1 ) then repincell = .true. nclines = nclines + 2 end if ! ! ........................................................ else if( cmd .eq. 'SECT' )then i7 = i7 + 1 if( summary ) write(7,'(i4,t20,a)')i7,'PRODUCTION' ! ! ........................................................ else if( cmd .eq. 'BEGS' )then if( summary ) write(7,'(t20,a)') 'BEGIN SECTION LOOPING' i7b = i7 ! ! ........................................................ else if( cmd .eq. 'ROTA' )then read(1,*,iostat=ioc) rotang,iaxis,iclass i7 = i7 + 1 if( summary ) then write(7,'(i4,t20,a,f12.4,a,i3)')i7,'ROTATION =',rotang, & ' AXIS = ',iaxis end if if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 ! ! ........................................................ else if( cmd .eq. 'TAPE' )then read(1,*,iostat=ioc) idum i7 = i7 + 1 if( summary ) then write(7,'(i4,t20,a,i4)')i7,'TAPER, cells = ',idum end if ! ........................................................ else if( cmd .eq. 'OUTP' )then i7 = i7 + 1 if( summary ) write(7,'(i4,t20,a)')i7,'OUTPUT' if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 1 end if if( rep1 ) nreplines = nreplines + 1 ! ! ........................................................ else if( cmd .eq. 'APER' )then read(1,*,iostat=ioc) iapertype,aperlim1,aperlim2,aperlim3 i7 = i7 + 1 if( summary ) write(7,'(i4,t20,a)')i7,'APERTURE' if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 ! ! ........................................................ else if( cmd .eq. 'BACK' )then read(1,*,iostat=ioc) bentbkg,prefbkg,ztotalbkg read(1,*,iostat=ioc) xlogr,dxgr,nxgr, & ylogr,dygr,nygr,slogr,dsgr,nsgr,interbkg i7 = i7 + 1 if( summary ) write(7,'(i4,t20,a)') i7,'BACKGROUND FIELD' bkgon = .true. saccbkg = 0. ! ! ........................................................ else if( cmd .eq. 'BFIE' )then ! background field read(1,*,iostat=ioc) zoffbkg,rmaxbkg,zminbkg,zmaxbkg read(1,'(a4)',iostat=ioc) bftag read(1,*,iostat=ioc)(bfparm(j),j=1,nfparm) i7 = i7 + 1 if( summary ) write(7,'(i4,t25,a)') i7,bftag ! ! ........................................................ else if( cmd .eq. 'ENDB' )then read(1,*,iostat=ioc) mfile i7 = i7 + 1 if( summary ) write(7,'(i4,t20,a)') i7,'END BACKGROUND FIELD' ! ! ........................................................ else if( cmd .eq. 'TRAN' )then read(1,*,iostat=ioc) ((trprt(i,j),j=1,6),i=1,6) i7 = i7 + 1 jsrg = i7 zhi = zlo if( summary ) write(7,402) i7,'USER TRANSPORT' 402 format(i4,t20,a) ! zlo = zhi if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 7 end if if( rep1 ) nreplines = nreplines + 7 ! ! ........................................................ else if( cmd .eq. 'REFP' )then ! read(1,*,iostat=ioc) refpar,pz0ref,t0ref,gradref,phmodref read(1,*) i7 = i7 + 1 if( summary ) write(7,'(i4,t20,a)')i7,'REF PARTICLE' ipstart = 0 ngoodp = -1 if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 ! ! ........................................................ else if( cmd .eq. 'DUMM' )then i7 = i7 + 1 if( summary ) write(7,'(i4,t20,a)')i7,'DUMMY' if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 1 end if if( rep1 ) nreplines = nreplines + 1 ! ! ........................................................ else if( cmd .eq. 'TILT' )then read(1,*,iostat=ioc) psi,phi i7 = i7 + 1 if( summary ) write(7,'(i4,t20,a,2f12.4)')i7,'TILT',psi,phi if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 ! ! ........................................................ else if( cmd .eq. 'DISP' )then read(1,*,iostat=ioc) d,phi i7 = i7 + 1 if( summary ) write(7,'(i4,t20,a,2f12.4)')i7,'DISP',d,phi if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 ! ! ........................................................ else if( cmd .eq. 'REF2' )then ! read(1,*,iostat=ioc) ref2par,pz2ref,t2ref,grad2ref read(1,*) i7 = i7 + 1 if( summary ) write(7,'(i4,t20,a)')i7,'2nd REF PARTICLE' ipstart = -1 ngoodp = -2 if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 ! ! ........................................................ else if( cmd .eq. 'CUTV' )then read(1,*,iostat=ioc) ivar,ivtyp,val i7 = i7 + 1 if( summary ) then write(7,'(i4,t6,a4,1x,2i4,1x,f9.5)')i7, & 'CUTV',ivar,ivtyp,val end if ! sr0 = zlo if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 ! ! ........................................................ else if( cmd .eq. 'EDGE' )then read(1,'(a4)',iostat=ioc) grtype read(1,*,iostat=ioc) model,p1,p2,p3,p4 i7 = i7 + 1 if( summary ) then write(7,'(i4,t20,a,a4,i4,4f11.5)')i7,'EDGE: ',grtype, & model,p1,p2,p3,p4 end if if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 3 end if if( rep1 ) nreplines = nreplines + 3 ! ! ........................................................ else if( cmd .eq. 'GRID' )then read(1,*,iostat=ioc) jgrid read(1,'(a4)',iostat=ioc) grtype read(1,*,iostat=ioc) fpar i7 = i7 + 1 if( summary ) then write(7,'(i4,t20,a,i4,5x,a4)')i7,'GRID: ',jgrid,grtype end if ! ! ........................................................ else if( cmd .eq. 'RESE' )then i7 = i7 + 1 if( summary ) write(7,'(i4,t20,a,i4,5x,a4)')i7,'RESET TIME ' if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 1 end if if( rep1 ) nreplines = nreplines + 1 ! ! ........................................................ else if( cmd .eq. 'DVAR' )then read(1,*,iostat=ioc) ivar,val,iclass i7 = i7 + 1 if( ivar .eq. 3 ) then slen = val else slen = 0. end if zhi = zlo + slen if( summary ) then write(7,'(i4,t6,a4,1x,i2,1x,f9.5,t23,3e13.5)')i7, & 'DVAR',ivar,val,slen,zlo,zhi end if zlo = zhi ! sr0 = zlo if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 ! ! ........................................................ else if( cmd .eq. 'DENS' )then read(1,*,iostat=ioc) mat,fact i7 = i7 + 1 if( summary ) then write(7,'(i4,t20,a,1x,a4,f12.4)')i7,'DENSITY ',mat,fact end if if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 ! ! --------------------------------------------------------- else if( cmd .eq. 'RKIC' )then read(1,*,iostat=ioc) mat,xm,sd,coup,azflg i7 = i7 + 1 icoup = coup ! if( is .eq. 1 ) then ! get random values for first section if( mat .eq. 'SCAL' ) then ! set scale factors rkscmn = xm rkscsd = sd else ! make or use kicks if( icoup .eq. 0 ) then ! generate random value d = GAUSS(xm*rkscmn,sd*rkscsd) else if( icoup .eq. 1 ) then ! generate random value and save it d = GAUSS(xm*rkscmn,sd*rkscsd) if( mat .eq. 'SOL' ) then rkickl = d else rkickt = d end if else if( icoup .eq. 2 ) then ! use saved random value if( mat .eq. 'SOL' ) then d = rkickl else d = rkickt end if else if( icoup .eq. 3 ) then ! use saved value with reverse sign if( mat .eq. 'SOL' ) then d = -rkickl else d = -rkickt end if end if ! if( azflg .lt. -0.5 ) then ! use random azimuthal kick direction if( icoup .eq. 0 ) then ! generate random value rkaz = RANV(0.d0,twopi) else if( icoup .eq. 1 ) then ! generate random value and save rkaz = RANV(0.d0,twopi) rphi = rkaz else if( icoup .eq. 2 ) then ! use saved random value rkaz = rphi else if( icoup .eq. 3 ) then ! use saved value with rev sign rkaz = rphi ! magnitude is reversed, so use same azimuth end if ! end of test on ICOUP else rkaz = 0.d0 end if ! end of test on AZFLG < -0.5 ! if( mat .eq. 'SOL' ) rkaz = azflg ! save solenoid length ! nrkval = nrkval + 1 if( nrkval .gt. 100 ) then write(2,*) 'Error: requested more than 100 RKICK commands' write(*,*) 'Error: requested more than 100 RKICK commands' stop else rkmat(nrkval) = mat rkval(nrkval) = d rkphi(nrkval) = azflg rkazm(nrkval) = rkaz end if ! end if ! end of test on MAT = SCAL ! end if ! end of test on IS = 1 ! if( summary ) then if( mat .eq. 'SCAL' ) then write(7,'(i4,t20,a,1x,a4,2f14.6)')i7,'RKICK ',mat,xm,sd else jrkval = jrkval + 1 if( jrkval .gt. nrkval ) jrkval = 1 v1 = rkval(jrkval) v2 = rkazm(jrkval) if( mat .ne. 'SOL' ) v2 = v2*180./pi write(7,'(i4,t20,a,1x,a4,f14.6,2x,f8.2)')i7, & 'RKICK ',mat,v1,v2 end if end if ! if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 ! ! ........................................................ else if( cmd .eq. 'DENP' )then read(1,*,iostat=ioc) mat,idir,p1,p2,p3,p4 i7 = i7 + 1 if( summary ) then write(7,'(i4,t20,a,1x,a4,i4)') & i7,'DENSITY PROFILE',mat,idir end if if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 ! ! --------------------------------------------------------- else if( cmd .eq. 'SREG' )then read(1,*,iostat=ioc) slen,nrreg,zstep i7 = i7 + 1 jsrg = i7 sr0 = zlo zhi = zlo + slen ! if( bkgon ) then if( saccbkg .le. slogr+ztotalbkg ) then bftag = 'BACK' saccbkg = saccbkg + slen else bftag = 'NONE' bkgon = .false. end if end if ! ctag = cftag if( cftag.eq.'BLOC' .and. cfparm(1).eq.3 ) then ig = cfparm(2) ctag(3:3) = 'g' write(ctag(4:4),'(i1)') ig end if ! if( cftag.eq.'COIL' .and. cfparm(1).eq.3 ) then ig = cfparm(2) ctag(3:3) = 'g' write(ctag(4:4),'(i1)') ig end if ! if( cftag .eq. 'SHEE' ) then if( cfparm(1) .eq. 3 ) then mfile = cfparm(2) if( mfile .gt. 0 ) then write(ctag(3:4),'(i2)') mfile else write(ctag(2:4),'(i3)') mfile end if else if( cfparm(1) .eq. 5 ) then ig = cfparm(2) ctag(3:3) = 'g' write(ctag(4:4),'(i1)') ig end if end if ! do 430 ir=1,nrreg read(1,*,iostat=ioc) idm,rlow(ir),rhigh(ir) read(1,'(a4)',iostat=ioc) ftag(ir) read(1,*,iostat=ioc)(fparm(j,ir),j=1,nfparm) read(1,'(3a4)',iostat=ioc) mtag(ir),mtag2(ir),mtag3(ir) read(1,'(a6)',iostat=ioc) mgeom(ir) read(1,*,iostat=ioc) (gparm(j,ir),j=1,ngparm) ! if( summary ) then write(7,405) i7,is,kcel,ctag,bftag,slen,zlo,zhi,ir,ftag(ir), & mtag(ir),mgeom(ir) 405 format(i5, t7, i2, t10, i3, t14, a4, t19, a4, t23, 3e13.5, + t63, i1, t65, a4, t70, a4, t75, a6) end if ! 430 continue ! ! Check number of curved fields present ! npres = 0 ! do i=1,ncrvfld ! do j=1,nrreg ! if( ftag(j) .eq. crvfld(i) ) npres = npres + 1 ! end do ! if( cftag .eq. crvfld(i) ) npres = npres + 1 ! end do ! ! do j=1,nrreg ! if( ftag(j).eq.'STUS' .and. fparm(3,j).eq.1 )npres = npres + 1 ! end do ! if( cftag.eq.'STUS' .and. cfparm(3).eq.1 ) npres = npres + 1 ! ! if( bftag.eq.'BACK' .and. bentbkg ) npres = npres + 1 ! ! if( npres .gt. 1 ) then ! write(2,*) 'Cant use more than 1 curved field type' ! iflag = -11 ! go to 900 ! end if ! if( neighbor ) then ! save region info for neighbor logic nregnb = nregnb + 1 if( nregnb .gt. 100 ) then write(2,*) 'Error: requested > 100 NEIGHBOR regions' write(*,*) 'Error: requested > 100 NEIGHBOR regions' stop else i7regnb(nregnb) = i7 sr0nb(nregnb) = zlo tagnb(nregnb) = ftag(1) do j=1,15 fparnb(j,nregnb) = fparm(j,1) end do end if end if ! zlo = zhi ! if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 6 * nrreg + 2 end if ! if( rep1 ) nreplines = nreplines + 6 * nrreg + 2 ! ! ........................................................ else if( cmd .eq. 'ENDR' )then if( rep1 ) nreplines = nreplines + 1 if( cell1 .and. rep1 ) nclines = nclines + 1 rep1 = .false. ! stop incrementing nreplines nrep = nrep - 1 if( nrep .gt. 0 ) then ! more repetitions left c do i=1,nreplines backspace (unit = 1) end do c else ! done until start of any new repetitions nreplines = 0 repincell = .false. end if ! ! ........................................................ else if( cmd .eq. 'ENDC' )then if( cell1 ) nclines = nclines + 1 cell1 = .false. ! stop incrementing nclines ncells = ncells - 1 ! if( ncells .gt. 0 ) then ! more cells left kcel = kcel + 1 sc0 = zlo ! reset SC0 for next cell if( cellflip ) cellslp = -cellslp c do i=1,nclines backspace (unit = 1) end do if( summary ) write(7,*) c else ! done until start of any new cells if( summary ) write(7,'(t20,a)') 'END CELL' nclines = 0 kcel = 0 cftag = 'NONE' end if ! end of test on ncells ! ! ........................................................ else if( cmd .eq. 'ENDS' )then ! write(2,*) go to 450 end if ! ........................................................ go to 400 ! get next region command c 450 continue ! end of loop over sections ! ! -------------------------------------------------------- ! write(2,'(78(''*''))') n7records = i7 ! 900 continue return end ! ******************************************************** subroutine RF_BUCKET(jj,kk) ! ! Sample initial longitudinal phase space from an rf bucket ! ! Ref. Weidemann, Vol 1, Chap. 8 ! implicit double precision (a-h,o-z) include 'icommon.inc' ! external fpsih ! data rad/57.29578d0/ ! emag = corr1(jj,kk) / 1.d3 ! field into GV/m psis = corr2(jj,kk) / rad ! psis into radians frf = corr3(jj,kk) * 1.d6 ! freq into Hz ic = corrtyp(jj,kk) ! correlation type ibd = bdistyp(kk) ! beam distribution type ! gamma = SQRT( 1.d0 + (pp(3)/pmass)**2 ) beta = SQRT( 1.d0 - (1.d0/gamma)**2 ) wl = cc / frf omsyn = cc/gamma * SQRT( twopi*beta*emag/pp(3)/wl ) ! if( ic .eq. 3 ) then ! small amplitude, elliptical amp = omsyn*SQRT(COS(psis))/twopi/frf*gamma**2 ! if( ibd .eq. 1 ) then ! matched,bi-gaussian sdph = x2bt(3,kk) / wl * twopi sdpz = amp * sdph * pp(3) pp(3) = GAUSS( p1bt(3,kk), sdpz ) ! else ! uniform ellipse phi = xp(3) / wl * twopi dphimx = (x2bt(3,kk) - x1bt(3,kk)) / 2.d0 / wl * twopi del = amp * SQRT( dphimx**2 - phi**2 ) pzlo = p1bt(3,kk) * (1.d0 - del) pzhi = p1bt(3,kk) * (1.d0 + del) pp(3) = RANV( pzlo,pzhi ) end if ! else if( ic .eq. 4 ) then ! small amplitude, separatrix amp = 1.414d0*omsyn*SQRT(COS(psis))/twopi/frf*gamma**2 phi = xp(3) / wl * twopi ! if( ibd .eq. 1 ) then ! gaussian distribution in z ! Divide allowed momentum range into 4 sigmas sdpz = amp * pp(3) * SQRT( 1.d0 + COS(phi) ) / 4.d0 ! 80 pp(3) = GAUSS( p1bt(3,kk), sdpz ) if( ABS(pp(3)-p1bt(3,kk)).gt.4.d0*sdpz ) go to 80 ! else ! uniform phi = RANV(-piov2,piov2) xp(3) = phi / twopi * wl del = amp * SQRT( 1.d0 + COS(phi) ) pzlo = p1bt(3,kk) * (1.d0 - del) pzhi = p1bt(3,kk) * (1.d0 + del) pp(3) = RANV( pzlo,pzhi ) end if ! else if( ic .eq. 5 ) then ! large amplitude, separatrix amp = 1.414d0*omsyn/twopi/frf*gamma**2 psit = pi - psis ! ! Get psih from solution of transcendental equation ! Use Newton-Raphson method if( psis .ge. 0. ) then ! accelerating psih = RTSAFE( FPSIH,-pi,psis,1d-4 ) else psih = RTSAFE( FPSIH,psit,3.d0*pi,1d-4 ) end if ! ! Sample the phase if( ibd .eq. 1 ) then ! "gaussian" distribution ! Use separate gaussians on each side of psis ! Divide allowed range into 4 sigmas rn = RAN1(irnarg) if( rn .lt. 0.5d0 ) then ! tail side sd = ABS(psis - psit) / 4.d0 ! 100 psi = GAUSS(psis,sd) if( psis .ge. 0.d0 ) then if( psi.lt.psis .or. ABS(psi-psis).gt.4.d0*sd )go to 100 else if( psi.gt.psis .or. ABS(psi-psis).gt.4.d0*sd )go to 100 end if ! else ! head side sd = ABS(psih - psis) / 4.d0 ! 110 psi = GAUSS(psis,sd) if( psis .ge. 0.d0 ) then if( psi.gt.psis .or. ABS(psi-psis).gt.4.d0*sd )go to 110 else if( psi.lt.psis .or. ABS(psi-psis).gt.4.d0*sd )go to 110 end if end if ! end of test on rn (head-tail) ! ......................................................... else ! uniform distribution psi = RANV(psit,psih) end if ! phi = psi - psis dph = phi - psis ! Get delta limit for this phase func = COS(psis)+COS(psis+phi)+(2.d0*psis+phi-pi)*SIN(psis) del = amp * SQRT( func ) ! ! Sample the momentum deviation if( ibd .eq. 1 ) then ! gaussian distribution sd = del / 4.d0 * pp(3) ! 150 dp = GAUSS(0.d0,sd) if( ABS(dp) .gt. 4.d0*del ) go to 150 ! else ! uniform distribution dpz = del * pp(3) dp = RANV( -dpz,dpz ) end if ! ! Adjust particle vector xp(3) = xp(3) + dph / twopi * wl pp(3) = pp(3) + dp ! end if ! end of test on correlation type ! return end c ********************************************************* subroutine RKSTEP(y,dydt,n,tt,h,yout) c c Propagate particle one step thru field (or drift) c Use 4th order Runge-Kutta integration c cf. Num Rec routine RK4 c implicit double precision (a-h,o-z) include 'icommon.inc' parameter (nmax = 15) dimension y(n),dydt(n),yout(n),yt(nmax),dyt(nmax),dym(nmax) c h2 = h / 2.d0 h6 = h / 6.d0 th = tt + h2 c c First step ......................................... c do i=1,n ! variables at midpoint yt(i) = y(i) + h2 * dydt(i) end do c c Second step ......................................... call DERIV(th,yt,dyt) if( iflag .ne. 0 ) go to 900 c do i=1,n ! variables at midpoint yt(i) = y(i) + h2 * dyt(i) end do c c Third step ......................................... call DERIV(th,yt,dym) if( iflag .ne. 0 ) go to 900 c do i=1,n ! variables at end yt(i) = y(i) + h * dym(i) dym(i) = dym(i) + dyt(i) end do c c Fourth step ......................................... call DERIV(tt+h,yt,dyt) if( iflag .ne. 0 ) go to 900 c do i=1,n ! variables at end yout(i) = y(i) + h6 * ( dydt(i)+2.d0*dym(i)+dyt(i) ) end do c 900 continue end c ********************************************************* subroutine ROTATE_COORD(ang,iaxis) ! ! Rotate particle coordinates by angle ANG around IAXIS ! Rotations around the x and y axes due to J.S. Berg (20 Oct 2005) ! implicit double precision (a-h,o-z) include 'icommon.inc' ! dimension x(3),p(3) ! ca = COS(ang) sa = SIN(ang) do i=1,3 x(i) = xp(i) p(i) = pp(i) end do ! if( iaxis .eq. 1 ) then pp(2) = p(2)*ca + p(3)*sa pp(3) = -p(2)*sa + p(3)*ca xp(2) = x(2)*p(3)/pp(3) xp(1) = x(1) + x(2)*p(1)/pp(3)*sa tp = tp + x(2)*sa*sqrt(p(1)**2+p(2)**2+p(3)**2+pmass**2)/ $ (pp(3)*cc) ! else if( iaxis .eq. 2) then pp(1) = p(1)*ca + p(3)*sa pp(3) = -p(1)*sa + p(3)*ca xp(1) = x(1)*p(3)/pp(3) xp(2) = x(2) + x(1)*p(2)/pp(3)*sa tp = tp + x(1)*sa*sqrt(p(1)**2+p(2)**2+p(3)**2+pmass**2)/ $ (pp(3)*cc) ! else if( iaxis .eq. 3) then xp(1) = x(1)*ca + x(2)*sa xp(2) = -x(1)*sa + x(2)*ca pp(1) = p(1)*ca + p(2)*sa pp(2) = -p(1)*sa + p(2)*ca end if ! return end c ********************************************************* subroutine SAVE_FILE(izp) c c Save particle information at plane izp ! ( or on error condition ) c implicit double precision (a-h,o-z) include 'icommon.inc' c if( fsav .and. (izp .eq. izfile) ) then if( izp .gt. 0 ) then ! want info at this plane if( fsavset ) then tt = tp - tpref zz = 0. else tt = tp zz = xp(3) end if ! write(4,20) ievt,ipnum,iptyp,ipflg,tt,evtwt, & xp(1),xp(2),zz,pp,pol 20 format(2i8,2i6,2e14.6,9e14.6) ! else ! want initial conditions (also for errors) write(4,20) ievt,ipnum,iptyp,ipflg,tpint,1., & xpint,ppint,polinit end if end if ! return end c ********************************************************* subroutine SETUP ! ! Set up parameters and data files for the simulation ! ! File useage 1 : input parameters ! 2 : simulation LOG file ! 3 : initial particles (input) ! 4 : particles info at specified s-planes ! 5 : (reserve for UNIX compatibility) ! 6 : (reserve for UNIX compatibility) ! 7 : region summary table ! 8 : particle info (internal) ! 9 : "n-tuple" file ! 10 : (reserve for future ICOOL use) ! 11 : beam envelope output ! 12-13 : ELMS ! 14-19 : (reserve for future ICOOL use) ! 20-99 : {coil,sheet,block} segment data; ! : {coil,sheet,block,background} grid data; ! : rf diagnostics;rf phases; neutrino production data; ! : etc... ! * : terminal output ! implicit double precision (a-h,o-z) include 'icommon.inc' c dimension fpar(15) integer*2 iyr,imon,ihr,imin,isec logical cloop,rloop,bgloop character*4 cmd,mat,ftagdat(21) character*4 cloopdat(20),rloopdat(18),bgloopdat(2),grtype character*6 mgeomdat(9) character*10 fname character*80 title ! namelist/BMT/nbeamtyp,bmalt,dtalt namelist/CONT/beep,betaperp,bgen,bunchcut,bzfldprd,dectrk, & diagref, & epsf,epsreq,epsstep,ffcr,forcerp,fsav,fsavset,f9dp, & goodtrack,izfile,magconf,mapdef, & neighbor,neutrino,nnudk,npart,nprnt,npskip,nuthmin, & nuthmax,nsections, & ntuple,output1,phantom,phasemodel, & prlevel,prnmax,pzmintrk,rfdiag,rfphase, & rnseed,rtuple,rtuplen,run_env,scalestep,spin,spinmatter, & spintrk,stepmin, & stepmax,steprk,summary,termout,timelim,varstep ! namelist/INTS/ldedx,lstrag,lscatter,ldray,ldecay, & linteract,lspace,lelms,lsamcs, & delev,straglev,scatlev,draylev,declev,intlev,spacelev, & dcute,dcutm,elmscor, & facfms,facmms,fastdecay, & frfbunsc,gfactsc,parbunsc,pdelev5,rwallsc,trainsc, & wanga,wangb,wangc,wangd,wangpmx,wangfmx namelist/NEM/nemit,discorr,ipzcor,pxycorr,edetail, & sigmacut,sig_cut namelist/NHS/nhist,hauto,hcprn namelist/NRH/nrhist,rhauto,rhprin namelist/NSC/nscat,sauto namelist/NZH/nzhist,zhauto,zhprin namelist/NCV/ncovar c data nftags/21/ data ftagdat/'NONE','ACCE','BLOC','BROD', & 'BSOL','COIL','DIP ','EFLD','FOFO','HDIP', & 'HELI','HORN','KICK','QUAD','ROD ','SEX ','SHEE','SOL ', & 'STUS','SQUA','WIG '/ ! data nmgeoms/9/ data mgeomdat/'CBLOCK','WEDGE','PWEDGE','ASPW','HWIN', & 'ASRW','RING','NIA','NONE '/ ! data ncloop/20/ ! allowed commands in CELL loops data cloopdat/'SREG','REPE','ENDR','OUTP','APER', & 'ROTA','TRAN','DUMM','TILT','EDGE','REF2', & 'DISP','DVAR','RESE','REFP','DENS','RKIC', & 'CUTV','DENP','ENDC'/ ! data nrloop/18/ ! allowed commands in REPEAT loops data rloopdat/'SREG','OUTP','APER','ROTA','TRAN', & 'DUMM','TILT','EDGE','DVAR','RESE','REFP', & 'REF2','DENS','DISP','RKIC','CUTV','DENP', & 'ENDR'/ ! data nbgloop/2/ ! allowed commands in BACKGROUND loops data bgloopdat/'BFIE','ENDB'/ ! data cloop/.false./,rloop/.false./,bgloop/.false./ ! data bzfldprd/0./ c Open required files ! ! Output log open(unit=2, & file=pathname(1:LASTNB(pathname)) // 'for002.dat', & status='unknown',iostat=ioc) if( ioc .ne. 0 ) then if(termout) write(*,*) 'Cant open log file FOR002.DAT' go to 800 end if ! write(2,'(a,a)')' Welcome to ICOOL - version ',version call ICOOLDATE(iyr,imon,iday0) call ICOOLTIME(ihr0,imin0,isec0) write(2,'(a,4i6)')' Date: ',iyr,imon,iday0 write(2,'(a,4i6)')' Start time: ',ihr0,imin0,isec0 ! ! Strip blank lines and comments out of problem command file call STRIPCOM ! ! Replace any named variables in the problem command file call SUBNAMES ! ! Input data open(unit=1, & file=pathname(1:LASTNB(pathname)) // 'for001.dsb', & status='unknown',iostat=ioc) if( ioc .ne. 0 ) then write(2,*) 'Cant find input file FOR001.DSB' go to 800 end if ! c File headers read(1,'(a79)',iostat=ioc) title if( ioc .ne. 0 ) then write(2,*) 'Cant read problem title' go to 800 end if write(2,'(a79)') title ! c Read simulation control parameters read(1,cont,iostat=ioc) if( ioc .ne. 0 ) then write(2,*) 'Cant read control variable namelist' go to 800 end if if(termout) then ! termout is a control variable write(*,'(a,a)')' Welcome to ICOOL - version ',version end if ! ! Output file header open(unit=9, & file=pathname(1:LASTNB(pathname)) // 'for009.dat', & status='unknown',iostat=ioc) if( ioc .ne. 0 ) then write(2,*) 'Cant open postprocessor file FOR009.DAT' go to 800 end if write(9,'(a1,a79)') '#',title write(9,'(a)')'# units = [s] [m] [GeV/c] [T] [V/m] ' ! if( f9dp .eq. 4 ) then write(9,24) 24 format('#',t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg', & t31,'time',t43,'x',t55,'y',t67,'z',t79,'Px',t91,'Py', & t103,'Pz',t115,'Bx',t127,'By',t139,'Bz',t151,'wt', & t163,'Ex',t175,'Ey',t187,'Ez',t199,'arclength',t211,'polX', & t223,'polY',t235,'polZ') else if( f9dp .eq. 6 ) then write(9,26) 26 format('#',t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg', & t31,'time',t45,'x',t59,'y',t73,'z',t87,'Px',t101,'Py', & t115,'Pz',t129,'Bx',t145,'By',t157,'Bz',t171,'wt', & t185,'Ex',t199,'Ey',t213,'Ez',t227,'arclength',t241,'polX', & t255,'polY',t269,'polZ') else if( f9dp .eq. 8 ) then write(9,28) 28 format('#',t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg', & t31,'time',t47,'x',t63,'y',t79,'z',t95,'Px',t111,'Py', & t127,'Pz',t143,'Bx',t159,'By',t175,'Bz',t191,'wt', & t207,'Ex',t223,'Ey',t239,'Ez',t255,'arclength',t271,'polX', & t287,'polY',t303,'polZ') else if( f9dp .eq. 10 ) then write(9,30) 30 format('#',t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg', & t31,'time',t49,'x',t67,'y',t85,'z',t103,'Px',t121,'Py', & t139,'Pz',t157,'Bx',t175,'By',t193,'Bz',t211,'wt', & t229,'Ex',t247,'Ey',t265,'Ez',t283,'arclength',t301,'polX', & t319,'polY',t337,'polZ') else if( f9dp .eq. 12 ) then write(9,32) 32 format('#',t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg', & t31,'time',t51,'x',t71,'y',t91,'z',t111,'Px',t131,'Py', & t151,'Pz',t171,'Bx',t191,'By',t211,'Bz',t231,'wt', & t251,'Ex',t271,'Ey',t291,'Ez',t311,'arclength',t331,'polX', & t351,'polY',t371,'polZ') else if( f9dp .eq. 14 ) then write(9,34) 34 format('#',t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg', & t31,'time',t53,'x',t75,'y',t97,'z',t119,'Px',t141,'Py', & t163,'Pz',t185,'Bx',t207,'By',t229,'Bz',t251,'wt', & t273,'Ex',t295,'Ey',t317,'Ez',t339,'arclength',t361,'polX', & t383,'polY',t405,'polZ') else if( f9dp .eq. 16 ) then write(9,36) 36 format('#',t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg', & t31,'time',t55,'x',t79,'y',t103,'z',t127,'Px',t151,'Py', & t175,'Pz',t199,'Bx',t223,'By',t247,'Bz',t271,'wt', & t295,'Ex',t319,'Ey',t343,'Ez',t367,'arclength',t391,'polX', & t415,'polY',t439,'polZ') else if( f9dp .eq. 17 ) then write(9,37) 37 format('#',t5,'evt',t9,'par',t13,'typ',t18,'flg',t24,'reg', & t31,'time',t56,'x',t81,'y',t106,'z',t131,'Px',t156,'Py', & t181,'Pz',t206,'Bx',t231,'By',t256,'Bz',t281,'wt', & t306,'Ex',t331,'Ey',t356,'Ez',t381,'arclength',t406,'polX', & t431,'polY',t456,'polZ') else write(9,24) ! default end if ! c ENVELOPE EQUATIONS ADDITION if (run_env) then call ENV_SETUP(title,ioc) if( ioc .ne. 0 ) go to 800 end if write(2,cont) c irnarg = rnseed ! ! The following tests must be made after the CONT namelist has been read if( .not. bgen ) then ! existing initial particle values open(unit=3, & file=pathname(1:LASTNB(pathname)) // 'for003.dat', & status='unknown',iostat=ioc) if( ioc .ne. 0 ) then write(2,*) 'Cant open beam data input file FOR003.DAT' go to 800 end if end if ! if( fsav ) then ! new particle data at specified plane open(unit=4, & file=pathname(1:LASTNB(pathname)) // 'for004.dat', & status='unknown',iostat=ioc) if( ioc .ne. 0 ) then write(2,*) 'Cant open beam data output file FOR004.DAT' go to 800 end if write(4,'(a80)') title end if ! ! Region summary table if( summary ) then open(unit=7, & file=pathname(1:LASTNB(pathname)) // 'for007.dat', & status='unknown',iostat=ioc) if( ioc .ne. 0 ) then write(2,*) 'Cant open file FOR007.DAT' go to 800 end if end if ! if( npart .gt. mxpar ) then ! particle info (internal) open(unit=8, & file=pathname(1:LASTNB(pathname)) // 'for008.dat', & access='direct', & form='unformatted',recl=200,status='unknown',iostat=ioc) if( ioc .ne. 0 ) then write(2,*) 'Cant open internal excess particle file FOR008.DAT' go to 800 end if end if ! ! Neutrino production data if( neutrino.gt.19 .and. neutrino.lt.100 ) then write(fname,40) neutrino 40 format('for0',i2,'.dat') open(unit=neutrino, & file=pathname(1:LASTNB(pathname)) // fname, & status='unknown',iostat=ioc) if( ioc .ne. 0 ) then write(2,*) 'Cant open neutrino production file' go to 800 end if write(neutrino,'(a)') 'Neutrino production data' write(neutrino,'(a80)') title end if ! ! Rf diagnostics if( rfdiag.gt.19 .and. rfdiag.lt.100 ) then write(fname,40) rfdiag open(unit=rfdiag, & file=pathname(1:LASTNB(pathname)) // fname, & status='unknown',iostat=ioc) if( ioc .ne. 0 ) then write(2,*) 'Cant open rf diagnostics output file' go to 800 end if write(rfdiag,'(a)') 'RF diagnostics' write(rfdiag,45) 45 format(t3,'JSRG',t10,'IEVT',t19,'DPHASE',t32,'DT',t43,'DPZ', & t56,'FREQ',t69,'GRAD',t82,'PHASE',t95,'EZMN',t108,'BPHIMN', & t121,'TREFMN',t134,'T2REFMN',t147,'ZMN') end if ! ! Rf phases if( rfphase.gt.19 .and. rfphase.lt.100 ) then write(fname,40) rfphase open(unit=rfphase, & file=pathname(1:LASTNB(pathname)) // fname, & status='unknown',iostat=ioc) if( ioc .ne. 0 ) then write(2,*) 'Cant open rf phases input file' go to 800 end if end if if( rtuple ) ntuple = .false. ! c Read beam distribution parameters if( bgen ) then read(1,bmt,iostat=ioc) if( ioc .ne. 0 ) then write(2,*) 'Cant read beam parameter namelist' go to 800 end if write(2,bmt) c if( nbeamtyp.lt.1 .or. nbeamtyp.gt.5 ) then if( ioc .ne. 0 ) then write(2,*) 'Must have 0 < NBEAMTYP < 6' go to 800 end if end if ! do i=1,nbeamtyp read(1,*,iostat=ioc)idm,bmtype(i),fracbt(i),bdistyp(i) if( ioc .ne. 0 ) then write(2,*) 'Error reading IDM,BMTYPE,FRACBT,BDISTYP' go to 800 end if write(2,50)idm,bmtype(i),bdistyp(i),fracbt(i) 50 format(' #,BMTYPE,bdistyp,fracbt: ',3i5,f12.4) ! if( bdistyp(i) .eq. 1 ) then ! gaussian read(1,*,iostat=ioc)(x1bt(j,i),j=1,3), + (p1bt(j,i),j=1,3) if( ioc .ne. 0 ) then write(2,*) 'Error reading mean positions and momenta' go to 800 end if write(2,53)(x1bt(j,i),j=1,3),(p1bt(j,i),j=1,3) 53 format(' Xmn,Pmn: ',6e11.4) read(1,*,iostat=ioc)(x2bt(j,i),j=1,3), + (p2bt(j,i),j=1,3) if( ioc .ne. 0 ) then write(2,*) 'Error reading positions and momenta sigmas' go to 800 end if write(2,56)(x2bt(j,i),j=1,3),(p2bt(j,i),j=1,3) 56 format(' Xsg,Psg: ',6e11.4) ! else if( bdistyp(i) .eq. 2 ) then ! uniform circle read(1,*,iostat=ioc)x1bt(1,i),x2bt(1,i),x1bt(2,i), + x2bt(2,i),x1bt(3,i),x2bt(3,i) if( ioc .ne. 0 ) then write(2,*) 'Error reading position ranges' go to 800 end if write(2,63)x1bt(1,i),x2bt(1,i),x1bt(2,i), + x2bt(2,i),x1bt(3,i),x2bt(3,i) 63 format(' R,PHI,Z: ',6e10.3) read(1,*,iostat=ioc)p1bt(1,i),p2bt(1,i),p1bt(2,i), + p2bt(2,i),p1bt(3,i),p2bt(3,i) if( ioc .ne. 0 ) then write(2,*) 'Error reading momentum ranges' go to 800 end if write(2,66)p1bt(1,i),p2bt(1,i),p1bt(2,i), + p2bt(2,i),p1bt(3,i),p2bt(3,i) 66 format(' Pr,Pphi,Pz: ',6e10.3) x1bt(2,i) = x1bt(2,i) * pi /180. x2bt(2,i) = x2bt(2,i) * pi /180. ! else write(2,*) ' Must have 0 < BDISTYP < 3' go to 800 end if ! bdistyp ! read(1,*,iostat=ioc) nbcorr(i) if( ioc .ne. 0 ) then write(2,*) 'Error reading NBCORR' go to 800 end if write(2,'(a,i5)') ' NBCORR = ',nbcorr(i) if( nbcorr(i).lt.0 .or. nbcorr(i).gt.6 ) then write(2,*) 'Must have -1 < NBCORR < 7' go to 800 end if ! do j=1,nbcorr(i) read(1,*,iostat=ioc) corrtyp(j,i),corr1(j,i),corr2(j,i), & corr3(j,i) if( ioc .ne. 0 ) then write(2,*) 'Error reading CORRTYP,CORR1..3' go to 800 end if write(2,78) j,corrtyp(j,i),corr1(j,i),corr2(j,i),corr3(j,i) 78 format(' I,Ctyp,CORR1,CORR2,CORR3: ',2i5,3f12.4) end do ! end do ! end of loop on nbeamtyp ! c ENVELOPE EQUATIONS ADDITION if (run_env) then if ((p2bt(1,1).eq.0.).or.(x2bt(1,1).eq.0.)) then write(2,*)' << ENVELOPE EQNS REQUIRE SOME SIGMA X,PX >>' go to 800 else if (p1bt(3,1).le.0.) then write(2,*)' <<<< ENVELOPE EQNS REQUIRE PZ>0 >>>>' go to 800 end if end if ! end if ! bgen c c Read physics interactions control parameters read(1,ints,iostat=ioc) if( ioc .ne. 0 ) then write(2,*) 'Error reading interactions namelist' go to 800 end if ! ! Cant do spin & space charge simultaneously if( spin .and. lspace ) then write(2,*) 'Cant do spin & space charge simultaneously' write(*,*) 'Cant do spin & space charge simultaneously' stop end if ! c ENVELOPE EQUATIONS ADDITION if (run_env) then env_decay = ldecay ldecay = .false. lspace = .false. end if ! write(2,ints) ! c Set up quantities to histogram read(1,nhs,iostat=ioc) if( ioc .ne. 0 ) then write(2,*) 'Error reading histograms namelist' go to 800 end if write(2,nhs) if(nhist .gt. 20) then write(2,*) 'Must have NHIST < 21' go to 800 end if nhword = 1 c if(nhist.gt.0) then ! Don't read if nhist=0 do 260 i=1,nhist read(1,*,iostat=ioc)hxmin(i),hdx(i),nhbins(i),ihvar(i),ihdes(i) if( ioc .ne. 0 ) then write(2,*) 'Cant read HXMIN,HDX,NHBINS,IHVAR,IHDES' go to 800 end if if(nhbins(i) .gt. 50) nhbins(i) = 50 write(2,*)hxmin(i),hdx(i),nhbins(i),ihvar(i),ihdes(i) lhpnt(i) = nhword nhword = nhword + nhbins(i) 260 continue end if ! Don't read if nhist=0 c c Set up quantities to scatter plot read(1,nsc,iostat=ioc) if( ioc .ne. 0 ) then write(2,*) 'Error reading scatterplot namelist' go to 800 end if write(2,nsc) if(nscat .gt. 20) then write(2,*) 'Must have NSCAT < 21' go to 800 end if nsword = 1 c if(nscat.gt.0) then ! Don't read if nscat=0 do 280 i=1,nscat read(1,*,iostat=ioc)sxmin(i),sdx(i),nsxbin(i),isxvar(i), + isxdes(i),symin(i),sdy(i),nsybin(i),isyvar(i),isydes(i) if( ioc .ne. 0 ) then write(2,*) 'Cant read SXMIN,SDX,NSXBIN,ISXVAR,ISXDES', & ' SYMIN,SDY,NSYBIN,ISYVAR,ISYDES' go to 800 end if if(nsxbin(i) .gt. 50) nsxbin(i) = 50 if(nsybin(i) .gt. 23) nsybin(i) = 23 write(2,*)sxmin(i),sdx(i),nsxbin(i),isxvar(i),isxdes(i) write(2,*)symin(i),sdy(i),nsybin(i),isyvar(i),isydes(i) lspnt(i) = nsword nsword = nsword + nsxbin(i)*nsybin(i) 280 continue end if ! Don't read if nscat=0 c c Set up Z-HISTORY plots read(1,nzh,iostat=ioc) if( ioc .ne. 0 ) then write(2,*) 'Error reading Z-history namelist' go to 800 end if write(2,nzh) if(nzhist .gt. 20) then write(2,*) 'Must have NZHIST < 21' go to 800 end if c if(nzhist.gt.0) then ! Don't read if nzhist=0 do i=1,nzhist read(1,*,iostat=ioc)nzhpar(i),zhxmin(i),zhdx(i),nzhxbin(i), + zhymin(i),zhymax(i),izhyvar(i) if( ioc .ne. 0 ) then write(2,*) 'Cant read NZHPAR,ZHXMIN,ZHDX,NZHXBIN', & ' ZHYMIN,ZHYMAX,IZHYVAR' go to 800 end if if( ABS(nzhpar(i)) .gt. 10 ) nzhpar(i) = 10 if(nzhxbin(i) .gt. 70) nzhxbin(i) = 70 write(2,*)nzhpar(i),zhxmin(i),zhdx(i),nzhxbin(i) write(2,*)zhymin(i),zhymax(i),izhyvar(i) zhxmax(i) = zhxmin(i) + nzhxbin(i) * zhdx(i) end do ! do i=1,nzhist ! Initialize z-history location pointer do j=0,10 zhxcur(j,i) = zhxmin(i) end do end do end if ! Don't read if nzhist=0 ! c Set up R-HISTORY plots read(1,nrh,iostat=ioc) if( ioc .ne. 0 ) then write(2,*) 'Error reading R-history namelist' go to 800 end if write(2,nrh) if(nrhist .gt. 10) then write(2,*) 'Must have NRHIST < 11' go to 800 end if c if(nrhist.gt.0) then ! Don't read if nrhist=0 do i=1,nrhist read(1,*,iostat=ioc) irhzmin(i),irhdz(i), + rhymin(i),rhymax(i),irhyvar(i) if( ioc .ne. 0 ) then write(2,*) 'Cant read IRHZMIN,IRHDZ,RHYMIN,RHYMAX,IRHYVAR' go to 800 end if if( irhdz(i) .eq. 0 ) then write(2,*) ' IRHDZ = 0' go to 800 end if write(2,*) irhzmin(i),irhdz(i),rhymin(i),rhymax(i), & irhyvar(i) irhzmax(i) = irhzmin(i) + 69*irhdz(i) end do end if ! Don't read if nrhist=0 ! c Set up EMITTANCE calculations read(1,nem,iostat=ioc) if( ioc .ne. 0 ) then write(2,*) 'Error reading NEM namelist' go to 800 end if write(2,nem) if(nemit .gt. 100) then write(2,*) 'NEMIT > 100' go to 800 end if c if(nemit.gt.0) then ! Don't read if nemit=0 read(1,*,iostat=ioc) (iemdes(j),j=1,nemit) if( ioc .ne. 0 ) then write(2,*) 'Cant read IEMDES' go to 800 end if write(2,*) (iemdes(j),j=1,nemit) end if ! Don't read if nemit=0 ! c Set up COVARIANCE calculations read(1,ncv,iostat=ioc) if( ioc .ne. 0 ) then write(2,*) 'Error reading NCV namelist' go to 800 end if write(2,ncv) if(ncovar .gt. 100) then write(2,*) 'NCOVAR > 100' go to 800 end if c if(ncovar.gt.0) then ! Don't read if ncovar=0 read(1,*,iostat=ioc) (icvdes(j),j=1,ncovar) if( ioc .ne. 0 ) then write(2,*) 'Cant read ICVDES' go to 800 end if write(2,*) (icvdes(j),j=1,ncovar) end if ! Don't read if ncovar=0 ! -------------------------------------------------------- c c Check validity of the region data c nrlines = 0 ! count number of lines of region data nplines = 0 ! number of pre-looping lines of region data c 300 continue ! loop over region commands c read(1,'(a4)',iostat=ioc)cmd if( ioc .ne. 0 ) then write(2,*) 'Cant read CMD' go to 800 end if nrlines = nrlines + 1 write(2,*) cmd ! ! Don't allow unsupported commands in parsing loops if( cloop ) then do i=1,ncloop if( cmd .eq. cloopdat(i) )go to 303 end do write(2,*) 'Command not supported in CELL loop' go to 800 end if ! 303 continue if( rloop ) then do i=1,nrloop if( cmd .eq. rloopdat(i) )go to 306 end do write(2,*) 'Command not supported in REPEAT loop' go to 800 end if ! 306 continue if( bgloop ) then do i=1,nbgloop if( cmd .eq. bgloopdat(i) )go to 309 end do write(2,*) 'Command not supported in BACKGROUND loop' go to 800 end if ! 309 continue c ........................................................ if( cmd .eq. 'SECT' )then go to 300 ! ! ......................................................... else if( cmd .eq. 'BEGS' ) then ! begin repeating part of section nplines = nrlines go to 300 ! ! ........................................................ else if( cmd .eq. 'CELL' )then read(1,*,iostat=ioc) ncells if( ioc .ne. 0 ) then write(2,*) 'Cant read NCELLS' go to 800 end if write(2,*) ncells if( ncells .le. 0 ) then write(2,*) 'Must have NCELLS > 0' go to 800 end if read(1,*,iostat=ioc) cellflip if( ioc .ne. 0 ) then write(2,*) 'Cant read CELLFLIP' go to 800 end if write(2,*) cellflip read(1,'(a4)',iostat=ioc) cftag if( ioc .ne. 0 ) then write(2,*) 'Cant read CFTAG' go to 800 end if write(2,*) cftag do i=1,nftags if( cftag .eq. ftagdat(i) ) go to 315 end do write(2,*) 'Dont recognize cell field tag' go to 800 c 315 continue read(1,*,iostat=ioc)(cfparm(j),j=1,nfparm) if( ioc .ne. 0 ) then write(2,*) 'Cant read CFPARM' go to 800 end if write(2,*) (cfparm(j),j=1,nfparm) cloop = .true. nrlines = nrlines + 4 c ! ........................................................ else if( cmd .eq. 'ENDC' )then cloop = .false. go to 300 ! ! ........................................................ else if( cmd .eq. 'SREG' )then read(1,*,iostat=ioc) slen,nrreg,zstep if( ioc .ne. 0 ) then write(2,*) 'Error reading SLEN,NRREG,ZSTEP' go to 800 end if write(2,*) slen,nrreg,zstep if( slen .le. 0. ) then write(2,*) 'Must have SLEN > 0' go to 800 end if if( nrreg .le. 0 .or. nrreg .gt. 4 ) then write(2,*) 'Must have 0 < NRREG < 5' go to 800 end if if( scalestep*zstep .le. 0. ) then write(2,*) 'Must have SCALESTEP*ZSTEP > 0' go to 800 end if c do 330 ir=1,nrreg read(1,*,iostat=ioc) idm,rlow(ir),rhigh(ir) if( ioc .ne. 0 ) then write(2,*) 'Error reading IDM,RLOW,RHIGH' go to 800 end if write(2,*) idm,rlow(ir),rhigh(ir) if( rlow(ir) .lt. 0. ) then write(2,*) 'Must have RLOW .ge. 0' go to 800 end if if( rhigh(ir) .lt. rlow(ir) ) then write(2,*) 'Must have RHIGH > RLOW' go to 800 end if read(1,'(a4)',iostat=ioc) ftag(ir) if( ioc .ne. 0 ) then write(2,*) 'Error reading FTAG' go to 800 end if write(2,*) ftag(ir) do i=1,nftags if( ftag(ir) .eq. ftagdat(i) ) go to 320 end do write(2,*) 'Dont recognize region field tag' go to 800 c 320 continue ! good ftag read(1,*,iostat=ioc)(fparm(j,ir),j=1,nfparm) if( ioc .ne. 0 ) then write(2,*) 'Error reading FPARM' go to 800 end if write(2,*) (fparm(j,ir),j=1,nfparm) ! read(1,'(3a4)',iostat=ioc) mtag(ir),mtag2(ir),mtag3(ir) if( ioc .ne. 0 ) then write(2,*) 'Error reading MTAG' go to 800 end if write(2,*) mtag(ir),mtag2(ir),mtag3(ir) ! if( mtag(ir) .eq. 'VAC' ) go to 322 do i=1,nmat if( mtag(ir) .eq. matid(i) ) go to 322 end do write(2,*) 'Dont recognize region material tag' go to 800 ! 322 continue ! good mtag read(1,'(a6)',iostat=ioc) mgeom(ir) if( ioc .ne. 0 ) then write(2,*) 'Error reading MGEOM' go to 800 end if write(2,*) mgeom(ir) ! do i=1,nmgeoms if( mgeom(ir) .eq. mgeomdat(i) ) go to 324 end do write(2,*) 'Dont recognize region geometry tag' go to 800 ! 324 continue ! good mgeom read(1,*,iostat=ioc) (gparm(j,ir),j=1,ngparm) if( ioc .ne. 0 ) then write(2,*) 'Error reading GPARM' go to 800 end if write(2,*) (gparm(j,ir),j=1,ngparm) ! 330 continue ! nrlines = nrlines + 6 * nrreg + 1 c ! ........................................................ else if( cmd .eq. 'DENS' )then ! change material density read(1,*,iostat=ioc) mat,fact if( ioc .ne. 0 ) then write(2,*) 'Error reading DENSITY modification' write(2,*) mat,fact go to 800 end if do i=1,nmat if( mat .eq. matid(i) ) go to 335 end do write(2,*) 'Dont recognize material type' go to 800 335 continue write(2,*) mat,fact nrlines = nrlines + 1 ! ! ........................................................ else if( cmd .eq. 'DENP' )then ! change density profile read(1,*,iostat=ioc) mat,idir,a,b,c,d if( ioc .ne. 0 ) then write(2,*) 'Error reading DENSITY profile parameters' write(2,*) mat,idir,a,b,c,d go to 800 end if do i=1,nmat if( mat .eq. matid(i) ) go to 337 end do write(2,*) 'Dont recognize material type' go to 800 337 continue write(2,*) mat,idir,a,b,c,d nrlines = nrlines + 1 ! ! ........................................................ else if( cmd .eq. 'REFP' )then read(1,*,iostat=ioc) refpar,v1ref,v2ref,v3ref,phmodref if( ioc .ne. 0 ) then write(2,*) 'Error reading REFP data' go to 800 end if write(2,*) refpar,v1ref,v2ref,v3ref,phmodref nrlines = nrlines + 1 ! ! ........................................................ else if( cmd .eq. 'REF2' )then read(1,*,iostat=ioc) ref2par,v1ref2,v2ref2,v3ref2 if( ioc .ne. 0 ) then write(2,*) 'Error reading REF2 data' go to 800 end if write(2,*) ref2par,v1ref2,v2ref2,v3ref2 nrlines = nrlines + 1 ! ! ........................................................ else if( cmd .eq. 'REPE' )then read(1,*,iostat=ioc) nrep if( ioc .ne. 0 ) then write(2,*) 'Error reading NREP' go to 800 end if write(2,*) nrep if( nrep .le. 0 ) then write(2,*) 'Must have NREP > 0' go to 800 end if rloop = .true. nrlines = nrlines + 1 ! ! ........................................................ else if( cmd .eq. 'RKIC' )then ! random field kick read(1,*,iostat=ioc) mat,xm,sd,coup,rkaz if( ioc .ne. 0 ) then write(2,*) 'Error reading RKICK parameters' write(2,*) mat,xm,sd,coup,rkaz go to 800 end if write(2,*) mat,xm,sd,coup,rkaz nrlines = nrlines + 1 ! ! ........................................................ else if( cmd .eq. 'ROTA' )then read(1,*,iostat=ioc) rotang,iaxis,iclass if( ioc .ne. 0 ) then write(2,*) 'Error reading ROTATION parameters' go to 800 end if write(2,*) rotang,iaxis,iclass nrlines = nrlines + 1 ! ! ........................................................ else if( cmd .eq. 'TILT' )then read(1,*,iostat=ioc) psi,phi if( ioc .ne. 0 ) then write(2,*) 'Error reading TILT parameters' go to 800 end if write(2,*) psi,phi nrlines = nrlines + 1 ! ! ........................................................ else if( cmd .eq. 'DISP' )then read(1,*,iostat=ioc) d,phi if( ioc .ne. 0 ) then write(2,*) 'Error reading DISP parameters' go to 800 end if write(2,*) d,phi nrlines = nrlines + 1 ! ! ........................................................ else if( cmd .eq. 'EDGE' )then read(1,'(a4)',iostat=ioc) grtype if( ioc .ne. 0 ) then write(2,*) 'Error reading EDGE type' go to 800 end if write(2,'(a4)') grtype read(1,*,iostat=ioc) model,p1,p2,p3,p4 if( ioc .ne. 0 ) then write(2,*) 'Error reading EDGE parameters' go to 800 end if write(2,*) model,p1,p2,p3,p4 nrlines = nrlines + 2 ! ! ........................................................ else if( cmd .eq. 'GRID' )then read(1,*,iostat=ioc) jgrid if( ioc .ne. 0 ) then write(2,*) 'Error reading grid number' go to 800 end if read(1,'(a4)',iostat=ioc) grtype if( ioc .ne. 0 ) then write(2,*) 'Error reading grid type' go to 800 end if if( grtype.ne.'SHEE' .and. grtype.ne.'BLOC' .and. & grtype.ne.'COIL' .and. grtype.ne.'SOL' .and. & grtype.ne.'BROD' .and. grtype.ne.'STUS' .and. & grtype.ne.'G4RZ' .and. grtype.ne.'G43D') then write(2,*)' Illegal grid type' go to 800 end if read(1,*,iostat=ioc) fpar if( ioc .ne. 0 ) then write(2,*) 'Error reading field parameters' go to 800 end if write(2,*) jgrid write(2,'(1x,a4)') grtype write(2,*) fpar nrlines = nrlines + 3 ! ! ........................................................ else if( cmd .eq. 'APER' )then read(1,*,iostat=ioc) iapertype,aperlim1,aperlim2,aperlim3 if( ioc .ne. 0 ) then write(2,*) 'Error reading IAPERTYPE,APERLIM1..3' go to 800 end if write(2,*) iapertype,aperlim1,aperlim2,aperlim3 nrlines = nrlines + 1 ! ! ........................................................ else if( cmd .eq. 'DVAR' )then ! change variable read(1,*,iostat=ioc) ivar,val,iclass if( ioc .ne. 0 ) then write(2,*) 'Error reading DVAR' go to 800 end if write(2,*) ivar,val,iclass nrlines = nrlines + 1 ! ! ........................................................ else if( cmd .eq. 'CUTV' )then ! cut on variable read(1,*,iostat=ioc) ivar,ivtyp,val if( ioc .ne. 0 ) then write(2,*) 'Error reading CUTV' go to 800 end if write(2,*) ivar,ivtyp,val nrlines = nrlines + 1 ! ! ........................................................ else if( cmd .eq. 'TAPE' )then ! taper read(1,*,iostat=ioc) ivar if( ioc .ne. 0 ) then write(2,*) 'Error reading NCELLTAP' go to 800 end if write(2,*) ivar nrlines = nrlines + 1 ! ! ........................................................ else if( cmd .eq. 'TRAN' )then if( refpar .eq. 0 ) then write(2,*)' Transport requires a reference particle' go to 800 end if do i=1,6 read(1,*,iostat=ioc) (trprt(i,j),j=1,6) if( ioc .ne. 0 ) then write(2,*) 'Error reading TRPRT' go to 800 end if end do write(2,'(6f12.4)') ((trprt(i,j),j=1,6),i=1,6) nrlines = nrlines + 6 ! ! ........................................................ else if( cmd .eq. 'OUTP' )then go to 300 ! ! ......................................................... else if( cmd .eq. 'DUMM' )then go to 300 ! ! ......................................................... else if( cmd .eq. 'RESE' )then go to 300 ! ! ........................................................ else if( cmd .eq. 'ENDR' )then rloop = .false. go to 300 ! ! ........................................................ else if( cmd .eq. 'BACK' )then read(1,*,iostat=ioc) bentbkg,prefbkg,ztotalbkg if( ioc .ne. 0 ) then write(2,*) 'Error reading BENTBKG,PREFBKG,ZTOTALBKG' go to 800 end if write(2,*) bentbkg,prefbkg,ztotalbkg if( prefbkg .le. 0. ) then write(2,*)' Must have PREFBKG > 0' go to 800 end if if( ztotalbkg .le. 0. ) then write(2,*)' Must have ZTOTALBKG > 0' go to 800 end if read(1,*,iostat=ioc) xlogr,dxgr,nxgr, & ylogr,dygr,nygr,slogr,dsgr,nsgr,interbkg if( ioc .ne. 0 ) then write(2,*) 'Error reading XLOBKG,DXBKG,NXBKG', & 'YLOBKG,DYBKG,NYBKG,ZLOBKG,DZBKG,NZBKG,INTERBKG' go to 800 end if write(2,*) xlogr,dxgr,nxgr, & ylogr,dygr,nygr,slogr,dsgr,nsgr,interbkg if( nxgr.le.0 .or. nxgr.gt.mxxg ) then write(2,*)' Must have 0 < NXBKG < 31' go to 800 end if if( nygr.le.0 .or. nygr.gt.mxyg ) then write(2,*)' Must have 0 < NYBKG < 31' go to 800 end if if( nsgr.le.0 .or. nsgr.gt.mxsg ) then write(2,*)' Must have 0 < NYBKG < 201' go to 800 end if if( dxgr .le. 0. ) then write(2,*)' Must have DXBKG > 0' go to 800 end if if( dygr .le. 0. ) then write(2,*)' Must have DYBKG > 0' go to 800 end if if( dsgr .le. 0. ) then write(2,*)' Must have DZBKG > 0' go to 800 end if if( interbkg.lt.0 .or. interbkg.gt.3 ) then write(2,*)' Must have -1 < INTERBKG < 4' go to 800 end if bgloop = .true. ! ncrvbkgfld = 0 nrlines = nrlines + 2 ! ! ........................................................ else if( cmd .eq. 'BFIE' )then ! background field read(1,*,iostat=ioc) zoffbkg,rmaxbkg,zminbkg,zmaxbkg if( ioc .ne. 0 ) then write(2,*) 'Error reading ZOFFBKG,RMAXBKG,ZMINBKG,ZMAXBKG' go to 800 end if write(2,*) zoffbkg,rmaxbkg,zminbkg,zmaxbkg read(1,'(a4)',iostat=ioc) bftag if( ioc .ne. 0 ) then write(2,*) 'Error reading BFTAG' go to 800 end if write(2,*) bftag do i=1,nftags if( bftag .eq. ftagdat(i) ) go to 340 end do write(2,*) 'Dont recognize background field tag' go to 800 c 340 continue ! read(1,*,iostat=ioc)(bfparm(j),j=1,nfparm) if( ioc .ne. 0 ) then write(2,*) 'Error reading BFPARM' go to 800 end if write(2,*) (bfparm(j),j=1,nfparm) nrlines = nrlines + 3 ! ! ........................................................ else if( cmd .eq. 'ENDB' )then read(1,*,iostat=ioc) idm if( ioc .ne. 0 ) then write(2,*) 'Error reading endb idm' go to 800 end if write(2,*) idm ! if( ncrvbkgfld .gt. 1 ) go to 820 bgloop = .false. nrlines = nrlines + 1 ! ! ........................................................ else if( cmd .eq. 'ENDS' )then go to 350 c ! ........................................................ else if(termout) write(*,*) & '<<< INPUT ERROR: see log file FOR002.DAT >>>' write(2,*) 'Bad region command in SETUP' go to 800 end if c go to 300 ! get next region command c 350 continue ! ! -------------------------------------------------------- ! c Print summary of region information ! Set up reference particle data ! --------------------------------------------------------- call REG_SUMMARY ! --------------------------------------------------------- if( iflag .ne. 0 ) go to 800 write(2,*) ! ! Define magnet configuration for the problem if( magconf.gt.19 .and. magconf.lt.100 ) then m = magconf write(fname,40) m write(2,'(a,a10)') 'Opening file ',fname open(unit=m, & file=pathname(1:LASTNB(pathname)) // fname, & status='unknown',iostat=ioc) if( ioc .ne. 0 ) then write(2,*) 'Error opening magnet configuration file' go to 800 end if read(m,'(a80)',iostat=ioc) title if( ioc .ne. 0 ) then write(2,*) 'Cant read magnet configuration title' go to 800 end if read(m,*,iostat=ioc) nmagmc if( ioc .ne. 0 ) then write(2,*) 'Error reading NMAGMC' go to 800 end if ! do i=1,nmagmc read(m,*,iostat=ioc) idum,z1mc(i),dzmc(i),r1mc(i),drmc(i), & curdmc(i),nshmc(i) if( ioc .ne. 0 ) then write(2,*) 'Error reading IDUM,Z1MC,DZMC,R1MC,DRMC', & 'CURDMC,NSHMC' go to 800 end if end do ! close(m) end if ! ! Define magnet map regions for the problem if( mapdef.gt.19 .and. mapdef.lt.100 ) then m = mapdef write(fname,40) m write(2,'(a,a10)') 'Opening file ',fname open(unit=m, & file=pathname(1:LASTNB(pathname)) // fname, & status='unknown',iostat=ioc) if( ioc .ne. 0 ) then write(2,*) 'Error opening magnet map definition file' go to 800 end if read(m,'(a80)',iostat=ioc) title if( ioc .ne. 0 ) then write(2,*) 'Error opening magnet map definition title' go to 800 end if read(m,*,iostat=ioc) nmapmd if( ioc .ne. 0 ) then write(2,*) 'Error reading NMAPMD' go to 800 end if ! do i=1,nmapmd read(m,*,iostat=ioc) idum,z1md(i),dzmd(i),nzmd(i),drmd(i), & nrmd(i),dzcutmd(i) if( ioc .ne. 0 ) then write(2,*) 'Error reading NMAPMD' go to 800 end if end do ! close(m) end if ! c ........................................................ c Initialize call CLEAR ! clear histograms, scatterplots, statistics call ICOOLTIME(ihr,imin,isec) write(2,'(a,4i6)')' Completed SETUP: ',ihr,imin,isec go to 900 c 800 continue c Error in input if(termout) write(*,*)' <<<<< INPUT ERROR >>>>>' write(2,*)' <<<<< INPUT ERROR >>>>>' iflag = -1 c 900 continue return end ! ********************************************************* subroutine SIMULATE c c Execute simulation of the cooling sections c implicit double precision (a-h,o-z) include 'icommon.inc' ! logical kstat,spcharge,curved logical lscatsav,lstragsav,ldraysav,ldecaysav,lspacesav logical vstepsav,steprksav,fsavfirst logical*1 cell1,rep1,repincell,bkgon character*80 probtitle character*1 key character*4 cmd,mat ! data kstat/.false./ ,fsavfirst/.true./,jsold/-1/ ! if(.not.bgen ) then read(3,'(a80)') probtitle ! skip title card on beam input file ! Read reference particle data read(3,*) zbref,pbref,tbref,zb2ref,pb2ref,tb2ref,typref end if ! lscatsav = lscatter ! save user control parameters lstragsav = lstrag ldraysav = ldray ldecaysav = ldecay lspacesav = lspace vstepsav = varstep steprksav = steprk ! ngood = ngoodp ! i7 = 0 cell1 = .false. ! set true first time thru cell rep1 = .false. ! set true first time thru repetition repincell = .false. ! set true if there is a repeat in a cell nclines = 0 nreplines = 0 jcel = 0 kcel = 0 cftag = 'NONE' bftag = 'NONE' bkgon = .false. jrkval = 0 zlo = zlostart ncrvcell = 0 ncrvback = 0 loutput = .false. ! -------------------------------------------------------- ! do 600 is=1,nsections jsec = is ss0 = zlo ! ! Rewind input to top of region data in the input file if( is .eq. 1 ) then nlines = nrlines else ! we have already read the list before if( is .eq. 2 ) nreg = i7 - i7b ! save # of regions in a section if( nplines .gt. 0 ) then nlines = nrlines - nplines else write(2,*) 'Error in REG_SUMMARY (missing BEGS command?)' if(termout) & write(*,*) 'Error in REG_SUMMARY (missing BEGS command?)' iflag = -1 go to 900 end if end if ! do i=1,nlines backspace (unit=1) end do c ......................................................... ! 100 continue ! loop over "region" commands ! read(1,'(a4)',iostat=ioc)cmd ! ........................................................ ft = 'NONE' jpsr = 0 ncrvreg = 0 ! if( cmd .eq. 'CELL' )then read(1,*,iostat=ioc) ncells read(1,*) cellflip read(1,'(a4)',iostat=ioc) cftag read(1,*,iostat=ioc)(cfparm(j),j=1,nfparm) jcel = jcel + 1 ! jcel counts new CELL structures kcel = 1 ! kcel counts individual cells inside CELL structure ! Thus if kcel .ne. 0, then we are currently inside a CELL structure sc0 = zlo ! initial position of this CELL cellslp = 1. ! cell field symmetry cell1 = .true. if( CURVED(cftag,2) ) ncrvcell = 1 go to 100 ! ! ........................................................ else if( cmd .eq. 'REPE' )then read(1,*,iostat=ioc) nrep rep1 = .true. if( cell1 ) then repincell = .true. nclines = nclines + 2 end if go to 100 ! ! ........................................................ else if( cmd .eq. 'SECT' )then i7 = i7 + 1 jsrg = i7 ! jsec = 0 ! ! ........................................................ else if( cmd .eq. 'BEGS' )then i7b = i7 go to 100 ! ! ........................................................ else if( cmd .eq. 'ROTA' )then ! read rotang,iaxis,iclass read(1,*,iostat=ioc) val1ps,ival1ps,ival2ps i7 = i7 + 1 jsrg = i7 val1ps = val1ps*pi/180.d0 if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 jpsr = -1 ! ! ........................................................ else if( cmd .eq. 'OUTP' )then i7 = i7 + 1 if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 1 end if if( rep1 ) nreplines = nreplines + 1 jpsr = -2 ! call PSEUDOREGION(i7,1) loutput = .true. go to 100 ! ! ........................................................ else if( cmd .eq. 'APER' )then ! read iapertype,aperlim1,aperlim2 read(1,*,iostat=ioc) ival1ps,val1ps,val2ps,val3ps i7 = i7 + 1 jsrg = i7 if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 jpsr = -3 ! ! ........................................................ else if( cmd .eq. 'BACK' )then read(1,*,iostat=ioc) bentbkg,prefbkg,ztotalbkg read(1,*,iostat=ioc) xlogr,dxgr,nxgr, & ylogr,dygr,nygr,slogr,dsgr,nsgr,interbkg i7 = i7 + 1 bkgon = .true. saccbkg = 0. jpsr = -4 call PSEUDOREGION(i7,1) go to 100 ! ! ........................................................ else if( cmd .eq. 'BFIE' )then ! background field read(1,*,iostat=ioc) zoffbkg,rmaxbkg,zminbkg,zmaxbkg read(1,'(a4)',iostat=ioc) bftag read(1,*,iostat=ioc)(bfparm(j),j=1,nfparm) i7 = i7 + 1 jpsr = -5 call PSEUDOREGION(i7,1) if( CURVED(bftag,3) ) ncrvback = ncrvback + 1 go to 100 ! ! ........................................................ else if( cmd .eq. 'ENDB' )then ! jpsr = -6 ! read mfile read(1,*,iostat=ioc) ival1ps i7 = i7 + 1 jpsr = -6 call PSEUDOREGION(i7,1) ncrvback = 0 go to 100 ! ! ........................................................ else if( cmd .eq. 'TRAN' )then ! jpsr = -7 read(1,*,iostat=ioc) ((trprt(i,j),j=1,6),i=1,6) i7 = i7 + 1 jsrg = i7 zhi = zlo ! zlo = zhi if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 7 end if if( rep1 ) nreplines = nreplines + 7 jpsr = -7 ! ! ........................................................ else if( cmd .eq. 'REFP' )then ! jpsr = -8 read(1,*) refpar,v1ref,v2ref,v3ref,phmodref i7 = i7 + 1 if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 jpsr = -8 call PSEUDOREGION(i7,0) if( phmodref .eq. 2 ) then varstep = .false. else varstep = vstepsav end if go to 100 ! ! ........................................................ else if( cmd .eq. 'DUMM' )then ! jpsr = -9 i7 = i7 + 1 if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 1 end if if( rep1 ) nreplines = nreplines + 1 go to 100 ! ! ........................................................ else if( cmd .eq. 'TILT' )then ! jpsr = -10 ! read psi,phi read(1,*,iostat=ioc) val1ps,val2ps i7 = i7 + 1 jsrg = i7 if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 jpsr = -10 ! ! ........................................................ else if( cmd .eq. 'DISP' )then ! jpsr = -11 ! read d,phi read(1,*,iostat=ioc) val1ps,val2ps i7 = i7 + 1 jsrg = i7 if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 jpsr = -11 ! ! ........................................................ else if( cmd .eq. 'REF2' )then ! jpsr = -12 read(1,*) ref2par,v1ref2,v2ref2,v3ref2 i7 = i7 + 1 if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 jpsr = -12 call PSEUDOREGION(i7,-1) go to 100 ! ! ........................................................ else if( cmd .eq. 'CUTV' )then ! jpsr = -13 ! read ivar,ivtyp,val read(1,*,iostat=ioc) ival1ps,ival2ps,val1ps i7 = i7 + 1 jsrg = i7 slen = 0. ! sr0 = zlo if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 jpsr = -13 ! ........................................................ else if( cmd .eq. 'EDGE' )then ! jpsr = -14 ! read grtype ! read model,p1,p2,p3,p4 read(1,'(a4)',iostat=ioc) matps read(1,*,iostat=ioc) ival1ps,val1ps,val2ps,val3ps,val4ps i7 = i7 + 1 jsrg = i7 if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 3 end if if( rep1 ) nreplines = nreplines + 3 jpsr = -14 ! ! ........................................................ else if( cmd .eq. 'GRID' )then ! jpsr = -15 ! read jgrid ! read grtype ! read fpar read(1,*,iostat=ioc) ival1ps read(1,'(a4)',iostat=ioc) matps read(1,*,iostat=ioc) fparps i7 = i7 + 1 jpsr = -15 call PSEUDOREGION(i7,1) go to 100 ! ! ........................................................ else if( cmd .eq. 'RESE' )then ! jpsr = -16 i7 = i7 + 1 jsrg = i7 if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 1 end if if( rep1 ) nreplines = nreplines + 1 jpsr = -16 ! ! ........................................................ else if( cmd .eq. 'DVAR' )then ! jpsr = -17 ! read ivar,val,iclass read(1,*,iostat=ioc) ival1ps,val1ps,ival2ps i7 = i7 + 1 jsrg = i7 if( ival1ps .eq. 3 ) then slen = val1ps else slen = 0. end if zhi = zlo + slen zlo = zhi ! sr0 = zlo if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 jpsr = -17 ! ! ........................................................ else if( cmd .eq. 'DENS' )then ! jpsr = -18 ! read mat,fact read(1,*,iostat=ioc) matdens,densfac i7 = i7 + 1 if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 jpsr = -18 ! call PSEUDOREGION(i7,1) go to 100 ! ! ........................................................ else if( cmd .eq. 'RKIC' )then ! jpsr = -19 read(1,*,iostat=ioc) mat,xm,sd,coup,azflg i7 = i7 + 1 jsrg = i7 jrkval = jrkval + 1 if( jrkval .gt. nrkval ) jrkval = 1 if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 jpsr = -19 ! ! ........................................................ else if( cmd .eq. 'DENP' )then ! jpsr = -20 read(1,*,iostat=ioc) matdenp,dirdenp,adenp,bdenp,cdenp,ddenp i7 = i7 + 1 jsrg = i7 if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 jpsr = -20 go to 100 ! ! ........................................................ else if( cmd .eq. 'TAPE' )then ! jpsr = -21 read(1,*,iostat=ioc) ncelltap i7 = i7 + 1 jpsr = -21 go to 100 ! ........................................................ else if( cmd .eq. 'SREG' )then read(1,*,iostat=ioc) slen,nrreg,zstep i7 = i7 + 1 jsrg = i7 sr0 = zlo ! if( bkgon ) then if( saccbkg .le. slogr+ztotalbkg ) then bftag = 'BACK' saccbkg = saccbkg + slen else bftag = 'NONE' bkgon = .false. end if end if ! do ir=1,nrreg read(1,*,iostat=ioc) idm,rlow(ir),rhigh(ir) read(1,'(a4)',iostat=ioc) ftag(ir) read(1,*,iostat=ioc)(fparm(j,ir),j=1,nfparm) read(1,'(3a4)',iostat=ioc) mtag(ir),mtag2(ir),mtag3(ir) read(1,'(a6)',iostat=ioc) mgeom(ir) read(1,*,iostat=ioc) (gparm(j,ir),j=1,ngparm) end do ! do j=1,nrreg if( CURVED(ftag(j),1) ) ncrvreg = ncrvreg + 1 end do ! ! Check number of curved fields present ncurved = ncrvreg + ncrvcell + ncrvback if( ncurved .gt. 1 ) then write(2,*) 'Cant use more than 1 curved field type at a time' iflag = -11 go to 900 end if ! ......................................................... ! Patch for scaling tapered channels with BSOL(5) cell field if(cftag.eq.'BSOL' .and. cfparm(1).eq.5 .and. kbcell.eq.1 & .and. jctap.le.ncelltap)then ! if( jsrg .ne. jsold ) then ! first time in this region slentap = slen apertap = rhigh(1) if( ftag(1) .eq. 'ACCE' ) then ift = fparm(1,1) if( ift.eq.2 .or. ift.eq.11 ) then freqtap = fparm(2,1) gradtap = fparm(3,1) end if end if if( mgeom(1) .eq. 'WEDGE' ) then uwtap = gparm(2,1) zwtap = gparm(3,1) widthtap = gparm(5,1) heighttap = gparm(6,1) end if jsold = jsrg end if ! slen = slentap / ftap(jctap) rhigh(1) = apertap / ftap(jctap) if( ftag(1) .eq. 'ACCE' ) then fparm(2,1) = freqtap * ftap(jctap) ! Leave RF gradient constant to match scaling of absorber thickness end if if( mgeom(1) .eq. 'WEDGE' ) then gparm(2,1) = uwtap / ftap(jctap) gparm(3,1) = zwtap / ftap(jctap) gparm(5,1) = widthtap / ftap(jctap) gparm(6,1) = heighttap / ftap(jctap) end if end if ! acctap = acctap + slen ! if( kdiagtap.ge.20 .and. kdiagtap.le.99 ) then if( ftag(1) .eq. 'ACCE' ) then freq = fparm(2,1) else freq = 0. end if if( mgeom(1) .eq. 'WEDGE' ) then uw = gparm(2,1) zw = gparm(3,1) w = gparm(5,1) h = gparm(6,1) else uw = 0. zw = 0. w = 0. h = 0. end if ! write(kdiagtap,140) jsrg,jsec,kcel,cftag,ftag(1),ftap(jctap), & slen,acctap,rhigh(1),bfld(3),pertap,freq,uw,zw,w,h 140 format(t2,i4,t7,i3,t11,i3,t17,a4,t24,a4,t28,f8.3,t37,f9.4, & t49,f9.4,t59,f9.3,t70,f9.3,t80,f9.3,t90,f9.3,t100,f9.3, & t110,f9.3,t120,f9.3,t130,f9.3) end if ! ......................................................... ! zhi = zlo + slen zlo = zhi ! if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 6 * nrreg + 2 end if ! if( rep1 ) nreplines = nreplines + 6 * nrreg + 2 ! jsec = 1 ! I do not use a goto here, so the program goes to the next statement ! after the test on cmd finishes ! ........................................................ else if( cmd .eq. 'ENDR' )then if( rep1 ) nreplines = nreplines + 1 if( cell1 .and. rep1 ) nclines = nclines + 1 rep1 = .false. ! stop incrementing nreplines nrep = nrep - 1 if( nrep .gt. 0 ) then ! more repetitions left c do i=1,nreplines backspace (unit = 1) end do c else ! done until start of any new repetitions nreplines = 0 repincell = .false. end if go to 100 ! ! ........................................................ else if( cmd .eq. 'ENDC' )then if( cell1 ) nclines = nclines + 1 cell1 = .false. ! stop incrementing nclines ncells = ncells - 1 ! if( ncells .gt. 0 ) then ! more cells left kcel = kcel + 1 sc0 = zlo ! reset SC0 for next cell if( cellflip ) cellslp = -cellslp c do i=1,nclines backspace (unit = 1) end do c else ! done until start of any new cells nclines = 0 kcel = 0 cftag = 'NONE' ncrvcell = 0 end if ! end of test on ncells ! go to 100 ! ! ........................................................ else if( cmd .eq. 'ENDS' )then go to 600 ! end if ! end of test on cmd ! ! -------------------------------------------------------- ! j7rec = i7 if( termout .and. jsec.gt.0 ) then write(*,*) ' Processing region: ',jsrg,' of ', n7records write(*,*) ' Good tracks: ',npart-nfail, ' of ',npart end if ! if( i7 .eq. 1 ) then ! first region => set up for new particle ! The first "region" is defined to be the act of particle production ! do i=1,3 ! Initialize some variables for OUTPUT bfld(i) = 0.d0 efld(i) = 0.d0 end do phaserf = 0.d0 ! if( fsav .and. ( (izfile.eq.1).or.(izfile.lt.0) ) ) then ! Write out reference particle data if( bgen ) then ! new data write(4,'(7e14.6)') 0.,0.,0., 0.,0.,0., 0. else ! pass along info from input file write(4,'(7e14.6)') zbref,pbref,tbref, & zb2ref,pb2ref,tb2ref,typref end if end if ! ! Load particle data into arrays / file ipstop = npart knt = ipstart - 1 knp = 0 ! ......................................................... do 200 ip=ipstart,ipstop ! loop over particles ! jpar = ip iflag = 0 ! call NEW_PARTICLE(ip,ifl) ! if( ip .gt. 0 ) then knp = knp + 1 ! count true particles if( knp .le. npskip ) ifl = -1 end if ! if( ifl .eq. 0 ) then ngood = ngood + 1 knt = knt + 1 call WRITE_EVT(knt) ! save production information else ! bad or missing or skipped particle npart = npart -1 go to 200 end if ! if( ip .gt. 0 ) then call FILL_HIST(1) call FILL_SCAT(1) if( spin ) call FILL_POL(1) call SAVE_FILE(1) end if ! 200 continue ! ......................................................... ! Initialize accumulated distance if( .not. bgen ) then ! beam restart if( zbref .ne. 0. ) then ! ref part defined saccum = zbref else ! take last real particle saccum = xp(3) end if else ! new beam saccum = 0.d0 end if ! go to 580 ! end if ! end of test on i7 = 1 ! --------------------------------------------------------- if( cmd .eq. 'SREG' ) then ! initialize SREGIONs ft = ftag(1) mgt = mgeom(1) ! ! Don't allow variable step sizes in wedge regions if( mgt.eq.'WEDGE' .or. mgt.eq.'PWEDGE' .or. & mgt.eq.'ASPW' .or. mgt.eq.'ASRW' .or. mgt.eq.'RING' ) then varstep = .false. else varstep = vstepsav end if ! ! Force variable stepping in hemispherical window regions if( mgt .eq. 'HWIN' ) then varstep = .true. else varstep = vstepsav end if ! ! Don't allow variable step sizes in multiple radial regions ! if( nrreg .gt. 1 ) then ! varstep = .false. ! else ! varstep = vstepsav ! end if ! ! Check if we need space-charge correction for this region ! n.b. the variable LSPACE allows space charge calculations to ! take place. The variable SPCHARGE determines if corrections are ! actually made in a given region. if( lspace ) then spcharge = .true. else spcharge = .false. end if ! ! Space charge kick for start of region if( spcharge .and. jsec.gt.0 ) then call SC_KICK(1) end if ! ! Initialize dipole logic edge1 = .false. edge2 = .false. edgekick = .false. nbffag = 0 ! DIP model 3 ! ! Set field-on flag jfieldon = 0 if( cftag .ne. 'NONE' ) jfieldon = 1 if( bftag .ne. 'NONE' ) jfieldon = 1 do i=1,nrreg if( ftag(i) .ne. 'NONE' ) jfieldon = 1 end do ! ! Check if this is a curved region if( ncurved .eq. 0 ) then ! reset to straight htrk = 0.d0 hptrk = 0.d0 gtrk = 0.d0 gptrk = 0.d0 tortion = 0.d0 end if ! ! Force RK stepping in curved regions if( ncurved .ne. 0 ) steprk = .true. ! ! Logic for OUTPUT of initial track values if( output1 .and. firstsreg ) then jsrgsav = jsrg jparsav = jpar jsrg = 1 ! this tells OUTPUT to use initial values do i=ipstart,npart jpar = i call READ_EVT(i) call FIELD(xp,tp) call OUTPUT end do jsrg = jsrgsav jpar = jparsav firstsreg = .false. end if ! end if ! end of test on cmd = sreg ! -------------------------------------------------------- ! ! Propagate each particle thru the cooling region ! do 500 ip=ipstart,npart ! loop over particles jpar = ip iflag = 0 firstrf = .true. ! rf diagnostics flag c c Check elapsed time call CHKTIME(timelim) if( iflag .ne. 0 ) go to 900 ! c Check if user pressed key to interrupt call CHECK_KEY(kstat) if( kstat ) then call IN_KEY(key) if( key .eq. 'x' ) go to 620 if( key .eq. 'p' ) then write(*,*) 'Paused, hit Enter' read(*,*) end if end if ! jp = MOD( ip,250 ) if( jsec.ge.0 .and. jp.eq.0 ) then if(termout) write(*,*) ' Current particle: ',ip end if ! call READ_EVT(ip) if( ipflg .ne. 0 ) go to 500 ! ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! ! Process "pseudo-regions" if( jpsr .ne. 0 )then ! pseudo-region ! call PSEUDOREGION(i7,ip) ! if( iflag .eq. -1 ) go to 900 ! fatal error ! if( jpsr .eq. -1 ) then ! rotate go to 380 else if( jpsr.eq.-3 ) then ! aperture go to 390 else if( jpsr.eq.-7 ) then ! transport go to 380 else if( jpsr .eq. -10 ) then ! tilt go to 380 else if( jpsr .eq. -11 ) then ! disp go to 380 else if( jpsr .eq. -13 ) then ! cut on variable go to 390 else if( jpsr .eq. -14 ) then ! fringe field edge kick go to 380 else if( jpsr .eq. -16 ) then ! reset particle times go to 380 else if( jpsr .eq. -17 ) then ! change variable go to 380 else if( jpsr .eq. -19 ) then ! random field kick go to 380 ! else ! logic error write(*,*) 'Fatal error in pseudoregion logic' stop end if ! end if ! end of test on jsec.lt.0 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! ! Physical region ! ! Get reference particle data for this region ! if( ip .lt. 1 ) then ! mean dE/dx only for ref particle lscatter = .false. lstrag = .false. ldray = .false. ldecay = .false. lspace = .false. ! call PHMODEL(iskip) if( iskip .ne. 0 ) go to 380 end if ! end of test on ip < 1 ! ! ENVELOPE EQUATIONS ADDITION if (run_env .and. ((ievt.eq.1) .and. (i7.eq.3))) then ! first "real" region in moments equations run call FIELD(xp,tp) B0 = bfld(3) env_angm = env_angm + & 0.5d0 * pcharge * B0 * x2bt(1,1) / p2bt(1,1) end if ! ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c call GO_REGION ! take particle to end of this s region c ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c 380 continue ! ! Print out current data if( ip .le. nprnt ) then write(2,382)jsrg,ip,iptyp,xp,pp,tp 382 format(' End of s-REGION:',i4,' , particle:',i4, & ' , type:',i4,/,7e11.4) write(2,*) ! if( spintrk.gt.0 .and. ip.gt.0 ) then if( ABS(iptyp) .eq. 2 ) then hlab = HELICITY( pol ) write(2,384)pol,hlab 384 format(' Spin,Hlab:',3f7.3,5x,f8.4) else write(2,386)pol 386 format(' Spin:',3f7.3) end if ! else ! write(2,'(/)') end if ! end of test on spin end if ! end of test on ip.le.nprnt ! ! Check position of particle in rf bucket if( ip .gt. 0 ) then dt = ABS( tp - tpref ) if( dt .ge. bunchcut ) iflag = -36 end if ! ! RF debugging information if( ip.gt.0 .and. ft.eq.'ACCE' .and. rfdiag.gt.19 ) then dpzmn = pzparmean - pzrefmean dtmn = tparmean - trefmean dphasemn = freqrf * dtmn * 360. dpz = pp(3) - ppref(3) dt = tp - tpref write(rfdiag,387) jsrg,ievt,dphasemn,dtmn,dpzmn, & freqrf,emagrf,phaserf,ezrf,bphirf,trefmean,t2refmean, & zparmean 387 format(t1,i5,t8,i5,t16,f10.3,t26,10e13.5) end if ! ........................................................ 390 continue ! Come in here for some pseudo-regions ! if( ip .eq. 0 ) then ! restore interaction flags lscatter = lscatsav lstrag = lstragsav ldray = ldraysav ldecay = ldecaysav lspace = lspacesav end if ! p-regions can also generate an iflag .ne. 0 condition) ipflg = iflag ! if( iflag .eq. -1 ) go to 900 ! fatal error ! ! Update particle arrays ! Pseudoregions need the particle arrays to be updated, e.g. ROTATE ! n.b. I have to do this, even if the particle fails, so the following ! regions know that this is a bad particle ! call WRITE_EVT(ip) ! ........................................................ if( iflag .ne. 0 ) then if(termout) write(*,395)iflag,ievt,ipnum,ipflg,jsrg write(2,395)iflag,ievt,ipnum,ipflg,jsrg 395 format(' FAILURE; IFLAG,IEVT,IPNUM,IPFLG,JSRG: ',5i8) end if ! ! Generate OUTPUT file information for NTUPLE control variable ! and OUTPUT pseudocommand. (RTUPLE output originates in subr ! GO_REGION for normal regions) if( (loutput.and..not.ntuple.and..not.rtuple) .or. & (ntuple.and..not.rtuple) .or. & (rtuple.and.jpsr.ne.0) ) then if( xp(3) .gt. xpint(3) ) then call OUTPUT end if end if ! ! Determine header information for beam restart case if( fsav .and. i7.eq.izfile .and. fsavfirst) then if( ipstart.eq.-1 .and. ip.eq.0 ) then if( fsavset ) then write(4,'(6e14.6,i5)') 0.,ppref(3),0., & 0.,pp2ref(3),tp2ref-tpref,ref2par else write(4,'(6e14.6,i5)') xpref(3),ppref(3),tpref, & xp2ref(3),pp2ref(3),tp2ref,ref2par end if fsavfirst = .false. end if ! if( ipstart.eq.0 .and. ip.eq.0 ) then if( fsavset ) then write(4,'(6e14.6,i5)') 0.,ppref(3),0., 0.,0.,0.,refpar else write(4,'(6e14.6,i5)') xpref(3),ppref(3),tpref, 0.,0.,0., & refpar end if fsavfirst = .false. end if ! if( ipstart .eq. 1 ) then write(4,'(6e14.6,i5)') 0.,0.,0., 0.,0.,0., 0 fsavfirst = .false. end if ! end if ! if( ip .gt. 0 ) then ! leave ref particle out of distributions if( iflag .eq. 0 ) then ! good track call FILL_HIST(i7) call FILL_SCAT(i7) if( spin ) call FILL_POL(i7) call SAVE_FILE(i7) c else ! bad track nfail = nfail + 1 if( nfail .ge. npart ) go to 900 ! call FILL_HIST( iflag ) call FILL_SCAT( iflag ) call SAVE_FILE( iflag ) end if ! end of test on iflag = 0 end if ! end of test on ip .gt. 0 ! c -------------------------------------------------------- ! 500 continue ! end of loop over particles ! -------------------------------------------------------- 505 continue ! ! Space charge kick for end of region if( spcharge .and. jsec.gt.0 ) then call SC_KICK(2) end if ! ! Test if we have passed region of validity for background field if( jsec .gt. 0 ) then saccum = saccum + slen ! diff = saccum - sbg0 - ztotalbkg ! if( diff.gt.0. .or. ABS(diff).le.1e-6 ) then ! bftag = 'NONE' ! tell FIELD to ignore background ! end if end if ! 580 continue ! ! Calculate diagnostics at end of this region ! if( jsec.ge.0 .or. jpsr.eq.-7 ) then loutput = .false. ! if( nemit .gt. 0 ) then call EMIT2D call EMITCOV end if ! if( ncovar .gt. 0 ) call COVARIANCE ! if( nrhist .gt. 0 ) call RHISTORY ! end if ! steprk = steprksav ! ........................................................ ! go to 100 ! get next "region" command ! 600 continue ! end of loop over sections ! -------------------------------------------------------- n7records = i7 620 continue ! if( fsav .and. izfile.eq.-1 ) then ! write(4,'(7e14.6)') 0.,0.,0., 0.,0.,0., 0. ! end if ! ! Save information on the production properties of good tracks do 650 ip=1,npart ! call READ_EVT(ip) ! if( ipflg .ne. 0 ) go to 650 call FILL_HIST(-1) call FILL_SCAT(-1) if( spin ) call FILL_POL(-1) call SAVE_FILE(-1) 650 continue ! 900 continue c if( iflag .eq. -1 ) then write(*,*) 'FATAL ERROR IN ICOOL: CHECK LOG FILE' end if ! if( termout ) then write(*,*)' Total number of initial good particles = ',ngood write(*,*)' Total number of final good particles = ',npart-nfail end if write(2,*)' Total number of initial good particles = ',ngood write(2,*)' Total number of final good particles = ',npart-nfail c return end ! ********************************************************* subroutine STEPSIZE(pm0) ! ! Determine step size to use for next step ! ! hdid is the spatial step taken along reference orbit ! sdid is the spatial step taken along the particle's trajectory ! implicit double precision (a-h,o-z) include 'icommon.inc' ! ! if( .not. varstep ) then ! use fixed step size ! hdid = fixstep ! else ! use adaptive step size ! sacc = xp0(3) - sr0 dsleft = slen - sacc step = MIN(dsleft,hregn,hnextf,hnextde,hnextms,stepmax) step = MAX(step,stepmin) ! if( step .lt. dsleft ) then hdid = step else hdid = dsleft end if ! end if ! if( ncurved .eq. 0 ) then cth = pp0(3) / pm0 sdid = hdid / cth else xpr = pp0(1) / pp0(3) ypr = pp0(2) / pp0(3) sdid = SQRT( xpr**2 + ypr**2 + (1.d0+htrk*xp0(1))**2) * hdid if( nbffag .eq. 1.d0 ) then ! DIP model 3 zpr = 1.d0 + htrk*xp0(1) + dsffag*xp0(1) sdid = SQRT( xpr**2 + ypr**2 + zpr**2) * hdid end if end if ! return end c ********************************************************* subroutine STRIPCOM c c Stripcom takes a commented ICOOL control file (for001.dat) which c may contain comment lines, blank lines, and end-of-line comments, c and outputs a standard ICOOL control file in which comments and c blank lines are removed. c ! From Steve Bracker 4 March 2002 ! c Definitions: c c A "line" consists of zero or more characters of text followed by c an end-of-line terminator (typically carriage return, line feed, c or both). In fortran formatted i/o, the terminator is invisible c and may be ignored; hereafter "line" refers only to the "visible" c characters preceding the terminator. Lines are limited to being c no more than 132 characters long, though that can be changed by c redefining parameter maxLineLength and recompiling. c c A "whitespace" is a character which is rendered as one or more c empty spaces. The only two whitespace characters recognized here c are space and horizontal tab. c c A "trailing whitespace" is a whitespace which occurs to the c right of any non-whitespace character. c c A "blank line" is a line which consists of (1) no visible c characters, or (2) one or more whitespaces but no other visible c characters. c c A "comment" is a block of visible text within a line which is c present only for the edification of human readers. Its presence c and content is at best irrelevant and at worst confusing to a c program reading and parsing the text. c c A "start-of-comment character" is a visible character which is c guaranteed to occur only to mark the leftmost character in a c comment string. For ICOOL, the start-of-comment character is c an exclamation point (!). c c A "comment line" is a line consisting only of a comment, c preceded by zero or more whitespaces. c c An "end-of-line comment" is a comment that follows non-whitespace c text. c c Subroutine stripcom processes the input text line by line. Each c line of input produces zero or one lines of output. The c processing rules are: c c 1. Blank lines are removed (no output line) c 2. Comment lines are removed (no output line) c 3. Other lines are copied to the output, but trailing c whitespaces and end-of-line comments are removed. c c c Files: c c Input file is unit 1, for001.dat in the default directory. c Opened for reading; unaltered; closed prior to return. c c Output file is unit 9, for001.dcr in the default directory. c Opened for writing with status 'unknown' (created if it doesn't c already exist and overwritten if it does); written with a copy c of the input file, with blank lines, comment lines, end-of-line c comments and trailing whitespace removed; closed prior to c return. (ICOOL should open for001.dcr on unit 1 just as it now c opens for001.dat; unit 9 must be available for use, but is c released prior to return.) c c Error reports are written to unit 2 (which must be open for c writing to the report file, for002.dat) and to standard c output *. c implicit none include 'icommon.inc' c c Maximum line length in characters integer maxLineLength parameter (maxLineLength = 1024) c c Character strings: space, tab, start-of-comment c It's simplest to define tabChar in ordinary code below. character*1 spaceChar, tabChar, socChar parameter (spaceChar = ' ') parameter (socChar = '!') c c status returned from i/o operations integer ios c c line buffer character*(maxLineLength) line c c number of input characters to be output integer lineLength c c index over character strings integer i c Declare function types integer LASTNB c c c Ascii character code for horizontal tab tabChar = char(9) c c>> Open the input file for001.dat open (unit=1, & file= pathname(1:LASTNB(pathname)) // 'for001.dat', & status='old', & iostat=ios) c c If there is a problem opening the file . . . if (ios .ne. 0) then c c Report a fatal error and quit write (2,800) ios if (termout) write (*,800) ios 800 format & (/' STRIPCOM: Fatal error opening file for001.dat'/ & ' Runtime error code = ', i6/ & ' Execution terminated') close (unit=2) stop end if c c>> Open the output file for001.dcr open (unit=9, & file= pathname(1:LASTNB(pathname)) // 'for001.dcr', & status='unknown', & iostat = ios) c c If there is a problem opening the file . . . if (ios .ne. 0) then c c Report a fatal error and quit write (2,801) ios if (termout) write (*,801) ios 801 format & (/' STRIPCOM: Fatal error opening file for001.dcr'/ & ' Runtime error code = ', i6/ & ' Execution terminated') close (unit=2) stop end if c c Line processing loop 100 continue c>> Read a line from for001.dat; set the line length to maximum. c End-of-file and error-on-read processing is deferred. read (1,802, end=500, err=510) line 802 format (a) lineLength = maxLineLength c c>> Find the leftmost start-of-comment, if any. If there is one, c reduce the line length accordingly (may become zero). do i = 1, lineLength if (line(i:i) .eq. socChar) then lineLength = i-1 c If this is a comment line, we're done with it if (lineLength .le. 0) go to 100 c We don't care what's beyond socChar; next processing step go to 105 end if end do c 105 continue c c>> Find the length of the remaining line, with trailing c whitespace removed. do i = lineLength, 1, -1 if (line(i:i) .ne. spaceChar .and. & line(i:i) .ne. tabChar) then lineLength = i go to 110 end if end do c c All that remains is whitespace, so we're done with this line go to 100 c 110 continue c c>> Output the text remaining after stripping write (9, 803, err=520) line(1:lineLength) 803 format (a) c c>> Next input line . . . go to 100 c c c>> End-of-file reading input; close files, report, and return. 500 close (unit = 1) close (unit = 9) write (2,805) ! if (termout) write (*,805) 805 format & (' STRIPCOM: for001.dat comments stripped successfully') return c c c>> Error reading input 510 write (2,806) ios if (termout) write (*,806) ios 806 format & (/' STRIPCOM: Fatal error reading file for001.dat'/ & ' Runtime error code = ', i6/ & ' Execution terminated') close (unit=2) stop c c>> Error writing output file 520 write (2,807) ios if (termout) write (*,807) ios 807 format & (/' STRIPCOM: Fatal error writing file for001.dcr'/ & ' Runtime error code = ', i6/ & ' Execution terminated') close (unit=2) stop end c ********************************************************* subroutine SUBNAMES ! ! Substitute named variables with a fixed or random value ! ! ICOOL uses the expression ! &SUB NAME VALUE ! to associate NAME with a fixed VALUE ! ! The name is used in the file by ! &NAME ! ! Previously defined substitution values can be scaled using ! &SCL char type factor ! implicit double precision (a-h,o-z) include 'icommon.inc' c c Maximum line length in characters parameter (maxLineLength = 1024) ! Maximum number of characters in names parameter (maxNameLength = 20 ) ! Maximum number of characters in values parameter (maxValueLength = 30 ) ! Maximum number of substitutions parameter (maxSubs = 100 ) c character*1 spaceChar, subChar,ch,ch0,sctyp parameter (spaceChar = ' ') parameter (subChar = '&' ) ! Line buffers character*(maxLineLength) line,line2 ! Name buffers ! character*(maxNameLength) sname(maxSubs),rname(maxSubs),nam character*(maxNameLength) sname(maxSubs),nam character*(maxNameLength) schr ! Value buffers character*(maxValueLength) svalue(maxSubs),scfac ! character*(maxValueLength) rtype(maxSubs) ! character*(maxValueLength) rv1(maxSubs),rv2(maxSubs) character*30 token(50) logical blank,done external blank ! c>> Open the input file for001.dcr open (unit=1, & file= pathname(1:LASTNB(pathname)) // 'for001.dcr', & status='old', & iostat=ios) c c If there is a problem opening the file . . . if (ios .ne. 0) then c c Report a fatal error and quit write (2,20) ios if (termout) write (*,20) ios 20 format & (/' SUBNAMES: Fatal error opening file for001.dcr'/ & ' Runtime error code = ', i6/ & ' Execution terminated') close (unit=2) stop end if c c>> Open the output file for001.dsb open (unit=9, & file= pathname(1:LASTNB(pathname)) // 'for001.dsb', & status='unknown', & iostat = ios) c c If there is a problem opening the file . . . if (ios .ne. 0) then c c Report a fatal error and quit write (2,21) ios if (termout) write (*,21) ios 21 format & (/' SUBNAMES: Fatal error opening file for001.dsb'/ & ' Runtime error code = ', i6/ & ' Execution terminated') close (unit=2) stop end if c nfsub = 0 done = .false. c Line processing loop ! --------------------------------------------------------- 100 continue c>> Read a line from for001.dcr if( done ) go to 285 read(1,'(a)',iostat=ioc) line ! if( ioc .lt. 0 ) then ! EOF go to 285 ! done else if( ioc .gt. 0 ) then ! read error write (2,116) ioc if (termout) write (*,116) ioc 116 format & (/' SUBNAMES: Fatal error reading file for001.dcr'/ & ' Runtime error code = ', i6/ & ' Execution terminated') close (unit=2) stop end if if( line(1:4) .eq. 'ENDS' ) done = .true. ! ! Look for a substitution definition. If found, save the ! name and value in a table lineLength = LASTNB(line) ! do i=1,lineLength if( line(i:i) .eq. subChar ) then ! found subs marker if( line(i+1:i+3) .eq. 'SUB' ) then ! new definition nfsub = nfsub + 1 ! Look for blanks to extract name and value do j=i+4,lineLength if( BLANK(line(j:j)) ) then jb = j ! last blank else go to 200 end if end do go to 500 ! error ! 200 continue do j=jb+1,lineLength if( BLANK(line(j:j)) ) then ! found next blank jb2 = j sname(nfsub) = line(jb+1:jb2-1) go to 210 end if end do go to 500 ! error ! 210 continue do j=jb2,lineLength if( BLANK(line(j:j)) ) then jb = j ! last blank else go to 220 end if end do go to 500 ! error ! 220 continue do j=jb+1,lineLength+1 ! allow trailing blank if( BLANK(line(j:j)) ) then svalue(nfsub) = line(jb+1:j-1) go to 100 ! get next line end if end do go to 500 ! error ! ......................................................... else if( line(i+1:i+3) .eq. 'SCL' ) then ! scale values ! Look for blanks to extract parameters do j=i+4,lineLength if( BLANK(line(j:j)) ) then jb = j ! last blank else go to 272 end if end do go to 500 ! error ! 272 continue do j=jb+1,lineLength if( BLANK(line(j:j)) ) then ! found next blank jb2 = j schr = line(jb+1:jb2-1) nscn = jb2 - jb - 1 go to 274 end if end do go to 500 ! error ! 274 continue do j=jb2,lineLength if( BLANK(line(j:j)) ) then jb = j ! last blank else go to 276 end if end do go to 500 ! error ! 276 continue do j=jb+1,lineLength+1 ! allow trailing blank if( BLANK(line(j:j)) ) then jb2 = j if( j-1 .ne. jb+1 ) go to 500 ! error sctyp = line(jb+1:j-1) go to 278 end if end do go to 500 ! error ! 278 continue do j=jb2,lineLength if( BLANK(line(j:j)) ) then jb = j ! last blank else go to 280 end if end do go to 500 ! error ! 280 continue do j=jb+1,lineLength+1 ! allow trailing blank if( BLANK(line(j:j)) ) then jb2 = j scfac = line(jb+1:j-1) go to 282 end if end do go to 500 ! error ! 282 continue do j=1,nfsub if( sname(j)(1:nscn) .eq. schr ) then read(svalue(j),*) r1 read(scfac,*) r2 if( sctyp .eq. '*' ) then r1 = r1 * r2 else r1 = r1 + r2 end if write(svalue(j),'(e14.6)') r1 end if end do go to 100 end if ! end of test on command end if ! end of test on substitution character end do ! end of loop on characters in this line ! go to 100 ! end of loop over lines ! --------------------------------------------------------- ! 285 continue rewind(1) ! Process the title line separately. Allow & in the title. read(1,'(a)',iostat=ioc) line write (9,'(a)') line(1:LASTNB(line)) done = .false. ! 290 continue if( nfsub .eq. 0 ) then ! just make a copy if( done ) go to 800 read(1,'(a)',iostat=ioc) line ! if( ioc .lt. 0 ) then ! EOF go to 800 ! done else if( ioc .gt. 0 ) then ! read error write (2,116) ioc if (termout) write (*,116) ioc close (unit=2) stop end if if( line(1:4) .eq. 'ENDS' ) done = .true. ! write (9,'(a)') line(1:LASTNB(line)) go to 290 end if ! ! --------------------------------------------------------- ! Line processing loop when substitutions are present done = .false. ! 300 continue c>> Read a line from for001.dcr if( done ) go to 800 read(1,'(a)',iostat=ioc) line ! if( ioc .lt. 0 ) then ! EOF go to 800 ! done else if( ioc .gt. 0 ) then ! read error write (2,116) ioc if (termout) write (*,116) ioc close (unit=2) stop end if if( line(1:4) .eq. 'ENDS' ) done = .true. ! lineLength = LASTNB(line) ! Skip substitution definition lines do i=1,lineLength if( line(i:i+3).eq.'&SUB' & .or. line(i:i+3).eq.'&SCL' ) go to 300 end do ! ......................................................... ! Break up the line into character tokens ntokens = 0 ch0 = line(1:1) if( .not.BLANK(ch0) ) kt1 = 1 ! do i=2,lineLength+1 ! allow trailing blank ch = line(i:i) if( BLANK(ch) ) then ! this character is blank if( .not.BLANK(ch0) ) then ! prev character was not blank ntokens = ntokens + 1 token(ntokens) = line(kt1:i) end if else ! this character is not blank if( BLANK(ch0) ) then ! previous character was blank kt1 = i end if end if ch0 = line(i:i) end do ! ......................................................... ! Check if any tokens are substitution names ! msub = 0 do 350 i=1,ntokens if( token(i)(1:1) .eq. subChar ) then ! msub = msub + 1 mtoken = LASTNB(token(i)) nam = token(i)(2:mtoken) ! desired name ! do k=1,nfsub if( nam .eq. sname(k) ) then ! found it token(i) = svalue(k) ! msub = msub - 1 go to 320 end if end do ! 320 continue end if 350 continue ! ......................................................... ! Put tokens back together into a new line do i=1,maxLineLength line2(i:i) = spaceChar end do ! mtt = LASTNB(token(1)) line2(1:mtt) = token(1)(1:mtt) j2 = mtt ! do i=2,ntokens line2(j2+1:j2+1) = ' ' j1 = j2 + 2 mtt = LASTNB(token(i)) j2 = j1 + mtt do j=j1,j2 line2(j1:j2) = token(i)(1:mtt) end do end do ! ! Copy the line write (9,'(a)') line2(1:LASTNB(line2)) go to 300 ! --------------------------------------------------------- ! 500 continue ! error write (2,516) if (termout) write (*,116) 516 format & (/' SUBNAMES: Fatal execution error'/ & ' Execution terminated') close (unit=2) stop ! 800 continue close(1) close(9) write (2,805) 805 format & (' SUBNAMES: for001.dat substitutions replaced successfully') ! end c ********************************************************* subroutine TILT_COORD(psi,fi) ! ! Rotate particle coordinates by angle ANG1 and ANG2 around some axis ! implicit double precision (a-h,o-z) include 'icommon.inc' ! dimension x(3),p(3) ! ! psi=psi*ran1(idum) ! tild angle ! fi=fi*ran1(idum) ! direction of rotation respect to axis x in plane (x-y) cospsi = COS(psi) sinpsi = SIN(psi) cosfi = COS(fi) sinfi = SIN(fi) do i=1,3 x(i) = xp(i) p(i) = pp(i) end do ! c need to find z z=x(3)-left end of the element c particle is at the left end of the lement hence z=0.d0 xp(1)=x(1)*cospsi+(1.d0-cospsi)*cosfi*(x(1)*cosfi+ + x(2)*sinfi)-z*sinpsi*sinfi xp(2)=x(2)*cospsi+(1.d0-cospsi)*sinfi*(x(1)*cosfi+ + x(2)*sinfi)+z*sinpsi*cosfi z=z*cospsi+sinpsi*(x(1)*sinfi-x(2)*cosfi) c pp(1)=p(1)*cospsi+(1.d0-cospsi)*cosfi*(p(1)*cosfi+ + p(2)*sinfi)-p(3)*sinpsi*sinfi pp(2)=p(2)*cospsi+(1.d0-cospsi)*sinfi*(p(1)*cosfi+ + p(2)*sinfi)+p(3)*sinpsi*cosfi pp(3)=p(3)*cospsi+sinpsi*(p(1)*sinfi-p(2)*cosfi) ! if( z .lt. 0.d0 ) z = 0.d0 xp(3)=x(3)+z ! return end ! ********************************************************* subroutine TRANSPORT ! ! Transport particle using user defined transport matrix ! implicit double precision (a-h,o-z) include 'icommon.inc' ! dimension u(6),v(6) ! ! Convert ICOOL variables into TRANSPORT form (x,x',y,y',dl,dp/p) pref = ppref(3) call VMAG(pp,pm) gamma = SQRT( 1.d0 + (pm/pmass)**2 ) vpar = pm / pmass / gamma * cc v(1) = xp(1) v(2) = pp(1) / pref v(3) = xp(2) v(4) = pp(2) / pref v(5) = (tp - tpref) * vpar v(6) = (pm - pref) / pref ! ! Transport particle to end of region: u = T v do i=1,6 u(i) = 0.d0 do j=1,6 u(i) = u(i) + trprt(i,j) * v(j) end do end do ! ! Convert back to ICOOL variables ! Don't change z xp(1) = u(1) xp(2) = u(3) tp = tpref + u(5)/vpar pp(1) = u(2) * pref pp(2) = u(4) * pref pm = pref * (1.d0 + u(6)) pp(3) = pm*pm - pp(1)**2 - pp(2)**2 if( pp(3) .lt. 0.d0 ) then ! unphysical momentum iflag = -87 else pp(3) = SQRT( pp(3) ) end if ! return end ! ******************************************************** subroutine TWISS(jj,kk) ! ! Sample initial phase space using Twiss transformed ellipse ! ! Ref. Weidemann, Vol 1, Eq. 5-91 and 5-98 ! implicit double precision (a-h,o-z) include 'icommon.inc' ! alp = corr1(jj,kk) ! alpha bet = corr2(jj,kk) ! beta eps = corr3(jj,kk) ! epsilon ic = corrtyp(jj,kk) ! correlation type ibd = bdistyp(kk) ! beam distribution type ! gam = (1.d0 + alp**2) / bet ! if( ic .eq. 6 ) then ! x Px ! if( ibd .eq. 1 ) then ! gaussian sd = SQRT( bet*eps ) ! rms x ! 80 dx = GAUSS(0.d0,sd) xp(1) = x1bt(1,kk) + dx ! arg = (alp*dx)**2-bet*(gam*dx*dx-eps) if( arg .lt. 0d0 ) go to 80 ! xpr1 = -alp*dx + SQRT(arg) xpr1 = xpr1 / bet xpr2 = -alp*dx - SQRT(arg) xpr2 = xpr2 / bet xprm = (xpr1 + xpr2) / 2.d0 pxm = xprm * pp(3) ! average px sdp = (xpr1 - xpr2)/2.d0 * pp(3) ! dp = GAUSS( pxm,sdp ) ! pp(1) = p1bt(1,kk) + dp ! else ! uniform ellipse xmax = SQRT( bet*eps ) dx = RANV(-xmax,xmax) xp(1) = x1bt(1,kk) + dx xpr1 = -alp*dx + SQRT((alp*dx)**2-bet*(gam*dx*dx-eps)) xpr1 = xpr1 / bet xpr2 = -alp*dx - SQRT((alp*dx)**2-bet*(gam*dx*dx-eps)) xpr2 = xpr2 / bet pp(1) = p1bt(1,kk) + RANV(xpr1,xpr2)*pp(3) end if ! else if( ic .eq. 7 ) then ! y Py if( ibd .eq. 1 ) then ! gaussian sd = SQRT( bet*eps ) ! 100 dx = GAUSS(0.d0,sd) xp(2) = x1bt(2,kk) + dx ! arg = (alp*dx)**2-bet*(gam*dx*dx-eps) if( arg .lt. 0d0 ) go to 100 ! xpr1 = -alp*dx + SQRT(arg) xpr1 = xpr1 / bet xpr2 = -alp*dx - SQRT(arg) xpr2 = xpr2 / bet xprm = (xpr1 + xpr2) / 2.d0 pxm = xprm * pp(3) sdp = (xpr1 - xpr2)/2.d0 * pp(3) ! dp = GAUSS( pxm,sdp ) ! pp(2) = p1bt(2,kk) + dp ! else ! uniform ellipse xmax = SQRT( bet*eps ) dx = RANV(-xmax,xmax) xp(2) = x1bt(2,kk) + dx xpr1 = -alp*dx + SQRT((alp*dx)**2-bet*(gam*dx*dx-eps)) xpr1 = xpr1 / bet xpr2 = -alp*dx - SQRT((alp*dx)**2-bet*(gam*dx*dx-eps)) xpr2 = xpr2 / bet pp(2) = p1bt(2,kk) + RANV(xpr1,xpr2)*pp(3) end if ! else if( ic .eq. 8 ) then ! z Pz if( ibd .eq. 1 ) then ! gaussian sd = SQRT( bet*eps ) ! 120 dx = GAUSS(0.d0,sd) xp(3) = x1bt(3,kk) + dx ! arg = (alp*dx)**2-bet*(gam*dx*dx-eps) if( arg .lt. 0d0 ) go to 120 ! xpr1 = -alp*dx + SQRT(arg) xpr1 = xpr1 / bet xpr2 = -alp*dx - SQRT(arg) xpr2 = xpr2 / bet xprm = (xpr1 + xpr2) / 2.d0 pxm = xprm * pp(3) sdp = (xpr1 - xpr2)/2.d0 * pp(3) ! dp = GAUSS( pxm,sdp ) ! pp(3) = p1bt(3,kk) + dp ! else ! uniform ellipse xmax = SQRT( bet*eps ) dx = RANV(-xmax,xmax) xp(3) = x1bt(3,kk) + dx xpr1 = -alp*dx + SQRT((alp*dx)**2-bet*(gam*dx*dx-eps)) xpr1 = xpr1 / bet xpr2 = -alp*dx - SQRT((alp*dx)**2-bet*(gam*dx*dx-eps)) xpr2 = xpr2 / bet pp(3) = p1bt(3,kk) + RANV(xpr1,xpr2)*pp(3) end if ! end if ! end of test on correlation type ! return end ! ********************************************************* subroutine VMAG(p,pm) c c Return magnitude of 3-vector p c implicit double precision (a-h,o-z) dimension p(3) c pm = sqrt( p(1)*p(1) + p(2)*p(2) + p(3)*p(3) ) c return end ! ********************************************************* subroutine WRITE_EVT(ip) ! ! Write particle data to common arrays or unit(8) ! implicit double precision (a-h,o-z) include 'icommon.inc' ! if( ip .le. mxpar ) then if( ip .eq. 0 ) then ! reference particle do i=1,3 xpref(i) = xp(i) ppref(i) = pp(i) end do tpref = tp ! ipnum = 1 iptypref = iptyp ipflgref = ipflg sarcref = sarc ! ......................................................... else if( ip .eq. -1 ) then ! 2nd reference particle do i=1,3 xp2ref(i) = xp(i) pp2ref(i) = pp(i) end do tp2ref = tp ! ipnum = 1 iptyp2ref = iptyp ipflg2ref = ipflg sarc2ref = sarc ! ......................................................... else ! normal particles do i=1,3 xpar(i,ip) = xp(i) ppar(i,ip) = pp(i) polpar(i,ip) = pol(i) xpintpar(i,ip) = xpint(i) ppintpar(i,ip) = ppint(i) polinitpar(i,ip) = polinit(i) end do ! tpar(ip) = tp tpintpar(ip) = tpint sarcpar(ip) = sarc bzfldpar(ip) = bfld(3) ievtpar(ip) = ievt ipnumpar(ip) = ipnum iptypar(ip) = iptyp ipflgpar(ip) = ipflg evtwtpar(ip) = evtwt end if ! end of test on ip = 0 ! else ! ip > mxpar ! write(8,rec=ip-mxpar) xp,pp,tp,pol,sarc,evtwt, & ievt,ipnum,iptyp,ipflg,bfld(3), & xpint,ppint,tpint,polinit end if ! return end