% This program corresponds to Section 4.4.1-4 and Appendix 4.5 of the thesis.
% All of the MATLAB-generated results in these sections were created from this program.
% Note that the values obtained in the finite element validation of Section 4.4.4 were also obtained from this program
clc, clear, close all

%% 2-bay example - lateral loads (Appendix A4.5)
n = 2;
H = 7.315;
E0 = 2E11;
fy0 = 3.5E8;

% Column properties
Ec = E0*ones(1,n+1);
Ic = [39.9e-6 63.8e-6 39.9e-6];
Lc = H*ones(1,n+1);
Ac = [10100 11200 10100]*1e-6;
rl = [1 0.5 1];
w = [1 2 1];

% Beam properties
Eb = E0*ones(1,n);
Ib = 245e-6*ones(1,n);
Lb = 7.315*ones(1,n);
zL = [1 0]; zR = [0 1];
vRL = 1*ones(1,n);
[RL,RR] = Get_R_190204(Eb,Ib,Lb,zL,zR,vRL);
ru0 = Get_ru(Ec,Ic,Lc,RL,RR);  
K = Get_K(ru0,rl);           
KLr = K.*Lc./sqrt(Ic./Ac)  
Lr = Lc./sqrt(Ic./Ac)
Py = fy0*Ac;

% Proportional loading
l_vec = [0 0.01 0.02 0.05 0.1];
P_vec = 0:1e3:8e5;
d0 = 0; % In case I want to activate imperfections
D0 = 0; % In case I want to activate imperfections
worst_case = 1; % In case I want to activate imperfections
Kb = 0; % In case I want to activate bracing
n_seg = 5; % Set to more if using out-of-plumbness
plot_data_prop = zeros(length(P_vec),2,length(l_vec));
for H_case = 1:2 % Whether or not to include Q effect on N
    for l_case = 1:length(l_vec)
        for P_case = 1:length(P_vec)
            P = P_vec(P_case);
            Loads = P.*w;
            Q = l_vec(l_case)*sum(Loads);
            for i=1:n+1 % Inelastic modify
                Ec(i) = Get_E_IE(Loads(i),Ac(i),E0,fy0);
            end
            HQ = Q.*[-0.1715 0.0391 0.1325];
            ru = Get_ru(Ec,Ic,Lc,RL,RR);  

            % With HQ
            if H_case == 1
                N = Loads + HQ;
            else % Without HQ
                N = Loads;
            end
            
            phi = Get_phi(N,Ec,Ic,Lc);
            beta = Get_beta(phi,ru,rl);
            chi = Get_chi(phi,ru,rl);
            [D,ST,C1,C2,yi,ymax] = Get_Deflection(N,Ec,Ic,Lc,ru,rl,beta,chi,d0,D0,worst_case,Q,Kb,n_seg);
            if isinf(D)
                D = 10000;
            end
            plot_data_prop(P_case,H_case,l_case) = D;
        end
    end
end
figure('Renderer', 'painters', 'Position', [10 10 500 250])
set(gca, 'FontName', 'Times New Roman');
legend_titles = {'\lambda = 0.00 (WH)';'\lambda = 0.01 (WH)';'\lambda = 0.02 (WH)';'\lambda = 0.05 (WH)';'\lambda = 0.10 (WH)';
                 '\lambda = 0.00 (NH)';'\lambda = 0.01 (NH)';'\lambda = 0.02 (NH)';'\lambda = 0.05 (NH)';'\lambda = 0.10 (NH)'};
for H_case=1:2
    for l_case = 1:length(l_vec)
        line_style = length(legend_titles);
        line_color = (1-(1-l_case/6))*[1 1 1];
        if H_case == 1
            str_style = '-';
        else
            str_style = '--';
        end
        line_width = 0.5*H_case^2;
        y_vec = plot_data_prop(:,H_case,l_case);
        plot(y_vec,P_vec,'Color',line_color,'LineStyle',str_style,'LineWidth',line_width)
        hold on
        grid on
    end
