%% This program corresponds to Section 5.6.3 of the thesis.
% The stochastic variable fire resistance of a four bay frame is assessed.
% First, the parametric study on the random variables is conducted (set parametric = 1) and the figures are generated from this program.
% Next, the Monte Carlo simulation is conducted (set parametric = 0) for 10^5 test cases (this may take awhile to run)
% Note that the resulting data from the Monte Carlo simulation should be saved for analysis (the histogram was plotted using a different program)
clc, clear, close all
parametric = 0; % Set to 0 to conduct the Monte Carlo simulation, 1 to generate parametric study curves
if parametric == 1
    dt = 3;         % Time step for analysis, seconds
else
    dt = 60;        % Time step for analysis, seconds
end

%% 1. Frame constants setup
n_bays = 4;             % Number of bays
n_cmp = 2;              % Number of compartments
n_cols = n_bays+1;      % Number of columns
E0 = 2e11;              % Ambient modulus, N/m2
fy0 = 350e6;            % Ambient yield stress, N/m2
Ic = [129 129 34.1 34.1 129]*1e-6;  % Column inertia, m4
Lc = 9.754*ones(1,n_cols);          % Column length, m
Ac = [7610 7610 4570 4570 7610]*1e-6;% Column area, m2
Wc = [60 60 36 36 60];              % Column unit weight, kg/m
Dc = [1.4 1.4 1.05 1.05 1.4];       % Column perimeter, m
rl = [1 1 0 0 0];                   % Column lower end fixity
HT = zeros(1,n_cols);               % Column thermal restraint load, N
Ib = 245e-6*ones(1,n_bays);         % Beam inertia, m4
Lb = 7.315*ones(1,n_bays);          % Beam length, m
Ab = [8580 8580 8580 8580]*1e-6;    % Beam area, m2
Wb = [67 67 67 67];                 % Beam unit weight, kg/m
Db = [1.52 1.52 1.52 1.52];         % Beam perimeter, m
vRL = ones(1,n_bays);               % Beam shape parameter
rL = ones(1,n_bays);                % Beam-column conn fixity
rR = ones(1,n_bays);                % Beam-column conn fixity
Kbr = 1e5;                          % Total semi-bracing in storey, N/m 
tc = [0.044229 0.044229 0.029568 0.029568 0.026537];    % Insul thickness, columns, m
tb = [0.043592 0.026155 0.026155 0.026155];             % Insul thickness, beams, m

%% 2. Thermal setup
Acmp = [282.8 848.4];       % Area in compartments, m2
Ao = [5 45];                % Opening area in compartments, m2
ho = [2 2];              	% Height of opening in compartments, m
fv = (Ao.*ho.^0.5)./Acmp;   % Ventilation factor, m^0.5
kfuel = 330;                % Fuel burning rate property, kg/hm^2.5
Rm = kfuel*Ao.*ho.^0.5/3600;% Mass burning rate, kg/s
H0 = 18800;                 % Enthalpy, kJ/kg
eff = 0.573;                % Burning efficiency
Hp = H0*eff;                % Effective enthalpy, kJ/kg
T0 = 20;                    % Ambient temperature, C
Tu = T0;                    % Outside temperature, C
sigma = 5.67e-8;            % Stefan-Boltzman constant W/m2K4
em_i = 0.9;                 % Radiative emissivity, insulation
em_o = 0.818;               % Emissivity between flame and outside 
rs = 7850;                  % Density steel, kg/m3
ri = 400;                   % Density insulation, kg/m3
cs = 600;                   % Heat capacity steel, J/kgK
ci = 1500;                  % Heat capacity insulation, J/kgK
ks = 40;                    % Conductivity steel, W/mK
ki = 0.12;                  % Conductivity insulation, W/mK
tmax = 600*60;              % Maximum time step, seconds
n_incr = tmax/dt;           % Number of time step increments
a0 = 25;                    % Convective transfer coefficient, members, W/m2K
au = 8.7+0.033*Tu;          % Convective transfer coefficient, outside ,W/m2K
lw = 0.2;                   % Wall thickness, m
nw = 10;                    % Number of wall segments   
dx = lw/nw*ones(1,nw);      % Wall segment thickness
cpgam = 1.7e6*ones(1,nw);   % Thermal inertia, wall, J/m3K
lambda = 0.8*ones(1,nw);    % Conduction parameter, wall, W/mK   
x = 2000/3600;              % Based on Magnuson et al (1976), kg m^-2.5 / s
Ic_max = Rm*3600.*Hp/3600;  % Maximum rate of combustion, kW


