close all; clc; clear;

X1=[];
X2=[];
Y1=[];
Y2=[];
Y1nl=[];
Y2nl=[];

for K=1:32
    if K == 1
        fc0 = 20;
        D = 150;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 5;
    end
    
    if K == 2
        fc0 = 20;
        D = 150;
        fultf = 600;
        Ef = 2.5e4;
        dst = 5;
    end
    
    if K == 3
        fc0 = 20;
        D = 300;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 8;
    end
    
    if K == 4
        fc0 = 20;
        D = 300;
        fultf = 600;
        Ef = 2.5e4;
        dst = 8;
    end
    
    if K == 5
        fc0 = 20;
        D = 600;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 13;
    end
    
    if K == 6
        fc0 = 20;
        D = 600;
        fultf = 600;
        Ef = 2.5e4;
        dst = 13;
    end
    
    if K == 7
        fc0 = 20;
        D = 1200;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 15;
    end
    
    if K == 8
        fc0 = 20;
        D = 1200;
        fultf = 600;
        Ef = 2.5e4;
        dst = 15;
    end
    
    if K == 9
        fc0 = 30;
        D = 150;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 5;
    end
    
    if K == 10
        fc0 = 30;
        D = 150;
        fultf = 600;
        Ef = 2.5e4;
        dst = 5;
    end
    
    if K == 11
        fc0 = 30;
        D = 300;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 8;
    end
    
    if K == 12
        fc0 = 30;
        D = 300;
        fultf = 600;
        Ef = 2.5e4;
        dst = 8;
    end
    
    if K == 13
        fc0 = 30;
        D = 600;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 13;
    end
    
    if K == 14
        fc0 = 30;
        D = 600;
        fultf = 600;
        Ef = 2.5e4;
        dst = 13;
    end
    
    if K == 15
        fc0 = 30;
        D = 1200;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 15;
    end
    
    if K == 16
        fc0 = 30;
        D = 1200;
        fultf = 600;
        Ef = 2.5e4;
        dst = 15;
    end
    
    if K == 17
        fc0 = 50;
        D = 150;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 5;
    end
    
    if K == 18
        fc0 = 50;
        D = 150;
        fultf = 600;
        Ef = 2.5e4;
        dst = 5;
    end
    
    if K == 19
        fc0 = 50;
        D = 300;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 8;
    end
    
    if K == 20
        fc0 = 50;
        D = 300;
        fultf = 600;
        Ef = 2.5e4;
        dst = 8;
    end
    
    if K == 21
        fc0 = 50;
        D = 600;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 13;
    end
    
    if K == 22
        fc0 = 50;
        D = 600;
        fultf = 600;
        Ef = 2.5e4;
        dst = 13;
    end
    
    if K == 23
        fc0 = 50;
        D = 1200;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 15;
    end
    
    if K == 24
        fc0 = 50;
        D = 1200;
        fultf = 600;
        Ef = 2.5e4;
        dst = 15;
    end
    
    if K == 25
        fc0 = 70;
        D = 150;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 5;
    end
    
    if K == 26
        fc0 = 70;
        D = 150;
        fultf = 600;
        Ef = 2.5e4;
        dst = 5;
    end
    
    if K == 27
        fc0 = 70;
        D = 300;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 8;
    end
    
    if K == 28
        fc0 = 70;
        D = 300;
        fultf = 600;
        Ef = 2.5e4;
        dst = 8;
    end
    
    if K == 29
        fc0 = 70;
        D = 600;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 13;
    end
    
    if K == 30
        fc0 = 70;
        D = 600;
        fultf = 600;
        Ef = 2.5e4;
        dst = 13;
    end
    
    if K == 31
        fc0 = 70;
        D = 1200;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 15;
    end
    
    if K == 32
        fc0 = 70;
        D = 1200;
        fultf = 600;
        Ef = 2.5e4;
        dst = 15;
    end
    
    load(sprintf('P_%d_limit.mat',K));
    Pf0 = P(:,:,1);
    Pf0 = repmat(Pf0(2:end,1),1,5);
    load(sprintf('P_%d_.mat',K));
    Pf0nl = P(:,:,1);
    Pf0nl = repmat(Pf0nl(2:end,1),1,5);
    load(sprintf('P_%d_limit.mat',K));
    PP = P(2:end,2:end,1);
    PP = (PP-Pf0)./Pf0;
    load(sprintf('P_%d_.mat',K));
    PPnl = P(2:end,2:end,1);
    PPnl = (PPnl-Pf0nl)./Pf0nl;
    ind1 = P(2:end,2:end,3)>=0.08;
    ind2 = P(2:end,2:end,3)<0.08;
    p1 = PP(ind1);
    p2 = PP(ind2);
    p1nl = PPnl(ind1);
    p2nl = PPnl(ind2);
    CF = P(2:end,2:end,2);
    cf1 = CF(ind1);
    cf2 = CF(ind2);
    
    c=25;
    rho_l = 0.01;
    nbars = 6;
    dsl = sqrt(rho_l*D^2/nbars);
    rho_s = [0.5 1 2 3 4];
    ss = 4*(0.25*pi*dst^2)./(0.01*rho_s*(D-2*c));
    ke = (1-0.5*(ss-dst*ones(1,5))./(D-2*c));
    ke(ke<0)=0;
    ke(ke>0)=1;
    ke = ke.*(1-0.5*(ss-dst*ones(1,5))./(D-2*c)).^2/(1-nbars*dsl^2/D^2);
    ke(ke>1)=1;
    if Ef == 1e5
        Ef_rho_f_2 = 0.01*0.5*1e5*[0.25 0.5 1 2 3]';
    else
        Ef_rho_f_2 = 0.01*0.5*2.5e4*[0.5 1 2 4 6]';
    end
    fls = 0.5*ke.*rho_s*450/100;
    
    fls = repmat(fls,5,1);
    Ef_rho_f_2 = repmat(Ef_rho_f_2,1,5);
    x = ((D-2*c)/D)^2*fls./Ef_rho_f_2 ;
    x1 = x(ind1);
    x2 = x(ind2);
    
    X1 = [X1;x1];
    X2 = [X2;x2];
    Y1 = [Y1;p1];
    Y2 = [Y2;p2];
    Y1nl = [Y1nl;p1nl];
    Y2nl = [Y2nl;p2nl];
