function [spectra_surface,spectra_bulk,error] = fit_spectralsf(spectra,thickness)
 
% A least square fitting method is applied in fit_spectralsf to estimate the surface and
% bulk spectra for a dataset including spectra acquired at regions with known thickness. 
% 
% Input
% spectra       pxm matrix, p_ the number of spectra including in a ¡®thickness¡¯ series,
%                   m_the number of energy channels
% thickness   1xp matrix, unit: unit cells
% 
% Output
% spectra_surface             1xm matrix, surface spectra
% spectra_bulk                 1xm matrix, bulk spectra
% error                             2xm matrix, the first row includes the uncertainties for 
%                                      spectra_bulk, and the second row contains the uncertainties for 
%                                      spectra_surface.
% Guozhen Zhu
% Oct 11st, 2011
 
% check inputs
if isnumeric(spectra)&&isnumeric(thickness)
    [n1,m] = size(spectra);
    [m2,n2] = size(thickness);
    if m2 ==1
        p = n2;
        if n2==m
            spectra = spectra';
            m = n1;
        end
    elseif n2 ==1
        p = m2;
        if m2 == m
            spectra = spectra';
            m = n1;
        end
    end
else
    fprintf('Input Error! Invalid Input!')
end
 
% least square fitting
thickness_matrix = zeros(2,p);
for i=1:1:p
    thickness_matrix(1,i) = thickness(i)/(4+thickness(i));
    thickness_matrix(2,i) = 2/(4+thickness(i));
end
temp = inv(thickness_matrix*thickness_matrix')*thickness_matrix;
results = temp*spectra;
spectra_bulk = results(1,:);
spectra_surface = results(2,:);
 
residual = spectra - thickness_matrix'*results; 
 
% error
m = size(spectra,2);
estimat = zeros(1,m); 
for i=1:1:m
    for j = 1:1:p
    estimat(i) = residual(j,i)^2+estimat(i);
    end
end
estimat2 = (1/(p-2))*sqrt(estimat);

error = [];
for i=1:1:m
    epsilon2 = residual(:,i)*residual(:,i)';
    error = [error, sqrt(abs(diag(temp*epsilon2*temp')))];;
end
 
% plot results
figure; 
hold on;
plot(results(1,:),'b');
plot(results(2,:),'r');
plot(estimat2,'k');
plot(error(1,:),'g');
plot(error(2,:),'m');