%------------------------------------------------------------
% Tony Hyun Kim
% Analyzes the output from the monochromator (Labview prog.)
%------------------------------------------------------------
source = '4861_1000ms_50i_50o_03st.txt';
hshift = 0;
push_index = 0; % If fit seems a bit off, we can push it 
tighten = 1; % Factor for the initial guess of std dev.
boost = 1; % Factor for boosting baseline

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

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

N = length(data(:,1));
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))/(2*a(3))).^2)+a(4);
dYda = {@(x,a) 1}; sgn = 0;

% Estimate a decent guess
baseline = min(data(:,2));
[M M_index] = max(data(:,2));
M_index = M_index + push_index;
M = data(M_index,2) - baseline;
half_max = 1/2*M + baseline;
low_hm_index = 1;
high_hm_index = N;
for i = 1:M_index
    if (data(i,2) > half_max)
        low_hm_index = i;
        break;
    end
end
for i = M_index:N
    if (data(i,2) < half_max)
        high_hm_index = i-1;
        break;
    end
end
fwhm = data(high_hm_index,1) - data(low_hm_index,1);
a0 = [M data(M_index,1)+hshift tighten*fwhm/2.355 baseline*boost];

% Second method for finding decent guess;
% [a1,aerr1,chisq1,yfit1,corr1] = levmar(data((M_index-1):(M_index+2),1),data((M_index-1):(M_index+2),2),data((M_index-1):(M_index+2),3),...
%     function_gaussian,a0,dYda,sgn);

[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);
z = function_gaussian(x,a1);
plot(x,z,'g','linewidth',2);
hold on;
plot(x,y,'r','linewidth',2);
legend('Data','Initial guess','Final fit');
chisqdof = chisq/(N-4);

% Create a good window
extraX = 0.1*(data(N,1)-data(1,1));
[M M_index] = max(data(:,2));
M = M + data(M_index,3);
[m m_index] = min(data(:,2));
m = m - data(m_index,3);
extraY = 0.1*(M - m);
xlim([data(1,1)-extraX data(N,1)+extraX]);
ylim([m-extraY M+extraY]);

% Output results
fprintf(1,'Reduced chisq: %.4f\n',chisqdof);
fprintf(1,'Location of peak: %.4f, err: %.6f\n',a(2),aerr(2));