end

figure (1)
X=10^4*X1;Y=Y1;

color1=[0,0.4470, 0.7410];
l = plot(X,Y,'o','Linewidth',1,'Markersize',2,'MarkerEdge',color1,'MarkerFace',color1);hold on;
blfit = @(coeff,x) (x<coeff(3)).*(coeff(1)*x) + (~(x<coeff(3))).*(coeff(1)*coeff(3)+coeff(2)*(x-coeff(3)));
lsqOpts = optimset('Display','off');
IC = [1 1 150]; %Initial guess
LB = [-inf -inf min(X)];
UB = [inf inf max(X)];
coeffEst = lsqcurvefit(blfit,IC,X,Y,LB,UB,lsqOpts);
xfit = [min(X) coeffEst(3) 800];
yfit = blfit(coeffEst,xfit);
x_cutoff = coeffEst(3);
y_cutoff = blfit(coeffEst,x_cutoff);
sss = sprintf('(%.1f,%.2f)',x_cutoff,y_cutoff);
plot(xfit, yfit, 'k-', 'LineWidth', 1.);

%R^2
yfit = blfit(coeffEst,X);
SStot = sum((Y-mean(Y)).^2);                            % Total Sum-Of-Squares
SSres = sum((Y(:)-yfit(:)).^2);                         % Residual Sum-Of-Squares
Rsq = 1-SSres/SStot;                                    % R^2
sRsq = sprintf('R^2 = %.3f',Rsq);

legend('Num. results','Bilinear fit','Fontsize',12,'FontName','Times New Roman','Location','southeast');
grid on;
xlabel ('{\itC_s} (%)','Fontsize',12,'FontName','Times New Roman');
h = ylabel('$({\it\bar{P}}_{max}-{\it\bar{P}}_{f0})/{\it\bar{P}}_{f0}$');
set(h,'interpreter','Latex','FontSize',12)
set(gcf, 'Position',  [00, 400, 380, 270])
set(gca,'YLim',[0 1.25])
set(gca,'Fontsize',12,'FontName','Times New Roman')
legend boxoff
ss = '(a)  {\itf_{lf}} /{\itf''_c} \geq 0.08';
annotation('textbox', [0.18, 0.82, 0.1, 0.1], 'String', ss,'Fontsize',12,'FontName','Times New Roman','EdgeColor','none');
annotation('textbox', [0.37, 0.37, 0.3, 0.1], 'String', sss,'Fontsize',12,'FontName','Times New Roman','EdgeColor','none');
annotation('textbox', [0.58, 0.82, 0.3, 0.1], 'String', sRsq,'Fontsize',12,'FontName','Times New Roman','EdgeColor','none');

