% timesteps in each random walk
timesteps=10000;
plotstep=10;
% diffusion coefficient for the random walk
D=0.5;
observation_radius=10;
observation_steps=100;
lim_factor=1.4;
plots=10;
plotevery=1;

%calculate grid size
dr=sqrt(2*D);
t=1:timesteps;

% simulate random walks on grid
x=zeros(plots,timesteps);
y=x;
for k=1:plots
    for ct=2:timesteps
        x(k,ct)=x(k,ct-1)+dr*(randi(2)-1.5)*2;
        y(k,ct)=y(k,ct-1)+dr*(randi(2)-1.5)*2;
    end
end

% calculate MSD as ensemble average on a logarihmic lag time scale
msdtau=round(logspace(0,floor(log10(timesteps)),50));
msd(1:length(msdtau))=0;
for tau=1:length(msdtau)
    msd(tau)=mean( (x(:,msdtau(tau))-x(:,1)).^2 + (y(:,msdtau(tau))-y(:,1)).^2 );
end

% plot random walk
figure(1)
subplot(2,1,1)
for k=1:plotevery:plots
    plot(x(k,1:plotstep:timesteps),y(k,1:plotstep:timesteps), 'Color', hsv2rgb([k/(plots+1),1,1]));
    if (k==1) 
        hold on;
    end
end
for k=5:5:length(msdtau)
    drr=sqrt(msd(k));
    rectangle('Position',[-drr,-drr,2*drr,2*drr], 'Curvature',[1,1], 'LineWidth', 2);
end
daspect([1,1,1])
hold off
xlim([-lim_factor*drr lim_factor*drr]);
ylim([-lim_factor*drr lim_factor*drr]);
xlabel('coordinate x');
ylabel('coordinate y');
title('random walk trajectories');

% plot MSD
subplot(2,1,2)
loglog(msdtau, msd, 'LineWidth', 2);
hold on
loglog(msdtau, 4*D*msdtau, 'r--', 'LineWidth', 2);
hold off
xlabel('lag time t')
ylabel('MSD \langle r^2(\tau)\rangle');
M{1}='measured MSD';
M{2}='\langle r^2(\tau)\rangle=4D\cdot\tau';
legend(M,'Location', 'NorthWest');
title('mean squared displacement');
xlim([1 1e4]);
ylim([1 1e4]);
daspect([1,1,1])