1 | function rp = ringpara(THERING,varargin) |
---|
2 | %rp = ringpara, use global THERING |
---|
3 | %rp = ringpara(THERING) |
---|
4 | %rp = ringpara(THERING,U0), supply total radiation loss in MeV |
---|
5 | %calculate various ring parameters |
---|
6 | %(1) The calculation of emittance, mcf, momentum spread, bunch length, damping time, etc |
---|
7 | %is more accurate than atsummary.m because detailed |
---|
8 | %calculation of dispersion function and curly H function inside dipoles is performed. |
---|
9 | %(2) calculate contribution of dispersion to vertical emittance. |
---|
10 | % |
---|
11 | %Author: Xiaobiao Huang |
---|
12 | %created on 12/17/2007 |
---|
13 | %Part of this code was modified from atsummary.m |
---|
14 | % |
---|
15 | |
---|
16 | if nargin==0 |
---|
17 | global THERING; |
---|
18 | end |
---|
19 | Cq = 3.83E-13; |
---|
20 | a = findcells(THERING,'Energy'); |
---|
21 | if isempty(a); |
---|
22 | gamma = 3000/.511; |
---|
23 | else |
---|
24 | gamma = THERING{a(1)}.Energy/.511E6; |
---|
25 | end |
---|
26 | |
---|
27 | % EntranceAngle instead of BendingAngle otherwise a multipole is taken as a |
---|
28 | % dipole |
---|
29 | %dpindex = findcells(THERING,'BendingAngle'); |
---|
30 | dpindex = findcells(THERING,'EntranceAngle'); |
---|
31 | [tw,tune,chrom] = twissring(THERING,0,dpindex,'chrom',0.00001); |
---|
32 | beta = cat(1, tw.beta); |
---|
33 | alpha = cat(1, tw.alpha); |
---|
34 | mu = cat(1, tw.mu); |
---|
35 | disper = cat(1, tw.Dispersion); |
---|
36 | Dx = disper(1:4:end); |
---|
37 | Dxp = disper(2:4:end); |
---|
38 | Dy = disper(3:4:end); |
---|
39 | Dyp = disper(4:4:end); |
---|
40 | |
---|
41 | [tmptw,tune,chrom] = twissring(THERING,0,1:length(THERING),'chrom',0.00001); |
---|
42 | |
---|
43 | lendp = getcellstruct(THERING,'Length',dpindex); %bending magnet length |
---|
44 | |
---|
45 | theta = getcellstruct(THERING,'BendingAngle',dpindex); %bending angle |
---|
46 | rho = lendp./theta; %THERING{dpindex(1)}.Length/THERING{dpindex(1)}.BendingAngle; |
---|
47 | |
---|
48 | I1 = 0; |
---|
49 | I2 = 0; |
---|
50 | I3 = 0; |
---|
51 | I4 = 0; |
---|
52 | I5 = 0; |
---|
53 | for ii=1:length(dpindex) |
---|
54 | if isfield(THERING{dpindex(ii)},'K') |
---|
55 | K = THERING{dpindex(ii)}.K; |
---|
56 | else |
---|
57 | K = 0; |
---|
58 | end |
---|
59 | th1 = THERING{dpindex(ii)}.EntranceAngle; |
---|
60 | th2 = THERING{dpindex(ii)}.ExitAngle; |
---|
61 | [dI1,dI2,dI3,dI4,dI5,curHavg1(ii)] = calcRadInt(rho(ii),theta(ii), alpha(ii,1),beta(ii,1),Dx(ii),Dxp(ii),K,th1,th2); |
---|
62 | I1 = I1 + dI1; |
---|
63 | I2 = I2 + dI2; |
---|
64 | I3 = I3 + dI3; |
---|
65 | I4 = I4 + dI4; |
---|
66 | I5 = I5 + dI5; |
---|
67 | end |
---|
68 | % curHavg = sum(curHavg1.*lendp./abs(rho))/sum(lendp); |
---|
69 | % %emittx = Cq*gamma^2*curHavg/Jx/rho*1e9; %nm-rad |
---|
70 | % emittx = Cq*gamma^2*curHavg/Jx*1e9; %nm-rad |
---|
71 | R = findspos(THERING, length(THERING)+1)/2/pi; |
---|
72 | alphac = I1/2/pi/R; |
---|
73 | U0 = 14.085*(gamma*.511/1000)^4*I2*1000.; %eV |
---|
74 | if nargin>=2 |
---|
75 | fprintf('dipole radiation loss: %4.5f keV\n', U0/1000.); |
---|
76 | U0 = varargin{1}*1e6; %convert MeV to eV |
---|
77 | end |
---|
78 | sigma_E = gamma*sqrt(Cq*I3/(2*I2+I4)); |
---|
79 | Jx = 1-I4/I2; |
---|
80 | Jy = 1.00; |
---|
81 | Je = 2+I4/I2; |
---|
82 | emittx = Cq*gamma^2*I5/(I2-I4); |
---|
83 | |
---|
84 | cspeed = 2.99792458e8; %m/s |
---|
85 | T0 = 2*pi*R/cspeed; |
---|
86 | alpha0 = U0/1.0e6/2/T0/(gamma*.511); |
---|
87 | alphax = Jx*alpha0; %horizontal damping rate, 1/s |
---|
88 | alphay = Jy*alpha0; |
---|
89 | alphaE = Je*alpha0; |
---|
90 | |
---|
91 | rp.E0 = gamma*0.511E6; |
---|
92 | rp.R = R; |
---|
93 | rp.alphac = alphac; |
---|
94 | rp.U0 = U0; %eV |
---|
95 | rp.sigma_E = sigma_E; |
---|
96 | rp.emittx = emittx; |
---|
97 | rp.T0 = T0; |
---|
98 | rp.integrals = [I1,I2,I3,I4,I5]; |
---|
99 | rp.dampingalpha = [alphax, alphay, alphaE]; |
---|
100 | rp.dampingtime = 1./[alphax, alphay, alphaE]; |
---|
101 | rp.dampingJ = [Jx,Jy,Je]; |
---|
102 | |
---|
103 | rp.tunes = tune; |
---|
104 | rp.chroms = chrom; |
---|
105 | rp.etac = 1/gamma^2-alphac; |
---|
106 | |
---|
107 | cavind = findcells(THERING,'HarmNumber'); |
---|
108 | if ~isempty(cavind) |
---|
109 | freq_rf = THERING{cavind(1)}.Frequency; |
---|
110 | harm = THERING{cavind(1)}.HarmNumber; |
---|
111 | Vrf = 0; |
---|
112 | for ii=1:length(cavind) |
---|
113 | Vrf = Vrf+THERING{cavind(ii)}.Voltage; |
---|
114 | end |
---|
115 | else |
---|
116 | % Default |
---|
117 | fprintf('warning: SPEAR3 rf parameters are assume\n'); |
---|
118 | freq_rf = 476.314e6; |
---|
119 | harm = 372; |
---|
120 | Vrf = 3.2e6; |
---|
121 | end |
---|
122 | |
---|
123 | phi_s = pi-asin(rp.U0/Vrf); |
---|
124 | nus = sqrt(harm*Vrf*abs(rp.etac*cos(phi_s))/2/pi/rp.E0); |
---|
125 | rp.nus = nus; |
---|
126 | rp.phi_s = phi_s; |
---|
127 | rp.harm = harm; |
---|
128 | rp.bunchlength = rp.sigma_E*rp.harm*abs(rp.etac)/rp.nus/2/pi/freq_rf*cspeed; % rms bunchlength in meter |
---|
129 | delta_max = sqrt(2*U0/pi/alphac/harm/rp.E0)*sqrt( sqrt((Vrf/U0).^2-1) - acos(U0./Vrf)); |
---|
130 | rp.delta_max = delta_max; |
---|
131 | |
---|
132 | %calculate vertical emittance |
---|
133 | %1. contribution of vertical dispersion |
---|
134 | curVavg1 = 1./beta(:,2).*(Dy.^2 + (beta(:,2).*Dyp + alpha(:,2).*Dy).^2); |
---|
135 | curVavg = sum(curVavg1.*lendp./abs(rho))/sum(lendp); |
---|
136 | emitty_d = Cq*gamma^2*curVavg/Jy; %m-rad |
---|
137 | |
---|
138 | % %2. contribution of linear coupling resonance |
---|
139 | % [G,Delta] = calc_lcG(THERING); |
---|
140 | % %emitty_c = emittx*abs(G)^2/(Delta^2+abs(G)^2); |
---|
141 | % emitty_c = emittx*abs(G)^2/Delta^2/2.0; |
---|
142 | % rp.emitty_c = emitty_c; |
---|
143 | |
---|
144 | rp.emitty_d = emitty_d; |
---|
145 | % rp.emitty = emitty_d + emitty_c; |
---|
146 | |
---|
147 | if nargout == 0 |
---|
148 | fprintf('\n'); |
---|
149 | fprintf(' ******** AT Ring Parmater Summary ********\n'); |
---|
150 | fprintf(' Energy: \t\t\t%4.5f [GeV]\n', rp.E0/1E9); |
---|
151 | fprintf(' Circumference: \t\t%4.5f [m]\n', rp.R*2*pi); |
---|
152 | fprintf(' Revolution time: \t\t%4.5f [ns] (%4.5f [MHz]) \n', rp.T0*1e9,1./rp.T0*1e-6); |
---|
153 | 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); |
---|
154 | fprintf(' V: \t\t%4.5f (%4.5f [kHz])\n', rp.tunes(2),(rp.tunes(2)-floor(rp.tunes(2)))/rp.T0*1e-3); |
---|
155 | fprintf(' Momentum Compaction Factor: \t%4.5f\n', rp.alphac); |
---|
156 | fprintf(' Chromaticity H: \t\t%+4.5f\n', rp.chroms(1)); |
---|
157 | fprintf(' V: \t\t%+4.5f\n', rp.chroms(2)); |
---|
158 | fprintf(' Synchrotron Integral 1: \t%4.5f [m]\n', rp.integrals(1)); |
---|
159 | fprintf(' 2: \t%4.5f [m^-1]\n', rp.integrals(2)); |
---|
160 | fprintf(' 3: \t%4.5f [m^-2]\n', rp.integrals(3)); |
---|
161 | fprintf(' 4: \t%4.5f [m^-1]\n', rp.integrals(4)); |
---|
162 | fprintf(' 5: \t%4.5f [m^-1]\n', rp.integrals(5)); |
---|
163 | fprintf(' Damping Partition H: \t%4.5f\n', rp.dampingJ(1)); |
---|
164 | fprintf(' V: \t%4.5f\n', rp.dampingJ(2)); |
---|
165 | fprintf(' E: \t%4.5f\n', rp.dampingJ(3)); |
---|
166 | fprintf(' Radiation Loss: \t\t%4.5f [keV]\n', rp.U0/1000.); |
---|
167 | fprintf(' Natural Energy Spread: \t%4.5e\n', rp.sigma_E); |
---|
168 | fprintf(' Natural Emittance: \t\t%4.5e [nm]\n', rp.emittx*1e9); |
---|
169 | fprintf(' Radiation Damping H: \t%4.5f [ms]\n', rp.dampingtime(1)*1e3); |
---|
170 | fprintf(' V: \t%4.5f [ms]\n', rp.dampingtime(2)*1e3); |
---|
171 | fprintf(' E: \t%4.5f [ms]\n', rp.dampingtime(3)*1e3); |
---|
172 | fprintf(' Slip factor : \t%4.5f\n', rp.etac); |
---|
173 | fprintf('\n'); |
---|
174 | fprintf(' Assuming cavities Voltage: %4.5f [kV]\n', Vrf/1e3); |
---|
175 | fprintf(' Frequency: %4.5f [MHz]\n', freq_rf/1e6); |
---|
176 | fprintf(' Harmonic Number: %5d\n', rp.harm); |
---|
177 | fprintf(' Synchronous Phase: %4.5f [rad] (%4.5f [deg])\n', rp.phi_s, rp.phi_s*180/pi); |
---|
178 | fprintf(' Linear Energy Acceptance: %4.5f %%\n', rp.delta_max*100); |
---|
179 | fprintf(' Synchrotron Tune: %4.5f (%4.5f kHz or %4.2f turns) \n', rp.nus, rp.nus/rp.T0*1e-3, 1/rp.nus); |
---|
180 | fprintf(' Bunch Length: %4.5f [mm], %4.5f [ps]\n', rp.bunchlength*1e3, rp.bunchlength/cspeed*1e12); |
---|
181 | fprintf('\n'); |
---|
182 | % fprintf(' Vertical Emittance: %4.5f [nm]\n', rp.emitty*1e9); |
---|
183 | % fprintf(' Emitty from Dy: %4.5f [nm], from linear coupling: %4.5f\n', rp.emitty_d*1e9,rp.emitty_c*1e9); |
---|
184 | fprintf(' Emitty from Dy: %4.5f [nm]\n', rp.emitty_d*1e9); |
---|
185 | end |
---|
186 | |
---|
187 | |
---|
188 | function [dI1,dI2,dI3,dI4,dI5,curHavg] = calcRadInt(rho,theta, a0,b0,D0,D0p,K1,th1,th2) |
---|
189 | %[dI1,dI2,dI3,dI4,dI5,curHavg] = calcRadInt(rho,theta, a0,b0,D0,D0p,K1) |
---|
190 | %calcuate the contribution to the radiation integrals of a dipole. |
---|
191 | %rho, theta, radius and angle of the dipole |
---|
192 | %a0, b0, are horizontal alpha and beta at the entrance of the dipole, |
---|
193 | %D0, D0p are dispersion at the entrace of the dipole |
---|
194 | %K1, the gradient parameter in AT's convention, i.e., positive for |
---|
195 | %horizontal focusing, K1=0 by default |
---|
196 | %th1, th2, the entrance and exit angle, respectively, th1=th2=theta/2 by |
---|
197 | %default. |
---|
198 | % |
---|
199 | |
---|
200 | if nargin==6 |
---|
201 | K1=0; |
---|
202 | end |
---|
203 | if nargin<9 |
---|
204 | th1 = theta/2.0; |
---|
205 | th2 = theta/2.0; |
---|
206 | end |
---|
207 | |
---|
208 | N = 100; |
---|
209 | th = (0:N)/N*theta; |
---|
210 | for ii=1:length(th) |
---|
211 | [Dx(ii), Dxp(ii)] = calcdisp(rho, th(ii), D0, D0p, K1); |
---|
212 | [ax, bx] = calctwiss(rho, th(ii), a0, b0, K1); |
---|
213 | curHavg1(ii) = (Dx(ii)^2+(ax*Dx(ii)+bx*Dxp(ii))^2)/bx; |
---|
214 | end |
---|
215 | curHavg = ((curHavg1(1)+curHavg1(end))/2.0 + sum(curHavg1(2:end-1)))/(length(th)-1); |
---|
216 | |
---|
217 | dI1 = ((Dx(1) + Dx(end))/2.0 + sum(Dx(2:end-1)))*theta/N; |
---|
218 | dI2 = abs(theta/rho); |
---|
219 | dI3 = abs(theta/rho^2); |
---|
220 | dI4 = (1/rho^2 + 2*K1)*dI1 + Dx(1)/rho^2*tan(th1) + Dx(end)/rho^2*tan(th2); |
---|
221 | dI5 = curHavg*abs(theta/rho^2); |
---|
222 | |
---|
223 | function [Dx, Dxp] = calcdisp(rho, theta, D0, D0p, K1) |
---|
224 | %calcualte dispersion function inside a combined-function dipole |
---|
225 | s = rho*theta; |
---|
226 | if K1>-1/rho^2; %horizontal focusing |
---|
227 | sqK = sqrt(1/rho^2+K1); |
---|
228 | Dx = D0*cos(sqK*s) + D0p/sqK*sin(sqK*s)+(1-cos(sqK*s))/rho/sqK^2; |
---|
229 | Dxp = -D0*sqK*sin(sqK*s)+D0p*cos(sqK*s)+sin(sqK*s)/rho/sqK; |
---|
230 | else %horizontal defocusing |
---|
231 | sqK = sqrt(-(1/rho^2+K1)); |
---|
232 | Dx = D0*cosh(sqK*s) + D0p/sqK*sinh(sqK*s)+(-1+cosh(sqK*s))/rho/sqK^2; |
---|
233 | Dxp = D0*sqK*sinh(sqK*s)+D0p*cosh(sqK*s)+sinh(sqK*s)/rho/sqK; |
---|
234 | |
---|
235 | end |
---|
236 | |
---|
237 | function [ax, bx] = calctwiss(rho, theta, a0, b0, K1) |
---|
238 | %calculate twiss function inside a combined-function dipole manget |
---|
239 | Mx = calcMx(rho, K1,theta); |
---|
240 | g0 = (1+a0^2)/b0; |
---|
241 | twx2 = [Mx(1,1)^2, -2*Mx(1,1)*Mx(1,2), Mx(1,2)^2; |
---|
242 | -Mx(1,1)*Mx(2,1), Mx(1,1)*Mx(2,2)+Mx(1,2)*Mx(2,1),-Mx(1,2)*Mx(2,2); |
---|
243 | Mx(2,1)^2, -2*Mx(2,1)*Mx(2,2),Mx(2,2)^2] * [b0, a0, g0]'; |
---|
244 | ax = twx2(2); |
---|
245 | bx = twx2(1); |
---|
246 | |
---|
247 | function Mx = calcMx(rho,K1,theta) |
---|
248 | s = rho*theta; |
---|
249 | if K1>-1/rho^2; %horizontal focusing |
---|
250 | sqK = sqrt(1/rho^2+K1); |
---|
251 | Mx = [cos(sqK*s), sin(sqK*s)/sqK; -sqK*sin(sqK*s), cos(sqK*s)]; |
---|
252 | else %horizontal defocusing |
---|
253 | sqK = sqrt(-(1/rho^2+K1)); |
---|
254 | Mx = [cosh(sqK*s), sinh(sqK*s)/sqK; sqK*sinh(sqK*s), cosh(sqK*s)]; |
---|
255 | end |
---|
256 | |
---|