% This program corresponds to Section 3.4.1 of the thesis.
% The capacity of a lean-on column is investigated with/without considering
%    the effects of shear deformations
clc, clear, close all

%% Lean-on Frame Analysis (Example 3.4.1)
% Page 169 of research journal, Part 2 - February 14, 2019
% Column 1 is a supporting column for an arbitrary system of lean-on columns
% As any number of lean-on columns of any length can be placed in the
% system with the supporting column, any lateral stiffness can technically
% be diminished to zero via a series of lean-on columns as long as the
% lean-on columns do not rotationally buckle
E0 = 2e11;
Fy = 3.5E8;
vRL = 0.3;
k = 0.44;
Ix = 445e-6;        % W460x97, strong axis bending
A = 12300e-6;       % W460x97
rx = sqrt(Ix/A);
Rl = 0; rl = 0;
Ru = 1/0; ru = 1;
%K = 0.699;          % Use fixed-pinned

% Parameters of the study and output data
Lr_vec = 5:1:100;
PPsw_vec = 0:0.01:1;
data_L = zeros(length(Lr_vec),1);
data_Psw_N = zeros(length(Lr_vec),1);
data_Psw_dy = zeros(length(Lr_vec),1);
data_Psw_phi = zeros(length(Lr_vec),1);
data_P = zeros(length(Lr_vec),length(PPsw_vec));
data_E = zeros(length(Lr_vec),length(PPsw_vec));
data_S1_N = zeros(length(Lr_vec),length(PPsw_vec));
data_S1_dy = zeros(length(Lr_vec),length(PPsw_vec));
data_S1_phi = zeros(length(Lr_vec),length(PPsw_vec));

% Populate data
for Lr_case = 1:length(Lr_vec) % for varying lengths
    Lr = Lr_vec(Lr_case);
    L = Lr*rx;%/K;
    Psw_N = Get_Pcr_IE_robust(E0,Fy,vRL,k,Ix,L,A,Rl,Ru,0,0); % no shear
    Psw_phi = Get_Pcr_IE_robust(E0,Fy,vRL,k,Ix,L,A,Rl,Ru,1,0); %phi (Engesser)
    Psw_dy = Get_Pcr_IE_robust(E0,Fy,vRL,k,Ix,L,A,Rl,Ru,1,1); %dydx (Haringx)
    data_L(Lr_case) = L;
    data_Psw_N(Lr_case) = Psw_N;
    data_Psw_phi(Lr_case) = Psw_phi;
    data_Psw_dy(Lr_case) = Psw_dy;
    for PPsw_case = 1:length(PPsw_vec)
        PPsw = PPsw_vec(PPsw_case); % between 0 and 100% of the Psw load
        P = PPsw*Psw_dy;
        E = Get_E_IE(P,A,E0,Fy);
        G = E0/(2*(1+vRL)+1.5*(E0/E-1)); % Assume buckling is occuring
        kAG = k*A*G;
        data_P(Lr_case,PPsw_case) = P;
        if (P == 0)
            P = 0.001;
        end
        [data_S1_N(Lr_case,PPsw_case),~] = Get_S_Shear(P,E,Ix,L,rl,ru,kAG,0,0);
        if PPsw == 1 % Ignore epsilons            
            data_S1_phi(Lr_case,PPsw_case) = 0;
            data_S1_dy(Lr_case,PPsw_case) = 0;
        else
            [data_S1_phi(Lr_case,PPsw_case),~] = Get_S_Shear(P,E,Ix,L,rl,ru,kAG,1,0);
            [data_S1_dy(Lr_case,PPsw_case),~] = Get_S_Shear(P,E,Ix,L,rl,ru,kAG,1,1);
        end
    end
end

