function varargout = Emit_Parabola_Fit_kl_XProjection(x,erx,y,ery,me,En,S12,erS12,NumPoints) %This function will fit parabole for given rms values, plot them and from %parabola parameter extract emitance. % Written by Sadiq Setiniyaz, Jan 2010 %x %y %ery order = 2; % second order fit ParNum = order + 1; % finding row matrix beta for k=1:ParNum for i=1:NumPoints f(k)=power(x(i),k-1); %fprintf(' i=%d f(%d)=%g \t',i,k,f(k)); beta_element(k,i)=y(i)*f(k)/(ery(i)*ery(i)); %fprintf('y(i=%d)=%g ery(i=%d)=%g \n',i,y(i),i,ery(i)); %fprintf('beta(k=%d,i=%d)=%g \n',k,i,beta_element(k,i)); beta(k)=sum(beta_element(k,:)); %fprintf('beta_element(%d,%d)=%g %beta(%d)=%g\n',k,i,beta_element(k,i),k,beta(k)); end end %beta % finding matrix alpha for k=1:ParNum for l=1:ParNum for i=1:NumPoints f(k)=power(x(i),k-1); f(l)=power(x(i),l-1); alpha_element(k,l,i)=f(k)*f(l)/(ery(i)*ery(i)); alpha(k,l)=sum(alpha_element(k,l,:)); %fprintf('alpha_element(%d,%d,%d)=%g alpha(%d,%d)=%g\n',k,l,i,alpha_element(k,l,i),k,l,alpha(k,l)); end end end %alpha alpha_invert = inv(alpha); %to find fit parameters, need to invert matrix alpha. %alpha_invert % the fit is y = a(1) + a(2)*x + a(3)*x*x %fit parameters by matrix inversion method parameter = beta*alpha_invert; %a(l)=parameter(l) %parameters fit_m_inver = parameter(1) + parameter(2)*x + parameter(3)*x.*x; % error on the parameters er_a = sqrt(diag(alpha_invert)); %fit parameters by MatLab par = polyfit(x,y,2); %par fit_MatLab = par(3) + par(2)*x + par(1)*x.*x; %Get Emittance from fit gamma = En/me; A=parameter(3); B=-parameter(2)/(2*A); C=parameter(1)-A*B^2; erA = er_a(3); erB = sqrt( (-er_a(2)/(2*A))^2 + (parameter(2)*erA/(2*A^2))^2 ); erC = sqrt((er_a(1))^2 + (-B^2*erA)^2 + (-2*A*B*erB)^2); %fprintf('A=%g+-%g B=%g+-%g C=%g+-%g\n',A,erA,B,erB,C,erC); emit = sqrt(A*C)/((S12)^2); er_emit = sqrt( (1/2*A^(-1/2)*A^(1/2)*erA/(S12)^2)^2 + (1/2*A^(1/2)*A^(-1/2)*erC/(S12)^2)^2 + ( 2*sqrt(A*C)*erS12/(S12)^3)^2 ); emit_n = emit*gamma; er_emit_n = er_emit*gamma; beta_func = sqrt(A/C); er_beta_func = sqrt( (1/2*sqrt(1/(A*C))*erA)^2 + (-1/2*sqrt(A/C^3)*erC)^2 ); alpha_func = beta_func*(B + 1/S12);%alpha_func = sqrt(A/C)*(B + 1/S12) er_alpha_func = sqrt(((B+1/S12)*er_beta_func)^2 + (beta_func*erB)^2 + (-beta_func*erS12/(S12)^2)^2 ); set(axes,'FontSize',16) pardat = errorbar(x,y,ery,'MarkerSize',30,'Marker','.','LineStyle','none',... 'LineWidth',2,... 'DisplayName','Data Points',... 'Color',[0 0 1]); hold on %set(axes,'FontSize',16); parfit_m_inver = plot(x,fit_m_inver,'r'); %plot fit done by matrix inversion method set(parfit_m_inver,'LineWidth',3); %parfit_MatLab = plot(x,fit_MatLab,'g'); %plot fit done Matlab %set(parfit_MatLab,'LineWidth',1); legend('Data Points','Fit by Matrix Inversion', 'Fit by MatLab'); set(legend,'Position',[0.5491 0.7591 0.2813 0.138]); set(legend,'FontSize',26); hx = xlabel('k_{1}*L (m^{-1})','FontSize',24); %hx = xlabel('I_{quad} (A)','FontSize',24); hy = ylabel('\sigma^2 (mm^2)','FontSize',24); %fit = axes('Parent',parfit_m_inver,'FontSize',16); fit=title('14 MeV, X-Projection: k_{1}*L vs \sigma^2'); %fit=title('14 MeV, X-Projection: I_{quad} vs \sigma^2'); set(fit,'FontSize',24,'FontWeight','bold'); drawnow %saveas(png,'par_fit_x-projection.png') %fprintf('Fit Parameters by MATLAB:\n y = %.8f + (%.8f)*x + %.8f*x*x\n\n', par(3), par(2), par(1)); %fprintf('Fit Parameters by Matrix Inversion:\n y = (%.8f+-%.8f) + (%.8f+-%.8f)*x + (%.8f+-%.8f)*x*x\n\n', parameter(1),er_a(1), parameter(2),er_a(2), parameter(3),er_a(3)); %fprintf('ABC Matrix Inversion:\n A=%.8f +- %.8f, B=%.2f +- %.2f, C=%.7f +- %.7f \n\n', A, erA, B, erB, C, erC); fprintf('x-projection:\n\nemit=%.3f +- %.3f mm*mrad, emit_norm=%.2f +- %.2f mm*mrad \n\n', emit*10^6, er_emit*10^6, emit_n*10^6, er_emit_n*10^6 ); fprintf('beta=%.3f +- %.3f, alpha=%.2f +-%.2f \n\n', beta_func, er_beta_func, alpha_func, er_alpha_func ); fprintf('parabola fit for y-projection (y in mm unit): \n y = (%.5f +-%.5f) + (%.5f+-%.5f)*x + (%.5f+-%.5f)*x.*x \n\n', parameter(1)*10^6,er_a(1)*10^6, parameter(2)*10^6,er_a(2)*10^6, parameter(3)*10^6,er_a(3)*10^6 ); return