Description of TouschekPiwinskiLifeTime

0001 function [Tl,contributionsTL]=TouschekPiwinskiLifeTime(r,dpp,Ib,varargin)
0002 % function [Tl,contributionsTL]=TouschekPiwinskiLifeTime(r,dpp,Ib)
0003 %
0004 % evaluates Touschek Lifetime using Piwinski formula
0005 %
0006 % TouschekPiwinskiLifeTime(latticeATring,momentumaperturecolumnvector,0.002)
0007 %  or
0008 % TouschekPiwinskiLifeTime(
0009 %  latticeATring,
0010 %  momentumaperturecolumnvector,  % column array (size of r or positions)
0011 %  current per bunch in A,                 % scalar
0012 %  positions where to evaluate,  %(default all elements with length>0 )  column array
0013 %  emittancex, %(default atx modemittance(1))   scalar
0014 %  emittancey, %(default 10 pm)               scalar
0015 %  integration method,  % (default quad, may be: 'integral', 'quad', 'trapz', 'elegantLike', 'Approximate')
0016 %  energy sperad,    % scalar
0017 %  bunch length,        % scalar
0018 %               )
0019 %
0020 %  contributionsTL 1/T contribution at each element
0021 %
0022 %  Tl  Lifetime in seconds 1/Tl=sum(contributionsTL.*L)/sum(L);
0023 %
0024 %
0025 % "The Touscheck Effect in strong focusing storage rings"
0026 % A.Piwinski, DESY 98-179, November 1998
0027 %
0028 % "Touscheck Effect calculation and its applications to a transport line"
0029 % A.Xiao M. Borland, Proceedings of PAC07, Albuquerque, New Mexico, USA
0030 %
0031 %
0032 % created 31-10-2012
0033 % updated 22-01-2013 corrected and compared with elegant
0034 
0035 %ensure a column lattice
0036 r=reshape(r,numel(r),1);
0037 
0038 e0 = 1.60217646e-19; %Coulomb
0039 r0 = 2.817940327e-15; %m %  classical electron radius
0040 spl=299792458; % speed of ligth
0041 
0042 naddvar=length(varargin);
0043 if naddvar>=1
0044     positions=varargin{1};
0045 else
0046     %    positions=1:length(r);
0047     
0048     % positions default= non zero length elements
0049     positions=findcells(r,'Length');
0050     L=getcellstruct(r,'Length',positions);
0051     positions=positions(L>0);
0052     size(positions);
0053 end
0054 
0055 % get optics
0056 [lo,pa]=atx(r,0,positions);
0057 
0058 emitx=pa.modemittance(1);
0059 emity=emitx./2;
0060 integrationmethod='quad';
0061 sigp=pa.espread; % relative momentum spread
0062 sigs=pa.blength; % bunch length
0063 
0064 if naddvar==2
0065     positions=varargin{1};
0066     emitx=varargin{2};
0067     
0068     disp(['set defaults: ey=ex/2'])
0069     disp([' integration method is quad,'])
0070     disp([' energy sperad, bunch length from ATX'])
0071     
0072 elseif naddvar==3
0073     positions=varargin{1};
0074     emitx=varargin{2};
0075     emity=varargin{3};
0076     disp(['set defaults: '])
0077     disp([' integration method is quad,'])
0078     disp([' energy sperad, bunch length from ATX'])
0079     
0080 elseif naddvar==4
0081     positions=varargin{1};
0082     emitx=varargin{2};
0083     emity=varargin{3};
0084     integrationmethod=varargin{4};
0085     disp(['set defaults: '])
0086     disp([' energy sperad, bunch length from ATX'])
0087     
0088 elseif naddvar==5
0089     positions=varargin{1};
0090     emitx=varargin{2};
0091     emity=varargin{3};
0092     integrationmethod=varargin{4};
0093     sigp=varargin{5};
0094     disp(['set defaults: '])
0095     disp(['bunch length from ATX'])
0096     
0097 elseif naddvar==6
0098     positions=varargin{1};
0099     emitx=varargin{2};
0100     emity=varargin{3};
0101     integrationmethod=varargin{4};
0102     sigp=varargin{5};
0103     sigs=varargin{6};
0104     
0105 else
0106     disp(['set defaults: ey=ex/2'])
0107     disp([' integration method is quad,'])
0108     disp([' energy sperad, bunch length, x emittance from ATX'])
0109     disp([' evaluation at all points with non zero length'])
0110 end
0111 
0112 % if dpp is a scalar assume constant momentum aperture.
0113 if numel(dpp)==1
0114     dpp=dpp*ones(size(positions'));
0115 end
0116 
0117 dppinput=dpp;
0118 Tlcol=zeros(1,size(dppinput,2));
0119 
0120 for dppcolnum=1:size(dppinput,2)
0121     
0122     dpp=dppinput(:,dppcolnum);
0123     
0124     
0125     Circumference=findspos(r,length(r)+1);
0126     
0127     E0=r{1}.('Energy');%
0128     
0129     Nb = Ib/(spl/Circumference)/e0; %Number of particle per bunch.
0130     
0131     relgamma = E0/0.510998928e6;
0132     relbeta=sqrt(1-1./relgamma.^2);
0133     
0134     aaa=cat(1,lo.alpha);
0135     bbb=cat(1,lo.beta);
0136     ddd=cat(2,lo.Dispersion)';
0137     
0138     bx=bbb(:,1); % betax
0139     by=bbb(:,2); % betay
0140     Dx=ddd(:,1);
0141     Dy=ddd(:,3);
0142     
0143     ax=aaa(:,1);
0144     ay=aaa(:,2);
0145     Dpx=ddd(:,2);
0146     Dpy=ddd(:,4);
0147     
0148     
0149     
0150     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0151     %%%%%%%%
0152     %%%%%%%% From here calculation takes place.
0153     %%%%%%%%
0154     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0155     
0156     
0157     sigxb=sqrt(emitx.*bx);
0158     sigyb=sqrt(emity.*by);
0159     
0160     sigx=sqrt(emitx.*bx+sigp.^2.*Dx.^2);
0161     sigy=sqrt(emity.*by+sigp.^2.*Dy.^2); %%%  was mistaken! it was Dx!!!!!!
0162     
0163     
0164     Dtx=Dx.*ax+Dpx.*bx;%  % alpha=-b'/2
0165     Dty=Dy.*ay+Dpy.*by;%
0166     
0167     sigtx=sqrt(sigx.^2+sigp.^2.*Dtx.^2);
0168     sigty=sqrt(sigy.^2+sigp.^2.*Dty.^2);
0169     
0170     sigtx2=sigx.^2+sigp.^2.*Dtx.^2;
0171     sigty2=sigy.^2+sigp.^2.*Dty.^2;
0172     
0173     sigp2=sigp.^2;
0174     Dx2=Dx.^2;
0175     Dy2=Dy.^2;
0176     Dtx2=Dtx.^2;
0177     Dty2=Dty.^2;
0178     sigxb2=sigxb.^2;
0179     sigyb2=sigyb.^2;
0180     
0181     
0182     sighinv2=1./(sigp2) +(Dx2+Dtx2)./(sigxb2) + (Dy2+Dty2)./(sigyb2);
0183     
0184     sigh=sqrt(1./sighinv2);
0185     
0186     um=relbeta.^2*dpp.^2;
0187     
0188     B1=1./(2.*(relbeta.^2).*(relgamma.^2)).*( (bx.^2./(sigxb.^2)).*(1-(sigh.^2.*Dtx.^2./(sigxb.^2))) + (by.^2./(sigyb.^2)).*(1-(sigh.^2.*Dty.^2./(sigyb.^2))));
0189     
0190     B2sq=1./(4.*(relbeta.^4).*(relgamma.^4)).*((bx.^2./(sigxb.^2)).*(1-(sigh.^2.*Dtx.^2./(sigxb.^2)))-(by.^2./(sigyb.^2)).*(1-(sigh.^2.*Dty.^2./(sigyb.^2)))).^2+(sigh.^4.*bx.^2.*by.^2.*Dtx.^2.*Dty.^2)./((relbeta.^4).*(relgamma.^4).*sigxb.^4.*sigyb.^4);
0191     
0192     B2=sqrt(B2sq);
0193     
0194     em=bx.^2.*sigx.^2./(relbeta.^2.*relgamma.^2.*sigxb.^2.*sigtx2).*um;
0195     
0196     val=zeros(size(B1));
0197     
0198     km=atan(sqrt(um));
0199     
0200     FpiWfact=sqrt(pi.*(B1.^2-B2.^2)).*um;
0201     
0202     for ii=1:length(positions)
0203         
0204         % choose integration method
0205         switch integrationmethod
0206             
0207             case 'infiniteintegral'
0208                 
0209                 val(ii)= integral(@(u)TLT_IntPiw(u,um(ii),B1(ii),B2(ii)),um(ii),Inf); %,...um(ii)*1e4
0210                 
0211             case 'integral'
0212                 
0213                 val(ii) = integral(@(k)TLT_IntPiw_k(k,km(ii),B1(ii),B2(ii)),km(ii),pi/2); %,...,'Waypoints',evalpoints um(ii)*1e4
0214                 
0215             case 'quad'
0216                 
0217                 val(ii)= quad(@(k)TLT_IntPiw_k(k,km(ii),B1(ii),B2(ii)),km(ii),pi/2); %,...,'Waypoints',evalpoints um(ii)*1e4
0218             case 'trapz'
0219                 
0220                 
0221                 k=linspace(km(ii),pi/2,10000);
0222                 val(ii)= trapz(k,TLT_IntPiw_k(k,km(ii),B1(ii),B2(ii))); %,...,'Waypoints',evalpoints um(ii)*1e4
0223                 
0224             case 'elegantLike'
0225                 
0226                 val(ii)=IntegrateLikeElegant(km(ii),B1(ii),B2(ii));
0227                 
0228             case 'Approximate'
0229                 
0230                 val(ii)=integral(@(e)Cfun(e,em(ii)),em(ii),Inf);
0231                 
0232             otherwise % use default method quad (backward compatible)
0233                 
0234                 
0235                 val(ii)=quad(@(k)TLT_IntPiw_k(k,km(ii),B1(ii),B2(ii)),km(ii),pi/2); %,...,'Waypoints',evalpoints um(ii)*1e4
0236                 
0237         end
0238     end
0239     
0240     
0241     
0242     switch integrationmethod
0243         case 'infiniteintegral'
0244             frontfact=(r0.^2.*spl.*Nb)./(8.*pi.*(relgamma.^2).*sigs.*sqrt(...
0245                 (sigx.^2).*(sigy.^2)-sigp.^4.*Dx.^2.*Dy.^2).*um).*FpiWfact;
0246             
0247         case {'integral' 'quad' } %'elegantLike'
0248             
0249             frontfact=(r0.^2.*spl.*Nb)./(8.*pi.*(relgamma.^2).*sigs.*sqrt(...
0250                 (sigx.^2).*(sigy.^2)-sigp.^4.*Dx.^2.*Dy.^2).*um).*2.*FpiWfact;
0251         case {'trapz' 'elegantLike'}
0252             
0253             frontfact=(r0.^2.*spl.*Nb.*sigh.*bx.*by)./(4.*sqrt(pi).*(relbeta.^2).*(relgamma.^4).*sigxb.^2.*sigyb.^2.*sigs.*sigp);
0254             
0255         case 'Approximate'
0256             
0257             frontfact=(r0.^2.*spl.*Nb.*bx)./(...
0258                 8.*pi.*(relbeta.^3).*(relgamma.^3).*...
0259                 sigxb.*sigyb.*sigs.*sqrt(sigtx2).*dpp.^2 ...
0260                 ).*em;
0261             
0262         otherwise
0263             
0264             frontfact=(r0.^2.*spl.*Nb)./(8.*pi.*(relgamma.^2).*sigs.*sqrt(...
0265                 (sigx.^2).*(sigy.^2)-sigp.^4.*Dx.^2.*Dy.^2).*um).*2.*FpiWfact;
0266             
0267     end
0268     contributionsTL=frontfact.*val;
0269     
0270     
0271     L=getcellstruct(r,'Length',positions);
0272     Tlcol(dppcolnum)=1/(1/sum(L)*sum(contributionsTL.*L));
0273     
0274 end
0275 
0276 Tl=length(Tlcol)/(sum(1./Tlcol));
0277 
0278 return
0279 
0280 
0281
TouschekPiwinskiLifeTime

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
