clearvars
close all

seed = 0;
lognc = 6;
sLL = 50;
Lvalues = [50,80];
meas = {'r','dpq','K'};
xname = {'$r$','$|\delta\rho_1|^2$','$K^{\beta}$'};

nLvals=size(Lvalues,2);
nmvals=size(meas,2);

ncbound=zeros(nLvals,nmvals);

for k=1:nmvals
for j=1:nLvals
L=Lvalues(j);
formatSpec = '%f';
nc = 10^lognc;

e970 = 0;
a970 = 0;

matfile=sprintf('%d%sDistL%dnc%dsLL%d.mat',seed,meas{k},L,lognc,sLL);
load(matfile,'A');

nb = 20;
h = (max(A)-min(A))/nb;

bins = [];
for i = 1:nb+1
    bins(i) = min(A)+i*h;
end

pend970 = A(1:nc);
pave970 = A(3*nc+1:4*nc);

e970 = e970 + hist(pend970,bins);
a970 = a970 + hist(pave970,bins);

fda970 = a970/(nc*h);
fde970 = e970/(nc*h);

figure((k-1)*2+j);
plot(bins,e970/(nc*h),'ko','DisplayName','$p_{\rm end}$','MarkerFaceColor','b','MarkerSize',10,'LineWidth',2.4);
hold on, plot(bins,a970/(nc*h),'kx','DisplayName','$p_{\rm ave}$','MarkerSize',15,'LineWidth',2.4);

alldata=[fde970,fda970];
ymax = 10^(ceil(log10(max(alldata))));
alldata(~alldata)=inf;
ymin = 10^(floor(log10(min(alldata))));
ylim([ymin,ymax]);
set(gca,'YTick', [10^(-8),10^(-6),10^(-4),10^(-2)]);

llim = min(0,max(A)*1.8); %Lower Limit
ulim = max(0,max(A)*1.8); %Upper Limit

xlim([llim,ulim]);

xlabel(xname{k},'interpreter','latex','FontSize',24);
ylabel('Probability Density','FontSize',24);

multi=0;
multi2=0;
for i=1:size(a970,2)
    if(a970(i)>0 && e970(i)>0)
       multi=multi+a970(i)/nc*log10(a970(i)/e970(i));
       multi2=multi2+a970(i)/nc*(a970(i)/e970(i));
    end
end
enum = sprintf('%d',10^(multi));
gnum2 = sprintf('%d',multi2);
outdata(1,1)=10^(multi);
outdata(1,2)=multi2;

meane=mean(pend970);
stde=std(pend970);

llim=round(llim);
ulim=round(ulim);
xlin=linspace(llim,ulim,ulim-llim+1);

efit = (1/(stde*sqrt(2*pi)))*exp((-1/2)*((xlin-meane)/stde).^2);
hold on, plot(xlin,efit,'k-','LineWidth',3.2,'DisplayName','$p_{\rm end}$ Fit');

meana=mean(pave970);
stda=std(pave970);

afit = (1/(stda*sqrt(2*pi)))*exp((-1/2)*((xlin-meana)/stda).^2);
hold on, plot(xlin,afit,'r-','LineWidth',3.2,'DisplayName','$p_{\rm ave}$ Fit');

set(gca, 'FontSize', 24);
set(gcf,'units','points','position',[10,10,400,300]);
set(gca,'box','off');
set(gca,'linewidth',2);

set(gca,'YScale','log');
lgd = legend('show');
lgd.FontSize = 24;
set(lgd,'Interpreter','Latex');
legend('boxoff','Location','Best');

plotname = sprintf('%sFits1,%d,%d,%d.png',meas{k},L,lognc,sLL);
saveas(gcf,plotname);

if meane>meana
    intlim=1;
    idx=round(bins(find(fda970~=0, 1, 'first'))-h/2);
else
    intlim=ulim;   
    idx=round(bins(find(fda970~=0, 1, 'last'))-h/2);
end
[maxf,argmaxf]=max(afit);
integral=0;
for i=min(intlim,argmaxf):max(intlim,argmaxf)-1
    integral=integral+(afit(i+1)+afit(i))/2;
end
outdata(2,1)=integral;

integral=0;
for i=min(intlim,idx):max(intlim,idx)-1
    integral=integral+(afit(i+1)+afit(i))/2;
end
outdata(2,2)=integral;

integral=0;
for i=min(intlim,argmaxf):max(intlim,argmaxf)-1
    integral=integral+(efit(i+1)+efit(i))/2;
end
outdata(2,3)=integral;

integral=0;
for i=min(intlim,idx):max(intlim,idx)-1
    integral=integral+(efit(i+1)+efit(i))/2;
end
outdata(2,4)=integral;

outfile=sprintf('%sIntsL%dnc%dsLL%d.mat',meas{k},L,lognc,sLL);
save(outfile,'outdata');

ncbound(j,k)=exp(((meane-meana)/(stde))^2);

end
end