%% 3. Variable setup
rng('default');
P_int_mean = 250e3;
P_ext_mean = 125e3;
P_int_std = 50e3;
P_ext_std = 100e3;
q1_mean = 807;
q1_std = 123;
q2_mean = 393;
q2_std = 132;
SR_mean = 45;
SR_std = 15;

% 3A. Monte Carlo Simulation
if parametric == 0
    num_sims = 100000;
    P = zeros(n_cols,num_sims);     % Gravity load, N
    orig = zeros(1,num_sims);       % Location of origin, cmp #
    fuel = zeros(n_cmp,num_sims);   % Fuel load density, MJ/m2
    SR = zeros(n_cmp-1,num_sims);     % Separation resistance, min
    for i_sim = 1:num_sims

        % Applied loads - normally distributed 
        P(1,i_sim) = max(0,randn*P_ext_std + P_ext_mean);
        P(2,i_sim) = max(0,randn*P_int_std + P_int_mean);
        P(3,i_sim) = max(0,randn*P_int_std + P_int_mean);
        P(4,i_sim) = max(0,randn*P_int_std + P_int_mean);
        P(5,i_sim) = max(0,randn*P_ext_std + P_int_mean); % Was P_int_mean -- leaving as is

        % Location of origin
        roll = rand();
        if roll < 0.85 % Chance of fire occuring in compartment 1
            orig(i_sim) = 1; 
        else
            orig(i_sim) = 2;
        end

        % Fuel loads - normally distributed
        fuel(1,i_sim) = max(0,randn*q1_std+q1_mean);
        fuel(2,i_sim) = max(0,randn*q2_std+q2_mean);

        % Separation resistance
        roll = rand();
        if roll < 0.20 % Chance of separation being left open
            SR(1,i_sim) = 0;
        else
            SR(1,i_sim) = round(max(0,randn*SR_std+SR_mean));
        end
    end
end

% 3B. Parametric analyses
if parametric == 1
    num_sims = 19;
    P = zeros(n_cols,num_sims);     % Gravity load, N
    orig = zeros(1,num_sims);       % Location of origin, cmp #
    fuel = zeros(n_cmp,num_sims);   % Fuel load density, MJ/m2
    SR = zeros(n_cmp-1,num_sims);     % Separation resistance, min
    
    % Plots 1, 2 and 4 - Tf and Ts curves
    P(1,1:6) = P_ext_std;
    P(2,1:6) = P_int_std;
    P(3,1:6) = P_int_std;
    P(4,1:6) = P_int_std;
    P(5,1:6) = P_ext_std;
    orig(1:3) = 1;
    orig(4:6) = 2;
    SR(1,1:6) = SR_mean;
    fuel(1,1:3) = q1_mean+q1_std*[-1 0 1];
    fuel(1,4:6) = fuel(1,1:3);
    fuel(2,1:3) = q2_mean+q2_std*[-1 0 1];
    fuel(2,4:6) = fuel(2,1:3);
    
    % Plot 3 - Vary SR
    P(1,7:12) = P_ext_std;
    P(2,7:12) = P_int_std;
    P(3,7:12) = P_int_std;
    P(4,7:12) = P_int_std;
    P(5,7:12) = P_ext_std;
    orig(7:9) = 1;
    orig(10:12) = 2;
    SR(1,7:9) = SR_mean+SR_std*[-1 0 1];
    SR(1,10:12) = SR(1,7:9);
    fuel(1,7:12) = q1_mean;
    fuel(2,7:12) = q2_mean;
    
    % Plot 5 - Vary P
    P(1,13:18) = P_ext_std + P_ext_std*[-1 0 1 -1 0 1];
    P(2,13:18) = P_int_std + P_int_std*[-1 0 1 -1 0 1];
    P(3,13:18) = P_int_std + P_int_std*[-1 0 1 -1 0 1];
    P(4,13:18) = P_int_std + P_int_std*[-1 0 1 -1 0 1];
    P(5,13:18) = P_ext_std + P_ext_std*[-1 0 1 -1 0 1];
    orig(13:15) = 1;
    orig(16:18) = 2;
    SR(1,13:18) = SR_mean;
    fuel(1,13:18) = q1_mean;
    fuel(2,13:18) = q2_mean;
    
    % Plot 6 - Worst case scenario
    P(1:5,19) = [0;281.5;291.6;427.5;123.3]*1e3;
    orig(19) = 2;
    SR(1,19) = 46;
    fuel(1:2,19) = [976.5;767.6];
