%% This program corresponds to Sections 6.3.1 and 6.3.2 of the thesis.
% It is used to analyze the fire resistance of frames containing three-segment members
% For the localized fire example (Section 6.3.1), set example = 2
% For the blast fire example (Section 6.3.2), set example = 1
% The tag (%VALIDATION) is used to denote places in the code which are modified during the corresponding validation studies
clc, clear, close all
E0 = 2E11;      % Ambient modulus (N/m2)
fy0 = 3.5E8;    % Yield stress (N/m2)
example = 2;    % Set to 1 for blast example, or 2 for localized fire example

% Problem setup
n = 2;                              % # bays
H = 4.877;                          % Height of storey, m
T0 = 20;                            % Ambient temperature, °C
Kb = 0;                             % Semi-bracing lateral stiffness, N/m

% Column properties
Ic = [1.29E-4 34.1E-6 1.29E-4];     % Inertia, columns, m4
Ac = [7610 4570 7610]/10^6;         % Area, columns, m2
Lc = H      *ones(1,n+1);           % Length of columns, m
Ll = Lc/3;  %changed later          % Lower segment length, m
Lm = Lc/3;  %changed later          % Middle segment length, m
Lu = Lc/3;  %changed later          % Upper segment length, m
rl0 = [0 0 1];                      % Ambient end fixity factors, lower end columns
m_l = zeros(1,n+1);                 % Reduction slope for ground connections, Nm/rad/°C
G = 7.5E5*[1 1 1.5];                % Prescribed gravity loads, N
bool_H = "Y";                       % Thermally restrained columns?
                                    
% Beam properties
Ib = 2.45E-4*ones(1,n);             % Inertia, beams, m4
Lb = 7.315  *ones(1,n);             % Length of beams, m
LL = Lb/3;  %changed later          % Left segment length, m
LM = Lb/3;  %changed later          % Middle segment length, m
LR = Lb/3;  %changed later          % Right segment length, m
zL0 = 0.493245689*[1 1];            % Ambient connection fixity factors, beam left connections
zR0 = 0.493245689*[1 1];            % Ambient connection fixity factors, beam right connections
m_L = 2.88e4*[1 1];                 % m-parameter of semi-rigid connection, left, Nm/rad/°C
m_R = 2.88e4*[1 1];                 % m-parameter of semi-rigid connection, right, Nm/rad/°C
vRL = [1 1];    

%% Blast damage example (Section 6.3.2)
if example == 1

% Blast Damage Scenarios
t_vec = 100:0.1:140;   % Fire duration vector % VALIDATION - reduce the increment to 0.01 to better detect regions of the plot with rotational buckling
D0 = 0;                % Column imperfections, m
Ld = 1.0;              % Delamination length at blast location, m
H = 4.877;             % Height of storey, m (**repeated since it was over-written)

% Additional problem parameters for blast analysis
Ac = [7610 4570 7610]/10^6;               	% Area, columns, m2
Dc = [1.4 1.05 1.4];                        % Heated perimeter, m
Lc = H      *ones(1,n+1);                   % Length of columns, m
Wc = [60 36 60];                            % Unit weight, kg/m
Ab = 8580/10^6*ones(1,n);                   % Area of beam, m2
Db = 1.52*ones(1,n);                        % Heated perimeter
Wb = 67*ones(1,n);                          % Unit weight, kg/m
sect_col = ["W310x60", "W200x36", "W310x60"];
sect_beam = ["W410x67", "W410x67"];