end
axis([0 1 0 max(P_vec)]);
lgnd = legend(legend_titles,'location','eastoutside');
ylabel('Base Gravity Load, \itP \rm[N]')
xlabel('Inter-Storey Displacement \Delta [m]')
temp = [lgnd; lgnd.ItemText];
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');
set(gca, 'FontName', 'Times New Roman');
diff_HQ = 1-plot_data_prop(:,1,:) ./ plot_data_prop(:,2,:);
figure('Renderer', 'painters', 'Position', [10 10 550 300])
set(gca, 'FontName', 'Times New Roman');
legend_titles = {'\lambda = 0.01';'\lambda = 0.02';'\lambda = 0.05';'\lambda = 0.10'};
Pu = 765792; % From Excel
for l_case = 2:length(l_vec)
    line_style = length(legend_titles);
    line_color = (1-(1-l_case/6))*[1 1 1];
    str_style = '--';
    line_width = 0.5;
    y_vec = diff_HQ(:,1,l_case);
    for i=1:length(y_vec)
        if i>1 && y_vec(i) == 0
            y_vec(i) = 0/0; % Ad hoc get-around for when D is set to 10000
        end
    end
    plot(P_vec/Pu,y_vec,'Color',line_color,'LineStyle',str_style,'LineWidth',line_width)
    hold on
    grid on
end
axis([0 1 -0.01 0.01]);
lgnd = legend(legend_titles,'location','nw');
xlabel('Base Gravity Load, \itP/P_{sw}')
ylabel('Error in \Delta due to neglecting \itH_{Q,i}')
temp = [lgnd; lgnd.ItemText];
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');
set(gca, 'FontName', 'Times New Roman');

% Plot Psw versus Kb
Kb = logspace(0,6,1000);
Psw_data = zeros(1,length(Kb));
Pb = [2467187.15815262 2285756.97893318 2467187.15815262]; % From Excel
Pu_max = min(Pb./w); %lean-on columns
for K_case = 1:length(Kb)
    Psw_data(K_case) = Get_Psw_IE(E0,fy0,Ac,Ic,Lc,ru,rl,Kb(K_case),w,1.142e6);
    if (Psw_data(K_case) > Pu_max)
        Psw_data(K_case) = Pu_max;
    elseif (Psw_data(K_case)<0)
        Psw_data(K_case) = 0/0;
    end
end
figure('Renderer', 'painters', 'Position', [10 10 450 215])
set(gca, 'FontName', 'Times New Roman');
hold on
grid on
plot(Kb,Psw_data/Pu_max,'k');
xlabel('Bracing stiffness, \it{K_{br}} \rm[N/m]')
ylabel('Sway Load, \it{P_{sw}/P_u}')
axis([0 1e6 0.6 1]);
% ----- end of lateral load example


%% 4-bay example - Proportional Loading with Imperfections (Section 4.4.1-4)
n = 4;                              % # bays
H = 7.315;                          % Height of storey, m
E0 = 2E11;                          % Ambient modulus, N/m2
fy0 = 3.5E8;                        % Ambient yield stress, N/m2   

% Column properties
Ec = E0     *ones(1,n+1);                       % Modulus, columns, N/m2
Ic = [129e-6 54.7e-6 54.7e-6 54.7e-6 129e-6];   % Inertia, columns, m4
Lc = H      *ones(1,n+1);                       % Length of columns, m
Ac = [7610 9280 9280 9280 7610]*1E-6;           % Area of column, m2
rl = [1 0 0 0 1];                               % End fixity factors, lower end columns

% Beam properties
Eb = E0     *ones(1,n);             % Modulus, beams, N/m2
Ib = 2.45E-4*ones(1,n);             % Inertia, beams, m4
Lb = 7.315  *ones(1,n);             % Length of beams, m
zL = [0.9 0 0 0];                         % End fixity factors, beam left connections
zR = [0 0 0 0.9];                         % End fixity factors, beam right connections