end


%% 4. Conduct analyses
time_fail = -1*ones(1,num_sims);
sumS_hist = zeros(n_incr,num_sims);
Tf_hist_sim = zeros(n_cmp,n_incr,num_sims);             % Track fire curves
Ts_hist_sim = zeros(n_cols+n_bays,n_incr,num_sims);     % Track steel temperatures (1:5 cols, 6:9 beams)
for i_sim = 1:num_sims
    P_case = P(:,i_sim);
    orig_case = orig(i_sim);
    fuel_case = fuel(:,i_sim);
    SR_case = SR(:,i_sim);
    
    % A. Compute fire curves, Tf(t)
    Tf_hist = zeros(n_cmp, tmax/dt);
    for cmp = 1:n_cmp
        
        max_iter_Tf = 100;          % Maximum iterations to solve for Tf
        tol_Tf = 0.1;               % Tolerance for solving Tf
        Tf_guess = 500;             % Initial guess for Tf
        Tw = T0*ones(nw,nw);        % Wall temperature gradient
        Tw_prev = Tw;           
        for t_step = 1:tmax/dt
            
            % Get fire temperature
            tf = dt*(t_step);
            Tf = Tf_guess;
            converge_Tf = 0; 
            n_iters_Tf =  0;
            while converge_Tf == 0 && n_iters_Tf < max_iter_Tf
                
                % Solve for wall gradient
                ai = 23 + 5.77 * em_o / (Tf - Tw(1, 1)) * ((Tf + 273) ^ 4 - (Tw(1, 1) + 273) ^ 4) / 100 ^ 4;
                hi = 1 / (dx(1) / (2 * lambda(1)) + 1 / ai);
                hu = 1 / (dx(nw) / (2 * lambda(nw)) + 1 / au);
                A_matrix = zeros(nw,nw);
                b_vector = zeros(nw,1);
                for i = 1:nw
                    b_vector(i, 1) = dx(i) * cpgam(i) * Tw_prev(i, 1) / dt;
                end
                b_vector(1, 1) = b_vector(1, 1) + Tf * hi;
                b_vector(nw, 1) = b_vector(nw, 1) + Tu * hu;
                A_matrix(1, 1) = dx(1) * cpgam(1) / dt + hij(1, 2, dx, lambda) + hi;
                A_matrix(1, 2) = -hij(1, 2, dx, lambda);
                for k = 2:nw - 1
                    A_matrix(k, k) = dx(k) * cpgam(k) / dt + hij(k, k - 1, dx, lambda) + hij(k, k + 1, dx, lambda);
                    A_matrix(k, k - 1) = -hij(k, k - 1, dx, lambda);
                    A_matrix(k, k + 1) = -hij(k, k + 1, dx, lambda);
                end
                A_matrix(nw, nw) = dx(nw) * cpgam(nw) / dt + hij(nw, nw - 1, dx, lambda) + hu;
                A_matrix(nw, nw - 1) = -hij(nw, nw - 1, dx, lambda);
                Tw = inv(A_matrix)*b_vector;

                % Compute Ic
                T_dur = fuel_case(cmp)*1000*Acmp(cmp) / (Rm(cmp)*3600*H0);
                Icom = getIc(Ic_max(cmp), T_dur, tf / 3600);

                % Calculate resulting Tf
                cp = 1 + 0.14 / 1000 * Tf;
                alpha_i_new = 23 + 5.77 * em_o / (Tf - Tw(1, 1)) * ((Tf + 273) ^ 4 - (Tw(1, 1) + 273) ^ 4) / 100 ^ 4; % W/m2K
                hi_new = 1 / (1 / alpha_i_new + dx(1) / (2 * lambda(1))); % W/m2K
                Ir = Ao(cmp) * 5.77 / 100 ^ 4 * ((Tf + 273) ^ 4 - (T0 + 273) ^ 4) / 1000; % kW
                Tf_res = (Icom + x * cp * Ao(cmp) * ho(cmp) ^ 0.5 * T0 + ...
                            (Acmp(cmp) - Ao(cmp)) * (hi_new / 1000) * Tw(1, 1) - Ir) / ...
                            (x * cp * Ao(cmp) * ho(cmp) ^ 0.5 + (Acmp(cmp) - Ao(cmp)) * (hi_new / 1000));

                % Check convergence
                err = Tf_res - Tf;
                if (abs(err) > tol_Tf)
                    Tf = Tf_res;
                    n_iters_Tf = n_iters_Tf + 1;
                else
                    converge_Tf = 1;
                end
            end
            
            % Process results of time step
            if converge_Tf == 1
                Tf_guess = Tf;
                Tw_prev = Tw;
            else
                Tf = -100;
                for i = 1:n
                    Tw(i, 1) = 0;
                end
                break;
            end
            Tf_hist(cmp, t_step) = Tf;
        end
    end
    Tf_hist_sim(:,:,i_sim) = Tf_hist;
    
    % B. Frame Analysis
    global_time_step = 0;
    global_time = 0;
    stability = 1;
    Tb = T0*ones(1,n_bays);  % Beam temperature
    Tc = T0*ones(1,n_cols);  % Column temperature
    Eb = E0*ones(1,n_bays);  % Beam elasticity
    Ec = E0*ones(1,n_cols);  % Column elasticity
    Tf_curr = zeros(1,n_cmp);
    while stability == 1 && global_time < tmax
        
        % I. Retrieve compartment fire temperature
        global_time_step = global_time_step + 1;
        global_time = global_time + dt;
        if orig_case == 1
            cmp_step_1 = global_time_step;
            cmp_step_2 = max(0,global_time_step - SR_case(1)*60/dt);
        else
            cmp_step_2 = global_time_step;
            cmp_step_1 = max(0,global_time_step - SR_case(1)*60/dt);
        end
        if cmp_step_1 > 0
            Tf_curr(1) = Tf_hist(1,cmp_step_1);
        else
            Tf_curr(1) = T0;
        end
        if cmp_step_2 > 0
            Tf_curr(2) = Tf_hist(2,cmp_step_2);
        else
            Tf_curr(2) = T0;
        end
        
        % II. Update steel temperature
        for beam = 1:n_bays
            if beam == 1
                Tf_local = Tf_curr(1);
            else
                Tf_local = Tf_curr(2);
            end
            Tb(beam) = Tb(beam) +... % Buchanan 2001; multiple methods exist - see VBA code
                ki/tb(beam)*(Tf_local-Tb(beam))*dt/(cs*Wb(beam)/Db(beam)+ci*ri*tb(beam)/2) - 0;
            Eb(beam) = Get_E(Tb(beam), 0, Ab(beam), E0, fy0);
            Ts_hist_sim(n_cols+beam,global_time_step,i_sim) = Tb(beam);
        end
        for col = 1:n_cols
            if col == 1
                Tf_local = Tf_curr(1);
            elseif col == 2
                Tf_local = max(Tf_curr(1),Tf_curr(2)); % Take maximum of fires on both sides
            else
                Tf_local = Tf_curr(2);
            end
            Tc(col) = Tc(col) + ... % Buchanan 2001; multiple methods exist - see VBA code
                ki/tc(col)*(Tf_local-Tc(col))*dt/(cs*Wc(col)/Dc(col)+ci*ri*tc(col)/2) - 0;
            Ec(col) = Get_E(Tc(col), P_case(col), Ac(col), E0, fy0);
            Ts_hist_sim(col,global_time_step,i_sim) = Tc(col);
        end
        
        % III. Stability
        RL = zeros(1,n_bays);
        RR = zeros(1,n_bays);
        for beam = 1:n_bays
            RL(beam) = 6*Eb(beam)*Ib(beam)/Lb(beam) * rL(beam)/(4-rL(beam)*rR(beam)) *(2+ vRL(beam)*rR(beam));
            RR(beam) = 6*Eb(beam)*Ib(beam)/Lb(beam) * rR(beam)/(4-rL(beam)*rR(beam)) *(2+ 1/vRL(beam)*rR(beam));
        end
        phi = zeros(1,n_cols);
        ru = zeros(1,n_cols);
        Si = zeros(1,n_cols);
        sumS = 0;
        for col = 1:n_cols
            phi(col) = Lc(col)*sqrt(P_case(col)/(Ec(col)*Ic(col)));
            Ru = 0;
            if col < n_cols
                Ru = Ru + RL(col);
            end
            if col > 1
                Ru = Ru + RR(col-1);
            end
            ru(col) = 1/(1+3*Ec(col)*Ic(col)/(Ru*Lc(col)));
            a1 = 3*(rl(col)*(1-ru(col)) + ru(col)*(1-rl(col)));
            a2 = 9*rl(col)*ru(col) - (1-rl(col))*(1-ru(col))*phi(col)^2;
            a3 = 18*rl(col)*ru(col) + a1*phi(col)^2;
            if phi(col) > 0
                beta = 1/12*phi(col)^3 * ...
                        (a1*phi(col)*cos(phi(col)) + a2*sin(phi(col))) /...
                        (18*rl(col)*ru(col) - a3*cos(phi(col)) + (a1-a2)*phi(col)*sin(phi(col)));
            else % phi(col) == 0; tension not supported and will be treated as zero load here 
                beta = (rl(col)+ru(col)+rl(col)*ru(col))/(4-rl(col)*ru(col));
            end
            Si(col) = 12*Ec(col)*Ic(col)/Lc(col)^3*beta;
        	sumS = sumS + Si(col);
        end
        sumS = sumS + Kbr;
        sumS_hist(global_time_step,i_sim) = sumS;
        if sumS < 0 && parametric == 0
            stability = 0;
        end
        
        % IV. Check for rotational buckling or E = 0
        K_app = Get_K(ru,rl);
        Pu = pi^2 .* Ec .* Ic./(K_app .* Lc).^2;
        if sum(P_case'>Pu) > 0  && parametric == 0
            stability = 0;
        end
        if sum(Ec==0) > 0   && parametric == 0
            stability = 0;
        end
        
    end
    
    % V. Instability has occurred - store time of failure
    if stability == 0
        time_fail(i_sim) = global_time/60;
    else
        time_fail(i_sim) = Inf;
    end
    
end


%% 5. Post-processing
if parametric == 1
    plot_size = [10 10 450 275];
    time_vec = 1/60*[0:dt:tmax];
    figure('Renderer', 'painters', 'Position', [10 10 450 350]) % Plot 1
    for i_sim =1:3 % Cmpt 1
        plot(time_vec,[T0 Tf_hist_sim(1,:,i_sim)],'Color','black','LineWidth',i_sim*0.5,'LineStyle','-')
        hold on
    end
    for i_sim =1:3 % Cmpt 2
        plot(time_vec,[T0 Tf_hist_sim(2,:,i_sim)],'Color',[0.5 0.5 0.5],'LineWidth',i_sim*0.5,'LineStyle','--')
        hold on
    end
    xlabel('Duration of fire in compartment (minutes)')
    ylabel('Temperature of fire, T_f (°C)')
    grid on
    lgnd_titles = [{"Cmpt. 1, Small"};{"Cmpt. 1, Mean"};{"Cmpt. 1, Large"};{"Cmpt. 2, Small"};{"Cmpt. 2, Mean"};{"Cmpt. 2, Large"}];
    lgnd = legend(lgnd_titles,'Location','best');
    temp = [lgnd; lgnd.ItemText];
    set(temp, 'FontSize', 10);
    set(temp, 'FontName', 'Times New Roman');
    set(gca, 'FontName', 'Times New Roman');
    
    figure('Renderer', 'painters', 'Position', plot_size) % Plot 2 - Col 2
    for i_sim =1:3 % O = 1
        plot(time_vec,[T0 Ts_hist_sim(2,:,i_sim)],'Color','black','LineWidth',i_sim*0.5,'LineStyle','-')
        hold on
    end
    for i_sim =4:6 % O = 2
        plot(time_vec,[T0 Ts_hist_sim(2,:,i_sim)],'Color',[0.5 0.5 0.5],'LineWidth',(i_sim-3)*0.5,'LineStyle','--')
        hold on
    end
    xlabel('Duration of fire event, t_O (minutes)')
    ylabel('Temperature of Column 2, T_{c,2} (°C)')
    grid on
    lgnd_titles = [{"O=1, Small"};{"O=1, Mean"};{"O=1, Large"};{"O=2, Small"};{"O=2, Mean"};{"O=2, Large"}];
    lgnd = legend(lgnd_titles,'Location','best');
    temp = [lgnd; lgnd.ItemText];
    set(temp, 'FontSize', 10);
    set(temp, 'FontName', 'Times New Roman');
    set(gca, 'FontName', 'Times New Roman');
    
    figure('Renderer', 'painters', 'Position', plot_size) % Plot 2 - Col 1
    for i_sim =1:3 % O = 1
        plot(time_vec,[T0 Ts_hist_sim(1,:,i_sim)],'Color','black','LineWidth',i_sim*0.5,'LineStyle','-')
        hold on
    end
    for i_sim =4:6 % O = 2
        plot(time_vec,[T0 Ts_hist_sim(1,:,i_sim)],'Color',[0.5 0.5 0.5],'LineWidth',(i_sim-3)*0.5,'LineStyle','--')
        hold on
    end
    xlabel('Duration of fire event, t_O (minutes)')
    ylabel('Temperature of Column 1, T_{c,1} (°C)')
    grid on
    lgnd_titles = [{"O=1, Small"};{"O=1, Mean"};{"O=1, Large"};{"O=2, Small"};{"O=2, Mean"};{"O=2, Large"}];
    lgnd = legend(lgnd_titles,'Location','best');
    temp = [lgnd; lgnd.ItemText];
    set(temp, 'FontSize', 10);
    set(temp, 'FontName', 'Times New Roman');
    set(gca, 'FontName', 'Times New Roman');
    
    figure('Renderer', 'painters', 'Position', plot_size) % Plot 2 - Col 4
    for i_sim =1:3 % O = 1
        plot(time_vec,[T0 Ts_hist_sim(4,:,i_sim)],'Color','black','LineWidth',i_sim*0.5,'LineStyle','-')
        hold on
    end
    for i_sim =4:6 % O = 2
        plot(time_vec,[T0 Ts_hist_sim(4,:,i_sim)],'Color',[0.5 0.5 0.5],'LineWidth',(i_sim-3)*0.5,'LineStyle','--')
        hold on
    end
    xlabel('Duration of fire event (minutes)')
    ylabel('Temperature of Column 4, T_{c,4} (°C)')
    grid on
    lgnd_titles = [{"O=1, Small"};{"O=1, Mean"};{"O=1, Large"};{"O=2, Small"};{"O=2, Mean"};{"O=2, Large"}];
    lgnd = legend(lgnd_titles,'Location','best');
    temp = [lgnd; lgnd.ItemText];
    set(temp, 'FontSize', 10);
    set(temp, 'FontName', 'Times New Roman');
    set(gca, 'FontName', 'Times New Roman');
    
    figure('Renderer', 'painters', 'Position', plot_size) % Plot 3
    for i_sim =7:9 % O = 1
        plot(time_vec,[T0 Ts_hist_sim(2,:,i_sim)],'Color','black','LineWidth',(i_sim-6)*0.5,'LineStyle','-')
        hold on
    end
    for i_sim =10:12 % O = 2
        plot(time_vec,[T0 Ts_hist_sim(2,:,i_sim)],'Color',[0.5 0.5 0.5],'LineWidth',(i_sim-9)*0.5,'LineStyle','--')
        hold on
    end
    xlabel('Duration of fire event, t_O (minutes)')
    ylabel('Temperature of Column 2, T_{c,2} (°C)')
    grid on
    lgnd_titles = [{"O=1, Fast spread"};{"O=1, Mean"};{"O=1, Slow spread"};{"O=2, Fast spread"};{"O=2, Mean"};{"O=2, Slow spread"}];
    lgnd = legend(lgnd_titles,'Location','best');
    temp = [lgnd; lgnd.ItemText];
    set(temp, 'FontSize', 10);
    set(temp, 'FontName', 'Times New Roman');
    set(gca, 'FontName', 'Times New Roman');
    
    figure('Renderer', 'painters', 'Position', plot_size) % Plot 4
    for i_sim =1:3 % O = 1
        plot(time_vec,[sumS_hist(1,i_sim); sumS_hist(:,i_sim)],'Color','black','LineWidth',i_sim*0.5,'LineStyle','-')
        hold on
    end
    for i_sim =4:6 % O = 2
        plot(time_vec,[sumS_hist(1,i_sim); sumS_hist(:,i_sim)],'Color',[0.5 0.5 0.5],'LineWidth',(i_sim-3)*0.5,'LineStyle','--')
        hold on
    end
    xlabel('Duration of fire event, t_O (minutes)')
    ylabel('Frame lateral stiffness, \SigmaS (N/m)')
    grid on
    lgnd_titles = [{"O=1, Small"};{"O=1, Mean"};{"O=1, Large"};{"O=2, Small"};{"O=2, Mean"};{"O=2, Large"}];
    lgnd = legend(lgnd_titles,'Location','best');
    temp = [lgnd; lgnd.ItemText];
    ylim([0 1e6]);
    set(temp, 'FontSize', 10);
    set(temp, 'FontName', 'Times New Roman');
    set(gca, 'FontName', 'Times New Roman');
    
    figure('Renderer', 'painters', 'Position', plot_size) % Plot 5
    for i_sim =13:15 % O = 1
        plot(time_vec,[sumS_hist(1,i_sim); sumS_hist(:,i_sim)],'Color','black','LineWidth',(i_sim-12)*0.5,'LineStyle','-')
        hold on
    end
    for i_sim =16:18 % O = 2
        plot(time_vec,[sumS_hist(1,i_sim); sumS_hist(:,i_sim)],'Color',[0.5 0.5 0.5],'LineWidth',(i_sim-15)*0.5,'LineStyle','--')
        hold on
    end
    xlabel('Duration of fire event, t_O (minutes)')
    ylabel('Frame lateral stiffness, \SigmaS (N/m)')
    grid on
    lgnd_titles = [{"O=1, Light"};{"O=1, Mean"};{"O=1, Heavy"};{"O=2, Light"};{"O=2, Mean"};{"O=2, Heavy"}];
    lgnd = legend(lgnd_titles,'Location','best');
    temp = [lgnd; lgnd.ItemText];
    ylim([0 1e6]);
    set(temp, 'FontSize', 10);
    set(temp, 'FontName', 'Times New Roman');
    set(gca, 'FontName', 'Times New Roman');
    
    figure('Renderer', 'painters', 'Position', plot_size) % Plot 6
    i_sim = 19;
    plot(time_vec,[sumS_hist(1,i_sim); sumS_hist(:,i_sim)],'Color','black','LineWidth',3,'LineStyle','-')
    xlabel('Duration of fire event, t_O (minutes)')
    ylabel('Frame lateral stiffness, \SigmaS (N/m)')
    grid on
    xlim([0 100]);
    ylim([0 1e6]);
    set(gca, 'FontName', 'Times New Roman');
else
    figure
    histogram(time_fail,tmax/60/20,'FaceColor','k')
    set(gca, 'FontName', 'Times New Roman');
    xlabel('Fire resistance of frame (minutes)')
    ylabel('Number of realizations')
    grid on
    set(gca, 'FontName', 'Times New Roman');
    xlim([0 tmax/60])
    avg = mean(time_fail(~isinf(time_fail)));
    stdev = std(time_fail(~isinf(time_fail)));
    minim = min(time_fail(~isinf(time_fail)));
    maxim = max(time_fail(~isinf(time_fail)));
    time_sorted = sort(time_fail);
end

%% Helper functions
function h = hij(i, j, dx, lambda)
    h = 1 / (dx(i) / (2 * lambda(i)) + dx(j) / (2 * lambda(j)));
end

function Ic = getIc(Ic_max, T_dur, t_hours)
    t1_start = 0.11;
    t1_end = (1 - T_dur / 8) * T_dur;
    t_half = 1.5 * T_dur;
    t_fifth = 2.3 * T_dur;
    t_zero = 4 * T_dur ^ (0.66667);
    
    if (t_hours < t1_start) 
        ratio = t_hours / t1_start;
    elseif (t_hours < t1_end) 
        ratio = 1;
    elseif (t_hours < t_half) 
        ratio = 1 - (t_hours - t1_end) / (t_half - t1_end) * 0.5;
    elseif (t_hours < t_fifth) 
        ratio = 0.5 - (t_hours - t_half) / (t_fifth - t_half) * 0.3;
    elseif (t_hours < t_zero) 
        ratio = 0.2 - (t_hours - t_fifth) / (t_zero - t_fifth) * 0.2;
    else
        ratio = 0;
    end
    
    if (ratio > 1) 
        ratio = 1;
    elseif (ratio < 0) 
        ratio = 0;
    end
    
    Ic = Ic_max * ratio;
end