% Run blast analysis for each blast scenario (note: #16 is undamaged)
S_data = zeros(length(t_vec),16);
for B_case = 1:16 % VALIDATION - change to the specific case to check as necessary
    % 1 = Protected, 0 = Unprotected for cc and cb
    cc = ones(3,n+1);
    cb = ones(3,n);
    Ll = Lc/3; Lm = Lc/3; Lu = Lc/3; %changed later
    LL = Lb/3; LM = Lb/3; LR = Lb/3; %changed later
    buckle = 0; yielded = 0; fail = 0; % for output
    if B_case == 1, cc(1,1)=0; Ll(1) = Ld; Lm(1) = 0.5*(Lc(1)-Ld); Lu(1) = Lm(1); end % Col 1, lower
    if B_case == 2, cc(2,1)=0; Lm(1) = Ld; Ll(1) = 0.5*(Lc(1)-Ld); Lu(1) = Ll(1); end % Col 1, middle
    if B_case == 3, cc(3,1)=0; Lu(1) = Ld; Ll(1) = 0.5*(Lc(1)-Ld); Lm(1) = Ll(1); end % Col 1, upper
    if B_case == 4, cb(1,1)=0; LL(1) = Ld; LM(1) = 0.5*(Lb(1)-Ld); LR(1) = LM(1); end % Beam 1, left
    if B_case == 5, cb(2,1)=0; LM(1) = Ld; LL(1) = 0.5*(Lb(1)-Ld); LR(1) = LL(1); end % Beam 1, middle
    if B_case == 6, cb(3,1)=0; LR(1) = Ld; LL(1) = 0.5*(Lb(1)-Ld); LM(1) = LL(1); end % Beam 1, right
    if B_case == 7, cc(1,2)=0; Ll(2) = Ld; Lm(2) = 0.5*(Lc(2)-Ld); Lu(2) = Lm(2); end % Col 2, lower
    if B_case == 8, cc(2,2)=0; Lm(2) = Ld; Ll(2) = 0.5*(Lc(2)-Ld); Lu(2) = Ll(2); end % Col 2, middle
    if B_case == 9, cc(3,2)=0; Lu(2) = Ld; Ll(2) = 0.5*(Lc(2)-Ld); Lm(2) = Ll(2); end % Col 2, upper
    if B_case == 10, cb(1,2)=0; LL(2) = Ld; LM(2) = 0.5*(Lb(2)-Ld); LR(2) = LM(2); end % Beam 2, left
    if B_case == 11, cb(2,2)=0; LM(2) = Ld; LL(2) = 0.5*(Lb(2)-Ld); LR(2) = LL(2); end % Beam 2, middle
    if B_case == 12, cb(3,2)=0; LR(2) = Ld; LL(2) = 0.5*(Lb(2)-Ld); LM(2) = LL(2); end % Beam 2, right
    if B_case == 13, cc(1,3)=0; Ll(3) = Ld; Lm(3) = 0.5*(Lc(3)-Ld); Lu(3) = Lm(3); end % Col 3, lower
    if B_case == 14, cc(2,3)=0; Lm(3) = Ld; Ll(3) = 0.5*(Lc(3)-Ld); Lu(3) = Ll(3); end % Col 3, middle
    if B_case == 15, cc(3,3)=0; Lu(3) = Ld; Ll(3) = 0.5*(Lc(3)-Ld); Lm(3) = Ll(3); end % Col 3, upper 
    
    % Calibrate v % VALIDATION - the calibrated values of vRL from FEA are typed here.
    % Comment out this block to perform uncalibrated analysis (default vRL = 1 was specified earlier)
    if B_case == 1, vRL = [-0.24 -6.80]; end
    if B_case == 2, vRL = [3.49 -17.1]; end
    if B_case == 3, vRL = [-0.71 -12.36]; end
    if B_case == 4, vRL = [-0.12 -9.25]; end
    if B_case == 5, vRL = [-0.24 -5.53]; end
    if B_case == 6, vRL = [-0.26 -5.29]; end
    if B_case == 7, vRL = [-1.41 -1.14]; end
    if B_case == 8, vRL = [-1.39 -1.16]; end
    if B_case == 9, vRL = [-1.37 -1.18]; end
    if B_case == 10,vRL = [-0.23 -5.71]; end
    if B_case == 11,vRL = [-0.25 -5.52]; end
    if B_case == 12,vRL = [-0.11 -16.77]; end
    if B_case == 13,vRL = [0.08 14.85]; end
    if B_case == 14,vRL = [-0.06 -25.79]; end
    if B_case == 15,vRL = [0.02 19.22]; end
    if B_case == 16,vRL = [-0.23 -5.98]; end
    
    % Calculate the lateral stiffness for each time step
    for t_case = 1:length(t_vec)
        tf = t_vec(t_case);
        Tc = 0*cc;
        Tb = 0*cb;
        
        % Read column segment temperatures
        for col = 1:n+1
            for seg = 1:3
                % FEA temperature-time relationships
                if cc(seg,col) == 1 % Protected
                    Tc(seg,col) = Get_Ts_FEA(tf,sect_col(col),"UD");
                else % Unprotected
                    Tc(seg,col) = Get_Ts_FEA(tf,sect_col(col),"DI");
                end
            end
        end
        
        % Read beam segment temperatures
        for beam = 1:n
            for seg = 1:3
                % FEA temperature-time relationships
                if cb(seg,beam) == 1 % Protected
                    Tb(seg,beam) = Get_Ts_FEA(tf,sect_beam(beam),"UD");
                else % Unprotected
                    Tb(seg,beam) = Get_Ts_FEA(tf,sect_beam(beam),"DI");
                end
            end
        end
        
        % Lateral stiffness
        [Eb,Ec,mb,mc,taul,tauu,tauL,tauR,RL,RR,rl,ru,...
            Sb,H_,P,Pu,Py,Pb,p,pl,pm,pu,psi,Si,S_data(t_case,B_case),Defl] = ...
            Analyze_3S(LL,LM,LR,Ll,Lm,Lu,Ib,Ic,Ac,Tb,Tc,...
                        vRL,E0,fy0,T0,G,bool_H,rl0,zL0,zR0,m_L,m_R,m_l,D0,Kb);

        % Error processing and reporting of failure times
        Pb_no_errors = Pb+1/0*(Pb<0);
        indices = find(P>Pb_no_errors,1);
        if ~isempty(indices) && buckle == 0
            fprintf('\nCase %.0f, Column %.0f buckled at t = %.2f',B_case,indices(1),tf);
            buckle = 1;
        end
        indices = find(P>Py,1);
        if ~isempty(indices) && yielded == 0
            fprintf('\nCase %.0f, Column %.0f yielded at t = %.2f',B_case,indices(1),tf);
            yielded = 1;
        end
        if S_data(t_case,B_case)< 0 && fail == 0
            fprintf('\nCase %.0f, unstable at t = %.2f, P/Pu = %.3f',B_case,tf,max(P./Pu));
            fail = 1;
        end
        if t_case == length(t_vec) && S_data(t_case,B_case) > 0 && S_data(t_case,B_case) < 1/0
            fprintf('\nCase %.0f, no failures',B_case)
        end
    end
end

% Plot sumS vs t per case
figure(6) 
set(gca, 'FontName', 'Times New Roman');
hold on
grid on
for B_case = 1:16
    line_color = (B_case/16/1.5)*[1 1 1];
    line_style = '-';
    if rem(B_case,3)==1
        line_style = '-.';
    elseif rem(B_case,3)==2
        line_style = '--';
    end
    plot(t_vec,S_data(:,B_case),'LineStyle',line_style,'Color',line_color);
end
lgnd = legend('Case 1','Case 2', 'Case 3', 'Case 4', 'Case 5', 'Case 6', 'Case 7',...
    'Case 8', 'Case 9', 'Case 10', 'Case 11', 'Case 12', 'Case 13', 'Case 14', 'Case 15', 'Undamaged',...
    'Location','eastoutside');
xlabel('Duration of ASTM E119 Fire, \it{t_f} \rm[min]');
ylabel('Storey-based lateral stiffness, \it{\Sigma S} \rm[N/m]');
axis([0 150 0 1.5E6]);
temp = [lgnd; lgnd.ItemText];
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');


%% Localized fire example (Section 6.3.1)
else % example == 2

% Localized Fire Scenarios
T_vec = [0:3:475 475.2:0.1:490 491:2:515 516:0.1:530 533:3:600 601:0.1:622];    % Base temperature load, C -> parameter

% VALIDATION - if necessary, override the values of v_RL here (shown for Case 5):
% Note that the validation for this example was originally done in Excel
% vRL = [0.0420197858608345 23.0661115855394];

% Case 1 - Left end fire
omb(1,:,:) = [1.0 0;
              0.9 0;
              0.7 0];         % Temperature load factors (r,c) = (segment, beam)
                                % segment = [Left; Middle; Right]
omc(1,:,:) = [1.0 0.4 0;
              1.0 0.6 0;
              1.0 0.7 0];    % Temperature load factors (r,c) = (segment, col)
                               % segment = [lower; middle; upper]     
% Case 2 - Left bay fire
omb(2,:,:) = [1.0 0;
              1.0 0;
              1.0 0];         % Temperature load factors (r,c) = (segment, beam)
                                % segment = [Left; Middle; Right]
omc(2,:,:) = [0.6 0.6 0;
              0.7 0.7 0;
              0.9 0.9 0];    % Temperature load factors (r,c) = (segment, col)
                               % segment = [lower; middle; upper]     
% Case 3 - Central fire
omb(3,:,:) = [0.7 1.0;
              0.9 0.9;
              1.0 0.7];         % Temperature load factors (r,c) = (segment, beam)
                                % segment = [Left; Middle; Right]
omc(3,:,:) = [0.4 1.0 0.4;
              0.6 1.0 0.6;
              0.7 1.0 0.7];    % Temperature load factors (r,c) = (segment, col)
                               % segment = [lower; middle; upper]      
% Case 4 - Right bay fire
omb(4,:,:) = [0 1.0;
              0 1.0;
              0 1.0];         % Temperature load factors (r,c) = (segment, beam)
                                % segment = [Left; Middle; Right]
omc(4,:,:) = [0 0.6 0.6;
              0 0.7 0.7;
              0 0.9 0.9];    % Temperature load factors (r,c) = (segment, col)
                               % segment = [lower; middle; upper]  
% Case 5 - Right end fire
omb(5,:,:) = [0 0.7;
              0 0.9;
              0 1.0];         % Temperature load factors (r,c) = (segment, beam)
                                % segment = [Left; Middle; Right]
omc(5,:,:) = [0 0.4 1.0;
              0 0.6 1.0;
              0 0.7 1.0];    % Temperature load factors (r,c) = (segment, col)
                               % segment = [lower; middle; upper]     
                               
D0 = 0;                      % Interstorey drift, m
S_vec = zeros(length(T_vec),5);
Si_vec = zeros(3,length(T_vec),5);

% Determine the failure temperature for each case
for om_case = 1:5
    buckle = 0; yielded = 0; fail = 0; % for output reporting of failure time and mode
    for T_case = 1:length(T_vec)
        T = T_vec(T_case);
        Tc(:,:) = T*omc(om_case,:,:); % Segment temperatures in columns, C
        Tb(:,:) = T*omb(om_case,:,:); % Segment temperatures in beams, C
        
        % Calculate lateral stiffness
        [Eb,Ec,mb,mc,taul,tauu,tauL,tauR,RL,RR,rl,ru,...
            Sb,H_,P,Pu,Py,Pb,p,pl,pm,pu,psi,Si(:,T_case,om_case),S_vec(T_case,om_case),Defl] = ...
            Analyze_3S(LL,LM,LR,Ll,Lm,Lu,Ib,Ic,Ac,Tb,Tc,...
                        vRL,E0,fy0,T0,G,bool_H,rl0,zL0,zR0,m_L,m_R,m_l,D0,Kb);
        
        % Error handling and failure time reporting            
        Pb_no_errors = Pb+1/0*(Pb<0);
        indices = find(P>Pb_no_errors,1);
        if ~isempty(indices) && buckle == 0
            fprintf('\nCase %.0f, Column %.0f buckled at T = %.1f C',om_case,indices(1),T);
            buckle = 1;
        end
        indices = find(P>Py,1);
        if ~isempty(indices) && yielded == 0
            fprintf('\nCase %.0f, Column %.0f yielded at T = %.1f C',om_case,indices(1),T);
            yielded = 1;
        end
        if S_vec(T_case,om_case)< 0 && fail == 0
            fprintf('\nCase %.0f, unstable at T = %.1f, P/Pu = %.3f',om_case,T,max(P./Pu));
            fail = 1;
        end
        if T_case == length(T_vec) && S_vec(T_case,om_case) > 0 && S_vec(T_case,om_case) < 1/0
            fprintf('\nCase %.0f, no failures',om_case)
        end
        if buckle == 1 || yielded == 1 || fail == 1
            S_vec(T_case,om_case) = -1e15;
        end
    end
end

% Plot sumS vs T per case
figure('Renderer', 'painters', 'Position', [10 10 550 330])
set(gca, 'FontName', 'Times New Roman');
hold on
grid on
for om_case=1:5
    line_color = [0 0 0];
    line_style = '-';
    if om_case==2
        line_style = '-.';
    elseif om_case==3
        line_style = ':';
    elseif om_case==4
        line_style = '--';
    elseif om_case==5
        line_color = [0.6 0.6 0.6];
    end
    plot(T_vec,S_vec(:,om_case),'LineStyle',line_style,'Color',line_color);
end
lgnd = legend('Case 1, Left End Fire', 'Case 2, Left Bay Fire','Case 3, Central Fire',...
              'Case 4, Right Bay Fire','Case 5, Right End Fire','Location','southwest');
xlabel('Base temperature load, \it{T} \rm[°C]');
ylabel('Storey-based lateral stiffness, \it{\Sigma S} \rm[N/m]');
axis([0 700 0 1.5E6]);
temp = [lgnd; lgnd.ItemText];
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');
end

%% Function - takes all member properties and temperatures -> stiffness
function [Eb,Ec,mb,mc,taul,tauu,tauL,tauR,RL,RR,rl,ru,Sb,H,P,Pu,Py,Pb,p,pl,pm,pu,psi,Si,sumS,Defl] = ...
    Analyze_3S(LL,LM,LR,Ll,Lm,Lu,Ib,Ic,Ac,Tb,Tc,v,E0,fy0,T0,G,bool_H,rl0,rL0,rR0,m_L,m_R,m_l,D0,Kb)
    
    %Initialize
    n = length(Ib);
    Ec = zeros(3,n+1); Eb = zeros(3,n);
    fyc = zeros(3,n+1);
    mc = zeros(3,n+1); mb = zeros(3,n);
    taul = zeros(1,n+1); tauu = zeros(1,n+1);
    tauL = zeros(1,n); tauR = zeros(1,n);
    rl = zeros(1,n+1); ru = zeros(1,n+1);
    rL = zeros(1,n); rR = zeros(1,n);
    RL = zeros(1,n); RR = zeros(1,n);
    R0_L = zeros(1,n); R0_R = zeros(1,n);
    R0_l = zeros(1,n+1);
    Ru = zeros(1,n+1);
    Sb = zeros(1,n); H = zeros(1,n+1);
    Lb = LL+LM+LR; Lc = Ll+Lm+Lu;
    Pu = zeros(1,n+1); Py = zeros(1,n+1); Pb = zeros(1,n+1);
    Tl = Tc(1,:); Tm = Tc(2,:); Tu = Tc(3,:);
    TL = Tb(1,:); TM = Tb(2,:); TR = Tb(3,:);
    
    % Beam elastic modulus
    for s = 1:3
        for b = 1:n
            Eb(s,b) = Get_EC3_value(Tb(s,b),"Ea",E0); % no axial loads on beams --> E = Ea(T)
            mb(s,b) = Eb(s,b)/E0;
            %[Eb(s,b), mb(s,b)] = Get_E_mu(Tb(s,b),E0); % old code, Lie (1992)
        end
    end
    mL = mb(1,:); mM = mb(2,:); mR = mb(3,:);
    
    % Calculate thermal restraint forces and beam end fixities
    for b = 1:n
        [tauL(b), tauR(b)] = Get_tau(LL(b),LM(b),LR(b),mL(b),mM(b),mR(b));
        R0_L(b) = 3*E0*Ib(b)/Lb(b)*rL0(b)/(1-rL0(b));
        R0_R(b) = 3*E0*Ib(b)/Lb(b)*rR0(b)/(1-rR0(b));
        % Temperatures of connections
        T_conn_L = TL(b);
        T_conn_R = TR(b);
        if LL(b) == 0 
            T_conn_L = TM(b);
        end
        if LR(b) == 0 
            T_conn_R = TM(b);
        end
        rL(b) = Get_r_adj_181219(rL0(b),tauL(b),R0_L(b),m_L(b),T_conn_L);
        rR(b) = Get_r_adj_181219(rR0(b),tauR(b),R0_R(b),m_R(b),T_conn_R);
        [RL(b), RR(b)] = Get_R_3S_190111(E0,Ib(b),rL(b),rR(b),tauL(b),tauR(b),LL(b),LM(b),LR(b),mL(b),mM(b),mR(b),v(b));
    end
    if (bool_H == "Y")
        for b = 1:n
           Sb(b) = GetBeamStiffness(E0,Ib(b),LL(b),LM(b),LR(b),mL(b),mM(b),mR(b),...
                rL(b),rR(b),tauL(b),tauR(b));
        end
        [ki, H] = GetAxialRestraint_181218(Ll,Lm,Lu,Tl,Tm,Tu,T0,E0,fy0,G,Sb,Ac);        
    end
    P = G+H;
    
    % Column elastic modulus and end fixities
    for s = 1:3 % column segments
        for c = 1:(n+1)
            Ec(s,c) = Get_E(Tc(s,c),P(c),Ac(c),E0,fy0);
            mc(s,c) = Ec(s,c)/E0;
            fyc(s,c) = Get_EC3_value(Tc(s,c),"fy",fy0);
        end
    end
    ml = mc(1,:); mm = mc(2,:); mu = mc(3,:);
    for c = 1:(n+1)
        [taul(c), tauu(c)] = Get_tau(Ll(c),Lm(c),Lu(c),ml(c),mm(c),mu(c));
        R0_l(c) = 3*E0*Ic(c)/Lc(c)*rl0(c)/(1-rl0(c));
        % Temperatures of connections
        T_conn_l = Tl(c);
        if Ll(c) == 0
            T_conn_l = Tm(c);
        end
        rl(c) = Get_r_adj_181219(rl0(c),taul(c),R0_l(c),m_l(c),T_conn_l);
    end
    ru = Get_ru_3S_190128(E0,Ic,Ll+Lm+Lu,RL,RR,tauu);
    
    % Calculate Pu
    bool_buckle = 0;
    for c = 1:(n+1)
        Py(c) = Ac(c)*min(fyc(:,c));
        % !!!!NOTE - con_l, con_u and ru0 are hard-coded here for the specific example
        if (c<n+1)
            Ru(c) = Ru(c) + RL(c);
        end
        if (c>1)
            Ru(c) = Ru(c) + RR(c-1);
        end
        Pb(c) = Solve_Pb_3S(E0, fy0, Ic(c), Ac(c), Tu(c), Tm(c), Tl(c), 0, 1, 0, Ru(c),...
                            rl0(c), 0, Lu(c), Lm(c), Ll(c));
        if (Pb(c) < 0)
            Pu(c) = 1/0;
        else
            Pu(c) = min(Py(c),Pb(c));
        end
        if (P(c) > Pu(c))
            %bool_buckle = 1; % temporary hotfix for cases 3 and 4 Pu with v
        end
    end
    
    % Calcluate S
    [p,pl,pm,pu,psi,Si,sumS,Defl] = Get_S_3S_190124(P,Ll,Lm,Lu,ml,mm,mu,E0,Ic,rl,ru,taul,tauu,D0,Kb);
    if (bool_buckle == 1)
        sumS = 1/0;
    end
    
end


