program EMITCALC c c Program EMITCALC, version 1.0 c dated 3/16/2004 c c calculates many beam parameters and amplitude cuts c main new feature: the 3 invariant emittances c note that the order in which they are listed may not c be as expected, especially if coupling is strong c c inputs control parameters from 'emitcalc.inp' c reads particle data file in ICOOL for009.dat format c calculations are performed on blocks of particles c having identical region # and Z c uses particle weights c most calculations assume axisymmetry c uses transverse position, canonical transverse momentum, energy and time c c options: c broad particle cuts (particle type, pzmin, pzmax) c several models for calculating vector potential in axisymmetric field c set up to be able to add your own analytic model for vector potential c 0: assume constant B0, taken from input c 1: use Bz at particle plus reference particle c if reference particle not available, use model 2 c 2: calculate using particle data, assuming axisymmetric c 3: algebraic expansion from input file c may remove linear dependence of transverse correlations on E, t c (though implemented in a very poor way) c may remove nonlinear dependence of E, t, on transverse amplitude c may cut tails of distribution c may overlap bunches based on RF period (2 options for doing this) c 0: no overlap c 1: center about reference particle, if available (else model 2) c 2: use weighting to find proper centered phase c c outputs some beam averages, plus: c number of particles within hard amplitude cuts c 3 invariant emittances, no axisymmetry assumption c note: it is possible for emittances to be out of order c transverse, longitudinal emittance c transverse and longitudinal beta, alpha parameters c x and y component of energy dispersion c nonlinear correlation parameters, E and t with transverse amplitude c sigma_E c asymmetry test (if axisymmetric, should be very small number) c c units: meters, GeV, GeV/c, Tesla c note given in GeV-m/c, multiply by 3.33 to get in eV-s c phase space area is (particle rest mass) * c * emittance c RF frequencies are in MHz c c note that e_6D=e_1*e_2*e_3 exactly, same as LBL definition except c for option to remove some nonlinear correlations c c ************************************************************ c IMPLICIT NONE c integer npmax parameter ( npmax = 100000 ) c common/CONST/cc,ee,pi,twopi,avagn,relcl,mass(5) real*8 cc,ee,pi,twopi,avagn,relcl,mass c common/CONTROL/sigma_cut,rffreq,transcut,longcut, & bzaxis,refcharge,refmass,reftime,refenergy,zloc, & pzmax,pzmin,curr_region,vp_model,rf_model,useptype, & ampl_corr,disp_corr real*8 sigma_cut,rffreq,transcut,longcut real*8 bzaxis,refcharge,refmass,reftime,refenergy,zloc real*8 pzmax,pzmin integer curr_region,vp_model,rf_model,useptype logical ampl_corr,disp_corr c common/REGION/xpar(3,npmax),ppar(3,npmax),tpar(npmax), & epar(npmax),bzpar(npmax),evtwtpar(npmax), & iptypar(npmax),usepar(npmax),npart real*8 xpar,ppar,tpar,epar,bzpar,evtwtpar integer iptypar,usepar,npart c integer i,ip,ioc real*8 Bp(3),Ep(3),phasep,polp(3),energy integer t_ievt,t_ipnum,t_iptyp,t_ipflg real*8 t_evtwt,t_xp(3),t_pp(3),count,t_tp,store_bz logical loopflag, flag, found_ref integer iregion, store_rf_model ! ! const data cc/2.9979246d8/,ee/0.29979246d0/,pi/3.1415926d0/, + avagn/6.0221367d23/,relcl/2.8179409d-13/, & twopi/6.2831853d0/ data mass/0.51099906d-3, 0.10565839d0, 0.13956995d0, 0.493677d0, + 0.93827231d0/ ! GeV/c^2 ! call SETUP ! store_bz = bzaxis store_rf_model = rf_model loopflag = .true. read(3,*,iostat=ioc)t_ievt,t_ipnum,t_iptyp,t_ipflg,iregion, & t_tp,t_xp,t_pp,Bp,t_evtwt,Ep,phasep,polp if (ioc .ne. 0) then loopflag = .false. end if do while (loopflag) ! loop over regions ! flag = .true. found_ref = .false. reftime = 0. refenergy = 0. if (useptype .eq. 0) then refcharge = ee refmass = mass(2) else refcharge = SIGN( ee, REAL(useptype) ) refmass = mass( ABS(useptype) ) end if bzaxis = store_bz rf_model = store_rf_model zloc = 0.0 ip = 0 curr_region = iregion if (.not. loopflag) flag = .false. do while (flag) ! loop over particles in region ! ip = ip + 1 do i=1,3 xpar(i,ip) = t_xp(i) ppar(i,ip) = t_pp(i) end do ! tpar(ip) = t_tp bzpar(ip) = Bp(3) evtwtpar(ip) = t_evtwt iptypar(ip) = t_iptyp if(t_ipflg .eq. 0) then usepar(ip) = 0 energy = DSQRT(t_pp(1)**2 + t_pp(2)**2 + & t_pp(3)**2 + mass(ABS(iptypar(ip)))**2) else usepar(ip) = -1 energy = 0. end if epar(ip) = energy if((useptype .ne. 0) .and. (t_iptyp .ne. useptype)) & usepar(ip) = -2 if((pzmax .gt. 0.0) .and. (t_pp(3) .gt. pzmax)) & usepar(ip) = -3 if(t_pp(3) .lt. pzmin) & usepar(ip) = -4 if(t_ievt .eq. 0) then found_ref = .true. usepar(ip) = -5 if (vp_model .ne. 0) then bzaxis = Bp(3) end if zloc = t_xp(3) reftime = t_tp refenergy = energy if (useptype .eq. 0) then refcharge = SIGN( ee, REAL(t_iptyp) ) refmass = mass( ABS(t_iptyp) ) end if end if if ((.not. found_ref) .and. (usepar(ip) .eq. 0)) then zloc = t_xp(3) reftime = t_tp if (useptype .eq. 0) then refcharge = SIGN( ee, REAL(t_iptyp) ) refmass = mass( ABS(t_iptyp) ) end if end if c read(3,*,iostat=ioc)t_ievt,t_ipnum,t_iptyp,t_ipflg,iregion, & t_tp,t_xp,t_pp,Bp,t_evtwt,Ep,phasep,polp if (ioc .ne. 0) loopflag = .false. if (iregion .ne. curr_region) flag = .false. if (.not. loopflag) flag = .false. ! end do ! end loop over particles npart = ip if ((.not. found_ref) .and. (rf_model .eq. 1)) rf_model = 2 c c Compute emittances and related quantities for this region call EMITTANCE c end do ! loop over regions c close(3) close(2) c end c c ********************************************************* subroutine EMITTANCE c IMPLICIT NONE c integer npmax,nloopmax real*8 min_conf,rf_scale parameter ( npmax = 100000 ) parameter ( nloopmax = 20 ) parameter ( min_conf = 0.3) parameter ( rf_scale = 1.E6) c common/CONST/cc,ee,pi,twopi,avagn,relcl,mass(5) real*8 cc,ee,pi,twopi,avagn,relcl,mass c common/CONTROL/sigma_cut,rffreq,transcut,longcut, & bzaxis,refcharge,refmass,reftime,refenergy,zloc, & pzmax,pzmin,curr_region,vp_model,rf_model,useptype, & ampl_corr,disp_corr real*8 sigma_cut,rffreq,transcut,longcut real*8 bzaxis,refcharge,refmass,reftime,refenergy,zloc real*8 pzmax,pzmin integer curr_region,vp_model,rf_model,useptype logical ampl_corr,disp_corr c common/REGION/xpar(3,npmax),ppar(3,npmax),tpar(npmax), & epar(npmax),bzpar(npmax),evtwtpar(npmax), & iptypar(npmax),usepar(npmax),npart real*8 xpar,ppar,tpar,epar,bzpar,evtwtpar integer iptypar,usepar,npart c integer nscut,i,j,ip,iptyp,knt real*8 phasespace(6,npmax),transverse(4,npmax) real*8 phaseavg(6),phasecorr(6,6),transavg(4),transcorr(4,4) real*8 vector3(3),vector4(4),longavg(3),longcorr(3,3) real*8 vpx,vpy,bzp,xp(3),pp(3),pcharge,Lcanon,avgpz real*8 ampl_trans(npmax),ampl_long(npmax) real*8 wtot,weight,save_wtot,disp_t,disp_E,avg_E real*8 e_perp,beta_perp,alpha_perp,Ldim,gamma1 real*8 e_long,beta_long,alpha_long,delta_long real*8 invar_emit(3),sigma_Aperp,avg_Aperp,corr_t,corr_E real*8 w_aperture,disp_x,disp_y,check_axisym real*8 store_corr_t,store_corr_E real*8 temp_t,temp_E,conf real*8 cut_perp,cut_long logical cut_flag,store_ampl_corr real*8 TIMESHIFT ! cut_flag = .true. ! test for iterations knt = 1 ! count iterations store_ampl_corr = ampl_corr ! do while (cut_flag) ! c initializations ampl_corr = store_ampl_corr if (vp_model .eq. 2) then call CALC_VP(bzaxis) end if if ((rf_model .gt. 0) .and. (rffreq .ne. 0.)) then call CALC_TC(rffreq*rf_scale,reftime,conf) if (conf .lt. min_conf) then write(*,*)' difficulty calculating bunch center' end if end if c c calculate averages, canonical momentum c shift time, energy by reference values c store new values in phasespace c do ip=1,npart if (usepar(ip) .eq. 0) then do i=1,3 xp(i) = xpar(i,ip) pp(i) = ppar(i,ip) end do bzp = bzpar(ip) ! call GET_VP(bzp,xp(1),xp(2),vpx,vpy) ! iptyp = iptypar(ip) pcharge = SIGN( ee, REAL(iptyp) ) c c ***** use canonical momentum, subtract reference time and energy ***** c phasespace(1,ip) = xp(1) phasespace(2,ip) = pp(1) + pcharge * vpx phasespace(3,ip) = xp(2) phasespace(4,ip) = pp(2) + pcharge * vpy if (rf_model .eq. 0) then phasespace(5,ip) = tpar(ip) - reftime else phasespace(5,ip) = & TIMESHIFT(tpar(ip)-reftime,rffreq*rf_scale) end if phasespace(6,ip) = epar(ip) - refenergy end if end do ! ! particle weights and some averages wtot = 0. avgpz = 0. Lcanon = 0. do ip=1,npart if (usepar(ip) .eq. 0) then weight = evtwtpar(ip) wtot = wtot + weight avgpz = avgpz + weight * ppar(3,ip) Lcanon = Lcanon + weight * & (phasespace(1,ip)*phasespace(4,ip)- & phasespace(3,ip)*phasespace(2,ip)) end if end do if (knt.eq.1) save_wtot=wtot ! total particle weight before cuts if (wtot>0.) then avgpz = avgpz / wtot Lcanon = Lcanon / wtot end if ! ! 6D correlations matrix do i=1,6 phaseavg(i)=0. do j=1,6 phasecorr(i,j)=0. end do end do do ip=1,npart if (usepar(ip) .eq. 0) then weight = evtwtpar(ip) do i=1,6 phaseavg(i)=phaseavg(i)+weight*phasespace(i,ip) do j=1,6 phasecorr(i,j)=phasecorr(i,j)+ & weight*phasespace(i,ip)*phasespace(j,ip) end do end do end if end do if (wtot>0.) then do i=1,6 phaseavg(i)=phaseavg(i)/wtot end do do i=1,6 do j=1,6 phasecorr(i,j)=phasecorr(i,j)/wtot - & phaseavg(i)*phaseavg(j) end do end do end if avg_E = phaseavg(6) + refenergy ! ! ***** for transverse parameters, ! ***** subtract out linear dispersion with energy, time do i=1,4 if (phasecorr(5,5) .gt. 0.) then disp_t = phasecorr(i,5)/phasecorr(5,5) else disp_t = 0. end if if (phasecorr(6,6) .gt. 0.) then disp_E = phasecorr(i,6)/phasecorr(6,6) else disp_E = 0. end if c save some parameters for output if (i .eq. 1) disp_x = disp_E * avg_E if (i .eq. 3) disp_y = disp_E * avg_E if (.not. disp_corr) then disp_t = 0. disp_E = 0. end if do ip=1,npart if (usepar(ip) .eq. 0) then transverse(i,ip)=phasespace(i,ip) & - disp_t*(phasespace(5,ip)-phaseavg(5)) & - disp_E*(phasespace(6,ip)-phaseavg(6)) end if end do end do ! ! averages, correlation matrix do i=1,4 transavg(i)=0. do j=1,4 transcorr(i,j)=0. end do end do do ip=1,npart if (usepar(ip) .eq. 0) then weight = evtwtpar(ip) do i=1,4 transavg(i)=transavg(i)+weight*transverse(i,ip) do j=1,4 transcorr(i,j)=transcorr(i,j)+ & weight*transverse(i,ip)*transverse(j,ip) end do end do end if end do if (wtot>0.) then do i=1,4 transavg(i)=transavg(i)/wtot end do do i=1,4 do j=1,4 transcorr(i,j)=transcorr(i,j)/wtot - & transavg(i)*transavg(j) end do end do end if ! ! calculate transverse emittance, beta, alpha, assuming axisymmetry call AXISYM(refmass,avgpz,transcorr,e_perp, & beta_perp,alpha_perp,Ldim,check_axisym) ! calculate single particle amplitudes if (beta_perp .gt. 0.) then gamma1 = (1.+alpha_perp**2+Ldim**2)/beta_perp else gamma1 = 0. end if do ip=1,npart if (usepar(ip) .eq. 0) then do i=1,4 vector4(i)=transverse(i,ip)-transavg(i) end do ampl_trans(ip) = (1/refmass)*((beta_perp/avgpz)* & (vector4(2)**2+vector4(4)**2) + & gamma1*avgpz*(vector4(1)**2+vector4(3)**2) + & 2.*alpha_perp*(vector4(1)*vector4(2) + & vector4(3)*vector4(4)) - 2.*Ldim* & (vector4(1)*vector4(4)-vector4(2)*vector4(3))) end if end do ! ! ***** now for longitudinal parameters, ! ***** subtract out nonlinear correlation with transverse amplitude do i=1,3 longavg(i)=0. do j=1,3 longcorr(i,j)=0. end do end do do ip=1,npart if (usepar(ip) .eq. 0) then weight = evtwtpar(ip) vector3(1) = phasespace(5,ip) vector3(2) = phasespace(6,ip) vector3(3) = ampl_trans(ip) do i=1,3 longavg(i)=longavg(i)+weight*vector3(i) do j=1,3 longcorr(i,j)=longcorr(i,j)+ & weight*vector3(i)*vector3(j) end do end do end if end do if (wtot>0.) then do i=1,3 longavg(i)=longavg(i)/wtot end do do i=1,3 do j=1,3 longcorr(i,j)=longcorr(i,j)/wtot - & longavg(i)*longavg(j) end do end do end if avg_Aperp = longavg(3) sigma_Aperp = longcorr(3,3) if (sigma_Aperp .le. 0.) ampl_corr = .false. ! if (sigma_Aperp .gt. 0.) then corr_t = longcorr(1,3)/sigma_Aperp corr_E = longcorr(2,3)/sigma_Aperp else corr_t = 0. corr_E = 0. end if c save correlation parameters for output store_corr_t = corr_t store_corr_E = corr_E if (.not. ampl_corr) then corr_t = 0. corr_E = 0. end if do ip=1,npart if (usepar(ip) .eq. 0) then phasespace(5,ip)=phasespace(5,ip) & - corr_t*(ampl_trans(ip)-longavg(3)) phasespace(6,ip)=phasespace(6,ip) & - corr_E*(ampl_trans(ip)-longavg(3)) end if end do ! ! create revised 6D correlation matrix after amplitude corr removed do i=1,6 phaseavg(i)=0. do j=1,6 phasecorr(i,j)=0. end do end do do ip=1,npart if (usepar(ip) .eq. 0) then weight = evtwtpar(ip) do i=1,6 phaseavg(i)=phaseavg(i)+weight*phasespace(i,ip) do j=1,6 phasecorr(i,j)=phasecorr(i,j)+ & weight*phasespace(i,ip)*phasespace(j,ip) end do end do end if end do if (wtot>0.) then do i=1,6 phaseavg(i)=phaseavg(i)/wtot end do do i=1,6 do j=1,6 phasecorr(i,j)=phasecorr(i,j)/wtot - & phaseavg(i)*phaseavg(j) end do end do end if ! ! calculate longitudinal emittance, beta, alpha e_long = phasecorr(5,5)*phasecorr(6,6)-phasecorr(5,6)**2 if (e_long .gt. 0.) then e_long = DSQRT(e_long) c normalize to c/(mc^2) c refmass is actually mc^2 in GeV e_long = e_long * cc / refmass delta_long = phasecorr(5,5) / (e_long * refmass / cc) alpha_long = phasecorr(5,6) / (e_long * refmass / cc) else if (e_long .lt. 0.) & write(*,*)' error in e_long= ',e_long, & ' region #',curr_region,' npart= ',npart delta_long = phasecorr(5,5) alpha_long = phasecorr(5,6) end if c define longitudinal beta roughly like transverse, sigma_z/sigma_pz beta_long = delta_long * cc * avgpz**3 / & (refmass**2 + avgpz**2) ! calculate single particle amplitudes do ip=1,npart if (usepar(ip) .eq. 0) then if ((e_long .le. 0.) .or. (delta_long .le. 0.)) then ampl_long(ip) = 0. else temp_t = phasespace(5,ip) - phaseavg(5) temp_E = phasespace(6,ip) - phaseavg(6) ampl_long(ip) = (cc/refmass) * & ( temp_t**2/delta_long + delta_long* & (temp_E-alpha_long*temp_t/delta_long)**2 ) end if end if end do c c calculate invariant emittances for modified correlation matrix c call EMIT3D(refmass,phasecorr,invar_emit) ! ! calculate number within amplitude cuts, and apply sigma cuts w_aperture = 0. nscut = 0 ! count # particles removed by sigma cuts this round do ip=1,npart if (usepar(ip) .eq. 0) then weight = evtwtpar(ip) if ((transcut .le. 0.) .or. & (ampl_trans(ip) .le. transcut)) then if ((longcut .le. 0.) .or. & (ampl_long(ip) .le. longcut)) then w_aperture = w_aperture + weight end if end if if (sigma_cut .gt. 0.) then cut_perp = 2*sigma_cut**2 * e_perp ! factor of two to make up for combining 4D at once cut_long = sigma_cut**2 * e_long if ((cut_perp .gt. 0.) .and. & (ampl_trans(ip) .gt. cut_perp)) then usepar(ip) = knt nscut = nscut + 1 else if ((cut_long .gt. 0.) .and. & (ampl_long(ip) .gt. cut_long)) then usepar(ip) = knt nscut = nscut + 1 end if end if end if end do ! ! check whether should repeat loop knt = knt + 1 if ((sigma_cut .le. 0.) .or. (knt .gt. nloopmax)) then cut_flag=.false. else ! Test if tail cutting process has converged, stop if nscut=0 if (nscut .eq. 0) cut_flag = .false. end if ! end do ! loop to test cut_flag ! ! OUTPUTS write(2,150)curr_region,zloc,bzaxis,avgpz,Lcanon, & transavg(1),transavg(3),phaseavg(5)+reftime,save_wtot, & w_aperture,invar_emit(1),invar_emit(2),invar_emit(3), & e_perp,e_long,beta_perp,alpha_perp,Ldim,beta_long,alpha_long, & disp_x,disp_y,store_corr_E,store_corr_t, & DSQRT(phasecorr(6,6)),check_axisym 150 format(i6,25e12.4) ! c return end c ********************************************************* subroutine LUDCMP(a,n,np,indx,d) c IMPLICIT NONE c ! ! Make LU decomposition of matrix ! INTEGER n,np,indx(n),NMAX ! --------------------------------------------------------- ! RCF modified REAL*8 d,a(np,np),TINY ! --------------------------------------------------------- PARAMETER (NMAX=500,TINY=1.0e-20) INTEGER i,imax,j,k ! --------------------------------------------------------- ! RCF modified REAL*8 aamax,dum,sum,vv(NMAX) ! --------------------------------------------------------- d=1. do 12 i=1,n aamax=0. do 11 j=1,n if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j)) 11 continue if (aamax.eq.0.) then ! singular matrix in ludcmp d = 0. go to 900 end if vv(i)=1./aamax 12 continue do 19 j=1,n do 14 i=1,j-1 sum=a(i,j) do 13 k=1,i-1 sum=sum-a(i,k)*a(k,j) 13 continue a(i,j)=sum 14 continue aamax=0. do 16 i=j,n sum=a(i,j) do 15 k=1,j-1 sum=sum-a(i,k)*a(k,j) 15 continue a(i,j)=sum dum=vv(i)*abs(sum) if (dum.ge.aamax) then imax=i aamax=dum endif 16 continue if (j.ne.imax)then do 17 k=1,n dum=a(imax,k) a(imax,k)=a(j,k) a(j,k)=dum 17 continue d=-d vv(imax)=vv(j) endif indx(j)=imax if(a(j,j).eq.0.)a(j,j)=TINY if(j.ne.n)then dum=1./a(j,j) do 18 i=j+1,n a(i,j)=a(i,j)*dum 18 continue endif 19 continue ! 900 continue return END c ********************************************************* subroutine AXISYM(refmass,avgpz,transcorr,e_perp, & beta_perp,alpha_perp,Ldim,check_axisym) c IMPLICIT NONE c real*8 refmass,avgpz,transcorr(4,4) real*8 e_perp,beta_perp,alpha_perp real*8 Ldim,check_axisym c integer i,j,indperp(4) real*8 matrix(4,4),denom,asymmetry do i=1,4 do j=1,4 matrix(i,j) = transcorr(i,j) end do end do ! ! get transverse emittance (square root of 4D emittance) call LUDCMP(matrix,4,4,indperp,e_perp) ! if (e_perp .ne. 0.) then ! else singular do i=1,4 e_perp = e_perp * matrix(i,i) end do end if ! if (e_perp .gt. 0.) then e_perp = DSQRT(e_perp) c take another square root to get a pseudo 2D emittance e_perp = DSQRT(e_perp) c normalize to 1/(mc) e_perp = e_perp / refmass beta_perp = (transcorr(1,1) + transcorr(3,3)) * avgpz & / (2. * refmass * e_perp) alpha_perp = - (transcorr(1,2) + transcorr(3,4)) & / (2. * refmass * e_perp) Ldim = (transcorr(1,4) - transcorr(3,2)) & / ( 2. * refmass * e_perp) else if (e_perp .lt. 0.) & write(*,*)' error in e_perp= ',e_perp beta_perp = (transcorr(1,1) + transcorr(3,3)) * avgpz / 2. alpha_perp = - (transcorr(1,2) + transcorr(3,4)) / 2. Ldim = (transcorr(1,4) - transcorr(3,2)) / 2. end if ! check_axisym = 0. denom = transcorr(1,1) + transcorr(3,3) denom = denom**2 if (denom .gt. 0.) then asymmetry = transcorr(1,2)**2 + & (transcorr(1,1)-transcorr(3,3))**2 check_axisym = check_axisym + asymmetry/denom end if denom = transcorr(2,2) + transcorr(4,4) denom = denom**2 if (denom .gt. 0.) then asymmetry = transcorr(2,4)**2 + & (transcorr(2,2)-transcorr(4,4))**2 check_axisym = check_axisym + asymmetry/denom end if denom = transcorr(1,1)*transcorr(2,2) + & transcorr(3,3)*transcorr(4,4) + transcorr(1,1)*transcorr(4,4) & + transcorr(3,3)*transcorr(2,2) if (denom .gt. 0.) then asymmetry = (transcorr(1,4)+transcorr(3,2))**2 + & (transcorr(1,2)-transcorr(3,4))**2 check_axisym = check_axisym + asymmetry/denom end if return end c ********************************************************* subroutine SETUP c IMPLICIT NONE c common/CONTROL/sigma_cut,rffreq,transcut,longcut, & bzaxis,refcharge,refmass,reftime,refenergy,zloc, & pzmax,pzmin,curr_region,vp_model,rf_model,useptype, & ampl_corr,disp_corr real*8 sigma_cut,rffreq,transcut,longcut real*8 bzaxis,refcharge,refmass,reftime,refenergy,zloc real*8 pzmax,pzmin integer curr_region,vp_model,rf_model,useptype logical ampl_corr,disp_corr c integer ioc character*80 title character*40 file_name_in, file_name_out, file_name_vp ! ! input file structure: ! description ! input file name ! output file name ! type of particle (0 = all types) ! Pzmin, Pzmax (if Pzmax<=0. then no cut) ! vector potential model, Bz on axis ! transverse, longitudinal emittance cutoffs ! subtract out linear correlations with energy and time? ! subtract out nonlinear correlations with transverse amplitude? ! time overlap model, RF frequency (in units of MHz) ! sigma cut (no cut if <=0.) ! c default values useptype=0 pzmin=0.0 pzmax=0.0 vp_model=1 bzaxis=0.0 transcut=0.0 longcut=0.0 disp_corr = .true. ampl_corr = .true. rf_model = 0 rffreq=0.0 sigma_cut=0.0 ! open(unit=1,file='emitcalc.inp',status='old',iostat=ioc) if( ioc .ne. 0 ) then write(*,*) 'Error opening emitcalc.inp' stop end if ! read(1,'(a80)',iostat=ioc) title read(1,*,iostat=ioc) file_name_in read(1,*,iostat=ioc) file_name_out read(1,*,iostat=ioc) useptype read(1,*,iostat=ioc) pzmin,pzmax read(1,*,iostat=ioc) vp_model,bzaxis,file_name_vp read(1,*,iostat=ioc) transcut,longcut read(1,*,iostat=ioc) disp_corr read(1,*,iostat=ioc) ampl_corr read(1,*,iostat=ioc) rf_model,rffreq read(1,*,iostat=ioc) sigma_cut ! close(1) ! if (vp_model .ge. 3) then open(unit=4,file=file_name_vp,status='unknown',iostat=ioc) if( ioc .ne. 0 ) then write(*,*) 'Error opening vector potential file ', & file_name_vp write(*,*) 'switching to vp model 2, numerical fit' vp_model = 2 else call INPUT_VP end if close(4) end if if (rffreq .eq. 0.) rf_model = 0 ! open(unit=3,file=file_name_in,status='old',iostat=ioc) if( ioc .ne. 0 ) then write(*,*) 'Error opening data file ',file_name_in stop end if read(3,'(a80)') title ! title card on beam output file read(3,*) read(3,*) c open(unit=2,file=file_name_out,status='unknown',iostat=ioc) write(2,*)title write(*,*)title write(2,*)' output from Emitcalc version 1.0' write(*,*)' Emitcalc version 1.0' write(2,*)' settings: ' write(*,*)' settings: ' write(2,*)' input file: ',file_name_in write(*,*)' input file: ',file_name_in write(2,*)' output file: ',file_name_out write(*,*)' output file: ',file_name_out write(2,*)' particle type: ',useptype write(*,*)' particle type: ',useptype write(2,*)' pzmin/pzmax: ',pzmin,'/',pzmax write(*,*)' pzmin/pzmax: ',pzmin,'/',pzmax if (vp_model .eq. 0) then write(2,*)' vp model: ',vp_model,' Bz on axis: ',bzaxis write(*,*)' vp model: ',vp_model,' Bz on axis: ',bzaxis else if (vp_model .le. 2) then write(2,*)' vp model: ',vp_model write(*,*)' vp model: ',vp_model else write(2,*)' vp model: ',vp_model, & ' input file: ',file_name_vp write(*,*)' vp model: ',vp_model, & ' input file: ',file_name_vp end if write(2,*)' trans/long cuts: ',transcut,longcut write(*,*)' trans/long cuts: ',transcut,longcut if (disp_corr) then write(2,*)' remove linear dispersions' write(*,*)' remove linear dispersions' else write(2,*)' do not remove linear dispersions' write(*,*)' do not remove linear dispersions' end if if (ampl_corr) then write(2,*)' subtract out amplitude correlation' write(*,*)' subtract out amplitude correlation' else write(2,*)' do not subtract out amplitude correlation' write(*,*)' do not subtract out amplitude correlation' end if if (rf_model .eq. 0) then write(2,*)' not periodic in time' write(*,*)' not periodic in time' else write(2,*)' RF frequency (MHz): ',rffreq, & ' RF model: ',rf_model write(*,*)' RF frequency (MHz): ',rffreq, & ' RF model: ',rf_model end if write(2,*)' sigma cut: ',sigma_cut write(*,*)' sigma cut: ',sigma_cut write(2,*) write(2,30) 30 format(t1,'regn #',t10,'Z',t22,'Bz0', & t34,'',t46,'',t58,'',t70,'', & t82,'',t94,'total #',t106,'# in cuts', & t118,'emit_1',t130,'emit_2',t142,'emit_3', & t154,'e_perp',t166,'e_long', & t178,'beta_perp',t190,'alpha_perp',t202,'Ldim', & t214,'beta_long',t226,'alpha_long', & t238,'D_x',t250,'D_y',t262,'Corr_E',t274,'Corr_t', & t286,'sigma_E',t298,'asymmetry') c return end c ********************************************************* real*8 function TIMESHIFT(time,freq) c c returns time shifted by integer number of periods so that c new time is between +/- 1/(2*freq) IMPLICIT NONE c real*8 time, freq c integer count c if (freq .ne. 0.) then count = 0.5 + time * freq if (count .lt. 0.) count = count - 1. TIMESHIFT = time - INT(count) / freq else TIMESHIFT = time end if c return end c ********************************************************* subroutine INPUT_VP c IMPLICIT NONE c common/VP/vp_axi(3),vp_x(5,6),vp_y(5,6) real*8 vp_axi,vp_x,vp_y c integer ioc c read(4,*,iostat=ioc)vp_x(1,1),vp_x(1,2) read(4,*,iostat=ioc)vp_x(2,1),vp_x(2,2),vp_x(2,3) read(4,*,iostat=ioc)vp_x(3,1),vp_x(3,2),vp_x(3,3),vp_x(3,4) read(4,*,iostat=ioc)vp_x(4,1),vp_x(4,2),vp_x(4,3),vp_x(4,4), & vp_x(4,5) read(4,*,iostat=ioc)vp_x(5,1),vp_x(5,2),vp_x(5,3),vp_x(5,4), & vp_x(5,5),vp_x(5,6) read(4,*,iostat=ioc) read(4,*,iostat=ioc)vp_y(1,1),vp_y(1,2) read(4,*,iostat=ioc)vp_y(2,1),vp_y(2,2),vp_y(2,3) read(4,*,iostat=ioc)vp_y(3,1),vp_y(3,2),vp_y(3,3),vp_y(3,4) read(4,*,iostat=ioc)vp_y(4,1),vp_y(4,2),vp_y(4,3),vp_y(4,4), & vp_y(4,5) read(4,*,iostat=ioc)vp_y(5,1),vp_y(5,2),vp_y(5,3),vp_y(5,4), & vp_y(5,5),vp_y(5,6) c return end c ********************************************************* subroutine CALC_VP(bzaxis) c c numerical fit to Bz = B0 + B1 r^2 + B2 r^4 c IMPLICIT NONE c integer npmax parameter ( npmax = 100000 ) c common/REGION/xpar(3,npmax),ppar(3,npmax),tpar(npmax), & epar(npmax),bzpar(npmax),evtwtpar(npmax), & iptypar(npmax),usepar(npmax),npart real*8 xpar,ppar,tpar,epar,bzpar,evtwtpar integer iptypar,usepar,npart c common/VP/vp_axi(3),vp_x(5,6),vp_y(5,6) real*8 vp_axi,vp_x,vp_y c real*8 bzaxis c integer i,j,ip real*8 corrmatrix(3,3),avg(3),bfit(3),vector3(3),rsq,weight,wtot real*8 denom c wtot = 0. do i=1,3 avg(i)=0. do j=1,3 corrmatrix(i,j)=0. end do end do do ip=1,npart if (usepar(ip) .eq. 0) then weight = evtwtpar(ip) wtot = wtot + weight rsq = xpar(1,ip)**2 + xpar(2,ip)**2 vector3(1) = bzpar(ip) vector3(2) = rsq vector3(3) = rsq*rsq do i=1,3 avg(i)=avg(i)+weight*vector3(i) do j=1,3 corrmatrix(i,j)=corrmatrix(i,j)+ & weight*vector3(i)*vector3(j) end do end do end if end do if (wtot>0.) then do i=1,3 avg(i)=avg(i)/wtot end do do i=1,3 do j=1,3 corrmatrix(i,j)=corrmatrix(i,j)/wtot - & avg(i)*avg(j) end do end do end if c denom = corrmatrix(2,2)*corrmatrix(3,3) - corrmatrix(2,3)**2 if (corrmatrix(2,2) .le. 0.) then c no variation bfit(3) = 0. bfit(2) = 0. bfit(1) = avg(1) else if (denom .le. 0.) then c use linear fit bfit(3) = 0. bfit(2) = corrmatrix(1,2) / corrmatrix(2,2) bfit(1) = avg(1) - bfit(2) * avg(2) else bfit(3) = (corrmatrix(1,3)*corrmatrix(2,2) - & corrmatrix(2,3)*corrmatrix(1,2)) / denom bfit(2) = (corrmatrix(1,2) - corrmatrix(2,3) * bfit(3)) / & corrmatrix(2,2) bfit(1) = avg(1) - avg(2) * bfit(2) - avg(3) * bfit(3) end if c vp_axi(1) = bfit(1) / 2. vp_axi(2) = bfit(2) / 4. vp_axi(3) = bfit(3) / 6. c bzaxis = bfit(1) c return end c ********************************************************* subroutine CALC_TC(freq,centertime,conf) c c for systems periodic in time, find central time c IMPLICIT NONE c integer npmax parameter ( npmax = 100000 ) c common/CONST/cc,ee,pi,twopi,avagn,relcl,mass(5) real*8 cc,ee,pi,twopi,avagn,relcl,mass c common/REGION/xpar(3,npmax),ppar(3,npmax),tpar(npmax), & epar(npmax),bzpar(npmax),evtwtpar(npmax), & iptypar(npmax),usepar(npmax),npart real*8 xpar,ppar,tpar,epar,bzpar,evtwtpar integer iptypar,usepar,npart c real*8 freq,centertime,conf c real*8 tcos(npmax),tsin(npmax),xbar,ybar,wtot,weight,avgtime integer ip c wtot = 0. avgtime = 0. xbar = 0. ybar = 0. do ip=1,npart if (usepar(ip) .eq. 0) then weight = evtwtpar(ip) wtot = wtot + weight if (freq .eq. 0.) then avgtime = avgtime + weight*tpar(ip) else xbar = xbar + weight*DCOS(2.*pi*tpar(ip)*freq) ybar = ybar + weight*DSIN(2.*pi*tpar(ip)*freq) end if end if end do if (wtot>0.) then avgtime = avgtime/wtot xbar = xbar/wtot ybar = ybar/wtot end if c if (freq .eq. 0.) then if (wtot .eq. 0.) then conf = 0. else conf = 1. end if else conf = DSQRT(xbar**2 + ybar**2) end if if (conf .eq. 0.) then centertime = 0. else centertime = DATAN2(ybar,xbar) / (2.*pi*freq) end if c return end c ********************************************************* subroutine EMITXYZ(esq,f,g,h) c c associate emittances with x, y, z planes, and reorder them c actually uses esq = (m c emit)^2 c f = sum of 2x2 determinants across a row c g = sum of codeterminants across a row c h = projected emittances (just in case other method fails) c IMPLICIT NONE c real*8 SMALL parameter (SMALL = 1.d-8) c real*8 esq(3),f(3),g(3),h(3) c integer i,isingle,ntuple,iflip real*8 detxy,detxz,detyz,ex,ey,ez,best real*8 invx,invy,esingle,epair real*8 factor1,factor2,factor3,factor4 real*8 savex,savey,savez,temp c ntuple=1 if (dabs(esq(1)-esq(2)).le.small*(esq(1)+esq(2))) then ntuple=ntuple+1 isingle = 3 end if if (dabs(esq(1)-esq(3)).le.small*(esq(1)+esq(3))) then ntuple=ntuple+1 isingle = 2 end if if (dabs(esq(2)-esq(3)).le.small*(esq(2)+esq(3))) then ntuple=ntuple+1 isingle = 1 end if if (ntuple.ge.3) then ! three emittances basically equal c rearrange to be in same order as projected emittances call REORDER(3,h,esq) return else if (ntuple.eq.2) then c define solitary and double emittance esingle=esq(isingle) epair=1.d0 do i=1,3 if (i.ne.isingle) epair=epair*esq(i) end do epair=dsqrt(epair) invx=(esingle-f(1))/(esingle-epair) invy=(esingle-f(2))/(esingle-epair) if (invx.ge.0.5d0) then ! either ex=ey or ex=ez if (invy.ge.1.d0-0.5d0*invx) then ! ex=ey esq(isingle)=esq(3) esq(3)=esingle else ! ex=ez esq(isingle)=esq(2) esq(2)=esingle end if else ! should have ey=ez if ((invy.ge.0.5d0).and.(invy.ge.2.d0*(1.d0-invx))) & then esq(isingle)=esq(1) esq(1)=esingle else ! nothing seems to work (large coupling?) write(*,*)'difficulty ordering eigenvalues' call REORDER(3,h,esq) end if end if return end if c c if made it this far, ntuple=1, all three emittances are distinct c don't have to worry about zeroes in denominators c first use detxz,detyz to choose ez; independent of order, ex and ey c then use detxy to choose ordering between ex and ey best=4.d0 do i=1,3 ez=esq(i) ex=esq(1) ey=esq(2) if (i.eq.1) ex=esq(3) if (i.eq.2) ey=esq(3) factor1=g(1)-ex*ey+ez*(f(1)-ez) factor2=g(2)-ex*ez+ez*(f(2)-ey) if ((factor1.ne.0.).and.((factor1+factor2).ne.0)) then detxz = - (ez-ex)*(ez-ey)/factor1 detyz = factor1/(factor1+factor2) c detxz-1. represents 'amount' of coupling in some sense, c magnitude must be <=1 else corresponds to flipping the beam if ((dabs(detxz-1.).le.1).and.(dabs(detyz-1.).le.1)) & then ! good values, check if detxy okay factor3=g(1)-ex*ey+ex*(f(1)-ez) if (factor3.ne.0.) then detxy = (ex-ey)*(ex-ez)/(detxz*factor3) if (dabs(detxy-1.).le.1) then ! valid option temp=dabs(detxy-1.)+dabs(detxz-1.) & +dabs(detyz-1.) if(temp.lt.best) then ! best so far if (best.le.3.d0) then write(*,*)'multiple orderings found' end if best=temp savez=ez savex=ex savey=ey end if end if end if c try exchanging ex and ey factor4=g(1)-ey*ex+ey*(f(1)-ez) if (factor4.ne.0.) then detxy = (ey-ex)*(ey-ez)/(detxz*factor4) if (dabs(detxy-1.).le.1) then ! valid option temp=dabs(detxy-1.)+dabs(detxz-1.) & +dabs(detyz-1.) if(temp.lt.best) then ! best so far if (best.le.3.d0) then write(*,*)'multiple orderings found' end if best=temp savez=ez savex=ey savey=ex end if end if end if end if ! end of checking detxy end if ! end of test on detxz and detyz end do ! loop over choosing ez plane c if (best.gt.3.d0) then ! no good solution found write(*,*)'difficulty ordering eigenvalues' call REORDER(3,h,esq) else ! did find a valid choice esq(1)=savex esq(2)=savey esq(3)=savez end if c return end c ********************************************************* subroutine REORDER(ndim,x,y) c c reorder vector y, not sorted but in same order as x c IMPLICIT NONE c integer ndim real*8 x(ndim),y(ndim) c integer NMAX parameter (NMAX=500) integer ishift(NMAX),i,j,itemp real*8 temp,new(NMAX) c do i=1,ndim ishift(i)=i end do do i=1,ndim-1 do j=i+1,ndim if (x(i) .gt. x(j)) then temp = x(i) x(i)=x(j) x(j)=temp itemp=ishift(i) ishift(i)=ishift(j) ishift(j)=itemp end if end do end do do i=1,ndim-1 do j=i+1,ndim if (y(i) .gt. y(j)) then temp = y(i) y(i)=y(j) y(j)=temp end if end do end do do i=1,ndim new(i) = y(ishift(i)) end do do i=1,ndim y(i) = new(i) end do c return end c ********************************************************* subroutine GET_VP(bz,x,y,vpx,vpy) c c get transverse components of vector potential for particle c IMPLICIT NONE c common/CONTROL/sigma_cut,rffreq,transcut,longcut, & bzaxis,refcharge,refmass,reftime,refenergy,zloc, & pzmax,pzmin,curr_region,vp_model,rf_model,useptype, & ampl_corr,disp_corr real*8 sigma_cut,rffreq,transcut,longcut real*8 bzaxis,refcharge,refmass,reftime,refenergy,zloc real*8 pzmax,pzmin integer curr_region,vp_model,rf_model,useptype logical ampl_corr,disp_corr c common/VP/vp_axi(3),vp_x(5,6),vp_y(5,6) real*8 vp_axi,vp_x,vp_y c integer i,j real*8 bz,x,y,vpx,vpy real*8 Aphi,rsq c rsq = x**2 + y**2 if (vp_model .eq. 0) then Aphi = bzaxis / 2. vpx = - Aphi * y vpy = Aphi * x else if (vp_model .eq. 1) then Aphi = (bzaxis + bz)/4. vpx = - Aphi * y vpy = Aphi * x else if (vp_model .eq. 2) then Aphi = vp_axi(1) + vp_axi(2)*rsq + vp_axi(3)*(rsq**2) vpx = - Aphi * y vpy = Aphi * x else if (vp_model .ge. 3) then vpx = 0. vpy = 0. do i=1,5 do j=1,i+1 vpx = vpx + vp_x(i,j)*(x**(i+1-j))*(y**(j-1)) vpy = vpy + vp_y(i,j)*(x**(i+1-j))*(y**(j-1)) end do end do end if c return end c ********************************************************* subroutine EMIT3D(refmass,matrix,emittances) c c calculate the 3 separate invariant emittances c IMPLICIT NONE c common/CONST/cc,ee,pi,twopi,avagn,relcl,mass(5) real*8 cc,ee,pi,twopi,avagn,relcl,mass c real*8 refmass,matrix(6,6),emittances(3) real*8 projected(3),sum,product,det,newmatrix(6,6) real*8 factor1,factor2,drow(3),crow(3) integer i,j,indx(6) real*8 SMALLDET,CODET c c need to rescale t to (ct), so dimensions work out properly c including in cross-terms do i=1,6 if (i .eq. 5) then factor1 = cc else factor1 = 1. end if do j=1,6 if (j .eq. 5) then factor2 = cc else factor2 = 1. end if newmatrix(i,j) = matrix(i,j) * factor1 * factor2 end do end do do i=1,3 projected(i) = SMALLDET(i,i,6,newmatrix) end do sum = 0. product = 0. do i=1,3 drow(i)=0. crow(i)=0. do j=1,3 drow(i) = drow(i) + SMALLDET(i,j,6,newmatrix) crow(i) = crow(i) + CODET(i,j,6,newmatrix) end do sum = sum + drow(i) product = product + crow(i) end do call LUDCMP(newmatrix,6,6,indx,det) if (det .ne. 0.) then ! else singular do i=1,6 det = det * newmatrix(i,i) end do end if c c solve x^3 - sum x^2 + product x - det = 0, c solutions are x=(m c emittance)^2 c call SOLVE_CUBIC(-sum,product,-det,emittances) c c emittances need to be associated with x, y, z planes: c this is a rough method, which tries to minimize the c necessary decoupling matrix c call EMITXYZ(emittances,drow,crow,projected) c c take square root, and divide by (mc) to get emittance in meters c do i=1,3 if (emittances(i) .gt. 0.) then emittances(i) = DSQRT(emittances(i))/refmass end if end do c return end c ********************************************************* real*8 function SMALLDET(i,j,ndim,matrix) c c 2x2 determinant of submatrix using conjugate coordinate pair c ndim = # dimensions should be even c IMPLICIT NONE c integer i,j,ndim real*8 matrix(ndim,ndim) c SMALLDET = matrix(2*i-1,2*j-1)*matrix(2*i,2*j) - & matrix(2*i,2*j-1)*matrix(2*i-1,2*j) c return end c ********************************************************* real*8 function CODET(i,j,ndim,matrix) c c codeterminant of submatrix, from conjugate coordinate pair c IMPLICIT NONE c integer i,j,ndim real*8 matrix(ndim,ndim) integer NMAX parameter (NMAX = 500) integer m,n,tm,tn,indx(NMAX) real*8 reduced(NMAX,NMAX),determ c do m=1,ndim if ((m .ne. 2*i-1) .and. (m .ne. 2*i)) then if (m .lt. 2*i-1) then tm = m else tm = m - 2 end if do n=1,ndim if ((n .ne. 2*j-1) .and. (n .ne. 2*j)) then if (n .lt. 2*j-1) then tn = n else tn = n - 2 end if reduced(tm,tn)=matrix(m,n) end if end do end if end do c call LUDCMP(reduced,ndim-2,NMAX,indx,determ) ! if (determ .ne. 0.) then ! else singular do m=1,ndim-2 determ = determ * reduced(m,m) end do end if CODET = determ c return end c ********************************************************* subroutine SOLVE_CUBIC(a,b,c,solutions) c c get nonnegative solutions to cubic y=x^3 + a x^2 + b x + c c there must be a better way of doing this, hard to get accurate though c IMPLICIT NONE c real*8 SMALL parameter (SMALL = 1.d-12) c integer i,j real*8 a,b,c,solutions(3) real*8 x0,y0,y1,sigma,quadr real*8 u,us,f,df,err,close,xc,yc,d,ds,xcs c x0 = -a/3.d0 y0 = x0**3 + a * x0**2 + b * x0 + c quadr = a**2-3.d0*b if (quadr.lt.0.d0) then if (dabs(quadr).gt.SMALL) then write(*,*)'complex roots ' else quadr = 0.d0 end if else quadr = DSQRT(quadr) end if y1 = quadr/3.d0 sigma = SIGN(1.,y0) if (c.eq.0.d0) then if ((a.gt.0.d0).or.(b.lt.0.d0)) then write(*,*)'negative roots' end if solutions(1)=0.d0 quadr=a**2-4.d0*b if (quadr.lt.0.d0) then if (dabs(quadr).gt.SMALL*a**2) then write(*,*)'complex roots ' else quadr = 0.d0 end if else quadr = DSQRT(quadr) end if solutions(2) = 0.5d0*(-a + quadr) if (solutions(2).eq.0.d0) then solutions(3) = 0.5d0*(-a - quadr) else solutions(3)=2.d0*b/solutions(2) end if else if ((a.ge.0.d0).or.(b.le.0.d0)) then write(*,*)'negative roots ' end if if (y0.le.0.d0) then u=2.d0*y1+x0 else u=-2.d0*y1+x0 end if do i=1,20 f = u**3 + a * u**2 + b * u + c df = 3.d0*u**2 + 2.d0* a * u + b us = u u = u - f/df end do if (u.eq.0.d0) then err = (u-us)/a else err = (u-us)/u end if solutions(1)=u xc=x0+sigma*y1 do i=1,20 xcs = xc xc = xc - (3.d0*xc**2+2.d0*a*xc+b)/(6.d0*xc+2.d0*a) end do if (xc.eq.0.d0) then err = (xc-xcs)/y1 else err = (xc-xcs)/xc end if yc = xc**3 + a * xc**2 + b * xc + c close=-sigma*yc/(3.d0*y1) if (close.lt.0.d0) then if (dabs(close).gt.SMALL*y1**2) then write(*,*)'complex roots ' else close = 0.d0 end if end if if (close.lt.SMALL*y1**2) then d=dsqrt(close) do i=1,20 ds = d d= dsqrt(close-sigma*d**3/(3.d0*y1)) end do if (d.eq.0.d0) then err = (d-ds)/y1 else err = (d-ds)/d end if solutions(2)=xc+d d=-dsqrt(close) do i=1,20 ds = d d= -dsqrt(close-sigma*d**3/(3.d0*y1)) end do if (d.eq.0.d0) then err = (d-ds)/y1 else err = (d-ds)/d end if solutions(3)=xc+d else if (u.eq.0.d0) then quadr=a**2-4.d0*b else quadr=(a+u)**2+4.d0*c/u end if if (quadr.lt.0.d0) then if (dabs(quadr).gt.SMALL) then write(*,*)'complex roots ' else quadr = 0.d0 end if else quadr = DSQRT(quadr) end if solutions(2) = 0.5d0*(-a-u+quadr) if (u.eq.0.d0) then solutions(3) = 0.5d0*(-a-u-quadr) else solutions(3) = -c/(u*solutions(2)) end if if (quadr.lt.SMALL) then if (u.eq.0.d0) then solutions(2)=dsqrt(b) solutions(3)=solutions(2) else solutions(2)=x0+sigma*y1 solutions(3)=solutions(2) end if else do j=2,3 u=solutions(j) do i=1,20 f = u**3 + a * u**2 + b * u + c df = 3.d0*u**2 + 2.d0* a * u + b us = u u = u - f/df end do if (u.eq.0.d0) then err = (u-us)/a else err = (u-us)/u end if solutions(j)=u end do end if end if ! end of test on near-double root end if ! end of test on c=0 c return end c *********************************************************