% Process data and plot figures
% Assuming that there is negligible difference between Haringx and Engesser, use the Haringx data
diff_newtons = data_S1_N-data_S1_dy; % N/m, change in lateral stiffness due to shear
diff_Pcr_pct = (data_S1_dy./data_S1_N-1)*100; % percent, reduction in critical load due to shear (since it's proportional to S1)
diff_Psw = (data_Psw_phi - data_Psw_dy)./data_Psw_phi*100; % percent, difference in critical sway load from shear assumptions
figure('Renderer', 'painters', 'Position', [10 10 550 330]) % ---- Figure 1
PPsw_vec_plot = [0 60 80 100]+1;
for i=1:length(PPsw_vec_plot)
    semilogy(Lr_vec,diff_newtons(:,PPsw_vec_plot(i)),'-','Color',1*[1 1 1]*(i-1)/length(PPsw_vec_plot))
    hold on
end
grid on
axis([0 max(Lr_vec) 0 1e8]);
set(gca, 'FontName', 'Times New Roman');
xlabel('Slenderness ratio, \itL/r')
ylabel('Reduction in \itS_{n+\rm1} \rmdue to Shear Deformations \rm[N/m]')
legend_entries = [{'\itP_{n+\rm1}/P_{sw} \rm= 0.0'}
                {'\itP_{n+\rm1}/P_{sw} \rm= 0.6'}
                {'\itP_{n+\rm1}/P_{sw} \rm= 0.8'}
                {'\itP_{n+\rm1}/P_{sw} \rm= 1.0'}];
lgnd = legend(legend_entries,'Location','ne');
temp = [lgnd; lgnd.ItemText];
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');
figure('Renderer', 'painters', 'Position', [20 20 550 330]) % ---- Figure 2
set(gca, 'FontName', 'Times New Roman');
hold on
Lr_vec_plot = [5 10 20 40 75 100]-Lr_vec(1)+1;
for i=length(Lr_vec_plot):-1:1
    colour_plot = 1*[0.85 0.85 0.85]*(i-1)/length(Lr_vec_plot);
    linewidth_plot = 0.5;%1+2*(i-1)/length(KLr_vec_plot);
    linetype_plot = '--';
    if i == 1
        linetype_plot = '-';
    elseif i ==2
        linetype_plot = '-.';
    elseif i ==3
        linetype_plot = '--';
    elseif i ==4
        linetype_plot = ':';
    elseif i == 5
        linetype_plot = '-.';
    end
    plot(PPsw_vec,diff_Pcr_pct(Lr_vec_plot(i),:),linetype_plot,'Color',colour_plot,'LineWidth',linewidth_plot)
end
%ytickformat(axes, 'percentage');
axis([0 1 -100 0]);
dummy = cellstr(num2str(get(gca,'ytick')'));
pct = char(ones(size(dummy,1),1)*'%');
new_yticks = [char(dummy),pct];
set(gca,'yticklabel',new_yticks);
grid on
xlabel('Axial Load on Supporting Column, \itP_{n+\rm1}/P_{sw}')
ylabel('Reduction in \itP_{cr} \rm due to shear deformations')
legend_entries = [{'\itL/r \rm= 100'}
                {'\itL/r \rm= 75'}
                {'\itL/r \rm= 40'}
                {'\itL/r \rm= 20'}
                {'\itL/r \rm= 10'}
                {'\itL/r \rm= 5'}];
lgnd = legend(legend_entries,'Location','southwest');
temp = [lgnd; lgnd.ItemText];
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');
figure('Renderer', 'painters', 'Position', [30 30 450 230]) % ---- Figure 3
set(gca, 'FontName', 'Times New Roman');
hold on
plot(Lr_vec,diff_Psw,'k-')
%ytickformat(axes, 'percentage');
axis([0 max(Lr_vec) 0 0.06]);
dummy = cellstr(num2str(get(gca,'ytick')'));
pct = char(ones(size(dummy,1),1)*'%');
new_yticks = [char(dummy),pct];
set(gca,'yticklabel',new_yticks);
grid on
xlabel('Slenderness ratio, \itL/r')
ylabel('Difference in \itP_{sw}')

%% Function for finding Pcr (no inelasticity consideration)
function [Pcr,success] = Get_Pcr(P_trial,Ec,Ic,Lc,rl,ru,kAG,bool_shear,shear_assumption)
    P_curr = P_trial;
    max_iter = 200;
    n_iter = 1;
    dP = 1;
    tol = 1e-6;
    success = 0;
    while n_iter < max_iter
        [~,S_curr] = Get_S_Shear(P_curr,Ec,Ic,Lc,rl,ru,kAG,bool_shear,shear_assumption);
        if (abs(S_curr) < tol)
            success = 1;
            break;
        end
        P_step = P_curr+dP;
        [~,S_delta] = Get_S_Shear(P_step,Ec,Ic,Lc,rl,ru,kAG,bool_shear,shear_assumption);
        dSdP = (S_delta - S_curr)/dP;
        P_curr = P_curr-S_curr/dSdP;
        n_iter = n_iter + 1;
    end
    Pcr = -1;
    if (success == 1)
        Pcr = P_curr;
    end
end


%% Function for finding Pcr (with inelasticity consideration)
function [Pcr] = Get_Pcr_IE_robust(E0,Fy,v,k,I,L,A,Rl,Ru,bool_shear,shear_assumption)
    Py = Fy*A;
    n_tries = 20;
    min_Pcr = Py*0.85;
    found = false;
    for i = 1:n_tries
        P_trial = 0.85*Py*(i/n_tries)^1.5;
        Pcr_i = Get_Pcr_IE(P_trial,E0,Fy,v,k,I,L,A,Rl,Ru,bool_shear,shear_assumption);
        if (Pcr_i > 0 && Pcr_i < min_Pcr)
            min_Pcr = Pcr_i;
            found = true;
        end
    end
    if (found == true)
        Pcr = min_Pcr;
    else
        Pcr = -1;
    end
end
function [Pcr,success] = Get_Pcr_IE(P_trial,E0,Fy,v,k,I,L,A,Rl,Ru,bool_shear,shear_assumption)
    Py = Fy*A;
    P_curr = P_trial;
    max_iter = 30;
    n_iter = 1;
    dP = 1;
    tol = 1e0;
    success = 0;
    while n_iter < max_iter

        % Update end fixity factor with new E
        if (P_curr/Py) > 0.3333
            E = E0 * (-7.38 * P_curr / Py * log10(P_curr / Py / 0.85));
        elseif (P_curr > 0.85*Py)
            E = 0;
        else
            E = E0;
        end
        ru = 1 / (1 + 3 * E * I / (L * Ru));
        rl = 1 / (1 + 3 * E * I / (L * Rl));
        G = E0/(2*(1+v)+1.5*(E0/E-1)); 
        kAG = k*A*G;
        [~,S_curr] = Get_S_Shear(P_curr,E,I,L,rl,ru,kAG,bool_shear,shear_assumption);
        if (abs(S_curr) < tol)
            success = 1;
            break;
        end
        P_step = P_curr+dP;
        if (P_step/Py) > 0.3333
            E = E0 * (-7.38 * P_step / Py * log10(P_step / Py / 0.85));
        elseif (P_curr > 0.85*Py)
            E = 0;
        else
            E = E0;
        end
        ru = 1 / (1 + 3 * E * I / (L * Ru));
        rl = 1 / (1 + 3 * E * I / (L * Rl));
        G = E0/(2*(1+v)+1.5*(E0/E-1)); 
        kAG = k*A*G;
        [~,S_delta] = Get_S_Shear(P_step,E,I,L,rl,ru,kAG,bool_shear,shear_assumption);
        dSdP = (S_delta - S_curr)/dP;
        P_curr = P_curr-S_curr/dSdP;
        n_iter = n_iter + 1;
    end
    Pcr = -1;
    if (success == 1)
        Pcr = P_curr;
    end
end
