clear all; close all; clc;

load meany_TinT_SCff0_2k0_5.mat; % Load a pre-saved workspace which results after running mean1Dden.m. Should contain the variables time, x, meany.
global smallit

xi = zeros(length(time),1);

xst = 10; % For your convenience, you can choose the segment of the density to be fitted by specifying start and end values of x
xnd = 25; % which are then converted to indices. Current conversion assumes Rf=30 and step=0.1 in tdepsims.m was used.
% Decide fitting segment by watching the movies generated in mean1Dden.m and looking for the linear section in the log plots.

st = 10*xst + 301; % start fit index in x
nd = 10*xnd + 301; % end fit index in x

for k = 51:101 % Time index. By watching the movies in mean1Dden.m, you should choose *when* to start fitting. Pick a time at which the density more or less settles down, becomes exponential, and only changes slowly afterwards. This will usually happen at t ~ 50t_0.
k

xuse = x(st:nd);
smallit = meany(k,st:nd);
err = meany(k,st:nd)/20; % error arbitrarily taken as 5%

xuse = xuse(:);
smallit = smallit(:);
err = err(:);

blah = polyfit(xuse,log(smallit),1);

guess = -blah(1);
plot(x,meany(k,:),'bo'); grid on; hold on;

[p,perr,nchi2] = nonlinft('exp_model',xuse,smallit,err,guess,1);

Mat = zeros(length(xuse),2);
Mat(:,1) = exp(-p*xuse);
Mat(:,2) = ones(length(xuse),1);
format long; C = Mat\smallit
Mat = zeros(length(xuse),2);
Mat(:,1) = exp(-p*xuse);
Mat(:,2) = ones(length(xuse),1);
fit_it = Mat*C;

plot(xuse,fit_it,'r-','linewidth',2); hold off; pause;

xi(k) = 2/p; % Record time-dependent \xi in a vector; it's in units of \ell.

end

plot(time,xi,'ro'); grid on;

