%% This program corresponds to Section 6.5 of the thesis.
% It is used to analyze the post-earthquake fire resistance of a frame via three-segment members
% The tag (%VALIDATION) is used to denote places in the code which are modified during the corresponding validation study in Section 6.5.2

clc, clear, close all
stochastic = 0; % 0 to run the parametric analyses (Sections 6.5.1, 6.5.3), 1 to run the stochastic analysis (Section 6.5.4)

% Problem setup
E0 = 2E11;      % Ambient modulus (N/m2)
fy0 = 3.5E8;    % Yield stress (N/m2)
n = 4;       % # bays
H = 6.096;   % Height of storey, m
T0 = 20;     % Ambient temperature, °C
Kb = 2e5;    % Bracing stiffness, N/m

% Column properties
Ic = [1.29E-4 34.1E-6 34.1E-6 34.1E-6 1.29E-4];    % Inertia, columns, m4
Ac = [7610 4570 4570 4570 7610]/10^6;              % Area, columns, m2
Lc = H      *ones(1,n+1);                          % Length of columns, m
rl0 = [0.3865 0.7045 0.7045 0.7045 0];             % Ambient end fixity factors, lower end columns
m_l = 0.012E6*ones(1,n+1);                         % Reduction slope for connection stiffness
G = 300E3*[1 2 2 2 1];                             % Prescribed gravity loads, N
bool_H = "N";                                      % Thermally restrained columns?
sect_col = ["W310x60", "W200x36", "W200x36","W200x36","W310x60"];

% Beam properties
Ib = 2.45E-4*ones(1,n);   % Inertia, beams, m4
Lb = 7.315  *ones(1,n);   % Length of beams, m
Ab = 8580/10^6*ones(1,n); % Area of beam, m2
rL0 = 0.284741*[1 1 1 1]; % Ambient end fixity factors, beam left connections
rR0 = 0.284741*[1 1 1 0]; % Ambient end fixity factors, beam right connections
m_L = 0.012e6*ones(1,n);
m_R = 0.012e6*ones(1,n);
vRL = 1*ones(1,n);        % VALIDATION - replace with the values obtained from FEA as necessary
sect_beam = ["W410x67", "W410x67","W410x67","W410x67"];

% Curves information: 1 = protected, 0 = damaged
Curves_c(1,:,:) = [0 0 0 0 1;
               1 1 1 1 1;
               0 0 0 0 1];   % Fire curves (r,c) = (segment, col)
                              % segment = [lower; middle; upper]    
Curves_b(1,:,:) = [0 0 0 0;
               1 1 1 1;
               0 0 0 1];     % Fire curves (r,c) = (segment, beam)
                              % segment = [Left; Middle; Right]
Curves_c(2,:,:) = [1 1 1 1 1;
               1 1 1 1 1;
               1 1 1 1 1];   % Fire curves (r,c) = (segment, col)
                              % segment = [lower; middle; upper]    
Curves_b(2,:,:) = [1 1 1 1;
               1 1 1 1;
               1 1 1 1];     % Fire curves (r,c) = (segment, beam)
                              % segment = [Left; Middle; Right]

%% Effect of interstorey drift and insulation damage
if stochastic == 0
t_vec = 0:1:180;  % VALIDATION - change to the desired time stamp to get material properties for input into FEA
D_vec = H*[0 0.01 0.02 0.03 0.04 0.05]; % VALIDATION - change to constant desired analysis case value as necessary
S_data = zeros(length(t_vec),length(D_vec),2);
D_data = zeros(length(t_vec),length(D_vec),2);

% Delamination
Ld = 0.2; % Length of delamination at connections, m
LL = Ld*ones(1,n);   LR = Ld*ones(1,n);   LM = Lb-LL-LR;
Ll = Ld*ones(1,n+1); Lu = Ld*ones(1,n+1); Lm = Lc-Ll-Lu; 

