[z lambda sigma mu] = HASSINSKYFIT(SIGMA0, R, L)
This function will generate a charge distribution for a given SIGMA0
(zero current bunch length), q (bunch total charge in Coloumbs and can be
a vector), using a generalised impedance model with R (OHms) and L
(Henry) according to the Haissinski equation.
SIGMA0 (meters, scalar) : zero current bunch length
q (coloumbs, vector): total charge in bunch
R (Ohms, scalar) : ring real effective resistance
L (Henry, scalar) : ring real effecitve impedance
frf (Hz, scalar) : RF frequency
h (integer, scalar) : harmonic number of the ring
V (Volts, scalar) : cavity voltage
phis (radians, scalar) : synchronous phase
Output:
z (meters, vector) : longitudinal position
lambda (normalised, vector) : charge distribution along z.
sigma (meters, vector) : bunch length for each q using fitgaussian.
mu (meters, vector) : mu of gaussian fit from fitgaussian.
WARNING: The gaussian fit uses the following command:
fitgaussian(I2,'plot',0,'bg_grad',0,'DC',0);
and will thus be a sigma in "points/pixels", no background
gradient and zero DC. Will fit the area, sigma, assymetry factor
and the centre of gravity for the gaussian.
See also FITGAUSSIAN
02/03/2010 Eugene0001 function varargout = hassinskyfit(sigma0, q, R, L, frf, h, V, phis, varargin) 0002 0003 % [z lambda sigma mu] = HASSINSKYFIT(SIGMA0, R, L) 0004 % This function will generate a charge distribution for a given SIGMA0 0005 % (zero current bunch length), q (bunch total charge in Coloumbs and can be 0006 % a vector), using a generalised impedance model with R (OHms) and L 0007 % (Henry) according to the Haissinski equation. 0008 % 0009 % SIGMA0 (meters, scalar) : zero current bunch length 0010 % q (coloumbs, vector): total charge in bunch 0011 % R (Ohms, scalar) : ring real effective resistance 0012 % L (Henry, scalar) : ring real effecitve impedance 0013 % frf (Hz, scalar) : RF frequency 0014 % h (integer, scalar) : harmonic number of the ring 0015 % V (Volts, scalar) : cavity voltage 0016 % phis (radians, scalar) : synchronous phase 0017 % 0018 % Output: 0019 % 0020 % z (meters, vector) : longitudinal position 0021 % lambda (normalised, vector) : charge distribution along z. 0022 % sigma (meters, vector) : bunch length for each q using fitgaussian. 0023 % mu (meters, vector) : mu of gaussian fit from fitgaussian. 0024 % 0025 % WARNING: The gaussian fit uses the following command: 0026 % fitgaussian(I2,'plot',0,'bg_grad',0,'DC',0); 0027 % and will thus be a sigma in "points/pixels", no background 0028 % gradient and zero DC. Will fit the area, sigma, assymetry factor 0029 % and the centre of gravity for the gaussian. 0030 % 0031 % See also FITGAUSSIAN 0032 % 0033 % 02/03/2010 Eugene 0034 0035 DEBUG = 0; 0036 TOL = 1e-10; 0037 MAXIT = 500; 0038 0039 if DEBUG 0040 % sigma0 = 20.3e-12*3e8; %natural bunch length in m (145.5/5.6)*e-12*3e8 0041 % R=1600; %resistance (in Ohm) 0042 % L=25e-9; %Inductance (in nanoHenry) 0043 % frf = 499.672e6; %RF frequency 0044 % V=3000000; %RF voltage 0045 % q = ([1:10]/1000)/1.388e6; 0046 % h = 360; 0047 % phis = 2.83964; 0048 end 0049 0050 c = 2.998e8; 0051 dz = sigma0/50; 0052 z = [-sigma0*6:dz:sigma0*6]; 0053 0054 % initialise 0055 sigmas = zeros(size(q)); 0056 mus = zeros(size(q)); 0057 lambda = zeros(length(z),length(q)); 0058 0059 for i=1:length(q) 0060 mu = 0; 0061 K = -c^2*q(i)/(2*pi*frf*V*cos(phis)*(sigma0^2)); 0062 0063 if 0 0064 % test solution to Hassinski equation based on Oide and Yokoya. 0065 % The method has a different method of finding the solution to the 0066 % Haissinski equation where you integrate from "infinity" at the 0067 % head the bunch towards the tail. This avoids the a fixed point 0068 % type algorithm to find the equilibrium solution. This bit has yet 0069 % to be tested and debugged fully 0070 % Eugene 2/3/2010 0071 warning('This algorithm not tested. Use with caution.'); 0072 I2 = zeros(size(z)); 0073 A = 1/(sigma0*sqrt(2*pi)); % Initial normalisation factor 0074 sumint = 0; 0075 while abs(sumint-1) > 0.0001 0076 for j=(length(z)-1):-1:1 0077 % Assume that the solution for I2 is zero for suffidciently 0078 % large z. 0079 intR = R*sum(I2(j:end))*dz; 0080 intL = c*L*sum(gradient(I2(j:end),dz))*dz; 0081 I2(j) = A*(exp(-(z(j)-mu).^2./(2*sigma0^2))).*(exp((-K.*(intR + intL)))); 0082 end 0083 sumint = (sum(I2)*dz); 0084 A = A/sumint; 0085 end 0086 converge = 1; 0087 else 0088 % Fixted point algorithm based on Rohan's old code that has been 0089 % put into functinal form. 0090 if i==1 0091 % initial distribution, subsequent starting distributions set at 0092 % the bottom of the loop. 0093 % kl = R/2/sqrt(pi)/sigma0; 0094 mu = 0; %c*kl*q(i)/(2*pi*frf*V*cos(phis)); 0095 I = (1/sqrt(2*pi)/sigma0)*exp(-(z-mu).^2/(2*sigma0^2)); 0096 I = I./(sum(I)*dz); % Normalise 0097 Ip = -((z-mu)/sigma0^2).*I; 0098 else 0099 I = I_stable; 0100 Ip = Ip_stable; 0101 end 0102 0103 converge = 0; 0104 itno = 1; 0105 while ~converge 0106 intR = R*cumsum(I)*dz; 0107 intL = c*L*cumsum(Ip)*dz; 0108 I2 = (1/(sigma0*sqrt(2*pi)))*(exp(-(z-mu).^2./(2*sigma0^2))).*(exp((-K.*(intR + intL)))); 0109 0110 %Renormalise the function 0111 I2 = I2./(sum(I2)*dz); 0112 if DEBUG 0113 fprintf('Iteration: %03d, STD: %e\n',itno,std(I2(:)-I(:))); 0114 end 0115 stderror = std(I2(:)-I(:)); 0116 if stderror < TOL; %convegence condition 0117 converge=1; 0118 % plot(I2,'r'); %plot final distribution in red 0119 else 0120 I = I2; 0121 Ip = gradient(I2,dz); 0122 end 0123 itno = itno + 1; 0124 0125 % if DEBUG && rem(itno,10) == 0 0126 % figure(98); 0127 % plot(I2); 0128 % end 0129 0130 if itno > MAXIT 0131 fprintf('WARNING: Did not converge during fit of Hassinsky''s equation!!\n'); 0132 fprintf('q: %g, std error: %g\n',q(i),stderror); 0133 figure(99); 0134 plot(z,I); 0135 break; 0136 end 0137 end 0138 end 0139 0140 lambda(:,i) = I2; 0141 0142 if nargout > 2 0143 if converge 0144 % units of meters 0145 temp = fitgaussian(I2,'plot',1,'bg_grad',0,'DC',0,'scale',dz); 0146 sigmas(i) = temp.sigma; 0147 mus(i) = temp.mu - z(end); 0148 0149 % Set I and Ip for the next sequence and only keep "stable" 0150 % solutions that converge. Otherwise small instabilities will 0151 % build up too quickly for subsequent solutions with higher 0152 % currents, q. 0153 % kl = R/2/sqrt(pi)/sigmas(i-1); 0154 mu = 0; %c*kl*q(i)/(2*pi*frf*V*cos(phis)); 0155 I_stable = I2; 0156 Ip_stable = gradient(I2,dz); 0157 else 0158 % Will try subset (plot envelope). There are "solutions" where 0159 % a sinusoidal pattern is observed at the head of the bunch. 0160 % This part of the code attempts to plot an envelope gaussian 0161 % function by looking for the peaks of the oscilating function 0162 % and only fitting to those points. 0163 ii = find(gradient((lambda(:,i))) > 0); 0164 ii = ii(abs(diff(ii)) > 1); 0165 if ~isempty(ii) && length(ii) > 1 0166 % it usually can't find where the gaussian tail of the 0167 % bunch starts. So assume the distance from the last peak 0168 % to the start of the tail of the gaussian bunch is the 0169 % same as the distance between the head of the gaussian 0170 % part of the bunch and the first peak. 0171 ii(end+1) = ii(end) + (ii(2)-ii(1))*2; % Generate last point as it can't seem to find the last turn around point. 0172 % datai = [1:ii(1), ii(2:end-1)', ii(end):length(z)]; 0173 datai = [1:ii(1), ii(end):length(z)]; 0174 else 0175 % sometimes convertance isn't met however the nongaussian 0176 % element is not strong so a standard gaussian fit is 0177 % sufficient. 0178 datai = [1:length(z)]'; 0179 end 0180 0181 z_subset = z(datai); 0182 I2_subset = lambda(datai,i); 0183 0184 temp = fitgaussian(I2_subset,'plot',0,'bg_grad',0,'DC',0,'scale',1,'x',z_subset); 0185 % pause; 0186 % must have reached iteration limit without convergence 0187 % therefore cannot trust fitdata 0188 sigmas(i) = temp.sigma; 0189 mus(i) = temp.mu - z(end); 0190 end 0191 else 0192 I_stable = I2; 0193 Ip_stable = gradient(I2,dz); 0194 end 0195 end 0196 0197 varargout{1} = z; 0198 varargout{2} = lambda; 0199 varargout{3} = sigmas; 0200 varargout{4} = mus; 0201 0202 0203 return 0204 0205 0206 0207 % Loop over currents 0208 current = []; 0209 sigmas = []; 0210 I2 = []; 0211 initstep = 0.5; 0212 itno = 0; 0213 for loop = 1:endloop 0214 if itno >= 600 0215 initstep = initstep/2; 0216 end 0217 % in mA 0218 if loop == 1 0219 current(loop) = 0.001; 0220 else 0221 current(loop) = current(loop-1) + initstep; 0222 end 0223 %current(loop) = initstep*(loop-1) +0.001; 0224 %I = zeros(1001,1); %initial charge distribution 0225 %I2 =zeros(1001,1); %final charge distribution 0226 %Ip = zeros(1001,1); %deriviative of charge distribution 0227 c=3e8; 0228 Frf = 499.672e6; %RF frequency 0229 V=3000000; %RF voltage 0230 Q=(current(end)/1000)/1.3880e6; %convert current to charge 0231 0232 %I2p/I2 = -z/sigma^2 + (c^2*Q*(R*I + c*L*Ip))/(2*pi*Frf*V*cosd(158)*sigma^2) 0233 0234 stop = 0; 0235 itno=0; 0236 %iterate! 0237 SIGMAx = sigma; 0238 K = -c^2*Q/(2*pi*Frf*V*cosd(158)*(SIGMAx^2)); 0239 %set the initial distribution as a gaussian 0240 if loop==1 0241 Ax = 1; 0242 MUx = 0; %centre it 0243 dz = SIGMAx/50; 0244 z = [-SIGMAx*6:dz:SIGMAx*6]; 0245 I = 1 * (1/(SIGMAx*sqrt(2*pi)))*exp(-(((z-MUx).^2)/(2*SIGMAx^2))); 0246 % Normalise the initial gaussian 0247 I = I./(sum(I)*dz); 0248 % figure(100); hold on; plot(I); %plot initial distribution 0249 % make derivative disribution 'Ip' 0250 Ip = -((z-MUx)/SIGMAx^2).*I; 0251 else 0252 I = gaussianY; 0253 Ip = gradient(I,dz); 0254 end 0255 0256 while stop==0; 0257 itno=itno+1; 0258 0259 for a=1:length(z) %do the I(z'-z)dz itergrations first 0260 intR = sum(R*I(1:a))*dz; 0261 intL = sum(c*L*Ip(1:a))*dz; 0262 %the Hassinski equation: 0263 I2(a) = (1/(SIGMAx*sqrt(2*pi)))*(exp(-(z(a)-MUx)^2/(2*SIGMAx^2)))*(exp((-K*(intR + intL)))); 0264 end 0265 %Renormalise the function 0266 I2 = I2./(sum(I2)*dz); 0267 0268 %figure(300); 0269 %plot(I2,'g'); %plot iterations in green 0270 %hold on; 0271 0272 fprintf('Iteration: %03d, STD: %e\n',itno,std(I2(:)-I(:))); 0273 if std(I2(:)-I(:)) < 1e-13; %convegence condition 0274 stop=1; 0275 % plot(I2,'r'); %plot final distribution in red 0276 itno; %print out the iteration number 0277 else 0278 I = I2; 0279 Ip = gradient(I2,dz); 0280 % Ip(1)=0; 0281 % for a=1:1001 0282 % I(a) = I2(a); 0283 % if a>1&&a<1001; %take trapezoidal derivative 0284 % Ip(a) = 0.5*((I2(a)-I2(a-1)) + (I2(a+1)-I2(a)) ); 0285 % else if a==1 %cant do this on the edges 0286 % Ip(a) = I2(a+1) - I2(a); 0287 % else 0288 % Ip(a) = I2(a) - I2(a-1); 0289 % end 0290 % end 0291 % end 0292 end 0293 if itno==500; %stop after 100 iterations. 0294 stop =1; 0295 end 0296 end 0297 % Fit a gaussian 0298 temp = fitgaussian(I2,'plot',0,'bg_grad',0,'DC',0); 0299 %gaussianY = ones(sizey,1); 0300 xx = [1:length(I2)]; 0301 gaussianY = temp.area * exp(-0.5*((xx-temp.mu)./((1+sign(xx-temp.mu)*temp.Assym_factor)*temp.sigma)).^2) / sqrt(2*pi*temp.sigma^2) + temp.bg_gradient*xx + temp.DC; 0302 sigmas(loop) = temp.sigma/c*1e12; % convert into ps 0303 %means(i) = Estimates(2); 0304 0305 figure(201); 0306 clf; 0307 plot(z,I2,'b-'); hold on; 0308 plot(z,gaussianY,'-r'); 0309 0310 end % for each current 0311 save([pref '_hass'],'fit*','cal*','ave*','measured*','current','sigmas'); 0312 0313 measuredsigma = avesigmas; 0314 vectormeas = [avecurrents; measuredsigma]; 0315 vectorsim = [current; sigmas]; 0316 vectorerror = [avecurrents; measurederrorsigma]; 0317 0318 figure(950); xlabel('Current (mA)'); ylabel('Bunch Sigma (arb. units)');hold on; 0319 plot(vectorsim(1,:),vectorsim(2,:),color); 0320 %errorbar(vectormeas(1,17:19),vectormeas(2,17:19),vectorerror(2,17:19),'.'); 0321 errorbar(vectormeas(1,:),vectormeas(2,:),vectorerror(2,:),['.',color]); 0322 %errorbar(vectormeas(1,17:19),vectormeas(2,17:19),vectorerror(2,17:19),'.'); 0323