source = '4046_500ms_10i_10o_03st.txt';

raw_data = load(source);
N = length(raw_data(:,1));

baseline = min(raw_data(:,3));

data = zeros(N,3);
for i = 1:N
    data(i,1) = raw_data(i,2);
    %data(i,2) = raw_data(i,3)-baseline; 
    data(i,2) = raw_data(i,3);
    data(i,3) = sqrt(raw_data(i,3));
end

errorbar(data(:,1),data(:,2),data(:,3),'o');
xlabel('Nominal wavelength from monochromator (A)');
ylabel('Photon count (cps)');
headline = sprintf('Data from %s',source);
title(headline);

% Fitting
function_gaussian = @(x,a) a(1)*exp(-((x-a(2))/a(3)).^2)+a(4);
dYda = {@(x,a) 1}; sgn = 0;
[M M_index] = max(data(:,2));
a0 = [M data(M_index,1) 0.15 baseline];
[a,aerr,chisq,yfit,corr] = levmar(data(:,1),data(:,2),data(:,3),...
    function_gaussian,a0,dYda,sgn);
hold on;
x = linspace(data(1,1),data(N,1));
y = function_gaussian(x,a);
plot(x,y,'r','linewidth',2);
chisqdof = chisq/(N-4);