% Inter-dependent properties
vRL = 1*ones(1,n);                              % Beam end rotation ratio
[RL,RR] = Get_R_190204(Eb,Ib,Lb,zL,zR,vRL);  % End rigidities, beams, Nm/rad
ru = Get_ru(Ec,Ic,Lc,RL,RR);        % Upper end fixity factors
K = Get_K(ru,rl);                   % Effective length factors
KLr = K.*Lc./sqrt(Ic./Ac)           % Slenderness ratio
Lr = Lc./sqrt(Ic./Ac)
Py = fy0*Ac                         % Squash load, elsatic, N
Pb = pi^2*Ec.*Ic./(K.*Lc).^2        % Buckling load, elastic, N
Pb./Py
Pu = min(Py,Pb);                    % Upper bound load, elastic, N
gamma = Get_gamma(ru,rl);           % Column gamma values
%S1 = Get_S1(Ec,Ic,Lc,ru,rl);       % Column first order stiffness, N/m

%% Parametric studies - D0 and d0
w = [1 2 2 2 1];                        % Gravity load factor (proportional)
Pmax = min(Pu./w);                  % Maximum base gravity load, N
G_vec = linspace(0,Pmax*0.999,10000); % Gravity loads, N

%% Effects of D0 and d0
Kb = 1E4;
D0 = [0 H/500 H/400 H/300 H/200 H/100];
Psw = Get_Psw(Ec,Ic,Lc,ru,rl,Kb,w);
Q = 0;
labels_D0 = [{'No Imperfections'}
            {'\Delta_0 = \it{H}\rm/500'}
            {'\Delta_0 = \it{H}\rm/400'}
            {'\Delta_0 = \it{H}\rm/300'}
            {'\Delta_0 = \it{H}\rm/200'}
            {'\Delta_0 = \it{H}\rm/100'}];
d0 = [0 -H/100 -H/500 H/500 H/100];
labels_d0 = [{'\delta_0 = 0'}
            {'\delta_0 = -\it{H}\rm/100'}
            {'\delta_0 = -\it{H}\rm/500'}
            {'\delta_0 = +\it{H}\rm/500'}
            {'\delta_0 = +\it{H}\rm/100'}];
plot_data_D0 = zeros(length(G_vec),length(D0));
plot_data_d0 = zeros(length(G_vec),length(d0));
plot_data_d0_ymax = zeros(length(G_vec),length(d0));
for G_case = 1:length(G_vec)
    
    % Gravity loading properties
    G = w.*G_vec(G_case);   % (nx1) Column gravity loads in case
    P = G;                  % (nx1) Column axial loads in case
    phi = Get_phi(P,Ec,Ic,Lc);
    beta = Get_beta(phi,ru,rl);
    chi = Get_chi(phi,ru,rl);
    
    %  Out-of-plumbness imperfection
    for D_case = 1:length(D0)
        [D,ST,C1,C2,yi,ymax] = Get_Deflection(P,Ec,Ic,Lc,ru,rl,beta,chi,0,D0(D_case),0,Q,Kb,5);
        if (D == 1/0)
            D = 100;
        end
        plot_data_D0(G_case,D_case) = D;  % sumS = f(P,D0)
    end
    
%  Out-of-straightness imperfection
    for d_case = 1:length(d0)
        [D,ST,C1,C2,yi,ymax] = Get_Deflection(P,Ec,Ic,Lc,ru,rl,beta,chi,d0(d_case),0,0,Q,Kb,500);
        if (D == 1/0)
            D = 100;
        end
        plot_data_d0(G_case,d_case) = D;  % sumS = f(P,d0)
        plot_data_d0_ymax(G_case,d_case) = ymax;  % sumS = f(P,d0)
    end
end

