% This program corresponds to the examples in Sections 3.3.5, 3.3.6 and 3.4.2 of the thesis
% It was used to obtain the results and figures for the corresponding examples
%   relating to axial beam deformations
clc, clear, close all

% Example 1: Section 3.3.5 (Axial example #1)
% Example 2: Section 3.4.2 (Lean-on frame)
% Example 3: Section 3.3.6 (Axial example #2)
example = 1; % Set to 1, 2 or 3

%% Example 1 (Section 3.3.5): Adopted from Xu (2003), A NLP approach for evaluating...
if example == 1
disp('*************** Example 1 ***************')
m = 4;      % Number of bays
n = m + 1;  % Number of columns
E0 = 2E11;  % Modulus of elasticity (Pa)
fy = 3.5E8; % Yield stress (Pa)
bool_inel = 1;  % Consider inelasticity

% Beams
Lb = 7.315*ones(1,m);   % Length of beam (m)
Ib = 245E-6*ones(1,m);  % Inertia of beam (m4)
Ab = 8580E-6*ones(1,m); % Area of beam (m2)
Eb = E0*ones(1,m);      % Elasticity of beam (Pa)
kb = E0.*Ab./Lb;        % Axial stiffness of beam (Pa)
zL = [0 0 0 0];         % Connection fixities, left ends
zR = [0 0 0 0];         % Connection fixities, right ends
vRL = [1 1 1 1];          % End connection rotation ratios
[RL,RR] = Get_R_190204(Eb,Ib,Lb,zL,zR,vRL); % Rotational rigidity of connection to column (Nm)

% Columns
Lc = 4.877*ones(1,n);                   % Length of column (m)
Ic = [129 34.1 34.1 34.1 129]*1e-6;     % Inertia of column (m4)
Ac = [7610 4570 4570 4570 7610]*1e-6;   % Area of column (m2)
ru0 = Get_ru(E0*ones(1,n),Ic,Lc,RL,RR); % End fixity factor, upper
rl0 = [0 1 1 1 0];                      % End fixity factor, lower
S1 = Get_S1(E0*ones(1,n),Ic,Lc,ru0,rl0);% First order lateral stiffness (N/m)
lambda = [1 2 2 2 1];                   % Proportional loading factors
for i=1:n
    Pu(i) = Iterate_Pu(E0,fy,0.3,Lc(i),Ic(i),Ac(i),ru0(i),rl0(i),0,1,"NA",0);
    %Pu(i) = Get_Pu_inel(E0,fy,Ac(i),Ic(i),Lc(i),ru0(i),rl0(i)); % Rotational buckling load
end
Gu = min(Pu./lambda);

Kb_vec = logspace(0,7,100);
data_2001 = zeros(1,length(Kb_vec));
data_axial= zeros(1,length(Kb_vec));
data_diff = zeros(1,length(Kb_vec));
for Kb_case = 1:length(Kb_vec)
    Kb = [0 1 0 0 1]*Kb_vec(Kb_case);                                           % Additional lateral bracing (N/m)
    Gcr_2001  = Get_Psw_2001(E0,fy,Ac,Ic,Lc,ru0,rl0,Kb,lambda,bool_inel);        % Sway load, without axial (Xu, 2001)
    Gcr_axial = Get_Psw_axial(E0,fy,Ac,Ic,Lc,ru0,rl0,Kb,lambda,kb,bool_inel);    % Sway load, with axial (Ma, Zhuang & Xu, 2019)

    % Bug fix - when it can't find the correct critical load because the
    % stiffness decreases too sharply near the rotational buckling load,
    % just take the rotational buckling load
    if (Gcr_2001 > Gu)
        Gcr_2001 = Gu;
    end
    if (Gcr_axial > Gu)
        Gcr_axial = Gu;
    end
    Diff = Gcr_axial/Gcr_2001-1;
    
    data_2001(Kb_case) = Gcr_2001;
    data_axial(Kb_case) = Gcr_axial;
    data_diff(Kb_case) = Diff;
end

% Bug fix - if Newton's method fails to converge, remove the data
while i <= length(data_2001)
    if data_2001(i) < 0 || data_axial(i) < 0
        for j=i:length(data_2001)-1
            data_2001(j) = data_2001(j+1);
            data_axial(j) = data_axial(j+1);
            data_diff(j) = data_diff(j+1);
            Kb_vec(j) = Kb_vec(j+1);
        end
        data_2001 = data_2001(1:length(data_2001)-1);
        data_axial = data_axial(1:length(data_axial)-1);
        data_diff = data_diff(1:length(data_diff)-1);
        Kb_vec = Kb_vec(1:length(Kb_vec)-1);
        i = i-1;
    end
    i = i+1;
end

% Plot the figure
figure('Renderer', 'painters', 'Position', [10 10 550 275]) 
labels_all= [{'Neglect Axial Deformations'}
            {'With Axial Deformations'}];
set(gca, 'FontName', 'Times New Roman');
set(gca, 'XScale', 'log');
hold on
grid on
semilogx(Kb_vec/1000,data_2001/1000,'-','Linewidth',2.5,'Color',[0.7 0.7 0.7])
semilogx(Kb_vec/1000,data_axial/1000,'k--')
xlabel('Bracing stiffness \itK_{br} \rm[kN/m]');
ylabel('Critical load, \itP_{cr} \rm[kN]');
axis([0 1e7/1000 0 7e5/1000])
lgnd = legend(labels_all,'Location','best');
temp = [lgnd; lgnd.ItemText];
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');

% Plot the percent difference
figure('Renderer', 'painters', 'Position', [20 20 550 330]) 
set(gca, 'FontName', 'Times New Roman');
set(gca, 'XScale', 'log');
hold on
grid on
semilogx(Kb_vec,data_diff*100,'k-')
xlabel('Bracing stiffness \itK_{br} \rm[N/m]');
ylabel('% Difference in \itP_{cr}');
axis([0 1e7 -0.05 0])


%% Example 2 (Section 3.4.2): n-bay lean-on frame
elseif example == 2
disp('*************** Example 2 ***************')
m = 1;      % Number of bays
n = m + 1;  % Number of columns
E0 = 2E11;  % Modulus of elasticity (Pa)
fy = 3.5E8; % Yield stress (Pa)
bool_inel = 1;  % Consider inelasticity
Kb = zeros(1,n);     % Additional lateral bracing (N/m)

% Beams
Lb = 8.534*ones(1,m);   % Length of beam (m)
Ib = 246E-6*ones(1,m);  % Inertia of beam (m4)
Ab = 8600E-6*ones(1,m); % Area of beam (m2)
Eb = E0*ones(1,m);      % Elasticity of beam (Pa)
kb = E0.*Ab./Lb;        % Axial stiffness of beam (Pa)
zL = [0];         % Connection fixities, left ends
zR = [1];         % Connection fixities, right ends
vRL = [1];          % End connection rotation ratios
[RL,RR] = Get_R_190204(Eb,Ib,Lb,zL,zR,vRL); % Rotational rigidity of connection to column (Nm)

% Columns
Lc = [2.439 7.315];             % Length of column (m)
Ic = [8620 8620]*1e-6;          % Inertia of column (m4)
Ac = [74300 74300]*1e-6;        % Area of column (m2)
ru0 = Get_ru(E0*ones(1,n),Ic,Lc,RL,RR);    % End fixity factor, upper
rl0 = [0 1];                     % End fixity factor, lower

% Loading
lambda = [1 0.8];  % Proportional loading factors, Case 1 (keep P2 = 15 MN)
Gcr_2001  = Get_Psw_2001(E0,fy,Ac,Ic,Lc,ru0,rl0,Kb,lambda,bool_inel)        % Sway load, without axial (Xu, 2001)

lambda = [1 0.8];  % Proportional loading factors, Case 2 (keep P2 = 15 MN)
Gcr_axial = Get_Psw_axial(E0,fy,Ac,Ic,Lc,ru0,rl0,Kb,lambda,kb,bool_inel)    % Sway load, with axial (Ma, Zhuang & Xu, 2019)
Diff = Gcr_axial/Gcr_2001-1

%% New Example 2, Lean-on Frame
disp('*************** New Example 2 ***************')
E0 = 2E11;  % Modulus of elasticity (Pa)
fy = 3.5E8; % Yield stress (Pa)
A_vec = 1e-6*[1630 8600 22200 50000 99800];
k_vec = zeros(1,length(A_vec));
max_bays = 15;
data = zeros(length(A_vec),max_bays);
for A_case = 1:length(A_vec)
    
    for m = 1:max_bays % Number of bays
        n = m + 1;  % Number of columns
        bool_inel = 1;  % Consider inelasticity
        Kb = zeros(1,n);     % Additional lateral bracing (N/m)

        % Beams
        Lb = 7.315*ones(1,m); % Length of beam (m)
        Ib = 99999E-6*ones(1,m);  % Inertia of beam (m4) ** doesn't affect results!!
        Ab = A_vec(A_case)*ones(1,m); % Area of beam (m2)
        Eb = E0*ones(1,m);      % Elasticity of beam (Pa)
        kb = E0.*Ab./Lb;        % Axial stiffness of beam (Pa)
        zL = zeros(1,m);        % Connection fixities, left ends
        zR = zeros(1,m);        % Connection fixities, right ends
        vRL = ones(1,m);          % End connection rotation ratios
        [RL,RR] = Get_R_190204(Eb,Ib,Lb,zL,zR,vRL); % Rotational rigidity of connection to column (Nm)

        % Columns
        Lc = 7.315*ones(1,n);           % Length of column (m)
        Ic = 8620e-6*ones(1,n);         % Inertia of column (m4)
        Ac = 74300e-6*ones(1,n);        % Area of column (m2)
        ru0 = zeros(1,n);               % End fixity factor, upper
        rl0 = zeros(1,n);               % End fixity factor, lower
        rl0(n) = 1;
        S1 = Get_S1(E0*ones(1,n),Ic,Lc,ru0,rl0);
        k_vec(A_case) = Eb(1)*Ab(1)/Lb(1)/S1(n);

        % Loading
        lambda = ones(1,n);
        lambda(n) = 0;
        Gcr_2001  = Get_Psw_2001(E0,fy,Ac,Ic,Lc,ru0,rl0,Kb,lambda,bool_inel);        % Sway load, without axial (Xu, 2001)
        Gcr_axial = Get_Psw_axial(E0,fy,Ac,Ic,Lc,ru0,rl0,Kb,lambda,kb,bool_inel);    % Sway load, with axial (Ma, Zhang & Xu, 2019)
        SP_2001 = Gcr_2001*m;
        SP_axial = Gcr_axial*m;
        Diff = SP_axial/SP_2001-1;
        if (SP_2001<0 || SP_axial<0)
            Diff = 0;
        end
        data(A_case,m)=Diff*100;
    end
end

figure('Renderer', 'painters', 'Position', [30 30 550 330])  % ---- Figure 1
for A_case=1:length(A_vec)
    plot(1:max_bays,data(A_case,:),'-x','Color',1*[1 1 1]*(A_case-1)/length(A_vec))
    hold on
end
grid on
axis([5 max_bays -100 0]);
set(gca, 'FontName', 'Times New Roman');
xlabel('Number of bays')
ylabel('Reduction in \itP_{cr} \rmdue to Beam Axial Deformations \rm')
legend_entries = [{'W100×19, \itA\rm = 1630 mm^2'}
                {'W410×67, \itA\rm = 8600 mm^2'}
                {'W610×174, \itA\rm = 22200 mm^2'}
                {'W840×392, \itA\rm = 50000 mm^2'}
                {'W920×784, \itA\rm = 99800 mm^2'}];
lgnd = legend(legend_entries,'Location','sw');
temp = [lgnd; lgnd.ItemText];
dummy = cellstr(num2str(get(gca,'ytick')'));
pct = char(ones(size(dummy,1),1)*'%');
new_yticks = [char(dummy),pct];
set(gca,'yticklabel',new_yticks);
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');

%% Example 3 (Section 3.3.6):  11-bay semi-rigid, semi-braced frame
elseif example == 3
disp('*************** Example 3 ***************')
m = 11;      % Number of bays
n = m + 1;  % Number of columns
E0 = 2E11;  % Modulus of elasticity (Pa)
fy = 3.5E8; % Yield stress (Pa)
bool_inel = 1;  % Consider inelasticity

% Beams
Lb = 14.631*ones(1,m);   % Length of beam (m)
Ib = 257E-6*ones(1,m);  % Inertia of beam (m4)
Ab = 4360E-6*ones(1,m); % Area of beam (m2)
Eb = E0*ones(1,m);      % Elasticity of beam (Pa)
kb = E0.*Ab./Lb;        % Axial stiffness of beam (Pa)
%Z_val = 4e6;
%z_val = 1/(1+3*E0*Ib(1)/Z_val/Lb(1)); 
z_val = 0.25;
zL = z_val*ones(1,m);   % Connection fixities, left ends
zR = z_val*ones(1,m);   % Connection fixities, right ends
vRL = [1 1 1 1 1 1 1 1 1 1 1];    % End connection rotation ratios
[RL,RR] = Get_R_190204(Eb,Ib,Lb,zL,zR,vRL); % Rotational rigidity of connection to column (Nm)

% Columns
Lc = 4.572*ones(1,n);                       % Length of column (m)
Ie = 61.5e-6; Ii = 5.94e-6;                 % W360x122, W250x39 (weak axis)
Ae = 15500e-6; Ai = 4920e-6;    
Ic = [Ie Ii Ii Ii Ii Ii Ii Ii Ii Ii Ii Ie]; % Inertia of column (m4)
Ac = [Ae Ai Ai Ai Ai Ai Ai Ai Ai Ai Ai Ae]; % Area of column (m2)
ru0 = Get_ru(E0*ones(1,n),Ic,Lc,RL,RR);     % End fixity factor, upper
rl0 = [1 0 0 0 0 0 0 0 0 0 0 1];            % End fixity factor, lower
S1 = Get_S1(E0*ones(1,n),Ic,Lc,ru0,rl0);    % First order lateral stiffness (N/m)
lambda = [1 2 2 2 2 2 2 2 2 2 2 1];         % Proportional loading factors
for i=1:n
    Pu(i) = Iterate_Pu(E0,fy,0.3,Lc(i),Ic(i),Ac(i),ru0(i),rl0(i),0,1,"NA",0);
    %Pu(i) = Get_Pu_inel(E0,fy,Ac(i),Ic(i),Lc(i),ru0(i),rl0(i)); % Rotational buckling load
end
Gu = min(Pu./lambda);

Kb_vec = logspace(0,7,100);
data_2001 = zeros(1,length(Kb_vec));
data_axial= zeros(1,length(Kb_vec));
data_diff = zeros(1,length(Kb_vec));
for Kb_case = 1:length(Kb_vec)
    Kb = [0 0 1 0 0 0 0 0 0 0 1 0]*Kb_vec(Kb_case);                              % Additional lateral bracing (N/m)
    Gcr_2001  = Get_Psw_2001(E0,fy,Ac,Ic,Lc,ru0,rl0,Kb,lambda,bool_inel);        % Sway load, without axial (Xu, 2001)
    Gcr_axial = Get_Psw_axial(E0,fy,Ac,Ic,Lc,ru0,rl0,Kb,lambda,kb,bool_inel);    % Sway load, with axial (Ma, Zhuang & Xu, 2019)

    % Bug fix - when it can't find the correct critical load because the
    % stiffness decreases too sharply near the rotational buckling load,
    % just take the rotational buckling load
    Gu = 442.7e3; % hotfix - remove these data points
    if (Gcr_2001 > Gu)
        %Gcr_2001 = Gu; 
        Gcr_2001 = -1; %edit - just flag the data for removal
    end
    if (Gcr_axial > Gu)
        %Gcr_axial = Gu;
        Gcr_axial = -1; %edit - just flag the data for removal
    end
    Diff = Gcr_axial/Gcr_2001-1;
    
    data_2001(Kb_case) = Gcr_2001;
    data_axial(Kb_case) = Gcr_axial;
    data_diff(Kb_case) = Diff;
end

% Bug fix - if Newton's method fails to converge, remove the data
while i <= length(data_2001)
    if data_2001(i) < 0 || data_axial(i) < 0 || data_diff(i)<data_diff(max(i-1,1)) % added third criteria as hotfix
        for j=i:length(data_2001)-1
            data_2001(j) = data_2001(j+1);
            data_axial(j) = data_axial(j+1);
            data_diff(j) = data_diff(j+1);
            Kb_vec(j) = Kb_vec(j+1);
        end
        data_2001 = data_2001(1:length(data_2001)-1);
        data_axial = data_axial(1:length(data_axial)-1);
        data_diff = data_diff(1:length(data_diff)-1);
        Kb_vec = Kb_vec(1:length(Kb_vec)-1);
        i = i-1;
    end
    i = i+1;
end

% Plot the figure
figure('Renderer', 'painters', 'Position', [10 10 450 275])
labels_all= [{'Neglect Axial Deformations'}
            {'With Axial Deformations'}];
set(gca, 'FontName', 'Times New Roman');
set(gca, 'XScale', 'log');
hold on
grid on
semilogx(Kb_vec/1000,data_2001/1000,'k-')
semilogx(Kb_vec/1000,data_axial/1000,'k--')
% xlabel('Bracing stiffness \itS_{L} \rm[kN/m]');
% ylabel('Critical gravity load, \itG_{cr} \rm[kN]');
% for thesis, enable this block (different variable names)
xlabel('Bracing stiffness \itK_{br} \rm[kN/m]');
ylabel('Critical load, \itP_{cr} \rm[kN]');
% end block
axis([0 1e7/1000 0 7e5/1000])
lgnd = legend(labels_all,'Location','nw');
temp = [lgnd; lgnd.ItemText];
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');

% Plot the percent difference
figure('Renderer', 'painters', 'Position', [10 10 400 220])
set(gca, 'FontName', 'Times New Roman');
set(gca, 'XScale', 'log');
hold on
grid on
semilogx(Kb_vec/1000,data_diff*100,'k-')
xlabel('Bracing stiffness \itS_{L} \rm[kN/m]');
ylabel('% Difference in \itG_{cr}');
% for thesis, enable this block (different variable names)
xlabel('Bracing stiffness \itK_{br} \rm[kN/m]');
ylabel('% Difference in \itP_{cr}');
% end block
axis([0 1e7/1000 -1.5 0])

end