ringpara

rp = ringpara, use global THERING
rp = ringpara(THERING)
rp = ringpara(THERING,U0), supply total radiation loss in MeV
calculate various ring parameters
(1) The calculation of emittance, mcf, momentum spread, bunch length, damping time, etc 
is more accurate than atsummary.m because detailed
calculation of dispersion function and curly H function inside dipoles is performed. 
(2) calculate contribution of dispersion to vertical emittance.

Author: Xiaobiao Huang
created on 12/17/2007
Part of this code was modified from atsummary.m

Modified by Peace Chang (check if theta(ii) ~= 0.0)
Modified by S.Liuzzo and B.Nash (Dipole gradient may be in PolynomB(2),
also coupled damping added) 7/24/2014

0001 function rp = ringpara(THERING,varargin)
0002 %rp = ringpara, use global THERING
0003 %rp = ringpara(THERING)
0004 %rp = ringpara(THERING,U0), supply total radiation loss in MeV
0005 %calculate various ring parameters
0006 %(1) The calculation of emittance, mcf, momentum spread, bunch length, damping time, etc
0007 %is more accurate than atsummary.m because detailed
0008 %calculation of dispersion function and curly H function inside dipoles is performed.
0009 %(2) calculate contribution of dispersion to vertical emittance.
0010 %
0011 %Author: Xiaobiao Huang
0012 %created on 12/17/2007
0013 %Part of this code was modified from atsummary.m
0014 %
0015 %Modified by Peace Chang (check if theta(ii) ~= 0.0)
0016 %Modified by S.Liuzzo and B.Nash (Dipole gradient may be in PolynomB(2),
0017 %also coupled damping added) 7/24/2014
0018 
0019 if nargin==0
0020     global THERING; %#ok<TLEV>
0021 end
0022 Cq = 3.8319E-13; 
0023 a = findcells(THERING,'Energy');
0024 if isempty(a)
0025    gamma = 3000/PhysConstant.electron_mass_energy_equivalent_in_MeV.value;
0026 else
0027    gamma = THERING{a(1)}.Energy/(PhysConstant.electron_mass_energy_equivalent_in_MeV.value*1e6); 
0028 end
0029 
0030 dpindex = findcells(THERING,'BendingAngle');
0031 [tw,tune,chrom] = twissring(THERING,0,dpindex,'chrom',0.00001);
0032 beta = cat(1, tw.beta);
0033 alpha = cat(1, tw.alpha);
0034 mu = cat(1, tw.mu);
0035 disper = cat(1, tw.Dispersion);
0036 Dx  = disper(1:4:end);
0037 Dxp = disper(2:4:end);
0038 Dy  = disper(3:4:end);
0039 Dyp = disper(4:4:end);
0040 
0041 [tmptw,tune,chrom] = twissring(THERING,0,1:length(THERING),'chrom',0.00001);
0042 
0043 lendp = getcellstruct(THERING,'Length',dpindex); %bending magnet length
0044 lendp(lendp==0)=1;
0045 theta = getcellstruct(THERING,'BendingAngle',dpindex); %bending angle
0046 rho = lendp./theta;%THERING{dpindex(1)}.Length/THERING{dpindex(1)}.BendingAngle;
0047 
0048 I1 = 0;
0049 I2 = 0;
0050 I3 = 0;
0051 I4 = 0;
0052 I5 = 0;
0053 len = length(dpindex);
0054 curHavg1 = 1:len;
0055 for ii=1:len
0056   if theta(ii) ~= 0.0
0057       K = 0;
0058       Kk = 0;
0059       Kp = 0;
0060       if isfield(THERING{dpindex(ii)},'K')
0061           Kk = THERING{dpindex(ii)}.K;
0062       end
0063       if isfield(THERING{dpindex(ii)},'PolynomB')
0064           Kp = THERING{dpindex(ii)}.PolynomB(2);
0065       end
0066       if Kk~=Kp && (Kk~=0 && Kp~=0)
0067           warning('Values in K and PolynomB(2) are different. Using larger absolute value'); 
0068       end
0069       Ks=[Kk,Kp];
0070       [~,i]=max(abs(Ks));
0071       K=Ks(i);
0072       
0073     th1 = THERING{dpindex(ii)}.EntranceAngle;
0074     th2 = THERING{dpindex(ii)}.ExitAngle;
0075     [dI1,dI2,dI3,dI4,dI5,curHavg1(ii)] = calcRadInt(rho(ii),theta(ii), ...
0076          alpha(ii,1),beta(ii,1),Dx(ii),Dxp(ii),K,th1,th2);
0077     I1 = I1 + dI1;
0078     I2 = I2 + dI2;
0079     I3 = I3 + dI3;
0080     I4 = I4 + dI4;
0081     I5 = I5 + dI5;
0082   end
0083 end
0084 % curHavg = sum(curHavg1.*lendp./abs(rho))/sum(lendp);
0085 % %emittx =  Cq*gamma^2*curHavg/Jx/rho*1e9; %nm-rad
0086 % emittx =  Cq*gamma^2*curHavg/Jx*1e9; %nm-rad
0087 R = findspos(THERING, length(THERING)+1)/2/pi;
0088 alphac = I1/2/pi/R;
0089 U0 = 14.085*(gamma*PhysConstant.electron_mass_energy_equivalent_in_MeV.value/1000)^4*I2*1000.; %eV
0090 if nargin>=2
0091     fprintf('dipole radiation loss:  %4.5f keV\n', U0/1000.);
0092     U0 = varargin{1}*1e6; %convert MeV to eV
0093 end
0094 sigma_E = gamma*sqrt(Cq*I3/(2*I2+I4));
0095 Jx = 1-I4/I2;
0096 Jy = 1.00;
0097 Je = 2+I4/I2;
0098 emittx = Cq*gamma^2*I5/(I2-I4);
0099 
0100 % minimum emittance due to radiation 1/gamma cone (Handbook, Chao, eq23, pag 211)
0101 spos=findspos(THERING,1:length(THERING)+1);
0102 [tmptw,tune,chrom] = twissring(THERING,0,1:length(THERING),'chrom');
0103 betaall = cat(1, tmptw.beta);
0104 meanbetayovers=sum(betaall(:,2)'.*diff(spos))/spos(end);%mean(betaall(:,2));%
0105 meaninvr2=mean((1./rho.^2));%sum((1./rhos.^2).*diff(spos)/spos(end));%
0106 meaninvr3=mean((1./abs(rho.^3)));%sum((1./abs(rhos.^3)).*diff(spos)/spos(end));%
0107 emitty_lim = Cq*meanbetayovers/2/Jy*meaninvr3/meaninvr2;
0108 
0109 
0110 cspeed = PhysConstant.speed_of_light_in_vacuum.value; %m/s
0111 T0 = 2*pi*R/cspeed;
0112 alpha0 = U0/1.0e6/2/T0/(gamma*PhysConstant.electron_mass_energy_equivalent_in_MeV.value);
0113 alphax = Jx*alpha0;  %horizontal damping rate, 1/s
0114 alphay = Jy*alpha0;
0115 alphaE = Je*alpha0;
0116 
0117 rp.E0 = gamma*PhysConstant.electron_mass_energy_equivalent_in_MeV.value*1E6;
0118 rp.R = R;
0119 rp.alphac = alphac;
0120 rp.U0 = U0; %eV
0121 rp.sigma_E = sigma_E;
0122 rp.emittx = emittx;
0123 rp.emitty_lim= emitty_lim;
0124 rp.T0 = T0;
0125 rp.integrals = [I1,I2,I3,I4,I5];
0126 rp.dampingalpha = [alphax, alphay, alphaE];
0127 rp.dampingtime = 1./[alphax, alphay, alphaE];
0128 
0129 
0130 %compute coupled damping times
0131 %[nu,chi]=atTunesAndDampingRatesFromMat(findm66(atradon(THERING)));
0132 try
0133     [nu,chi]=atTunesAndDampingRatesFromMat(findm66((THERING)));
0134 catch exc
0135     warning('failed coupled damping times computation');
0136     nu=[NaN,NaN,NaN];
0137     chi=[NaN,NaN,NaN];
0138 end
0139 
0140 rp.coupleddampingtime=T0./chi;
0141 
0142 rp.dampingJ = [Jx,Jy,Je];
0143 
0144 %rp.tw = tw;
0145 %rp.tmptw = tmptw;
0146 rp.tunes = tune;
0147 rp.chroms = chrom;
0148 rp.etac = 1/gamma^2-alphac;
0149 
0150 cavind = findcells(THERING,'HarmNumber');
0151 if ~isempty(cavind)
0152     freq_rf = THERING{cavind(1)}.Frequency;
0153     harm = THERING{cavind(1)}.HarmNumber;
0154     Vrf = 0;
0155     for ii=1:length(cavind)
0156         Vrf = Vrf+THERING{cavind(ii)}.Voltage;
0157     end
0158 else
0159     % Default
0160     fprintf('warning: SPEAR3 rf parameters are assume\n');
0161     freq_rf = 476.314e6;
0162     harm = 372;
0163     Vrf = 3.2e6;
0164 end
0165 
0166 phi_s = pi-asin(rp.U0/Vrf);
0167 nus = sqrt(harm*Vrf*abs(rp.etac*cos(phi_s))/2/pi/rp.E0);
0168 rp.nus = nus;
0169 rp.phi_s = phi_s;
0170 rp.harm = harm;
0171 rp.bunchlength = rp.sigma_E*rp.harm*abs(rp.etac)/rp.nus/2/pi/freq_rf*cspeed; % rms bunchlength in meter
0172 delta_max = sqrt(2*U0/pi/alphac/harm/rp.E0)*sqrt( sqrt((Vrf/U0).^2-1) - acos(U0./Vrf));
0173 rp.delta_max = delta_max;
0174 
0175 %calculate vertical emittance
0176 %1. contribution of vertical dispersion
0177 curVavg1 = 1./beta(:,2).*(Dy.^2 + (beta(:,2).*Dyp + alpha(:,2).*Dy).^2);
0178 curVavg = sum(curVavg1.*lendp./abs(rho))/sum(lendp);
0179 emitty_d = Cq*gamma^2*curVavg/Jy; %m-rad
0180 
0181 % %2. contribution of linear coupling resonance
0182 % [G,Delta] = calc_lcG(THERING);
0183 % %emitty_c = emittx*abs(G)^2/(Delta^2+abs(G)^2);
0184 % emitty_c = emittx*abs(G)^2/Delta^2/2.0;
0185 % rp.emitty_c = emitty_c;
0186 
0187 rp.emitty_d = emitty_d;
0188 % rp.emitty = emitty_d + emitty_c;
0189 
0190 if nargout == 0
0191     fprintf('\n');
0192     fprintf('   ******** AT Ring Parmater Summary ********\n');
0193     fprintf('   Energy: \t\t\t%4.5f [GeV]\n', rp.E0/1E9);
0194     fprintf('   Circumference: \t\t%4.5f [m]\n', rp.R*2*pi);
0195     fprintf('   Revolution time: \t\t%4.5f [ns] (%4.5f [MHz]) \n', rp.T0*1e9,1./rp.T0*1e-6);
0196     fprintf('   Betatron tune H: \t\t%4.5f (%4.5f [kHz])\n', rp.tunes(1),(rp.tunes(1)-floor(rp.tunes(1)))/rp.T0*1e-3);
0197     fprintf('                 V: \t\t%4.5f (%4.5f [kHz])\n', rp.tunes(2),(rp.tunes(2)-floor(rp.tunes(2)))/rp.T0*1e-3);
0198     fprintf('   Momentum Compaction Factor: \t%4.5f\n', rp.alphac);
0199     fprintf('   Chromaticity H: \t\t%+4.5f\n', rp.chroms(1));
0200     fprintf('                V: \t\t%+4.5f\n', rp.chroms(2));
0201     fprintf('   Synchrotron Integral 1: \t%4.5f [m]\n', rp.integrals(1));
0202     fprintf('                        2: \t%4.5f [m^-1]\n', rp.integrals(2));
0203     fprintf('                        3: \t%4.5f [m^-2]\n', rp.integrals(3));
0204     fprintf('                        4: \t%4.5f [m^-1]\n', rp.integrals(4));
0205     fprintf('                        5: \t%4.5f [m^-1]\n', rp.integrals(5));
0206     fprintf('   Damping Partition H: \t%4.5f\n', rp.dampingJ(1));
0207     fprintf('                     V: \t%4.5f\n', rp.dampingJ(2));
0208     fprintf('                     E: \t%4.5f\n', rp.dampingJ(3));
0209     fprintf('   Radiation Loss: \t\t%4.5f [keV]\n', rp.U0/1000.);
0210     fprintf('   Natural Energy Spread: \t%4.5e\n', rp.sigma_E);
0211     fprintf('   Natural Emittance: \t\t%4.5e [nm]\n', rp.emittx*1e9);
0212     fprintf('   Radiation Damping H: \t%4.5f [ms]\n', rp.dampingtime(1)*1e3);
0213     fprintf('                     V: \t%4.5f [ms]\n', rp.dampingtime(2)*1e3);
0214     fprintf('                     E: \t%4.5f [ms]\n', rp.dampingtime(3)*1e3);
0215     fprintf('   Slip factor : \t%4.5f\n', rp.etac);
0216     fprintf('\n');
0217     fprintf('   Assuming cavities Voltage: %4.5f [kV]\n', Vrf/1e3);
0218     fprintf('                   Frequency: %4.5f [MHz]\n', freq_rf/1e6);
0219     fprintf('             Harmonic Number: %5d\n', rp.harm);
0220     fprintf('   Synchronous Phase:  %4.5f [rad] (%4.5f [deg])\n', rp.phi_s, rp.phi_s*180/pi);
0221     fprintf('   Linear Energy Acceptance:  %4.5f %%\n', rp.delta_max*100);
0222     fprintf('   Synchrotron Tune:   %4.5f (%4.5f kHz or %4.2f turns) \n', rp.nus, rp.nus/rp.T0*1e-3, 1/rp.nus);
0223     fprintf('   Bunch Length:       %4.5f [mm], %4.5f [ps]\n', rp.bunchlength*1e3, rp.bunchlength/cspeed*1e12);
0224     fprintf('\n');
0225 %     fprintf('   Vertical Emittance:  %4.5f [nm]\n', rp.emitty*1e9);
0226 %     fprintf('   Emitty from Dy:  %4.5f [nm], from linear coupling: %4.5f\n', rp.emitty_d*1e9,rp.emitty_c*1e9);
0227     fprintf('   Emitty from Dy:  %4.5f [nm]\n', rp.emitty_d*1e9);
0228     fprintf('   Emitty 1/gamma cone limit:  %4.5f [pm]\n', rp.emitty_lim*1e12);
0229 end
0230 
0231 
0232 function [dI1,dI2,dI3,dI4,dI5,curHavg] = calcRadInt(rho,theta, a0,b0,D0,D0p,K1,th1,th2)
0233 %[dI1,dI2,dI3,dI4,dI5,curHavg] = calcRadInt(rho,theta, a0,b0,D0,D0p,K1)
0234 %calcuate the contribution to the radiation integrals of a dipole.
0235 %rho, theta, radius and angle of the dipole
0236 %a0, b0, are horizontal alpha and beta at the entrance of the dipole,
0237 %D0, D0p are dispersion at the entrace of the dipole
0238 %K1, the gradient parameter in AT's convention, i.e., positive for
0239 %horizontal focusing, K1=0 by default
0240 %th1, th2, the entrance and exit angle, respectively, th1=th2=theta/2 by
0241 %default.
0242 %
0243 
0244 if nargin==6
0245    K1=0; 
0246 end
0247 if nargin<9
0248    th1 = 0; %theta/2.0;
0249    th2 = 0; %theta/2.0;
0250 end
0251 
0252 M21 = tan(th1)/rho;
0253 D0p = M21*D0+D0p;
0254 a0 = -M21*b0+a0;
0255 
0256 N = 100;
0257 th = (0:N)/N*theta;
0258 len = length(th);
0259 Dx = zeros(len,1); Dxp = zeros(len,1); curHavg1 = zeros(len,1);
0260 for ii=1:len
0261        [Dx(ii), Dxp(ii)] = calcdisp(rho, th(ii), D0, D0p, K1);
0262        [ax, bx] = calctwiss(rho, th(ii), a0, b0, K1);
0263        curHavg1(ii) = (Dx(ii)^2+(ax*Dx(ii)+bx*Dxp(ii))^2)/bx;
0264 end
0265 curHavg = ((curHavg1(1)+curHavg1(end))/2.0 + sum(curHavg1(2:end-1)))/(length(th)-1);
0266 
0267 dI1 = ((Dx(1) + Dx(end))/2.0 + sum(Dx(2:end-1)))*theta/N;
0268 dI2 = abs(theta/rho);
0269 dI3 = abs(theta/rho^2);
0270 dI4 = (1/rho^2 + 2*K1)*dI1  - (Dx(1)/rho^2*tan(th1) + Dx(end)/rho^2*tan(th2));
0271 dI5 = curHavg*abs(theta/rho^2);
0272 
0273 function [Dx, Dxp] = calcdisp(rho, theta, D0, D0p, K1)
0274 %calcualte dispersion function inside a combined-function dipole
0275 s = rho*theta;
0276 if K1>-1/rho^2 %horizontal focusing
0277     sqK = sqrt(1/rho^2+K1);
0278     Dx =  D0*cos(sqK*s) + D0p/sqK*sin(sqK*s)+(1-cos(sqK*s))/rho/sqK^2;
0279     Dxp = -D0*sqK*sin(sqK*s)+D0p*cos(sqK*s)+sin(sqK*s)/rho/sqK;
0280 else %horizontal defocusing
0281     sqK = sqrt(-(1/rho^2+K1));
0282     Dx =  D0*cosh(sqK*s) + D0p/sqK*sinh(sqK*s)+(-1+cosh(sqK*s))/rho/sqK^2;
0283     Dxp = D0*sqK*sinh(sqK*s)+D0p*cosh(sqK*s)+sinh(sqK*s)/rho/sqK;
0284 
0285 end
0286 
0287 function [ax, bx] = calctwiss(rho, theta, a0, b0, K1)
0288 %calculate twiss function inside a combined-function dipole manget
0289 Mx = calcMx(rho, K1,theta);
0290 g0 = (1+a0^2)/b0;
0291 twx2 = [Mx(1,1)^2, -2*Mx(1,1)*Mx(1,2), Mx(1,2)^2; 
0292     -Mx(1,1)*Mx(2,1), Mx(1,1)*Mx(2,2)+Mx(1,2)*Mx(2,1),-Mx(1,2)*Mx(2,2);
0293     Mx(2,1)^2, -2*Mx(2,1)*Mx(2,2),Mx(2,2)^2] * [b0, a0, g0]';
0294 ax = twx2(2);
0295 bx = twx2(1);
0296 
0297 function Mx = calcMx(rho,K1,theta)
0298 s = rho*theta;
0299 if K1>-1/rho^2 %horizontal focusing
0300     sqK = sqrt(1/rho^2+K1);
0301     Mx = [cos(sqK*s), sin(sqK*s)/sqK; -sqK*sin(sqK*s), cos(sqK*s)];
0302 else %horizontal defocusing
0303     sqK = sqrt(-(1/rho^2+K1));
0304     Mx = [cosh(sqK*s), sinh(sqK*s)/sqK; sqK*sinh(sqK*s), cosh(sqK*s)];
0305 end
0306