clear all:
close all;

numOfBoxes=5;
nParticles=16*(numOfBoxes^3)
nRuns=112;

Pd=[0.01,0.02,0.03];

PdLabel={'0.01','0.02','0.03'};
PdLabel1={'0.01','0.01','0.02','0.02','0.03','0.03'};
cM1={'r-o','g-o','b-o','m-o','c-o'};
cM2={'r-s','g-s','b-s','m-s','c-s'};
cM={'r-','g-','b-','m-','c-'};
%PdLineY={[8e-4 2e-3],[2e-3 5e-3],[5e-3 9e-3]}
PdLineY={[1.2e-3 2.3e-3],[3e-3 6e-3],[8e-3 1.7e-2]}

topDir = './gardnerHS1/';

Npds = numel(Pd);

data1 = [];
for i =1:Npds

    datFile1 = [ topDir 'd_dab_betaP_phi_N' num2str(nParticles) '_s' num2str(Pd(i)) '.dat']
    data1{i} = load(datFile1);
    
end


dim1=500;
dim2=400;
fig1=figure('position',[0,0,dim1,dim2]);

for i=1:Npds
    semilogy(data1{i}(:,2),data1{i}(:,4),cM1{i});
    hold on;
end
for i=1:Npds
    semilogy(data1{i}(:,2),data1{i}(:,3),cM2{i});
    semilogy( [1/Pd(i) 1/Pd(i)], PdLineY{i} ,cM{i});
end
title([ 'HS1   MSD vs betaP   nParticles = ' num2str(nParticles)]);
ylabel('$$(\Delta,\Delta_{AB})/\bar{\sigma}^2$$','Interpreter','Latex')
xlabel('$${\beta P \bar{\sigma}^3}$$','Interpreter','Latex');
legend(['Pd=' PdLabel{1}],...
    ['Pd=' PdLabel{2}],...
    ['Pd=' PdLabel{3}],...
    'Location', 'NorthEast');
%axis equal
xlim([0,300])
%ylim([1e-5,2e+0])