clear all;
close all;

cm = [1,2]; %Cloning Methods
Lvalues = [50,80];
colours = {'ro','rx';'bo','bx'};
fitcol = {'r-','k-';'b-','g-'};

load('K.mat');
analytical=K(end,1); %K at sL2=50
figure(2), hold on, plot([1,10^6],[analytical,analytical],'--','Color',[0 0.5 0],'LineWidth',3.2);
figure(3), hold on, plot([1,10^6],[analytical,analytical],'--','Color',[0 0.5 0],'LineWidth',3.2);

x=zeros(2,size(Lvalues,2),size(cm,2)); %%Variables, system sizes, cloning methods
x3=zeros(2,size(Lvalues,2),size(cm,2)); %%Variables, system sizes, cloning methods
relerr=zeros(7,size(Lvalues,2),size(cm,2)); %%Relative error of each data point
nc6err=zeros(size(Lvalues,2),size(cm,2)); %%Error in 10^4 - 10^5 fit at 10^6

for j=1:size(Lvalues,2)
for n=1:size(cm,2)
for m = 0:6
seed = 0;

logt = 4;
L = Lvalues(j);
lognc = m;
dt = 10;

repspf=10; %Repeats per file

FILENAME = sprintf('%dKTrack%d,%d,%d,%d,50,%d.txt',0,cm(n),logt,L,lognc,dt);
fileID = fopen(FILENAME,'r');
formatSpec = '%f';

A = fscanf(fileID,formatSpec);

for i=1:repspf
    k(i) = A(i*5 - 4);
end

kave(m+1) = mean(k)/L;
diff(m+1) = kave(m+1)-analytical;
ksterr(n,m+1) = std(k/L)/sqrt(repspf);
end

y = [10^(-4),1;10^(-5),1];
b = [kave(5);kave(6)];
x(:,j,n) = inv(y'*y)*(y'*b);

nc=10.^linspace(0,6,7);
figure(1), hold on, plot(nc,kave,colours{j,n},'MarkerSize',10,'LineWidth',2);
set(gca,'xscale','log')
xlim([min(nc),max(nc)]);
ylim([0,0.6]);
xlabel('$n_c$','interpreter','latex');
yname = '$\hat{k}(\lambda)_{n_c,t_{obs}}$';
ylabel(yname,'interpreter','latex');
set(gca, 'FontSize', 24);
set(gca,'linewidth',2);
set(gca,'XTick', [1,100,10000,1000000]);
set(gca,'YTick', 0:0.2:0.4);
set(gcf,'units','points','position',[10,10,400,300]);

ncfit=10.^linspace(2,6,100);
figure(2), hold on, plot(nc,kave,colours{j,n},'MarkerSize',10,'LineWidth',2);
figure(2), hold on, plot(ncfit,x(1,j,n)./ncfit+x(2,j,n),fitcol{j,n});
set(gca,'xscale','log');
xlim([min(nc),max(nc)]);
ylim([0,0.6]);
xlabel('$n_c$','interpreter','latex');
yname = '$\hat{k}(\lambda)_{n_c,t_{obs}}$';
ylabel(yname,'interpreter','latex');
set(gca, 'FontSize', 20);
set(gca,'linewidth',2);
set(gca,'XTick', [1,100,10000,1000000]);
set(gca,'YTick', 0:0.2:0.4);


nc6err(j,n)=abs(x(1,j,n)/(10^6)+x(2,j,n)-kave(7));

y3 = [10^(-4),1;10^(-5),1;10^(-6),1];
b3 = [kave(5);kave(6);kave(7)];
x3(:,j,n) = inv(y3'*y3)*(y3'*b3);
relerr(:,j,n) = abs(x3(2,j,n)-kave(:))./kave(:);

nc3=10.^linspace(0,6,7);
ncfit3=10.^linspace(2,6,100);
figure(3), hold on, plot(nc3,kave,colours{j,n},'MarkerSize',10,'LineWidth',2);
figure(3), hold on, plot(ncfit3,x3(1,j,n)./ncfit3+x3(2,j,n),fitcol{j,n});
set(gca,'xscale','log');
xlim([min(nc3),max(nc3)]);
ylim([0,0.6]);
xlabel('$n_c$','interpreter','latex');
yname = '$\hat{k}(\lambda)_{n_c,t_{obs}}$';
ylabel(yname,'interpreter','latex');
set(gca, 'FontSize', 20);
set(gca,'linewidth',2);
set(gca,'XTick', [1,100,10000,1000000]);
set(gca,'YTick', 0:0.2:0.4);


figure(4), hold on, plot(nc,diff,colours{j,n},'MarkerSize',10,'LineWidth',2);
set(gca,'xscale','log');
set(gca,'yscale','log');

xlim([min(nc),max(nc)]);
ylim([0,5*max(max(diff))]);
xlabel('$n_c$','interpreter','latex');
ylabel('D');
set(gca, 'linewidth', 2);
set(gca,'XTick', [1,100,10000,1000000]);
set(gca,'YTick', [0.002,0.02,0.2]);

plotname = sprintf('ncDiffsLL50.png');
saveas(gcf,plotname);

end
end

figure(1), hold on, plot([1,10^6],[analytical,analytical],'--','Color',[0 0.5 0],'LineWidth',3.2);
figure(1), hold on, hold on, plot(ncfit,x(1,2,2)./ncfit+x(2,2,2),'-','Color',[0 0 0],'LineWidth',1.5);
figure(1), hold on, plot([10000,27000],[0.08,0.08],'-','Color',[0 0 0],'LineWidth',1.5);
figure(1), text(30000,0.08,'Fit to (28)','FontSize',15,'interpreter','latex');
legend({'$L=50$ [iid]','$L=50$ [eq]','$L=80$ [iid]','$L=80$ [eq]','Theory'},'FontSize',12,'location','northeast','interpreter','latex');
legend('boxoff');

pl2name = sprintf('nc45fit.png');
figure(2),saveas(gcf,pl2name);
pl3name = sprintf('nc456fit.png');
figure(3),saveas(gcf,pl3name);

plotname = sprintf('ncConvsLL50.png');
figure(1),saveas(gcf,plotname);

plotname = sprintf('ncDiffsLL50.png');
figure(4),saveas(gcf,plotname);