%------------------------------------------------------------
% EXAMPLE OUTPUT
%------------------------------------------------------------
% Marquardt Gradient-Expansion Algorithm 
% i 	 Chisqr 	 a1 	 a2 	 a3 	 a4 	 a5 	 lambda
%     1   661.6   0.032  -0.032  20.000 1.000e-003
%     2    16.5   0.029  -0.030  19.454 1.000e-004
%     3    16.5   0.029  -0.030  19.371 1.000e-005
% Normalized chi sq is 1.499373
% Estimate of T1 is 19.372128 pm 0.793155
%------------------------------------------------------------
files = dir('*.txt');
Nfiles = length(files);

time_digits = 2;
window = 3;
data = zeros(Nfiles,3);
subplot(211);
fs = 14;
xlabel('Time (in ms)','fontsize',fs);
ylabel('FID amplitude (arb. scale)','fontsize',fs);
title('90-90 pulse sequence','fontsize',fs);

for j = 1:Nfiles
    filename = getfield(files(j),'name');
    tau = str2num(filename([1:time_digits]));
    raw_data = load(filename);
    plot(1000*raw_data(:,1),raw_data(:,2));
    hold on;
    
    [M max_index] = max(raw_data(:,2));
    samples = raw_data((max_index-window):(max_index+window),2);
    
    data(j,1) = tau;
    data(j,2) = mean(samples);
    data(j,3) = std(samples);
end

subplot(212);
baseline = data(1,2);
N = length(data(:,1));
data = data(2:N,:); % Omit the baseline
N = N-1;
errorbar(data(:,1),data(:,2),data(:,3),'o','markersize',4); hold on;

% Fitting
function_exp = @(x,a) a(1) + a(2)*exp(-x/a(3));
a0 = [baseline -baseline 20];
dYda = {@(x,a) 1}; sgn = 0;
[a,aerr,chisq,yfit,corr] = levmar(data(:,1),data(:,2),data(:,3),...
    function_exp,a0,dYda,sgn);

M = 1000;
x = linspace(data(1,1),data(N,1),M);
y = function_exp(x,a);
plot(x,y,'r','linewidth',2);
legend('Measured amplitudes',...
       'Fitted exponential decay',4);
xlabel('Time (in ms)','fontsize',fs);
ylabel('Pulse amplitudes','fontsize',fs);

% Results
chisqdof = chisq/(N-length(a));
fprintf(1,'Normalized chi sq is %f\n', chisqdof);
fprintf(1,'Estimate of T1 is %f pm %f\n', a(3), aerr(3));