% EExtra_1.m
% Solution to problem 1 of ConvExtra. JMS, Jun.2024

clear

% Load data
load('EExtraDatos.mat');
[~,ns,np] = size(posR);  % number of steps and particles
ns = ns-1;

% Plot trajectory of three particles
figure(11)
s = 1:ns;
plot(posR(1,:,25), posR(2,:,25), ...
     posR(1,:,50), posR(2,:,50), ...
     posR(1,:,75), posR(2,:,75),'LineWidth',2)
legend('Particle 25','Particle 50','Particle 75','Location','NorthWest')
xlabel('x (Angstrom)')
xlabel('y (Angstrom)')
set(gca,'FontSize',15)
grid on

% Find jump lengths in each direction
dx = diff(posR(1,:,:));  % horizontal jumps
dy = diff(posR(2,:,:));  % vertical jumps

% Find jump probabilities
Pr = sum(dx(:)==+2.5)/ns/np
Pl = sum(dx(:)==-2.5)/ns/np
Pu = sum(dy(:)==+2.5)/ns/np
Pd = sum(dy(:)==-2.5)/ns/np
P0 = sum(dx(:)==0 & dy(:)==0)/ns/np
Ptot = Pr+Pl+Pu+Pd+P0

% Plot histogram of final distances
figure(12)
d = sqrt(sum((posR(:,101,:)-posR(:,1,:)).^2));
hist(d(:),20)
xlabel('d_{100} (Angstrom)')
ylabel('number of particles')
set(gca,'FontSize',15)

% Plot distance as a function of number of steps
figure(13)
d = sqrt(sum((posR-posR(:,1,:)).^2));
davg = sum(d,3)/np;
dstd = sqrt( sum((d-davg).^2,3)/np );
errorbar(0:ns,davg,dstd,'o')
xlabel('number of steps')
ylabel('distance (Angstrom)')
set(gca,'FontSize',15)
grid on
