function [spectra_surface,spectra_bulk,error] = fit_spectrawlsf(spectra,spectraerror, thickness)
 
% A weighted 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
% spectraeeor   pxm matrix, the uncertainties for each value in spectra
% 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 uncertainties for 
%                                      spectra_bulk, and the second row is uncertainties for 
%                                      spectra_surface.
%  
% Guozhen Zhu
% Oct 15th, 2012
 
% input
if isnumeric(spectra)&&isnumeric(spectraerror)&&isnumeric(thickness)
    [n1,m] = size(spectra);
    [n3,m3] = size(spectraerror);
    [m2,n2] = size(thickness);
    if n1==n3&& m==m3
        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
else
    fprintf('Input Error! Invalid Input!')
end
 
% weighted 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
 
spectra_bulk = zeros(1,m);
spectra_spectra = zeros(1,m);
estimat = zeros(1,m); 
error = [];
for i = 1:1:m
    matrix_weight = zeros(p,p);
    for k =1:1:p
        matrix_weight(k,k) = spectraerror(k,i);
    end
    temp = inv(thickness_matrix*matrix_weight*thickness_matrix')*thickness_matrix;
    results = temp*matrix_weight*spectra(:,i);
    spectra_bulk(1,i) = results(1);
    spectra_surface(1,i) = results(2);
 
    residual = spectra(:,i) - thickness_matrix'*results; 
    for k = 1:1:p
    estimat(i) = residual(k)^2+estimat(i);
    end
 
    epsilon2 = residual*residual';
    error = [error, sqrt(abs(diag(temp*epsilon2*temp')))];
end
 
estimat2 = (1/(p-2))*sqrt(estimat);
 
% plot results
figure; 
hold on;
plot(spectra_bulk,'b');
plot(spectra_surface,'r');
plot(estimat2,'k');
plot(error(1,:),'g');
plot(error(2,:),'m');