!MS$pagesize:59 program ICOOL c c Monte carlo simulation of muon ionization cooling c include 'icommon.inc' c version = 1.82 ! call SETUP if( iflag .ne. 0 ) go to 900 c call SIMULATE c call FINISH c 900 continue end !MS$page c ********************************************************* BLOCK DATA ! ! Initialize common block data in F77-compatible manner ! include 'icommon.inc' ! ! backgrnd data bftag/'NONE'/,gridbkg/.false./ data ztotalbkg/0./,sbg0/0./ ! ! const data cc/2.9979246e8/,ee/0.29979246/,pi/3.1415926/, + avag/6.0221367e23/ & ,twopi/6.2831853/ data mass/0.51099906e-3, 0.10565839, 0.13956995, 0.493677, + 0.93827231/ ! GeV/c^2 ! ! control data bgen/.true./,eps/1e-3/,fsav/.false./,nprnt/0/,nsections/1/, & prlevel/1/,goodtrack/.true./,refpar/.true./,t0ref/0./ & phasemodel/1/,rnseed/-1/,neighbor/.false./,ntuple/.false./, & bunchcut/1.e6/,rtuple/.false./,rtuplen/5/,fsavset/.false./, & rfdiag/0/,rfphase/0/,spin/.false./,timelim/120./, & varstep/.true./ ! ! dvars data labvar/' X',' Y',' Z',' Px',' Py',' Pz',' t',' Pm', & ' E',' KE', ' X''',' Y''',' TH',' R','PHI', & ' Pr','Pph',' Lz',' L2',' S', & ' X0',' Y0',' Z0','Px0','Py0','Pz0',' t0',' A2', & ' R2',' ',' Bx',' By',' Bz',' Ex',' Ey',' Ez', & ' Sx',' Sy',' Sz','phs'/ ! ! emits data emevt/30*0/,emevcut/30*0/ data sigmacut/.true./,pxycorr/.false./,pzcorr/.false./ data sig_cut/4./ ! ! hist data ihcon/20*0/,ihund/20*0/,ihovf/20*0/,hval/1000*0/ ! ! intdat data ldedx/.false./,lscatter/.false./,lbrem/.false./, + linteract/.false./,ldecay/.false./,lspace/.false./, & delev/1/,scatlev/1/,spacelev/1/, + bremlev/1/,intlev/1/,declev/1/ ! ! matdata data nmat/9/,matid/'LH','LI','BE','B','C','AL','TI','FE','CU', & 'W','LIH'/ ! data order: Z A rho IP Lr Xo X1 Sa Sc Sm ! : nelcomp { zelcomp(i), felcomp(i) } data matdat/ & 1.,1.008,0.071,21.8,890.,0.4759,1.9215,0.1348,-3.2632, & 5.6249, 1,1.,1.,0.,0., ! Liq H & 3.,6.94,0.534,40.0,155., 0.1304,1.6397,0.9514,-3.1221, & 2.4993, 1,3.,1.,0.,0., ! Li & 4.,9.01,1.848,63.9,35.3, 0.0592,1.6922,0.8039,-2.7847, & 2.4339, 1,4.,1.,0.,0., ! Be & 5.,10.811,2.370,76.0,22.2,0.0305,1.9688,0.5622,-2.8477, & 2.4512, 1,5.,1.,0.,0., ! B & 6.,12.011,2.265,78.0,18.8,-0.0178,2.3415,0.2614,-2.8680, & 2.8697, 1,6.,1.,0.,0., ! C & 13.,26.982,2.699,166.,8.9,0.1708,3.0127,0.0802,-4.2395, & 3.6345, 1,13.,1.,0.,0., ! Al & 22.,47.88,4.54,233.,3.56,0.0957,3.0386,0.1566,-4.4450, & 3.0302, 1,22.,1.,0.,0., ! Ti & 26.,55.845,7.87,286.,1.76,-0.0012,3.1531,0.1468,-4.2911, & 2.9632, 1,26.,1.,0.,0., ! Fe & 29.,63.546,8.96,322.,1.43,-0.0254,3.2792,0.1434,-4.4190, & 2.9044, 1,29.,1.,0.,0., ! Cu & 74.,183.85,19.3,722.,0.35, 0.2167,3.4960,0.1551,-5.406, & 2.845, 1,74.,1.,0.,0., ! W & 4.,7.949,0.820,36.5,97.1, -0.0988,1.4515,0.9057,-2.3580, & 2.5849, 2,3.,1.,1.,1./ ! LiH ! ! polar data polx/30*0./,poly/30*0./,polz/30*0./ ! ! regndat data nfparm/15/,ngparm/10/ ! ! scat data iscon/20*0/,isxund/20*0/,isyund/20*0/,isxovf/20*0/, & isyovf/20*0/,sval/23000*0/ ! ! simul data krad/1/ ! ! zhist data izhknt/130*0/ ! end !MS$page c ********************************************************* subroutine BACKGET(x,t,e,b) ! ! Get background field ! include 'icommon.inc' ! real*8 x(3),t,e(3),b(3),xloc(3),bb(3) ! dum = t ! defeat annoying compiler warnings dum = e(1) ! ! --------------------------------------------------------- if( gridbkg ) then ! Interpolate field from background grid x4 = x(1) y4 = x(2) z4 = x(3) call HUNT(xbkg,nxbkg,x4,jx) ! returns lower grid corner call HUNT(ybkg,nybkg,y4,jy) call HUNT(zbkg,nzbkg,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-xbkg(jx)) / dxbkg v = (y4-ybkg(jy)) / dybkg w = (z4-zbkg(jz)) / dzbkg ! call TRILINEAR(u,v,w,jx,jy,jz,bxbkg,mxxg,mxyg,mxzg,b(1)) call TRILINEAR(u,v,w,jx,jy,jz,bybkg,mxxg,mxyg,mxzg,b(2)) call TRILINEAR(u,v,w,jx,jy,jz,bzbkg,mxxg,mxyg,mxzg,b(3)) ! else if( interbkg.eq.2 .or. interbkg.eq.3 ) then ! call TRIPOLY(jx,jy,jz,xbkg,ybkg,zbkg,bxbkg,mxxg,mxyg,mxzg, & interbkg,x4,y4,z4,b(1),df) call TRIPOLY(jx,jy,jz,xbkg,ybkg,zbkg,bybkg,mxxg,mxyg,mxzg, & interbkg,x4,y4,z4,b(2),df) call TRIPOLY(jx,jy,jz,xbkg,ybkg,zbkg,bzbkg,mxxg,mxyg,mxzg, & interbkg,x4,y4,z4,b(3),df) ! end if ! ! ......................................................... else ! Compute field analytically ! do i=1,3 b(i) = 0. end do ! do i=reclobkg,rechibkg read(7,rec=i) dum,dum,zoffbkg,rmaxbkg,zcutbkg,bftag,bfparm xloc(1) = x(1) xloc(2) = x(2) xloc(3) = x(3) - sbg0 - zoffbkg call FIELDFNC(bftag,bfparm,xloc,t,e,bb) do j=1,3 b(j) = b(j) + bb(j) end do end do ! end if ! --------------------------------------------------------- htrk = ee / prefbkg * b(2) hptrk = 0. ordtrk = 1. ! return end !MS$page c ********************************************************* subroutine BACKSET(iarg) ! ! Generate background field grid ! include 'icommon.inc' ! real*8 x(3),t,e(3),b(3),xloc(3) ! character*10 fname ! if( iarg .eq. 0 ) then ! initialize ! do ix=1,nxbkg xbkg(ix) = (ix-1) * dxbkg + xlobkg end do do iy=1,nybkg ybkg(iy) = (iy-1) * dybkg + ylobkg end do do iz=1,nzbkg zbkg(iz) = (iz-1) * dzbkg + zlobkg end do ! do ix=1,nxbkg do iy=1,nybkg do iz=1,nzbkg bxbkg(ix,iy,iz) = 0. ! clear background array bybkg(ix,iy,iz) = 0. bzbkg(ix,iy,iz) = 0. ! exbkg(ix,iy,iz) = 0. ! eybkg(ix,iy,iz) = 0. ! ezbkg(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 = zoffbkg ! write(fname,40) mfile !40 format('for0',i2,'.dat') ! open(unit=mfile,file=fname,iostat=ioc) ! if( ioc .ne. 0 ) then ! iflag = -48 ! go to 900 ! end if ! ! FILE FORMAT NEEDS TO BE SPECIFIED AND READ IN HERE ! close(mfile) ! ! else ! use predefined field types ! do 250 ix=1,nxbkg x(1) = xbkg(ix) ! do 230 iy=1,nybkg x(2) = ybkg(iy) r = SQRT(x(1)**2 + x(2)**2) if( r .gt. rmaxbkg ) go to 230 ! do 210 iz=1,nzbkg x(3) = zbkg(iz) if( ABS(x(3)-zoffbkg) .gt. zcutbkg ) go to 210 ! ! Translate into magnet local coordinate system xloc(1) = x(1) xloc(2) = x(2) xloc(3) = x(3) - sbg0 - zoffbkg ! call FIELDFNC(bftag,bfparm,xloc,t,e,b) ! bxbkg(ix,iy,iz) = bxbkg(ix,iy,iz) + b(1) bybkg(ix,iy,iz) = bybkg(ix,iy,iz) + b(2) bzbkg(ix,iy,iz) = bzbkg(ix,iy,iz) + b(3) ! exbkg(ix,iy,iy) = exbkg(ix,iy,iz) + e(1) ! eybkg(ix,iy,iy) = eybkg(ix,iy,iz) + e(2) ! ezbkg(ix,iy,iy) = ezbkg(ix,iy,iz) + e(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 !MS$page ! ********************************************************* 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} ! real*4 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 !MS$page c ********************************************************* subroutine CHKTIME(tlim) c c Check elapsed time from start of execution c time calculated in minutes c 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) iflag = -10 c return end !MS$page c ********************************************************* SUBROUTINE DCROSS(A,B,C) c c Vector cross product real*8 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 !MS$page c ********************************************************* real*8 function DDOT(a,b) ! ! Double precision vector dot product ! real*8 a(3),b(3) ! ddot = 0. ! do i=1,3 ddot = ddot + a(i)*b(i) end do ! return end !MS$page c ********************************************************* subroutine DERIV(x,y,dydx) c c Compute values of derivatives DYDX for the variables Y at X c include 'icommon.inc' c real*8 x,y(9),dydx(9) real*8 gamma,mpgam,vm,ddot,v(3),fac(3),exv(3) real*8 px,py,ps,pm real*8 u,term,xterms,yterms,vv real*8 cxbx,cxby,cxbz,cybx,cyby,cybz data g/1.166e-3/ ! anomalous g-factor of muon ! ! Use space stepping ! X = s ! Y = {x,y,t,x',y',Ps...} ! dydx(1) = y(4) ! dx/ds = x' dydx(2) = y(5) ! dy/ds = y' ! ps =y(6) px = ps * y(4) py = ps * y(5) pm = SQRT( px**2 + py**2 + ps**2 ) gamma = SQRT(1. + pm*pm/mp/mp) mpgam = mp * gamma v(1) = px / mpgam * cc v(2) = py / mpgam * cc v(3) = ps / mpgam * cc dydx(3) = 1. / v(3) ! dt/ds = 1 / vs ! xp(1) = y(1) xp(2) = y(2) xp(3) = x ! call FIELD( xp(1),y(3) ) ! h = htrk hp = hptrk ! hx = h * y(1) u = 1. + hx ! ........................................................ if( ordtrk .eq. 1 ) then xterms = h * u yterms = 0. cxbx = 0. cxby = 1. + 2.*hx cxbz = -y(5) cybx = 1. + 2.*hx cyby = 0. cybz = -y(4) ! else if ( ordtrk .eq. 2 ) then term = 2.*h*y(4) + hp * y(1) xterms = h*u + y(4)*term yterms = y(5)*term cxbx = -y(4)*y(5) cxby = 1. + 2.*hx +1.5*hx*hx + 1.5*y(4)**2 + 0.5*y(5)**2 cxbz = -y(5)*u cybx = 1. + 2.*hx + 1.5*hx*hx + 0.5*y(4)**2 + 1.5*y(5)**2 cyby = -y(4)*y(5) cybz = -y(4)*u ! else if ( ordtrk .eq. 3 ) then vv = 1. - hx term = (2.*h*y(4) + hp * y(1)) * vv xterms = h*u + y(4)*term yterms = y(5)*term cxbx = -y(4)*y(5) cxby = 1. + 2.*hx + hx*hx - 0.5*hx**3 + 1.5*y(4)**2 + & 0.5*y(5)**2 cxbz = -y(5) * (1. - hx + 0.5*y(5)**2 + 0.5*y(4)**2) cybx = 1. + 2.*hx + hx*hx + 0.5*y(4)**2 + 1.5*y(5)**2 - & 0.5*hx**3 cyby = -y(4)*y(5) cybz = -y(4)*(1. + hx + 0.5*y(4)**2 + 0.5*y(5)**2) ! else ! use exact equations vv = SQRT( y(4)**2 + y(5)**2 + u*u ) term = (2.*h*y(4) + hp * y(1)) / u xterms = h*u + y(4)*term yterms = y(5)*term cxbx = -y(4)*y(5) * vv/u cxby = (1. + hx + y(4)**2/u) * vv cxbz = -y(5) * vv cybx = (1. + hx + y(5)**2/u) * vv cyby = -y(4)*y(5) * vv/u cybz = -y(4) * vv ! end if ! end of test on tracking order ! ........................................................ term = ee * efld(3) / ps / ps / v(3) ! dydx(4) = xterms - ee / pm * & (cxbx*bfld(1)+cxby*bfld(2)+cxbz*bfld(3)) & + ee*efld(1)/ps/v(3) - px*term dydx(5) = yterms + ee / pm * & (cybx*bfld(1)+cyby*bfld(2)+cybz*bfld(3)) & + ee*efld(2)/ps/v(3) - py*term ! term = ee / ps / v(3) ! dydx(6) = ( term* (ps*efld(3)+px*efld(1)+py*efld(2)) & - px*dydx(4) - py*dydx(5) ) / (1.+y(4)**2+y(5)**2) ! if( spin ) then hh = 1. / SQRT(1.+y(4)**2+y(5)**2) call VMAG(v,vm) call DCROSS(efld,v,exv) ! do i=1,3 fac(i) = (g*gamma+1.)*bfld(i) & - (gamma-1.)*g*DDOT(bfld,v)*v(i)/vm/vm & + gamma*( g+1./(gamma+1.) )*exv(i)/cc/cc fac(i) = -ee*hh/pm * fac(i) end do ! call DCROSS( fac,pol,dydx(7) ) end if ! return end !MS$page c ********************************************************* subroutine FIELD(xyz,tt) c c Return EM field at location xyz and time tt c include 'icommon.inc' c real*8 xyz(3),x(3),tt,e(3),b(3) real*4 fprm(15,4),fpar(15) integer*4 skip5(5),skip28(28) character*4 tag,ftg(4) ! do i=1,3 ! initialize; also zero field for drift space efld(i) = 0. bfld(i) = 0. x(i) = xyz(i) end do ! ! Region field if( ft .ne. 'NONE' ) then kbcell = 0 offset = sr0 x(3) = xyz(3) - offset ! transform to magnet coordinate system ! call FIELDFNC(ft,fparm(1,krad),x,tt,e,b) ! do i=1,3 efld(i) = efld(i) + e(i) bfld(i) = bfld(i) + b(i) end do end if ! ! Neighbor region fields if( neighbor ) then do 50 j=j7rec-1,1,-1 read(7,rec=j) js,skip5,sr,skip28,ftg,fprm if( js .lt. 0 ) go to 50 ! watch for pseudo-regions tag = ftg(1) if( tag .ne. 'NONE' ) then do i=1,nfparm fpar(i) = fprm(i,1) end do offset = sr x(3) = xyz(3) - offset ! transform to magnet system ! call FIELDFNC(tag,fpar,x,tt,e,b) ! do i=1,3 efld(i) = efld(i) + e(i) bfld(i) = bfld(i) + b(i) end do end if ! end of test on tag go to 55 50 continue ! 55 continue ! do 60 j=j7rec+1,n7records read(7,rec=j) js,skip5,sr,skip28,ftg,fprm if( js .lt. 0 ) go to 60 tag = ftg(1) if( tag .ne. 'NONE' ) then do i=1,nfparm fpar(i) = fprm(i,1) end do offset = sr x(3) = xyz(3) - offset ! transform to magnet system ! call FIELDFNC(tag,fpar,x,tt,e,b) ! do i=1,3 efld(i) = efld(i) + e(i) bfld(i) = bfld(i) + b(i) end do end if ! end of test on tag ! go to 65 60 continue ! 65 continue end if ! end of test on neighbor ! ! Cell field if( cftag .ne. 'NONE' ) then kbcell = 1 ! inform field routines that this is a cell field offset = sc0 x(3) = xyz(3) - offset ! transform to magnet coordinate system ! 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' ) then htrk = cellslp * htrk hptrk = cellslp * hptrk end if ! do i=1,3 efld(i) = efld(i) + e(i) bfld(i) = bfld(i) + b(i) end do end if ! ! Background field if( bftag.ne.'NONE' .and. cftag.ne.'BACK') then ! call BACKGET(xyz,tt,e,b) ! do i=1,3 efld(i) = efld(i) + e(i) bfld(i) = bfld(i) + b(i) end do end if ! return end !MS$page 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 ! real*8 x(3),t,e(3),b(3) real*4 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. '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. 'DEFL' ) then ! deflection cavity call DEFLECTION(fpar,x,tt,e,b) c else if( tag .eq. 'DIP' ) then ! sector dipole call DIPOLE(fpar,x,e,b) c else if( tag .eq. 'FOFO' ) then ! solenoid FOFO lattice call FOFO(fpar,x,e,b) ! 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. 'TROD' ) then ! tapered current rod call TAPERED_ROD(fpar,x,e,b) c else if( tag .eq. 'TSOL' ) then ! tapered solenoid call TAPSOL(fpar,x,e,b) c end if ! end of test for field tags ! return end !MS$page c ********************************************************* subroutine FINISH c c Take care of end of simulation tasks c include 'icommon.inc' c integer*2 iyr,imon,iday,ihr,imin,isec ! character*1 ff,cr ! data ff/#C/,cr/#D/ c ! write(2,'(2a1)')ff,cr call COVARIANCE 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 c 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 write(*,85)nh,nm,ns 85 format(' Total elapsed time = ',i4,' hours,',i4,' minutes,', + i4,' seconds') write(2,85)nh,nm,ns c do i=1,8 close(i) end do if( rfdiag .gt. 10 ) close(rfdiag) if( rfphase .gt. 10 ) close(rfphase) c if( nprnt.lt.0 ) then call SOUND(400,20) call SOUND(400,20) call SOUND(400,20) call SOUND(200,40) end if c return end !MS$page c ********************************************************* subroutine FSTEP(y,neqn,x,htry,hnext) c c Propagate particle one step thru field (or drift) c For adaptive step size control; cf Num Rec routine RKQC c include 'icommon.inc' c parameter ( grow = -0.20, shrink = -0.25, fcor = 1./15., + safety = 0.9, errcon = 1.89e-4 , nmax = 10) ! real*8 y(neqn),dydx(nmax),x,xsav,h,hh,htry,hnext,hold,yscal real*8 yt(nmax),ysav(nmax),dysav(nmax),errmax data yscal/1./ ! constant errors in all variables y data maxknt/1000/ ! h = htry call DERIV(x,y,dydx) ! xsav = x ! save initial values ! do i=1,neqn ysav(i) = y(i) dysav(i) = dydx(i) end do ! if( varstep ) then knt = 0 ! ........................................................ 10 continue knt = knt + 1 if( knt .gt. maxknt ) then iflag = -21 go to 900 end if c c Take two half-steps .................................. hh = 0.5 * h call RKSTEP(ysav,dysav,neqn,xsav,hh,yt) x = xsav + hh call DERIV(x,yt,dydx) call RKSTEP(yt,dydx,neqn,x,hh,y) x = xsav + h c c Take one large step .................................. call RKSTEP(ysav,dysav,neqn,xsav,h,yt) c ......................................................... errmax = 0. c do i=1,neqn yt(i) = y(i) - yt(i) ! yt holds error errmax = MAX( errmax, abs( yt(i)/yscal ) ) end do c errmax = errmax / eps ! scale relative to required tolerance if( errmax .gt. 1. ) then hold = h h = safety * h * (errmax**shrink) ! try smaller steps if( h .lt. 0.1*hold ) then h = 0.1 * hold ! no smaller than factor of 10 end if go to 10 ! ........................................................ else hdid = h ! this step OK; try bigger one next time if( errmax .gt. errcon ) then hnext = safety * h * (errmax**grow) else hnext = 5. * 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 hdid = h hnext = h ! end if ! 900 continue return end !MS$page 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 gauss = GASDEV(idum) c c renormalize to requested mean and standard deviation gauss = sd*gauss + xm c return end !MS$page c ********************************************************* subroutine GO_REGION c c Propagate particle thru specific s region c include 'icommon.inc' c real*8 x,y(9),htry,hnext,gamma,vz, + maxstepsize,rp,pm,pm0,xp0(3),ds c 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 knt = 0 call VMAG(pp,pm) gamma = SQRT( 1. + pm*pm/mp/mp ) vz = pp(3) / mp / gamma * cc if( vz .le. 0. ) then iflag = -32 go to 900 end if ! ! It defeats the step size adjustment if you put ! this inside the fine step loop below! ! ft = ftag(1) htry = zstep maxstepsize = slen / 5. maxknt = 100. * slen / htry ! htrk = 0. ! initialize for a straight region hptrk = 0. ordtrk = 1. ! sft = cftag krold = -1 ! flag when we change radial regions ! c --------------------------------------------------------- 100 continue ! loop over fine steps c ......................................................... c do 450 ir=1,nrreg ! loop over radial regions 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 krad = ir mt = mtag(ir) ft = ftag(ir) mgt = mgeom(ir) ! if( mt .ne. 'VAC') call LOAD_MATERIAL ! krold = ir end if ! 150 continue knt = knt + 1 if( knt .gt. maxknt ) then iflag = -33 go to 900 end if ! ! Compute initial values and derivatives for FSTEP ! x = xp(3) ! s y(1) = xp(1) ! x y(2) = xp(2) ! y y(3) = tp ! t y(4) = pp(1) / pp(3) ! x' y(5) = pp(2) / pp(3) ! y' y(6) = pp(3) ! Ps ! if( spin ) then do i=1,3 y(i+6) = pol(i) end do neqn = 9 else neqn = 6 end if ! call VMAG( pp,pm0 ) ! save initial magnitude pz0 = pp(3) do i=1,3 xp0(i) = xp(i) end do ! ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c Take a step thru appropriate field c call FSTEP(y,neqn,x,htry,hnext) c if( iflag .ne. 0 ) go to 900 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! c Check if we overshot the end of the region sacc = x - sr0 ! accumulated length along s ds = sacc - slen if( ds .gt. eps ) then ! overshoot htry = hdid - ds go to 150 end if ! ! 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) ! if( spin ) then do i=1,3 pol(i) = y(i+6) end do end if ! sarc = sarc + SQRT( (xp(1)-xp0(1))**2 + & (xp(2)-xp0(2))**2 + (xp(3)-xp0(3))**2 ) ! Ensure conservation of momentum call VMAG( pp,pm ) if( ft.ne.'ACCE' .and. ft.ne.'DEFL' ) then do i=1,3 pp(i) = pp(i) * pm0 / pm end do end if ! Use average p from the field step for material calculation pm = 0.5 * (pm0 + pm) ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c Take a step thru material c if ( mt .ne. 'VAC' ) then ! call MSTEP(pm) ! if( iflag .ne. 0 ) go to 900 end if ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c Check for other processes c if ( ldecay ) call DECAY ! particle decay ! if( iflag .ne. 0 ) go to 900 ! ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c if( jpar .le. ABS(nprnt) ) then if( knt .le. 300 ) then if( prlevel .gt. 1 ) then write(*,425)knt,tp,htry,hdid,hnext,xp,pp write(2,425)knt,tp,htry,hdid,hnext,xp,pp 425 format(' KNT,TP,HTRY,HDID,HNEXT:',i5,4e12.4,/, & ' XP,PP: ',6e12.4) end if c if( prlevel .gt. 2 ) then if( spin) then write(*,428) pol write(2,428) pol 428 format('POL:',3e12.4) end if write(*,430)efld,bfld write(2,430)efld,bfld 430 format(' E,B: ',6e12.4) end if c if( prlevel.gt.3 .and. xp(1)*xp(2).ne.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(*,435) rp,phi*180./pi,pr,pph,br,bph 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 c 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 ! ! Check if z-history is requested if( nzhist .gt. 0 ) call FILL_ZHIST ! if( rtuple .and. MOD(knt,rtuplen).eq.1 ) call OUTPUT ! ........................................................ ! Check if we are at end of region sacc = xp(3) - sr0 ! accumulated length along s ds = slen - sacc c if( ABS(ds) .le. eps ) go to 800 ! done c ........................................................ ! if( pp(3) .le. 0. ) then iflag = -32 go to 900 end if ! if( .not. varstep ) then ! use fixed step size if( zstep .lt. ds ) then htry = zstep else htry = ds end if ! else ! use adaptive step size if( hnext .lt. ds ) then htry = MIN( hnext, maxstepsize ) else htry = MIN( ds, maxstepsize ) end if end if ! go to 100 c c End of loop over fine steps c --------------------------------------------------------- 800 continue ! 900 continue if( jpar .le. abs(nprnt) ) then write(2,'(a,i5)')' Total steps = ',knt end if ! return end !MS$page c ********************************************************* subroutine LOAD_MATERIAL ! ! Load parameters for this material ! include 'icommon.inc' ! ! Load material parameters do i=1,nmat if( mt .eq. matid(i) ) then ! do j=1,10 mparm(j,krad) = matdat(j,i) end do ! nelcomp(krad) = matdat(11,i) do j=1,nelcomp(krad) zelcomp(j,krad) = matdat(10+2*j,i) felcomp(j,krad) = matdat(11+2*j,i) end do ! end if end do ! ! Load some material dependent quantities for Moliere scattering zsmcs = 0. zemcs = 0. ! 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 ! chiccmcs = 0.156915 * zsmcs * mparm(3,krad) / mparm(2,krad) ! return end !MS$page c ********************************************************* subroutine MAKE_BEAM c c Generate incident beam particles c include 'icommon.inc' c real*8 pm c ievt = jpar ipnum = 1 ipflg = 0 evtwt = 1. c c Assign particle mass randomly according to specified beam fractions rn = RAN1(idum) xlim = 0. 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 c -------------------------------------------------------- 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 c c Correlated gaussian ellipse ? ! do i=1,nbcorr(kk) if( corrtyp(i,kk) .eq. 1 ) then ix = ibcx(i,kk) ip = ibcp(i,kk) c = corr1(i,kk) if( ip .lt. 3 ) then th = pp(ip) / pp(3) else th = (pp(3)-p1bt(3,kk)) / p2bt(3,kk) end if xp(ix) = xp(ix) + th * c * x2bt(ix,kk) ! else ! 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 end if end do ! ........................................................ else if( bdistyp(kk) .eq. 2 ) then ! uniform circle r = RANV(x1bt(1,kk),x2bt(1,kk)) phi = RANV(x1bt(2,kk),x2bt(2,kk)) xp(3) = RANV(x1bt(3,kk),x2bt(3,kk)) pr = RANV(p1bt(1,kk),p2bt(1,kk)) pphi = RANV(p1bt(2,kk),p2bt(2,kk)) pp(3) = RANV(p1bt(3,kk),p2bt(3,kk)) ! cph = cos(phi) sph = sin(phi) xp(1) = r * cph xp(2) = r * sph pp(1) = pr*cph - pphi*sph pp(2) = pr*sph + pphi*cph end if c --------------------------------------------------------- ! c Convert initial bunch length into time spread call VMAG(pp,pm) mp = mass(iptyp) gamma = sqrt(1. + pm*pm/mp/mp) vz = pp(3) / mp / gamma * cc tp = xp(3) / vz xp(3) = 0. c if( spin ) then ! assume spin parallel to momentum do i=1,3 pol(i) = pp(i) / pm end do end if ! return end !MS$page c ********************************************************* subroutine MSTEP(pm) ! ! Take a step thru material medium ! include 'icommon.inc' ! real*8 pm,pm0,pavg logical in_wedge,in_pwedge,in_aspwedge ! pm0 = pm ! original magnitude of momentum cth = pp(3) / pm0 ! cosine of particle's angle ! hdid is the spatial step taken along z ! hstp is the spatial step taken along the particle's trajectory hstp = hdid / cth ! ! Normal absorber if( mgt .eq. 'CBLOCK') then ! 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 ! Use average p from energy loss step to get scattering ANGLE ! Use pm to scale momentum vector after scattering pavg = 0.5 * (pm + pm0) ! if( lscatter ) call SCATTER(hstp,pavg,pm) ! Multiple scattering ! end if ! ! Wedge absorber ! ! This is similar to block absober above, except that the stepsize ! is adjusted for part of step that is inside the wedge ! if( (mgt.eq.'WEDGE' .and. IN_WEDGE(hstp)) & .or. (mgt.eq.'PWEDGE' .and. IN_PWEDGE(hstp)) & .or. (mgt.eq.'ASPW' .and. IN_ASPWEDGE(hstp)) ) then ! if( ldedx ) call DEDX(hstp,pm) if( iflag .ne. 0 ) go to 900 ! pavg = 0.5 * (pm + pm0) if( lscatter ) call SCATTER(hstp,pavg,pm) if( iflag .ne. 0 ) go to 900 ! end if ! 900 continue return end !MS$page c ********************************************************* subroutine NEW_PARTICLE(ip,ifl) ! ! Set up for new particle ! include 'icommon.inc' ! ifl = 0 if( ip .eq. 0 ) then ! initialize REFERENCE particle ievt = 0 ipnum = 0 iptyp = 2 ipflg = 0 do i=1,3 xp(i) = 0. pp(i) = 0. end do mp = mass(iptyp) iptypref = 2 evtwt = 1. ! if( bgen ) then ! use CONTROL values pp(3) = pz0ref tp = t0ref else if( tbref .eq. 0. ) then ! use CONTROL values pp(3) = pz0ref tp = t0ref else ! use FOR003 values do i=1,3 xp(i) = xbref(i) pp(i) = pbref(i) end do tp = tbref end if end if ! ppref(3) = pp(3) ! define these for production plots tpref = tp ! and bunch cut ! else ! initialize NORMAL particles ! if( bgen ) then ! generate new particle call MAKE_BEAM ! else ! read in an old one read(3,*,iostat=ioc)ievt,ipnum,iptyp,ipflg,tp,evtwt,xp,pp if( ioc .ne. 0 ) then iflag = -49 go to 900 end if if( goodtrack .and. ipflg.ne.0 ) then ifl = -1 go to 900 end if mp = mass(iptyp) end if ! end if ! end of test on ip=0 ! if( ip .le. ABS(nprnt) ) then write(*,105)ip,iptyp,xp,pp,tp write(2,105)ip,iptyp,xp,pp,tp 105 format(' INITIAL:',2i4,/,7e11.4) end if ! do i=1,3 ! save initial particle properties xpint(i) = xp(i) ppint(i) = pp(i) polinit(i) = pol(i) end do tpint = tp ! 900 continue return end !MS$page c ********************************************************* subroutine OUTPUT ! ! Write particle information at this location to "n-tuple" file ! include 'icommon.inc' ! write(16,100)ievt,ipnum,iptyp,ipflg,jsrg,tp,xp,pp,bfld,evtwt, & efld,phaserf 100 format(1x,i6,2i3,2i4,15e14.5) ! return end !MS$page ! ********************************************************* subroutine PKICK(e,b) ! ! Kick particle momentum due to EM fields ! include 'icommon.inc' ! real*8 pm real*4 e(3),b(3),bet(3),dp(3),mmu ! mmu = mass(2) call VMAG(pp,pm) gamma = SQRT( 1. + (pm/mmu)**2 ) do i=1,3 bet(i) = pp(i) / mmu / gamma end do ! term = slen / 2. / bet(3) dp(1) = ee * ( e(1)/cc + bet(2)*b(3) - bet(3)*b(2) ) * term dp(2) = ee * ( e(2)/cc - bet(1)*b(3) + bet(3)*b(1) ) * term dp(3) = ee * ( 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 ! include 'icommon.inc' ! character*10 fname ! ........................................................ if( jsec .eq. -1 ) then ! rotation plane read(7,rec=i7) jsec,jsrg,ang ! call ROTATE_COORD(ang) ! if( ip .le. abs(nprnt) ) then write(2,75)ss0*180./pi,xp,pp 75 format(' ROT=',f7.2,6e11.4) end if ! ! ........................................................ ! Save output plane information else if( jsec.eq.-2 .and. .not.(ntuple.or.rtuple) ) then ! call OUTPUT ! ! ........................................................ else if( jsec .eq. -3 ) then ! aperture plane read(7,rec=i7) jsec,jsrg,iapertype,aperlim1,aperlim2 ! if( iapertype .eq. 1 ) then ! elliptical cut if( xp(1) .ge. aperlim1 ) then iflag = -28 go to 900 end if yell = aperlim2 * SQRT(1. - (xp(1)/aperlim1)**2) if( xp(2) .ge. yell ) then iflag = -28 go to 900 end if ! else ! rectangular cut if( ABS(xp(1)).ge.aperlim1 .or. & ABS(xp(2)).ge.aperlim2 ) then iflag = -28 go to 900 end if end if ! ! ........................................................ ! Background field needed else if( jsec.eq.-4 .and. & (jpar.eq.0 .or. (jpar.eq.1 .and. .not.refpar)) ) then read(7,rec=i7) jsec,jsrg,gridbkg,prefbkg,ztotalbkg,interbkg, & xlobkg,dxbkg,nxbkg,ylobkg,dybkg,nybkg,zlobkg,dzbkg,nzbkg if( gridbkg) then write(*,*)' Resetting background field grid' call BACKSET(0) ! reset background field map end if sbg0 = saccum ! ! ........................................................ ! Background field region else if( jsec.eq.-5 .and. & (jpar.eq.0 .or. (jpar.eq.1 .and. .not.refpar)) ) then read(7,rec=i7) jsec,jsrg,zoffbkg,rmaxbkg,zcutbkg,bftag,bfparm write(*,*)' Computing background field: ',bftag ! if( gridbkg) then call BACKSET(1) ! add field to background if( iflag .ne. 0 ) then write(2,*)' Error in computing background field:',bftag write(*,*)' Error in computing background field:',bftag go to 900 end if end if ! ! ........................................................ ! End background field region else if( jsec.eq.-6 .and. & (jpar.eq.0 .or. (jpar.eq.1 .and. .not.refpar)) ) then read(7,rec=i7) jsec,jsrg,mfile rechibkg = i7 - 1 ! rec num of last BFIELD reclobkg = i7 - 1 ! do i=i7-2,0,-1 read(7,rec=i) idum if( idum .eq. -4 ) go to 100 ! look for BACKGROUND reclobkg = i end do ! 100 continue if( mfile.gt.10 .and. gridbkg ) then ! Save field grid to file write(fname,120) mfile 120 format('for0',i2,'.dat') open(unit=mfile,file=fname,status='unknown',iostat=ioc) if( ioc .ne. 0 ) then iflag = -50 go to 900 end if write(mfile,'(i8)') nxbkg*nybkg*nzbkg do k=1,nzbkg do j=1,nybkg do i=1,nxbkg write(mfile,'(3i4,6e11.3)') i,j,k,xbkg(i),ybkg(j), & zbkg(k),bxbkg(i,j,k),bybkg(i,j,k),bzbkg(i,j,k) end do end do end do close(mfile) end if ! ........................................................ ! User transport matrix else if( jsec.eq.-7 ) then read(7,rec=i7) jsec,jsrg,slentrp,((trprt(i,j),j=1,6),i=1,6) firsttrp = .true. call TRANSPORT ! end if ! end of test on jsec ! ! ........................................................ 900 continue return end !MS$page c ********************************************************* function RANV(xlo,xhi) ! ! Return uniform distribution of variable over specified range ! rn = RAN1(idum) RANV = (xhi-xlo)*rn + xlo ! return end !MS$page ! ********************************************************* function RAN1(idum) ! ! Return uniformly distributed random variable in (0,1) ! ! Adapted from Numerical Recipes, 2nd ed., p. 271 ! include 'icommon.inc' ! INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV REAL ran1,AM,EPSs,RNMX PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836, *NTAB=32,NDIV=1+(IM-1)/NTAB,EPSs=1.2e-7,RNMX=1.-EPSs) INTEGER j,k,iv(NTAB),iy SAVE iv,iy DATA iv /NTAB*0/, iy /0/ ! if (idum.le.0.or.iy.eq.0) then idum = rnseed ! take seed from common/CONTROL/ idum=max(-idum,1) do 11 j=NTAB+8,1,-1 k=idum/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM if (j.le.NTAB) iv(j)=idum 11 continue iy=iv(1) endif ! k=idum/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum=idum+IM j=1+iy/NDIV iy=iv(j) iv(j)=idum ran1=min(AM*iy,RNMX) ! return END !MS$page ! ********************************************************* subroutine READ_EVT(ip) ! ! Read particle data from common arrays or unit(8) ! 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) end do tp = tpref ievt = 0 iptyp = iptypref ipflg = ipflgref evtwt = 1. ! 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) bzfld = bzfldpar(ip) ievt = ievtpar(ip) ipnum = ipnumpar(ip) iptyp = iptypar(ip) ipflg = ipflgpar(ip) evtwt = evtwtpar(ip) 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,bzfld, & xpint,ppint,tpint,polinit end if ! return end !MS$page c ********************************************************* subroutine REG_SUMMARY(nrlines) ! ! Prepare ordered list of region data ! Print summary of regions ! include 'icommon.inc' ! real*4 mmu logical*1 cell1,rep1,repincell character*80 probtitle character*4 cmd ! mmu = mass(2) write(2,*) write(2,358) 358 format(t1,20('*'),t30,'REGION SUMMARY',t58,20('*')) write(2,360) 360 format(t4,'I',t6,'SRG',t10,'SEC',t18,'SLEN',t29,'ZLO', & t41,'ZHI',t50,'IR',t53,'FTAG',t59,'MTAG',t66,'MGEOM') ! jsrg = 0 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 = 'xxxx' ! shift = 0. if( .not.bgen ) then ! determine shift from beam restart read(3,'(a80)') probtitle read(3,*) xbref,pbref,tbref read(3,*,iostat=ioc)ievt,ipnum,iptyp,ipflg,tp,evtwt,xp,pp if( ioc .ne. 0 ) then iflag = -49 end if shift = xp(3) rewind(3) end if zlo = shift ! -------------------------------------------------------- do 450 is=1,nsections jsec = is ss0 = zlo ! FIX THIS IF WE USE IT c ! rewind input to top of region data in the input file do i=1,nrlines 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,'(a4)',iostat=ioc) cftag read(1,*) cellflip 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. write(2,'(t20,a,i5,5x,a4)')' CELLS = ',ncells,cftag ! ! ........................................................ 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 write(2,'(i4,t20,a)')i7,'PRODUCTION' write(7,rec=i7) 0 ! ! ........................................................ else if( cmd .eq. 'ROTA' )then read(1,*,iostat=ioc) rotang i7 = i7 + 1 write(2,'(i4,t20,a,f12.4)')i7,'ROTATION =',rotang write(7,rec=i7) -1,jsrg,rotang*pi/180. if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 2 end if if( rep1 ) nreplines = nreplines + 2 ! ! ........................................................ else if( cmd .eq. 'OUTP' )then i7 = i7 + 1 write(2,'(i4,t20,a)')i7,'OUTPUT' write(7,rec=i7) -2,0,jsrg 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 i7 = i7 + 1 write(2,'(i4,t20,a)')i7,'APERTURE' write(7,rec=i7) -3,jsrg,iapertype,aperlim1,aperlim2 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) gridbkg,prefbkg,ztotalbkg read(1,*,iostat=ioc) xlobkg,dxbkg,nxbkg, & ylobkg,dybkg,nybkg,zlobkg,dzbkg,nzbkg,interbkg i7 = i7 + 1 write(2,'(i4,t20,a)') i7,'BACKGROUND FIELD' write(7,rec=i7) -4,jsrg,gridbkg,prefbkg,ztotalbkg,interbkg, & xlobkg,dxbkg,nxbkg,ylobkg,dybkg,nybkg,zlobkg,dzbkg,nzbkg ! ! ........................................................ else if( cmd .eq. 'BFIE' )then ! background field read(1,*,iostat=ioc) zoffbkg,rmaxbkg,zcutbkg read(1,'(a4)',iostat=ioc) bftag read(1,*,iostat=ioc)(bfparm(j),j=1,nfparm) i7 = i7 + 1 write(2,'(i4,t25,a)') i7,bftag write(7,rec=i7) -5,jsrg,zoffbkg,rmaxbkg,zcutbkg,bftag,bfparm ! ! ........................................................ else if( cmd .eq. 'ENDB' )then read(1,*,iostat=ioc) mfile i7 = i7 + 1 write(7,rec=i7) -6,jsrg,mfile write(2,'(i4,t20,a)') i7,'END BACKGROUND FIELD' ! ! ........................................................ else if( cmd .eq. 'TRAN' )then read(1,*,iostat=ioc) slentrp read(1,*,iostat=ioc) ((trprt(i,j),j=1,6),i=1,6) i7 = i7 + 1 jsrg = jsrg + 1 write(7,rec=i7) -7,jsrg,slentrp,((trprt(i,j),j=1,6),i=1,6) zhi = zlo + slentrp write(2,402) i7,jsrg,slentrp,zlo,zhi,'USER TRANSPORT' 402 format(i4,t5,i4,t14,e10.3,t25,e11.3,t37,e11.3,t50,a) zlo = zhi if( cell1.and.(.not.repincell.or.(repincell.and.rep1)))then nclines = nclines + 8 end if if( rep1 ) nreplines = nreplines + 8 ! ! ........................................................ else if( cmd .eq. 'SREG' )then read(1,*,iostat=ioc) slen,nrreg,zstep i7 = i7 + 1 jsrg = jsrg + 1 sr0 = zlo zhi = zlo + slen ! do 430 ir=1,nrreg read(1,*,iostat=ioc) idum,rlow(ir),rhigh(ir) read(1,'(a4)',iostat=ioc) ftag(ir) read(1,*,iostat=ioc)(fparm(j,ir),j=1,nfparm) read(1,'(a4)',iostat=ioc) mtag(ir) read(1,'(a6)',iostat=ioc) mgeom(ir) read(1,*,iostat=ioc) (gparm(j,ir),j=1,ngparm) write(2,405) i7,jsrg,is,slen,zlo,zhi,ir,ftag(ir), & mtag(ir),mgeom(ir) 405 format(i4,t5,i4,t9,i4,t14,e10.4,t25,e11.4,t37,e11.4,t49,i3, + t53,a4,t59,a4,t66,a6) ! 430 continue ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! write(7,rec=i7) jsec,jcel,jsrg,kcel,ss0,sc0,sr0,cellslp, & cftag,cfparm,slen,nrreg,zstep, & rlow,rhigh,ftag,fparm,mtag,mgeom,gparm ! ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ zlo = zhi c 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 write(2,*) c else ! done until start of any new cells nclines = 0 kcel = 0 cftag = 'NONE' ! clength = 0. end if ! end of test on ncells ! ! ........................................................ else if( cmd .eq. 'ENDS' )then 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 ! return end !MS$page 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 ! include 'icommon.inc' c parameter (nmax = 10) real*8 y(n),dydt(n),tt,th,h,yout(n), + yt(nmax),dyt(nmax),dym(nmax) c 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] h2 = h / 2. h6 = h / 6. 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) 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) 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) c do i=1,n ! variables at end yout(i) = y(i) + h6 * ( dydt(i)+2.*dym(i)+dyt(i) ) end do c return end !MS$page c ********************************************************* subroutine ROTATE_COORD(ang) ! ! Rotate particle coordinates by angle ANG around s-axis ! include 'icommon.inc' ! real*8 x(3),p(3) ! ca = COS(ang) sa = SIN(ang) do i=1,3 x(i) = xp(i) p(i) = pp(i) end do ! 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 ! return end !MS$page c ********************************************************* subroutine SAVE_FILE(izp) c c Save particle information at plane izp ! ( or on error condition ) c 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,'(4i6,2e14.6)') ievt,ipnum,iptyp,ipflg,tt,evtwt write(4,'(6e14.6)') xp(1),xp(2),zz,pp ! else ! want initial conditions (also for errors) write(4,'(4i6,2e14.6)') ievt,ipnum,iptyp,ipflg,tpint,1. write(4,'(6e14.6)') xpint,ppint end if end if ! return end !MS$page 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 ! 6 : "n-tuple" file ! 7 : region data file (internal) ! 8 : particle info (internal) ! * : terminal output ! ## : {coil,sheet} segment data; ! : {coil,sheet,background} grid data; ! : rf diagnostics;rf phases ! include 'icommon.inc' c integer*2 iyr,imon,ihr,imin,isec logical*4 cloop,rloop,bgloop character*4 cmd,ftagdat(16) character*4 cloopdat(8),rloopdat(6),bgloopdat(2) character*6 mgeomdat(5) character*10 fname character*79 probtitle ! namelist/BMT/nbeamtyp namelist/CONT/betaperp,bgen,bunchcut,eps,fsav,fsavset, & goodtrack,gradref,izfile, & neighbor,npart,nprnt,nsections,ntuple, & phasemodel,prlevel,pz0ref,refpar,rfdiag,rfphase, + rnseed,rtuple,rtuplen,spin,timelim,t0ref,varstep namelist/INTS/ldedx,lscatter,lbrem,linteract,ldecay,lspace, + delev,scatlev,bremlev,intlev,declev,spacelev namelist/NEM/nemit,sigmacut,pxycorr,pzcorr,sig_cut namelist/NHS/nhist,hauto,hcprn namelist/NSC/nscat,sauto namelist/NZH/nzhist,zhauto,zhprin c data nftags/16/ data ftagdat/'NONE','ACCE','BACK','BSOL','COIL','DEFL','DIP ', & 'FOFO','QUAD','ROD ','SEX ','SHEE','SOL ','STUS', & 'TROD','TSOL'/ data nmgeoms/5/ data mgeomdat/'CBLOCK','WEDGE','PWEDGE','ASPW', & 'NONE '/ data ncloop/8/ data cloopdat/'SREG','REPE','ENDR','OUTP','APER', & 'ROTA','TRAN','ENDC'/ data nrloop/6/ data rloopdat/'SREG','OUTP','APER','ROTA','TRAN','ENDR'/ data nbgloop/2/ data bgloopdat/'BFIE','ENDB'/ data cloop/.false./,rloop/.false./,bgloop/.false./ ! c Open required files ! Input data open(unit=1,file='for001.dat',status='unknown',iostat=ioc) if( ioc .ne. 0 ) go to 800 ! ! Output log open(unit=2,file='for002.dat',status='unknown',iostat=ioc) if( ioc .ne. 0 ) go to 800 ! ! Region data open(unit=7,file='for007.dat',access='direct', & form='unformatted',recl=596,status='unknown',iostat=ioc) if( ioc .ne. 0 ) go to 800 c c File headers write(*,'(a,f6.2)')' Welcome to ICOOL - version ',version write(2,'(a,f6.2)')' Welcome to ICOOL - version ',version call ICOOLDATE(iyr,imon,iday0) call ICOOLTIME(ihr0,imin0,isec0) write(2,'(a,4i6)')' Date: ',iyr,imon,iday0 write(*,'(a,4i6)')' Start time: ',ihr0,imin0,isec0 write(2,'(a,4i6)')' Start time: ',ihr0,imin0,isec0 c read(1,'(a79)',iostat=ioc) probtitle if( ioc .ne. 0 ) go to 810 write(2,'(a79)') probtitle ! ! Output file header open(unit=16,file='for016.dat',status='unknown',iostat=ioc) if( ioc .ne. 0 ) go to 800 write(16,'(a1,a79)') '#',probtitle write(16,'(a)')'# units = [s] [m] [GeV/c] [T]' write(16,'(a,a)')'>event partnum parttype flag region', & ' time x y z px py pz Bx By Bz weight Ex Ey Ez phase' ! c Read simulation control parameters read(1,cont,iostat=ioc) if( ioc .ne. 0 ) go to 810 write(2,cont) c if( .not. bgen ) then ! existing initial particle values open(unit=3,file='for003.dat',status='unknown',iostat=ioc) if( ioc .ne. 0 ) go to 800 end if ! if( fsav ) then ! new particle data at specified plane open(unit=4,file='for004.dat',status='unknown',iostat=ioc) if( ioc .ne. 0 ) go to 800 write(4,'(a79)') probtitle end if ! if( npart .gt. mxpar ) then ! particle info (internal) open(unit=8,file='for008.dat',access='direct', & form='unformatted',recl=200,status='unknown',iostat=ioc) if( ioc .ne. 0 ) go to 800 end if ! ! Rf diagnostics if( rfdiag.gt.10 .and. rfdiag.lt.100 ) then write(fname,40) rfdiag 40 format('for0',i2,'.dat') open(unit=rfdiag,file=fname,status='unknown',iostat=ioc) if( ioc .ne. 0 ) go to 800 write(rfdiag,'(a,a)') 'IEVT JSRG XP PP TP PZREF TPREF ', & ' PHASE EZ BPHI DPZ DT DPHASE' end if ! ! Rf phases if( rfphase.gt.10 .and. rfphase.lt.100 ) then write(fname,40) rfphase open(unit=rfphase,file=fname,status='unknown',iostat=ioc) if( ioc .ne. 0 ) go to 800 end if if( rtuple ) ntuple = .false. ! c Read beam distribution parameters if( bgen ) then read(1,bmt,iostat=ioc) if( ioc .ne. 0 ) go to 810 write(2,bmt) c if( nbeamtyp .gt. 0 ) then do i=1,nbeamtyp read(1,*,iostat=ioc)idum,bmtype(i),fracbt(i),bdistyp(i) if( ioc .ne. 0 ) go to 810 write(2,50)idum,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 ) go to 810 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 ) go to 810 write(2,56)(x2bt(j,i),j=1,3),(p2bt(j,i),j=1,3) 56 format(' Xsg,Psg: ',6e11.4) ! read(1,*,iostat=ioc) nbcorr(i) if( ioc .ne. 0 ) go to 810 write(2,57) nbcorr(i) 57 format(' NBCORR = ',i5) ! do j=1,nbcorr(i) read(1,*,iostat=ioc) corrtyp(j,i),ibcx(j,i),ibcp(j,i), & corr1(j,i),corr2(j,i) if( ioc .ne. 0 ) go to 810 write(2,58) j,corrtyp(j,i),ibcx(j,i),ibcp(j,i), & corr1(j,i),corr2(j,i) 58 format(' I,Ctyp,IBCX,IBCP,CORR1,CORR2: ',4i5,2f12.4) end do ! 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 ) go to 810 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 ) go to 810 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 go to 810 end if ! bdistyp end do end if ! nbeamtyp > 0 end if ! bgen c c Read physics interactions control parameters read(1,ints,iostat=ioc) if( ioc .ne. 0 ) go to 810 write(2,ints) c c Set up quantities to histogram read(1,nhs,iostat=ioc) if( ioc .ne. 0 ) go to 810 write(2,nhs) if(nhist .gt. 20) go to 820 nhword = 1 c do 260 i=1,nhist read(1,*,iostat=ioc)hxmin(i),hdx(i),nhbins(i),ihvar(i),ihdes(i) if( ioc .ne. 0 ) go to 810 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 c c Set up quantities to scatter plot read(1,nsc,iostat=ioc) if( ioc .ne. 0 ) go to 810 write(2,nsc) if(nscat .gt. 20) go to 820 nsword = 1 c 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(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 c c Set up Z-HISTORY plots read(1,nzh,iostat=ioc) if( ioc .ne. 0 ) go to 810 write(2,nzh) if(nzhist .gt. 10) go to 820 c do i=1,nzhist read(1,*,iostat=ioc)nzhpar(i),zhxmin(i),zhdx(i),nzhxbin(i), + zhymin(i),zhdy(i),nzhybin(i),izhyvar(i) if( ioc .ne. 0 ) go to 810 if( ABS(nzhpar(i)) .gt. 12 ) nzhpar(i) = 12 if(nzhxbin(i) .gt. 70) nzhxbin(i) = 70 if(nzhybin(i) .gt. 23) nzhybin(i) = 23 write(2,*)nzhpar(i),zhxmin(i),zhdx(i),nzhxbin(i) write(2,*)zhymin(i),zhdy(i),nzhybin(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,12 zhxcur(j,i) = zhxmin(i) end do end do ! c Set up EMITTANCE calculations read(1,nem,iostat=ioc) if( ioc .ne. 0 ) go to 810 write(2,nem) if(nemit .gt. 30) go to 820 c read(1,*,iostat=ioc) (iemdes(j),j=1,nemit) if( ioc .ne. 0 ) go to 810 write(2,*) (iemdes(j),j=1,nemit) ! -------------------------------------------------------- c c Check validity of the region data c nrlines = 0 ! count number of lines of region data c 300 continue ! loop over region commands c read(1,'(a4)',iostat=ioc)cmd if( ioc .ne. 0 ) go to 810 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 go to 820 end if ! 303 continue if( rloop ) then do i=1,nrloop if( cmd .eq. rloopdat(i) )go to 306 end do go to 820 end if ! 306 continue if( bgloop ) then do i=1,nbgloop if( cmd .eq. bgloopdat(i) )go to 309 end do go to 820 end if ! 309 continue c ........................................................ if( cmd .eq. 'SECT' )then go to 300 c ! ........................................................ else if( cmd .eq. 'CELL' )then read(1,*,iostat=ioc) ncells if( ioc .ne. 0 ) go to 810 write(2,*) ncells if( ncells .le. 0 ) go to 820 read(1,'(a4)',iostat=ioc) cftag if( ioc .ne. 0 ) go to 810 write(2,*) cftag do i=1,nftags if( cftag .eq. ftagdat(i) ) go to 315 end do go to 820 c 315 continue read(1,*,iostat=ioc) cellflip if( ioc .ne. 0 ) go to 810 write(2,*) cellflip read(1,*,iostat=ioc)(cfparm(j),j=1,nfparm) if( ioc .ne. 0 ) go to 810 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 ) go to 810 write(2,*) slen,nrreg,zstep if( slen .le. 0. ) go to 820 if( nrreg .le. 0 .or. nrreg .gt. 4 ) go to 820 if( zstep.le.0. .or. zstep.gt.slen ) go to 820 c do 330 ir=1,nrreg read(1,*,iostat=ioc) idum,rlow(ir),rhigh(ir) if( ioc .ne. 0 ) go to 810 write(2,*) idum,rlow(ir),rhigh(ir) if( rlow(ir) .lt. 0. ) go to 820 if( rhigh(ir) .lt. rlow(ir) ) go to 820 read(1,'(a4)',iostat=ioc) ftag(ir) if( ioc .ne. 0 ) go to 810 write(2,*) ftag(ir) do i=1,nftags if( ftag(ir) .eq. ftagdat(i) ) go to 320 end do go to 820 c 320 continue ! good ftag read(1,*,iostat=ioc)(fparm(j,ir),j=1,nfparm) if( ioc .ne. 0 ) go to 810 write(2,*) (fparm(j,ir),j=1,nfparm) ! read(1,'(a4)',iostat=ioc) mtag(ir) if( ioc .ne. 0 ) go to 810 write(2,*) mtag(ir) ! if( mtag(ir) .eq. 'VAC' ) go to 322 do i=1,nmat if( mtag(ir) .eq. matid(i) ) go to 322 end do go to 820 ! 322 continue ! good mtag read(1,'(a6)',iostat=ioc) mgeom(ir) if( ioc .ne. 0 ) go to 810 write(2,*) mgeom(ir) ! do i=1,nmgeoms if( mgeom(ir) .eq. mgeomdat(i) ) go to 324 end do go to 820 ! 324 continue ! good mgeom read(1,*,iostat=ioc) (gparm(j,ir),j=1,ngparm) if( ioc .ne. 0 ) go to 810 write(2,*) (gparm(j,ir),j=1,ngparm) ! 330 continue ! nrlines = nrlines + 6 * nrreg + 1 c ! ........................................................ else if( cmd .eq. 'REPE' )then read(1,*,iostat=ioc) nrep if( ioc .ne. 0 ) go to 810 write(2,*) nrep if( nrep .le. 0 ) go to 820 rloop = .true. nrlines = nrlines + 1 ! ! ........................................................ else if( cmd .eq. 'ROTA' )then read(1,*,iostat=ioc) rotang if( ioc .ne. 0 ) go to 810 write(2,*) rotang nrlines = nrlines + 1 ! ! ........................................................ else if( cmd .eq. 'APER' )then read(1,*,iostat=ioc) iapertype,aperlim1,aperlim2 if( ioc .ne. 0 ) go to 810 write(2,*) iapertype,aperlim1,aperlim2 nrlines = nrlines + 1 ! ! ........................................................ else if( cmd .eq. 'TRAN' )then read(1,*,iostat=ioc) slentrp if( ioc .ne. 0 ) go to 810 write(2,*) slentrp if( slentrp .le. 0. ) go to 820 do i=1,6 read(1,*,iostat=ioc) (trprt(i,j),j=1,6) if( ioc .ne. 0 ) go to 810 end do write(2,'(6f12.4)') ((trprt(i,j),j=1,6),i=1,6) if( .not. refpar ) go to 820 nrlines = nrlines + 7 ! ! ........................................................ else if( cmd .eq. 'OUTP' )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) gridbkg,prefbkg,ztotalbkg if( ioc .ne. 0 ) go to 810 write(2,*) gridbkg,prefbkg,ztotalbkg if( prefbkg .le. 0. ) go to 820 if( ztotalbkg .le. 0. ) go to 820 read(1,*,iostat=ioc) xlobkg,dxbkg,nxbkg, & ylobkg,dybkg,nybkg,zlobkg,dzbkg,nzbkg,interbkg if( ioc .ne. 0 ) go to 810 write(2,*) xlobkg,dxbkg,nxbkg, & ylobkg,dybkg,nybkg,zlobkg,dzbkg,nzbkg,interbkg if( gridbkg ) then if( nxbkg.le.0 .or. nxbkg.gt.mxxg ) go to 820 if( nybkg.le.0 .or. nybkg.gt.mxyg ) go to 820 if( nzbkg.le.0 .or. nzbkg.gt.mxzg ) go to 820 if( dxbkg .le. 0. ) go to 820 if( dybkg .le. 0. ) go to 820 if( dzbkg .le. 0. ) go to 820 if( interbkg.lt.0 .or. interbkg.gt.3 ) go to 820 end if bgloop = .true. nrlines = nrlines + 2 ! ! ........................................................ else if( cmd .eq. 'BFIE' )then ! background field read(1,*,iostat=ioc) zoffbkg,rmaxbkg,zcutbkg if( ioc .ne. 0 ) go to 810 write(2,*) zoffbkg,rmaxbkg,zcutbkg ! if( zcutbkg .le. 0. ) go to 820 read(1,'(a4)',iostat=ioc) bftag if( ioc .ne. 0 ) go to 810 write(2,*) bftag do i=1,nftags if( bftag .eq. ftagdat(i) ) go to 340 end do go to 820 c 340 continue read(1,*,iostat=ioc)(bfparm(j),j=1,nfparm) if( ioc .ne. 0 ) go to 810 write(2,*) (bfparm(j),j=1,nfparm) nrlines = nrlines + 3 ! ! ........................................................ else if( cmd .eq. 'ENDB' )then read(1,*,iostat=ioc) idum if( ioc .ne. 0 ) go to 810 write(2,*) idum bgloop = .false. nrlines = nrlines + 1 ! ! ........................................................ else if( cmd .eq. 'ENDS' )then go to 350 c ! ........................................................ else write(*,*) 'Bad region command in SETUP' write(2,*) 'Bad region command in SETUP' go to 820 end if c go to 300 ! get next region command c 350 continue ! ! -------------------------------------------------------- ! c Print summary of region information; set up region data file; ! Set up reference particle data ! call REG_SUMMARY(nrlines) c write(2,*) c ........................................................ c Initialize call CLEAR ! clear histograms, scatterplots, statistics call ICOOLTIME(ihr,imin,isec) write(*,'(a,4i6)')' Completed SETUP: ',ihr,imin,isec write(2,'(a,4i6)')' Completed SETUP: ',ihr,imin,isec go to 900 c 800 continue c Error in input write(*,*)' <<<<< ERROR OPENING FILE >>>>>' write(2,*)' <<<<< ERROR OPENING FILE >>>>>' iflag = -2 go to 900 c 810 continue c Error reading input write(*,*)' <<<<< ERROR READING INPUT DATA >>>>>' write(2,*)' <<<<< ERROR READING INPUT DATA >>>>>' iflag = -2 go to 900 c 820 continue c Consistency check error write(*,*)' <<<<< CONSISTENCY CHECK ERROR >>>>>' write(2,*)' <<<<< CONSISTENCY CHECK ERROR >>>>>' iflag = -11 c 900 continue return end !MS$page ! ********************************************************* subroutine SIMULATE c c Execute simulation of the cooling sections c include 'icommon.inc' ! real*8 pm,stem(2,2) real*4 esc(3),bsc(3),ddsc(3),mmu logical*4 kstat,spcharge logical*1 lscatsav,vstepsav character*1 key character*4 tag character*80 probtitle ! data nfail/0/ ! lscatsav = lscatter ! save user control parameters delevsav = delev vstepsav = varstep ! if(.not.bgen ) then ! skip title card on beam input file read(3,'(a80)') probtitle read(3,*) xbref,pbref,tbref end if ! saccum = 0. ! initialize background field ztotalbkg = 0. bftag = 'NONE' ! ! -------------------------------------------------------- do 600 i7=1,n7records j7rec = i7 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! read(7,rec=i7) jsec,jcel,jsrg,kcel,ss0,sc0,sr0,cellslp, & cftag,cfparm,slen,nrreg,zstep, & rlow,rhigh,ftag,fparm,mtag,mgeom,gparm ! ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ if( jsec .gt. 0 ) write(*,*) & ' Processing region: ',jsrg,', good tracks: ',npart-nfail ! ! Don't allow variable step sizes in wedge regions mgt = mgeom(1) if( mgt.eq.'WEDGE' .or. mgt.eq.'PWEDGE' .or. & mgt.eq.'ASPW' ) then varstep = .false. 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 ! Initialize space-charge for this region spcharge = .false. tag = ftag(1) ! tag = mtag(1) if( lspace .and. & tag .eq. 'ACCE' ) then ! & .not.(tag.eq.'LI' .or. tag.eq.'BE' .or. tag.eq.'AL' ! & .or. tag.eq.'TI' .or. tag.eq.'FE' .or. ! & tag.eq.'CU') ) then spcharge = .true. ! do i=1,2 do j=1,2 stem(i,j) = 0. end do end do stev = 0. 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 ! Load particle data into arrays / file do 200 ip=0,npart ! loop over particles if( ip.eq.0 .and. .not.refpar ) go to 200 ! jpar = ip iflag = 0 ! call NEW_PARTICLE(ip,ifl) ! if( iflag .ne. 0 ) then ! eof reading FOR003 beam input npart = npart - 1 go to 200 end if ! if( ifl .ne. 0 ) go to 200 ! old, bad track ! ! Save production information if( ntuple .or. rtuple ) call OUTPUT ! if( ip .gt. 0 ) then ! skip reference particle call FILL_HIST(1) call FILL_SCAT(1) if( spin ) call FILL_POL(1) end if ! call WRITE_EVT(ip) ! 200 continue ! go to 580 ! end if ! end of test on i7 = 1 ! --------------------------------------------------------- if( spcharge ) then do 220 ip=1,npart ! accumulate time statistics call READ_EVT(ip) if( ipflg .ne. 0 ) go to 220 stem(1,1) = stem(1,1) + tp stev = stev + 1. 220 continue ! if( stev .gt. 0. ) tavg = stem(1,1) / stev stem(1,1) = 0. mmu = mass(2) ! ! Transform z-plane values into spatial bunch do 250 ip=1,npart call READ_EVT(ip) if( ipflg .ne. 0 ) go to 250 dt = tp - tavg call VMAG(pp,pm) gamma = SQRT( 1. + (pm/mmu)**2 ) term = cc * dt / mmu / gamma do i=1,3 ddsc(i) = pp(i) * term xp(i) = xp(i) + ddsc(i) call WRITE_EVT(ip) end do ! ! Accumulate length and radius statistics stem(1,1) = stem(1,1) + ddsc(3) stem(2,1) = stem(2,1) + ddsc(3) * ddsc(3) r = SQRT( xp(1)**2 + xp(2)**2 ) stem(1,2) = stem(1,2) + r stem(2,2) = stem(2,2) + r * r 250 continue ! ......................................................... if( stev.gt.0. ) then zavg = stem(1,1) / stev sigzb = SQRT( stem(2,1)/stev - zavg*zavg ) ravg = stem(1,2) / stev sigrb = SQRT( stem(2,2)/stev - ravg*ravg ) end if ! end if ! ! Propagate each particle thru the cooling region ! -------------------------------------------------------- do 500 ip=0,npart ! loop over particles ! check for reference particle if( ip.eq.0 .and. .not.refpar ) go to 500 ! jpar = ip iflag = 0 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 900 end if ! jp = MOD( ip,20 ) if( jsec.gt.0 .and. jp.eq.0 ) then write(*,*) ' Current particle: ',ip end if ! call READ_EVT(ip) if( ipflg .ne. 0 ) go to 500 ! if( spcharge .and. ip.gt.0 ) then call VMAG(pp,pm) gamma = SQRT( 1. + (pm/mmu)**2 ) loc = 1 ! space-charge kick at start ! ! WRITE(2,*) I7,IP,LOC call SPACE_CHARGE(loc,gamma,sigzb,sigrb,esc,bsc) call PKICK(esc,bsc) ! ! call VMAG(pp,pm) ! gamma = SQRT( 1. + (pm/mmu)**2 ) dt = tp - tavg term = cc * dt / mmu / gamma do i=1,3 xp(i) = xp(i) - pp(i)*term ! back to z-plane values end do end if ! ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! ! Watch for "pseudo-regions" if( jsec .lt. 0 )then ! pseudo-region ! call PSEUDOREGION(i7,ip) ! if( jsec .eq. -7 ) then ! transport go to 380 else go to 390 end if end if ! end of test on jsec.lt.0 ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! ! Physical region ! ! Get reference particle data for this region ! if( ip .eq. 0 ) then lscatter = .false. ! no scattering delev = 1 ! no straggling else lscatter = lscatsav delev = delevsav end if ! ........................................................ if( ip.eq.0 .and. phasemodel.eq.3 ) then ! assume constant ref particle velocity ! do this for every region for plotting purposes u = pz0ref / mp ! use constant reference momentum betaz = u / SQRT(1.+u*u) tnow = tp dt = slen / (betaz*cc) tp = tnow + dt trefmean = tnow + 0.5*dt xp(3) = xp(3) + slen ft = ftag(1) call FIELD(xp,tp) ! get field value for an OUTPUT go to 380 ! skip tracking end if ! if( ip.eq.0 .and. phasemodel.eq.4 .and. & ftag(1).eq.'ACCE' ) then ! assume constant reacceleration in cavities call VMAG(pp,pm) ! use actual particle momentum u = pm / mp gamma = SQRT( 1. + u*u ) betaz = pp(3) / gamma / mp pp(3) = pp(3) + (gradref/1.e3)*slen/betaz betazend = pp(3) / gamma / mp betaz = 0.5 * (betaz + betazend) tnow = tp dt = slen / (betaz*cc) tp = tnow + dt trefmean = tnow + 0.5*dt xp(3) = xp(3) + slen ft = ftag(1) call FIELD(xp,tp) ! get field value for an OUTPUT go to 380 ! skip tracking end if c ........................................................ c call GO_REGION ! take particle to end of this s region c c ........................................................ c 380 continue ! ! Print out current data if( ip .le. ABS(nprnt) ) then write(2,385)jsrg,ip,xp,pp,tp 385 format(' End of s-REGION:',i4,' , particle:',i4,/,7e11.4) end if ! ! 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.10 ) then dpzrf = pp(3) - ppref(3) dtrf = tp - tpref dphaserf = omrf * dtrf * 360. / twopi write(rfdiag,387) ievt,jsrg,xp,pp,tp,ppref(3),tpref, & phaserf,ezrf,bphirf,dpzrf,dtrf,dphaserf 387 format(1x,2i5,15e12.4) end if ! ........................................................ 390 continue ! Come in here from a pseudo-region ! 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 write(*,395)iflag,ievt,ipnum,ipflg,jsrg write(2,395)iflag,ievt,ipnum,ipflg,jsrg 395 format(' FAILURE; IFLAG,IEVT,IPNUM,IPFLG,JSRG: ',5i6) end if ! if( jsec.gt.0 .and. (ntuple .or. (rtuple.and.iflag.ne.0)))then call OUTPUT end if ! if( ip.eq.0 .and. fsav .and. i7.eq.izfile ) then if( fsavset ) then write(4,'(7e14.6)') 0.,0.,0., 0.,0.,0., 0. else write(4,'(7e14.6)') xpref,ppref,tpref 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 ! -------------------------------------------------------- if( spcharge ) then do i=1,2 do j=1,2 stem(i,j) = 0. end do end do stev = 0. ! do 520 ip=1,npart ! accumulate time statistics call READ_EVT(ip) if( ipflg .ne. 0 ) go to 520 stev = stev + 1. stem(1,1) = stem(1,1) + tp 520 continue ! if( stev .gt. 0. ) tavg = stem(1,1) / stev ! stem(1,1) = 0. ! do 540 ip=1,npart ! accumulate bunch statistics call READ_EVT(ip) if( ipflg .ne. 0 ) go to 540 ! Transform z-plane values into spatial bunch dt = tp - tavg call VMAG(pp,pm) gamma = SQRT( 1. + (pm/mmu)**2 ) term = cc * dt / mmu / gamma do i=1,3 ddsc(i) = pp(i) * term xp(i) = xp(i) + ddsc(i) end do call WRITE_EVT(ip) ! ! Accumulate length and radius statistics stem(1,1) = stem(1,1) + ddsc(3) stem(2,1) = stem(2,1) + ddsc(3) * ddsc(3) r = SQRT( xp(1)**2 + xp(2)**2 ) stem(1,2) = stem(1,2) + r stem(2,2) = stem(2,2) + r * r 540 continue ! if( stev .gt. 0. ) then zavg = stem(1,1) / stev sigzb = SQRT( stem(2,1)/stev - zavg*zavg ) ravg = stem(1,2) / stev sigrb = SQRT( stem(2,2)/stev - ravg*ravg ) end if ! do 560 ip=1,npart call READ_EVT(ip) if( ipflg .ne. 0 ) go to 560 call VMAG(pp,pm) gamma = SQRT( 1. + (pm/mmu)**2 ) loc = 2 ! space-charge kick at end ! WRITE(2,*) I7,IP,LOC call SPACE_CHARGE(loc,gamma,sigzb,sigrb,esc,bsc) call PKICK(esc,bsc) ! dt = tp - tavg term = cc * dt / mmu / gamma do i=1,3 xp(i) = xp(i) - pp(i)*term ! back to z-plane values end do call WRITE_EVT(ip) ! 560 continue end if ! end of test on spcharge ! ! Test if we have passed region of validity for background field if( jsec .gt. 0 ) saccum = saccum + slen diff = saccum - sbg0 - ztotalbkg if( ABS(diff) .le. 1.0e-6 ) then bftag = 'NONE' ! tell FIELD to ignore background ! gridbkg = .false. end if ! 580 continue ! ! Calculate emittance at end of this region ! including possible cuts on tails and Bz corrections ! if( jsec .ge. 0 ) call EMITTANCE ! 600 continue ! end of loop over records ! -------------------------------------------------------- ! 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 write(*,*)' Total number of failed particles = ',nfail write(2,*)' Total number of failed particles = ',nfail c return end !MS$page ! ********************************************************* subroutine TRANSPORT ! ! Transport particle using user defined transport matrix ! include 'icommon.inc' ! real*8 pmref real*4 u(6),v(6) save pzref,tref,vref ! ! Save starting values of reference particle variables if( firsttrp ) then call VMAG( ppref,pmref ) gamma = SQRT( 1. + (pmref/mp)**2 ) pzref = ppref(3) vref = pzref / mp / gamma * cc tref = tpref firsttrp = .false. end if ! ! Convert ICOOL variables into TRANSPORT form (x,x',y,y',dl,dp/p) v(1) = xp(1) v(2) = pp(1) / pp(3) v(3) = xp(2) v(4) = pp(2) / pp(3) v(5) = (tp - tref) * vref v(6) = (pp(3) - pzref) / pzref ! ! Transport particle to end of region: u = T v do i=1,6 u(i) = 0. do j=1,6 u(i) = u(i) + trprt(i,j) * v(j) end do end do ! ! Convert back to ICOOL variables xp(1) = u(1) xp(2) = u(3) xp(3) = xp(3) + slentrp tp = tpref + slentrp/vref + u(5)/vref pp(3) = pzref + u(6)*pzref pp(1) = u(2) * pp(3) pp(2) = u(4) * pp(3) ! return end !MS$page ! ********************************************************* subroutine TRILINEAR(u,v,w,jx,jy,jz,fg,mx,my,mz,f) ! ! Do tri-linear interpolation ! real*8 f real*4 fg(mx,my,mz) ! f = (1.-u)*(1.-v)*(1.-w) * fg(jx,jy,jz) & + (1.-u)*(1.-v)*w * fg(jx,jy,jz+1) & + (1.-u)*v*w * fg(jx,jy+1,jz+1) & + (1.-u)*v*(1.-w) * fg(jx,jy+1,jz) & + u*(1.-v)*(1.-w) * fg(jx+1,jy,jz) & + u*(1.-v)*w * fg(jx+1,jy,jz+1) & + u*v*w * fg(jx+1,jy+1,jz+1) & + u*v*(1.-w) * fg(jx+1,jy+1,jz) ! return end !MS$page ! ********************************************************* subroutine TRIPOLY(jx,jy,jz,xg,yg,zg,fg,mp,np,pp,kk, & x,y,z,f,df) ! ! Polynomial interpolation on 3-dimensional grid ! ! cf. Num. Rec., 2nd ed, p.103,118; uses POLINT ! ! kk: polynomial order {1,2,3} ! integer*4 pp real*8 f real*4 fg(mp,np,pp),xg(mp),yg(np),zg(pp) real*4 f1tmp(4),f2tmp(4),f3tmp(4) ! kx = MIN( MAX(jx-kk/2,1), mp-kk ) ky = MIN( MAX(jy-kk/2,1), np-kk ) kz = MIN( MAX(jz-kk/2,1), pp-kk ) ! do 300 k=1,kk+1 ! loop over z grid points do 200 j=1,kk+1 ! loop over y grid points ! ! Load temporary array with x grid values for fixed y and z do 100 i=1,kk+1 f1tmp(i) = fg(kx+i-1,ky+j-1,kz+k-1) 100 continue ! call POLINT(xg(kx),f1tmp,kk+1,x,f2tmp(j),df) ! 200 continue ! call POLINT(yg(ky),f2tmp,kk+1,y,f3tmp(k),df) ! 300 continue ! call POLINT(zg(kz),f3tmp,kk+1,z,f4,df) ! f = f4 return end !MS$page ! ********************************************************* subroutine VMAG(p,pm) c c Return magnitude of 3-vector p c real*8 p(3),pm c pm = sqrt( p(1)*p(1) + p(2)*p(2) + p(3)*p(3) ) c return end !MS$page ! ********************************************************* subroutine WRITE_EVT(ip) ! ! Write particle data to common arrays or unit(8) ! include 'icommon.inc' ! if( ip .le. mxpar ) then if( ip .eq. 0 ) then do i=1,3 xpref(i) = xp(i) ppref(i) = pp(i) end do tpref = tp iptypref = iptyp ipflgref = ipflg ! else ! 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