% Compute cases - plot lateral stiffness vs. time for varying interstorey drift
for curves_case = 1:2 % VALIDATION - change to the desired case as necessary
    cc(:,:) = Curves_c(curves_case,:,:);
    cb(:,:) = Curves_b(curves_case,:,:);
    for D_case = 1:length(D_vec)
        D0 = D_vec(D_case);
        % Time steps
        for t_case = 1:length(t_vec)
            tf = t_vec(t_case);
            Tc = 0*cc;
            Tb = 0*cb;
            
            % Column segment temperatures
            for col = 1:n+1
                for seg = 1:3
                    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
            
            % Beam segment temperatures
            for beam = 1:n
                for seg = 1:3
                    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
            
            % Bracing stiffness, Kbr
            T_br = ASTM(tf);
            Kb_T = Kb*Get_EC3_value(T_br, "Ea", E0)/E0;
            
            % Lateral stiffness of frame
            [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,D_case,curves_case),D_data(t_case,D_case,curves_case)] = ...
                Analyze_3S(LL,LM,LR,Ll,Lm,Lu,Ib,Ic,Ac,Tb,Tc,...
                            vRL,E0,fy0,T0,G,bool_H,rl0,rL0,rR0,m_L,m_R,m_l,D0,Kb_T,0);
                        
            % VALIDATION Internal output of properties to input into FEA or compare against FEA
            format long
            tf
            Eb
            Ec
            Tb
            Kb_T
            Defl = D_data(t_case,D_case,curves_case)
        end
    end
end
% Process the deflection data
D_data(S_data<0) = max(max(max(D_data)))+1;

% Plot the figure
figure(1)
labels_all= [{'No Drift (DI)'}
            {'1% Drift (DI)'}
            {'2% Drift (DI)'}
            {'3% Drift (DI)'}
            {'4% Drift (DI)'}
            {'5% Drift (DI)'}
            {'No Drift (UI)'}
            {'1% Drift (UI)'}
            {'2% Drift (UI)'}
            {'3% Drift (UI)'}
            {'4% Drift (UI)'}
            {'5% Drift (UI)'}];
set(gca, 'FontName', 'Times New Roman');
hold on
grid on
for curves_case = 1:2
    for D_case = 1:length(D_vec)
        line_color = D_case/length(D_vec)/1.5*[1 1 1];
        if (curves_case == 1)
            line_style = '-.';
        else
            line_style = '-';
        end
        plot(D_data(:,D_case,curves_case),t_vec,line_style,'Color',line_color)
    end
end
xlabel('Deflection, \rm{\Delta} \rm[m]');
ylabel('Duration of ASTM E119 fire, \it{t_f} \rm[min]');
axis([0 1.0 0 160])
lgnd = legend(labels_all,'Location','eastoutside');
temp = [lgnd; lgnd.ItemText];
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');


%% Effect of delamination
D0 = 0;
Ld_vec = [0,25,100:100:500]*0.001;
% Compute cases - plot lateral stiffness versus fire duration for varying interstorey-drift
S_data = zeros(length(t_vec),length(Ld_vec));
cc(:,:) = Curves_c(1,:,:); % Delaminated case
cb(:,:) = Curves_b(1,:,:); % Delaminated case
for L_case = 1:length(Ld_vec)
    Ld = Ld_vec(L_case);
    LL = Ld*ones(1,n);   LR = Ld*ones(1,n);   LM = Lb-LL-LR;
    Ll = Ld*ones(1,n+1); Lu = Ld*ones(1,n+1); Lm = Lc-Ll-Lu; 
    % Time steps
    for t_case = 1:length(t_vec)
        tf = t_vec(t_case);
        Tc = 0*cc;
        Tb = 0*cb;
        % Column segment temperatures
        for col = 1:n+1
            for seg = 1:3
                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
        % Beam segment temperatures
        for beam = 1:n
            for seg = 1:3
                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
        
        % Bug fix relating to zero-length segments
        if (Ld == 0)
            Tb(1,:) = 0;
            Tb(3,:) = 0;
            Tc(1,:) = 0;
            Tc(3,:) = 0;
        end

        % Bracing stiffness, Kbr
        T_br = ASTM(tf);
        Kb_T = Kb*Get_EC3_value(T_br, "Ea", E0)/E0;

        % Lateral stiffness of frame
        [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,L_case),Defl] = ...
            Analyze_3S(LL,LM,LR,Ll,Lm,Lu,Ib,Ic,Ac,Tb,Tc,...
                        vRL,E0,fy0,T0,G,bool_H,rl0,rL0,rR0,m_L,m_R,m_l,D0,Kb_T,0);
    end
end
% Plot the figure
figure(2)
labels_ld= [{'Undamaged (UD)'}
            {'L_d = 25 mm (DI)'}
            {'L_d = 100 mm (DI)'}
            {'L_d = 200 mm (DI)'}
            {'L_d = 300 mm (DI)'}
            {'L_d = 400 mm (DI)'}
            {'L_d = 500 mm (DI)'}];
