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