% first exit time distributions as a function of discretization
%Remark:  This simulation may take a few minutes
clear all

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);

N = 2000;  %number of particles
L = 1; %reflect at x=0, absorb at x=L
D = 1; % diffusion coefficient

dt_list = [.004,.002,.001,.0005,.00025,0.000125];
for jex = 1:length(dt_list)
dt = dt_list(jex);  %time step
diff = sqrt(2*D*dt); % diffusion step size
 
jst = 20; % number of initial position increments
dx = L/jst;
for js = 1:jst
    x0 = (js-1)*dx;  % starting positions
for j = 1:N
    clear x X
    x = x0;
    k = 1;
    
    X(1) = 0;
    while (x<L) % step until x reaches L
        x = abs(x+diff*randn(1,1)); % absolute value gives reflection at the origin
        k = k + 1;
        X(k) = x;
        
    end
    Ts(j) = dt*(k-1);
    %if(js==1)
        %if (j<6)% make a plot of 5 sample paths
        %figure(1)
        %plot(dt*[0:k-1],X,'linewidth',2)
        %hold on
        %xlabel('t','fontsize',20)
        %ylabel('x','fontsize',20)
        %axis([0 max(Ts(1:j)) 0 L])
        %end
    %end
    
end 
hold off
   
Tm(js) = mean(Ts);
xs(js) = x0;
end

x = [0:.1:L];

%fx = 1-xs.^2/L^2;  % fit data to a parabolic curve
 %a = sum(fx.*Tm)/sum(fx.^2)
 
 Tx = (L^2-x.^2)/(2*D);

% plot the mean first exit times
figure(4)
plot(xs,Tm,'*',x,Tx,'--','linewidth',2)
xlabel('initial position','fontsize',20)
ylabel('Mean Exit Time','fontsize',20)
axis([0 L 0 1])

% fit the data to the parabola y = a(x^2-L^2)

al = sum((xs.^2-L^2).^2);
bt = sum(Tm.*(xs.^2-L^2));
a = bt/al;
aj(jex) = a;
dtj(jex) = dt;
figure(5)
plot(xs,Tm,'*',x,Tx,'k--',x,a*(x.^2-L^2),'r--','linewidth',2)
xlabel('initial position','fontsize',20)
ylabel('Mean Exit Time','fontsize',20)
axis([0 L 0 1])
 
end

% estimate the power law power
aj
p=polyfit(log(dtj),log(-0.5-aj),1)
yj=log(-aj-0.5);
xj=log(dtj);
A = sum(yj-xj/2)/length(yj);
y=exp(A)*dtj.^(1/2);
figure(6)
loglog(dtj,-aj-0.5,'*',dtj,y,'--','linewidth',2)
xlabel('dt','fontsize',20)
ylabel('-a-0.5','fontsize',20)

 