ax = gca;
% ax.XTick = [0 0.5 1 2 3 4];
ax.YTick = [0 0.25 0.5 0.75 1 1.25];

figure (2)
X=10^4*X2;Y=Y2;
color1=[0,0.4470, 0.7410];
l = plot(X,Y,'o','Linewidth',1,'Markersize',2,'MarkerEdge',color1,'MarkerFace',color1);hold on;
blfit = @(coeff,x) (x<coeff(3)).*(coeff(1)*x) + (~(x<coeff(3))).*(coeff(1)*coeff(3)+coeff(2)*(x-coeff(3)));
lsqOpts = optimset('Display','off');
IC = [1 1 200]; %Initial guess
LB = [-inf -inf min(X)];
UB = [inf inf max(X)];
coeffEst = lsqcurvefit(blfit,IC,X,Y,LB,UB,lsqOpts);
xfit = [min(X) coeffEst(3) 1500];
yfit = blfit(coeffEst,xfit);
x_cutoff = coeffEst(3);
y_cutoff = blfit(coeffEst,x_cutoff);
sss = sprintf('(%.1f,%.2f)',x_cutoff,y_cutoff);
plot(xfit, yfit, 'k-', 'LineWidth', 1.);

%R^2
yfit = blfit(coeffEst,X);
SStot = sum((Y-mean(Y)).^2);                            % Total Sum-Of-Squares
SSres = sum((Y(:)-yfit(:)).^2);                         % Residual Sum-Of-Squares
Rsq = 1-SSres/SStot;                                    % R^2
sRsq = sprintf('R^2 = %.3f',Rsq);

legend({'Num. results','Bilinear fit'},'Location','southeast','Fontsize',12,'FontName','Times New Roman','Location','southeast');
grid on;
xlabel ('{\itC_s} (%)','Fontsize',12,'FontName','Times New Roman');
h = ylabel('$({\it\bar{P}}_{max}-{\it\bar{P}}_{f0})/{\it\bar{P}}_{f0}$');
set(h,'interpreter','Latex','FontSize',12)
set(gcf, 'Position',  [00, 400, 380, 270])
set(gca,'YLim',[0 1.25])
set(gca,'Fontsize',12,'FontName','Times New Roman')
legend boxoff
ss = '(b)  {\itf_{lf}} /{\itf''_c} < 0.08';
annotation('textbox', [0.18, 0.82, 0.1, 0.1], 'String', ss,'Fontsize',12,'FontName','Times New Roman','EdgeColor','none');
annotation('textbox', [0.6, 0.38, 0.3, 0.1], 'String', sss,'Fontsize',12,'FontName','Times New Roman','EdgeColor','none');
annotation('textbox', [0.65, 0.82, 0.3, 0.1], 'String', sRsq,'Fontsize',12,'FontName','Times New Roman','EdgeColor','none');

ax = gca;
% ax.XTick = [0 0.5 1 2 3 4];
ax.YTick = [0 0.25 0.5 0.75 1 1.25];

figure (3)
X=10^4*X1;Y=Y1nl;
color1=[0,0.4470, 0.7410];
l = plot(X,Y,'o','Linewidth',1,'Markersize',2,'MarkerEdge',color1,'MarkerFace',color1);hold on;
blfit = @(coeff,x) (x<coeff(3)).*(coeff(1)*x) + (~(x<coeff(3))).*(coeff(1)*coeff(3)+coeff(2)*(x-coeff(3)));
lsqOpts = optimset('Display','off');
IC = [1 1 150]; %Initial guess
LB = [-inf -inf min(X)];
UB = [inf inf max(X)];
coeffEst = lsqcurvefit(blfit,IC,X,Y,LB,UB,lsqOpts);
xfit = [min(X) coeffEst(3) 800];
yfit = blfit(coeffEst,xfit);
x_cutoff = coeffEst(3);
y_cutoff = blfit(coeffEst,x_cutoff);
sss = sprintf('(%.1f,%.2f)',x_cutoff,y_cutoff);
plot(xfit, yfit, 'k-', 'LineWidth', 1.);

%R^2
yfit = blfit(coeffEst,X);
SStot = sum((Y-mean(Y)).^2);                            % Total Sum-Of-Squares
SSres = sum((Y(:)-yfit(:)).^2);                         % Residual Sum-Of-Squares
Rsq = 1-SSres/SStot;                                    % R^2
sRsq = sprintf('R^2 = %.3f',Rsq);