% Figure 1 - Plot P vs D+D0 for different cases of D0 
figure('Renderer', 'painters', 'Position', [10 10 525 250])
set(gca, 'FontName', 'Times New Roman');
hold on
for D_case=1:length(D0)
    line_color = D_case/length(D0)/1.5*[1 1 1];
    line_style = '-.';
    if (D_case==1)
        line_style = '-';
    end
    plot(plot_data_D0(:,D_case)+D0(D_case),G_vec/Psw,line_style,'Color',line_color);
end
ylabel('Base Gravity Load, \it{P/P_{sw}}');
xlabel('Total Inter-storey Displacement, {\Delta}+{\Delta_0} \rm[m]');
lgnd = legend(labels_D0,'Location','best');
axis([0 0.5 0 1]);
temp = [lgnd; lgnd.ItemText];
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');

% Figure 2 - Plot P vs D for different cases of d0 
figure('Renderer', 'painters', 'Position', [10 10 475 270])
set(gca, 'FontName', 'Times New Roman');
hold on
for d_case=1:length(d0)
    line_color = [0 0 0;
                  0.5 0.5 0.5;
                  0 0 0
                  0 0 0
                  0.5 0.5 0.5];
    line_style = '-.';
    if (d_case==1)
        line_style = '-';
    elseif (d_case < 4)
        line_style = ':';
    end
    plot(plot_data_d0(:,d_case),G_vec/Psw,line_style,'Color',line_color(d_case,:));
end
ylabel('Base gravity load, \it{P/P_{sw}}');
xlabel('Inter-storey displacement, {\Delta} \rm[m]');
lgnd = legend(labels_d0,'Location','sw');
axis([-0.3 0.300001 0 1]);
temp = [lgnd; lgnd.ItemText];
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');

% Figure 3 - Plot P vs ymax for difference cases of D0
figure('Renderer', 'painters', 'Position', [10 10 600 350])
set(gca, 'FontName', 'Times New Roman');
hold on
for d_case=1:5%3
    line_color = [0 0 0;
                  0.5 0.5 0.5;
                  0 0 0
                  0 0 0
                  0.5 0.5 0.5];
    line_style = '-.';
    if (d_case==1)
        line_style = '-';
    elseif(d_case<4)
        line_style = ':';
    end
    % Plot only the maximum deflection from both the negative and positive
    % cases since it's symmetrical
    %if (d_case == 1)
    %    max_abs = abs(plot_data_d0_ymax(:,d_case));
    %else
    %    max_abs = abs(max(abs(plot_data_d0_ymax(:,d_case)),abs(plot_data_d0_ymax(:,length(d0)-d_case+2))));
    %end
    %plot(max_abs,G_vec/Psw,line_style,'Color',line_color(d_case,:));
    plot(plot_data_d0_ymax(:,d_case),G_vec/Psw,line_style,'Color',line_color(d_case,:));
end
ylabel('Base gravity load, \it{P/P_{sw}}');
xlabel('Maximum deflection, \it\delta_{max} \rm[m]');
lgnd = legend(labels_d0,'Location','sw');
axis([-0.3 0.300001 0 1]);
temp = [lgnd; lgnd.ItemText];
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');

% Figure 4 - Plot Psw versus Kb
Kb = logspace(0,5.5,100);
Psw_data = zeros(1,length(Kb));
Pu_max = min(Pb./w); %lean-on columns
for K_case = 1:length(Kb)
    Psw_data(K_case) = Get_Psw(Ec,Ic,Lc,ru,rl,Kb(K_case),w);
    if (Psw_data(K_case) > Pu_max)
        Psw_data(K_case) = Pu_max;
    end
end
figure('Renderer', 'painters', 'Position', [10 10 400 215])
set(gca, 'FontName', 'Times New Roman');
hold on
grid on
plot(Kb,Psw_data/Pu_max,'k');
xlabel('Bracing stiffness, \it{K_{br}} \rm[N/m]')
ylabel('Elastic Sway Load, \it{P_{sw}/P_u}')
axis([0 max(Kb) 0.9 1]);
