% This program corresponds to Section 5.2.5 of the thesis.
% It plots the tangent shear modulus Gt vs. T for varying levels of stress
%  1. Gt calculated using direct application of tangent modulus, Gt = Et/2/(1+v)
%  2. Use inelastic buckling tangent modulus via Becque 2016 method
%  3. Use non-conservative incremental theory (e=0)
clc, clear, close all
E0 = 2e11;
fy0 = 3.5e8;
T_vec = 0:1:1200;
v = Get_nu(T_vec); % Note: caps at 725°C
ffy_vec = [0.25 0.5 0.75]; % f/fy stress level
data_Gt = zeros(length(ffy_vec),length(T_vec)); % Direct application of tangent modulus, Gt = Et/2/(1+v)
data_Gi = zeros(length(ffy_vec),length(T_vec)); % Use inelastic buckling tangent modulus, Becque (2016)
data_G0 = zeros(length(ffy_vec),length(T_vec)); % Use non-conservative incremental theory (e=0)
for ffy_case = 1:length(ffy_vec)
    ffy = ffy_vec(ffy_case);
    stress = ffy*fy0;
    for T_case = 1:length(T_vec)
        T = T_vec(T_case);
        Ea = Get_EC3_value(T, 'Ea', E0);
        Et = Get_E(T, stress, 1, E0, fy0); % Note: P/A = stress/1
        if Et == 0
            e = 99999999; % just make it not infinity
        else
            e = 0.5*(Ea/Et-1);
        end
        data_Gt(ffy_case,T_case) = Et/2/(1+v(T_case))       /10^9;
        data_Gi(ffy_case,T_case) = Ea/(2*(1+v(T_case))+3*e) /10^9;
        data_G0(ffy_case,T_case) = Ea/(2*(1+v(T_case)))     /10^9;
    end
end

% Plot Poisson's ratio assumed in the study
figure('Renderer', 'painters', 'Position', [10 10 425 225])
set(gca, 'FontName', 'Times New Roman');
plot(T_vec(1:726),v(1:726),'k-')
hold on
plot(T_vec(727:length(T_vec)),v(727:length(T_vec)),'k--')
grid on
xlabel('Temperature, \itT \rm[°C]')
ylabel('Poisson''s ratio, \nu')
set(gca, 'FontName', 'Times New Roman');

% Plot the tangent shear modulus under the varying assumptions
figure('Renderer', 'painters', 'Position', [10 10 550 260])
set(gca, 'FontName', 'Times New Roman');
legend_titles = {'With \ite = \rm0';
                    'Becque, \itf/f_{y0} = \rm 0.25'; 'Direct,   \itf/f_{y0} = \rm 0.25';
                    'Becque, \itf/f_{y0} = \rm 0.50'; 'Direct,   \itf/f_{y0} = \rm 0.50';
                    'Becque, \itf/f_{y0} = \rm 0.75'; 'Direct,   \itf/f_{y0} = \rm 0.75';};
plot(T_vec,data_G0(1,:),'k','LineWidth',2.5);
hold on
grid on
for ffy_case=1:length(ffy_vec)
    line_color = (1-(1-ffy_case/6))*[1 1 1];
    if ffy_case == 1
        str_style = '-';
    else
        str_style = '--';
    end
    %line_width = 0.5*ffy_case^2;
    plot(T_vec,data_Gi(ffy_case,:),'Color',line_color,'LineStyle','-' )%,'LineWidth',line_width)
    plot(T_vec,data_Gt(ffy_case,:),'Color',line_color,'LineStyle','--')%,'LineWidth',line_width)
end
axis([0 max(T_vec) 0 80]);
lgnd = legend(legend_titles,'location','northeast'); % Alternative: eastoutside
xlabel('Temperature, \itT \rm[°C]')
ylabel('Tangent Shear Modulus, \itG_t \rm[GPa]')
temp = [lgnd; lgnd.ItemText];
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');
set(gca, 'FontName', 'Times New Roman');