legend('Num. results','Bilinear fit','Location','southeast','Fontsize',12,'FontName','Times New Roman','Location','southeast');
grid on;
xlabel ('{\itC_s} (%)','Fontsize',12,'FontName','Times New Roman');
h = ylabel('$({\it{P}}_{max}-{\it{P}}_{f0})/{\it{P}}_{f0}$');
set(h,'interpreter','Latex','FontSize',12)
set(gcf, 'Position',  [00, 400, 380, 270])
set(gca,'YLim',[0 1.25])
set(gca,'Fontsize',12,'FontName','Times New Roman')
legend boxoff
ss = '(c)  {\itf_{lf}} /{\itf''_c} \geq 0.08';
annotation('textbox', [0.18, 0.82, 0.1, 0.1], 'String', ss,'Fontsize',12,'FontName','Times New Roman','EdgeColor','none');
annotation('textbox', [0.43, 0.38, 0.3, 0.1], 'String', sss,'Fontsize',12,'FontName','Times New Roman','EdgeColor','none');
annotation('textbox', [0.65, 0.82, 0.3, 0.1], 'String', sRsq,'Fontsize',12,'FontName','Times New Roman','EdgeColor','none');

ax = gca;
% ax.XTick = [0 0.5 1 2 3 4];
ax.YTick = [0 0.25 0.5 0.75 1 1.25];

figure (4)
X=10^4*X2;Y=Y2nl;
color1=[0,0.4470, 0.7410];
l = plot(X,Y,'o','Linewidth',1,'Markersize',2,'MarkerEdge',color1,'MarkerFace',color1);hold on;
blfit = @(coeff,x) (x<coeff(3)).*(coeff(1)*x) + (~(x<coeff(3))).*(coeff(1)*coeff(3)+coeff(2)*(x-coeff(3)));
lsqOpts = optimset('Display','off');
IC = [1 1 250]; %Initial guess
LB = [-inf -inf min(X)];
UB = [inf inf max(X)];
coeffEst = lsqcurvefit(blfit,IC,X,Y,LB,UB,lsqOpts);
xfit = [min(X) coeffEst(3) 1500];
yfit = blfit(coeffEst,xfit);
x_cutoff = coeffEst(3);
y_cutoff = blfit(coeffEst,x_cutoff);
sss = sprintf('(%.1f,%.2f)',x_cutoff,y_cutoff);
plot(xfit, yfit, 'k-', 'LineWidth', 1.);

%R^2
yfit = blfit(coeffEst,X);
SStot = sum((Y-mean(Y)).^2);                            % Total Sum-Of-Squares
SSres = sum((Y(:)-yfit(:)).^2);                         % Residual Sum-Of-Squares
Rsq = 1-SSres/SStot;                                    % R^2
sRsq = sprintf('R^2 = %.3f',Rsq);

legend('Num. results','Bilinear fit','Location','southeast','Fontsize',12,'FontName','Times New Roman','Location','southeast');
grid on;
xlabel ('{\itC_s} (%)','Fontsize',12,'FontName','Times New Roman');
% h = ylabel('$({\it\bar{P}}_{max}-{\it\bar{P}}_{f0})/{\it\bar{P}}_{f0}$');
h = ylabel('$({\it{P}}_{max}-{\it{P}}_{f0})/{\it{P}}_{f0}$');
%h = ylabel('${\it\bar{P}}_{max}/{\it\bar{P}}_{0}$');
set(h,'interpreter','Latex','FontSize',12)
% ylabel ('{\itP}_{max,lim}/{\itP_{f }}_{0,lim}','Fontsize',12,'FontName','Times New Roman');
set(gcf, 'Position',  [00, 400, 380, 270])
set(gca,'YLim',[0 1.25])
set(gca,'Fontsize',12,'FontName','Times New Roman')
legend boxoff
ss = '(d)  {\itf_{lf}} /{\itf''_c} < 0.08';
annotation('textbox', [0.18, 0.82, 0.1, 0.1], 'String', ss,'Fontsize',12,'FontName','Times New Roman','EdgeColor','none');
annotation('textbox', [0.5, 0.4, 0.3, 0.1], 'String', sss,'Fontsize',12,'FontName','Times New Roman','EdgeColor','none');
annotation('textbox', [0.65, 0.82, 0.3, 0.1], 'String', sRsq,'Fontsize',12,'FontName','Times New Roman','EdgeColor','none');

ax = gca;
% ax.XTick = [0 0.5 1 2 3 4];
ax.YTick = [0 0.25 0.5 0.75 1 1.25];