set(gca, 'FontName', 'Times New Roman');
hold on
grid on
for L_case = 1:length(Ld_vec)
    line_color = L_case/length(Ld_vec)/1.5*[1 1 1];
    if (Ld_vec(L_case) == 0)
        line_style = '-';
    else
        line_style = '-.';
    end
    plot(t_vec,S_data(:,L_case),line_style,'Color',line_color)
end
xlabel('Duration of ASTM E119 fire, \it{t_f} \rm[min]');
ylabel('Storey-based lateral stiffness, \it{\Sigma S} \rm[N/m]');
axis([0 160 0 8E5])
lgnd = legend(labels_ld,'Location','northeast');
temp = [lgnd; lgnd.ItemText];
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');


%% Stochastic analyses - for revision of PE paper
else %stochastic == 1
N_inst = 10^3;
D_crit = 0.05*H;
t_vec = 0:1:150;
cc(:,:) = Curves_c(1,:,:); % Delaminated case
cb(:,:) = Curves_b(1,:,:); % Delaminated case
num_analyses = 4;
tcr_results = Inf*ones(num_analyses,N_inst);
rng('default');
G_data = zeros(num_analyses,N_inst,5);
D0_data = zeros(num_analyses,N_inst);
Ld_data = zeros(num_analyses,N_inst);
Q_data = zeros(num_analyses,N_inst);
P_ext_mean = 150e3; P_ext_std = 50e3;
P_int_mean = 300e3; P_int_std = 100e3;
Q_const = 30e3;%Q_mean = 0; Q_std = 30e3; 
wbl_a = 1.8;
wbl_b = 1.5;
S_data_stoch = zeros(num_analyses,length(t_vec));
D_data_stoch = zeros(num_analyses,length(t_vec));

for analysis_num = 2:num_analyses % skip #1 because Ld doesn't affect the results much
    fprintf('\nStarting Analysis #%.0d',analysis_num)
    for inst = 1:N_inst
        
        %% I. Randomize variables
        G = zeros(1,5);
        G(1) = max(0,randn*P_ext_std + P_ext_mean);
        G(2) = max(0,randn*P_int_std + P_int_mean);
        G(3) = max(0,randn*P_int_std + P_int_mean);
        G(4) = max(0,randn*P_int_std + P_int_mean);
        G(5) = max(0,randn*P_ext_std + P_int_mean);
        Q = Q_const;%Q = randn*Q_std + Q_mean;
        Ld = 0.001+rand()*0.499;
        D0 = H*0.01*( wbl_b*(-log(1-rand()))^(1/wbl_a) ); % Weibull distribution based on Yazgan & Dazio (2011)
        
        % Analysis #1. Randomize Ld only
        if analysis_num == 1
            G = 150e3*[1 2 2 2 1];
            Q = Q_const;%Q = Q_mean;
            D0 = H*0.01;
        % Analysis #2. Randomize G & Q only
        elseif analysis_num == 2
            D0 = H*0.01;
            Ld = 100e-3;
        % Analysis #3. Randomize D0 via Weibull distribution only
        elseif analysis_num == 3
            Ld = 100e-3;
            Q = Q_const;%Q = Q_mean;
            G = 150e3*[1 2 2 2 1];
            
        % Analysis #5. Randomize everything
        elseif analysis_num == 4
            Ld = 100e-3; % EDIT: Make Ld constant!
        end
        
        % Save the data
        G_data(analysis_num,inst,1) = G(1);
        G_data(analysis_num,inst,2) = G(2);
        G_data(analysis_num,inst,3) = G(3);
        G_data(analysis_num,inst,4) = G(4);
        G_data(analysis_num,inst,5) = G(5);
        Q_data(analysis_num,inst) = Q;
        D0_data(analysis_num,inst) = D0;
        Ld_data(analysis_num,inst) = Ld;
        
        % Apply delamination
        LL = Ld*ones(1,n);   LR = Ld*ones(1,n);   LM = Lb-LL-LR;
        Ll = Ld*ones(1,n+1); Lu = Ld*ones(1,n+1); Lm = Lc-Ll-Lu; 
        
        
        %% II. Get failure time
        for t_case = 1:length(t_vec)
            tf = t_vec(t_case);
            Tc = 0*cc;
            Tb = 0*cb;
            % Column segment temperatures
            for col = 1:n+1
                for seg = 1:3
                    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
            % Beam segment temperatures
            for beam = 1:n
                for seg = 1:3
                    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

            % Bug fix for zero-length segments
            if (Ld == 0)
                Tb(1,:) = 0;
                Tb(3,:) = 0;
                Tc(1,:) = 0;
                Tc(3,:) = 0;
            end

            % Bracing stiffness
            T_br = ASTM(tf);
            Kb_T = Kb*Get_EC3_value(T_br, "Ea", E0)/E0;

            % 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,sumS,Defl] = ...
                Analyze_3S(LL,LM,LR,Ll,Lm,Lu,Ib,Ic,Ac,Tb,Tc,...
                            vRL,E0,fy0,T0,G,bool_H,rl0,rL0,rR0,m_L,m_R,m_l,D0,Kb_T,1);
            S_data_stoch(analysis_num,t_case) = sumS;
            Defl = D0 + (Q+sum(G)*D0/H)/sumS;
            D_data_stoch(analysis_num,t_case) = Defl;
            
            % Failure criterion
            if abs(Defl) > 0.05*H
                fprintf('\nFailure criterion reached at %.0d minutes',tf)
                tcr_results(analysis_num,inst) = tf;
                break;
            end
        end
        
    end
    
    % Plot histogram for analysis case
    figure('Renderer', 'painters', 'Position', [10 10 350 250])
    time_fail = tcr_results(analysis_num,:);
    histogram(time_fail,t_vec(length(t_vec))/10,'FaceColor','k')
    set(gca, 'FontName', 'Times New Roman');
    xlabel('Fire resistance based on Performance Level CP1 (minutes)')
    ylabel('Number of realizations')
    grid on
    set(gca, 'FontName', 'Times New Roman');
    xlim([0 t_vec(length(t_vec))])
    tcr_avg(analysis_num) = mean(time_fail(~isinf(time_fail)));
    tcr_std(analysis_num) = std(time_fail(~isinf(time_fail)));
    tcr_min(analysis_num) = min(time_fail(~isinf(time_fail)));
    tcr_max(analysis_num) = max(time_fail(~isinf(time_fail)));
    time_sorted = sort(time_fail);
