function varargout = Emit_Parabola_Fit_kl_XProjection(mean_sigy,er_sigy)
%This function will fit parabole for given rms values, plot them and from
%parabola parameter extract emitance.
% Written by Sadiq Setiniyaz, Jan 2010

%B = Sig(:,2);
%Er_B = Sig(:,3);
%r_quad = 0.0254; %distance from center of dipole to the pole face is 1 inche
%g = B/r_quad; %quadrupole gradiant
En = 10; % e- energy in Mev
me = 0.511; % electron rest mass in MeV
p = sqrt(En*En - me*me)/1000;% e- momentum in GeV/c
%k1=0.2998*(g/p);%quadruople strength

for n = 1:41
    k1(n) = (-6.5 + (n - 1) * 0.65); % magnetic field strength
end


k1L = 0.3*k1(:);
%kl = 0.15*Sig(:,2); %k1*L
%kl = 0.165*Sig(:,2);
%kl = 0.18*Sig(:,2);
%kl = 0.20*Sig(:,2);
%kl = 0.22*Sig(:,2);

%S12=2.07500;%m thin lense aproximation.
S12=2;%m (2.075 - 0.15/2)
%S12=1.9925;%m (2.075 - 0.165/2)
%S12=1.985;%m (2.075 - 0.18/2)
%S12=1.975;%m (2.075 - 0.20/2)
%S12=1.965;%m (2.075 - 0.22/2)

[col,NumPoints] = size(mean_sigy);
x = k1L;
y = mean_sigy'.*mean_sigy';
ery = er_sigy'.*er_sigy';

%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)
%parameter
fit_m_inver = parameter(1) + parameter(2)*x + parameter(3)*x.*x;

%fit parameters by MatLab
par = polyfit(x,y,2);
%par
fit_MatLab = par(3) + par(2)*x + par(1)*x.*x;

er_a = diag(alpha_invert);

%Get Emittance from fit
En = 10; % e- energy in Mev
me = 0.511; % electron rest mass in MeV
gamma = En/me;
S12=0.64; % from magnet cneter to the screen, 56.5 + 15/2
erS12 = 0.001; % assume 1 mm error.

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 );


pardat = errorbar(x,y,ery,'.');
hold on

parfit_m_inver = plot(x,fit_m_inver,'r'); %plot fit done by matrix inversion method
set(parfit_m_inver,'LineWidth',1); 

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');
hx = xlabel('k_{1}*L (m^{-1})');
hy = ylabel('\sigma^2 (m^2)');
fit=title('Y-Projection: k_{1}*L vs \sigma^2');
set(fit,'FontSize',12,'FontWeight','bold');
drawnow

fprintf('y-projection:\n emit=%.4f +- %.4f mm*mrad, emit_norm=%.4f +- %.4f mm*mrad \n', emit*10^6, er_emit*10^6, emit_n*10^6, er_emit_n*10^6 );
fprintf('beta=%12.15f +- %12.15f, alpha=%12.15f +-%12.15f \n', beta_func, er_beta_func, alpha_func, er_alpha_func );

return