end

% Plot CDF's
figure('Renderer', 'painters', 'Position', [10 10 500 250])
set(gca, 'FontName', 'Times New Roman');
for analysis_num =2:num_analyses
    hold on
    time_sorted = sort(tcr_results(analysis_num,:));
    CDF_values = [1:N_inst]/N_inst;
    if analysis_num==2
        linestyle = '--k';
    elseif analysis_num ==3
        linestyle = '-.k';
    else % analysis_num==4
        linestyle = '-k';
    end
    x_vec = [0 time_sorted];
    y_vec = [0 CDF_values];
    ii=1;
    while ii < length(x_vec) % remove kinks in the CDF
        if x_vec(ii+1) == x_vec(ii)
            for jj=ii:1:length(x_vec)-1
                x_vec(jj) = x_vec(jj+1);
                y_vec(jj) = y_vec(jj+1);
            end
            x_vec = x_vec(1:length(x_vec)-1);
            y_vec = y_vec(1:length(y_vec)-1);
        else
            ii=ii+1;
        end
    end
    plot(x_vec,y_vec,linestyle)
end
set(gca, 'FontName', 'Times New Roman');
xlabel('Fire resistance based on Performance Level CP1 (minutes)')
ylabel('CDF')
legend_titles = {'Loads Random','IDR Random','Loads & IDR Random'};
lgnd = legend(legend_titles,'location','northwest')
grid on
set(gca, 'FontName', 'Times New Roman');
xlim([0 t_vec(length(t_vec))])

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,skip_Pu)
    
    %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;
    if skip_Pu == 0
        for c = 1:(n+1)
            Py(c) = Ac(c)*min(fyc(:,c));
            % !!!!NOTE - con_l, con_u and ru0 are hard-coded here
            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));
            %Pb(c) = Solve_Pu_3S(ru(c), rl(c), tauu(c), taul(c),E0, ...
            %    Ic(c),Lu(c), Lm(c), Ll(c),...
            %                mu(c), mm(c), ml(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
    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

function T = ASTM(tf) % tf in minutes
    T = 345*log10(8*tf+1)+20;
end