% This program corresponds to Sections 2.1.2, 3.2.2, 3.2.5, A3.2.3 and A3.2.5
% The behaviour of the beta function is plotted to produce the figures in the thesis
clc, clear, close all

%% 1. Plot beta(phi,rl,ru) - Section 2.1.2
phi_vec = [0:pi/100:pi/2 pi/2+pi/1000:pi/1000:2*pi];
rl_vec = 1; % Can change this to zero for a second figure
ru_vec = 0:0.2:1;
beta_data = zeros(length(phi_vec),length(rl_vec),length(ru_vec));

for rl_case = 1:length(rl_vec)
    rl = rl_vec(rl_case);
    for ru_case = 1:length(ru_vec)
        ru = ru_vec(ru_case);
        K_app = Get_K(ru,rl);
        %phiu = pi/K_app;
        phiu = Get_Phiu(ru,rl,K_app,0,0);
        for phi_case = 1:length(phi_vec)
            phi = phi_vec(phi_case);
            beta = Get_beta(phi,rl,ru);
            
            % Remove past Nu
            if phi_case > 1 && (beta > beta_data(phi_case-1,rl_case,ru_case) ||phi > phiu)
                beta = -Inf;
                beta_data(phi_case:length(phi_vec),rl_case,ru_case) = -Inf;
                continue;
            end
            beta_data(phi_case,rl_case,ru_case) = beta;
        end
    end
end
figure('Renderer', 'painters', 'Position', [10 10 500 280])
set(gca, 'FontName', 'Times New Roman');
legend_titles = {};
for rl_case=1:length(rl_vec)
    for ru_case = rl_case:length(ru_vec)
        legend_titles(length(legend_titles)+1) = ...
            {strcat('\it r_l =\rm',num2str(rl_vec(rl_case)),...
                    ', \it r_u =\rm',num2str(ru_vec(ru_case)))};
        line_style = length(legend_titles);
        line_color = (1-ru_vec(ru_case))/1.3*[1 1 1];
        str_style = '-';
%         if rem(line_style,3) == 0
%             str_style = '-';
%         elseif rem(line_style,3) == 1
%             str_style = '--';
%         else
%             str_style = '-.';
%         end
        line_width = 0.5;
        %line_width = rem(line_style,2)+1;
        plot(phi_vec,beta_data(:,rl_case,ru_case),'Color',line_color,'LineStyle',str_style,'LineWidth',line_width)
        hold on
        grid on
    end
end
if rl_vec ==0
    lgnd = legend(legend_titles,'location','ne');
else
    lgnd = legend(legend_titles,'location','sw');
end
xlabel('\phi')
ylabel('\beta')
if rl_vec(1) == 0
    ylim([-4 1])
else
    ylim([-4 1])
end
xlim([0 7])
temp = [lgnd; lgnd.ItemText];
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');
set(gca, 'FontName', 'Times New Roman');

% Get Taylor series error 
r = 0:0.01:1;
beta_exact = zeros(length(phi_vec),length(r));
beta_tay = zeros(length(phi_vec),length(r));
for phi_case = 1:length(phi_vec)
    phi = phi_vec(phi_case);
    for r_case = 1:length(r)
        ru = r(r_case); rl = r(r_case);
        beta_exact(phi_case,r_case) = Get_beta(phi,rl,ru);
        beta_tay(phi_case,r_case) = (rl+ru+rl*ru)/(4-rl*ru) - phi^2* ...
            (40+8*(rl^2+ru^2)+rl*ru*(rl+ru+3*rl*ru-34))/(30*(4-rl*ru)^2);
    end
end
err = abs(beta_tay-beta_exact)./max(abs(beta_tay),abs(beta_exact));

%% 2. Plot beta versus shear flexibility - Section 3.2.5
n_vec = [0 1e-3 5e-3 1e-2 5e-2];
rl_vec = 0.75;
ru_vec = 0.5;
beta_shear_data = zeros(length(n_vec),length(phi_vec),length(rl_vec),length(ru_vec));
shear_assumption = 0; % Can change to 1 for a second plot
for n_case = 1:length(n_vec)
    neta = n_vec(n_case);
    for rl_case = 1:length(rl_vec)
        rl = rl_vec(rl_case);
        for ru_case = 1:length(ru_vec)
            ru = ru_vec(ru_case);
            K_app = Get_K(ru,rl);
            phiu = Get_Phiu(ru,rl,K_app,neta,shear_assumption);
            if phiu < 0
                phiu = Inf;
            end
            for phi_case = 1:length(phi_vec)
                phi = phi_vec(phi_case);
                [beta,zeta] = Get_beta_shear(phi,rl,ru,neta,shear_assumption);
                
                % Remove past Nu
                if phi_case > 1 && (beta/(1+zeta) > beta_shear_data(n_case,phi_case-1,rl_case,ru_case) || phi>phiu)
                    beta = -Inf;
                    beta_shear_data(n_case,phi_case:length(phi_vec),rl_case,ru_case) = -Inf;
                    continue;
                end
                beta_shear_data(n_case,phi_case,rl_case,ru_case) = beta/(1+zeta);
            end
        end
    end
end
figure('Renderer', 'painters', 'Position', [20 20 500 260])
set(gca, 'FontName', 'Times New Roman');
legend_titles = {};
for n_case = 1:length(n_vec)
    for rl_case=1:length(rl_vec)
        for ru_case = rl_case:length(ru_vec)
            legend_titles(length(legend_titles)+1) = ...
                {strcat('\it r_l =\rm',num2str(rl_vec(rl_case)),...
                        ', \it r_u =\rm',num2str(ru_vec(ru_case)),...
                        ', \it \eta =\rm',num2str(n_vec(n_case)))};
            line_style = length(legend_titles);
            line_color = n_case/length(n_vec)/1.3*[1 1 1];
            str_style = '-';
%             if rem(line_style,3) == 0
%                 str_style = '-.';
%             elseif rem(line_style,3) == 1
%                 str_style = '-';
%             else
%                 str_style = '--';
%             end
            line_width = 0.5;
            %line_width = rem(line_style,2)+1;
            plot(phi_vec,beta_shear_data(n_case,:,rl_case,ru_case),'Color',line_color,'LineStyle',str_style,'LineWidth',line_width)
            hold on
            grid on
        end
    end
end
lgnd = legend(legend_titles,'location','sw');
xlabel('\phi')
ylabel('\beta/(1+\zeta)')
ylim([-4 1])
xlim([0 5])
temp = [lgnd; lgnd.ItemText];
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');
set(gca, 'FontName', 'Times New Roman');

%% 3. Plot beta tension - no shear (Appendix A3.2.3)
phi_vec = 0:0.01:17;
rl_vec = 1; % Can change this to zero for a second figure
ru_vec = 0:0.5:1;
beta_tens_data = zeros(length(phi_vec),length(rl_vec),length(ru_vec));

for rl_case = 1:length(rl_vec)
    rl = rl_vec(rl_case);
    for ru_case = 1:length(ru_vec)
        ru = ru_vec(ru_case);
        for phi_case = 1:length(phi_vec)
            phi = phi_vec(phi_case);
            beta = Get_beta_tens_shear(phi,rl,ru,0,0);
            beta_tens_data(phi_case,rl_case,ru_case) = beta;
        end
    end
end
figure('Renderer', 'painters', 'Position', [30 30 500 260])
set(gca, 'FontName', 'Times New Roman');
legend_titles = {};
for rl_case=1:length(rl_vec)
    for ru_case = rl_case:length(ru_vec)
        legend_titles(length(legend_titles)+1) = ...
            {strcat('\it r_l =\rm',num2str(rl_vec(rl_case)),...
                    ', \it r_u =\rm',num2str(ru_vec(ru_case)))};
        line_style = length(legend_titles);
        line_color = (1-ru_vec(ru_case))/1.3*[1 1 1];
        str_style = '-';
%         if rem(line_style,3) == 0
%             str_style = '-';
%         elseif rem(line_style,3) == 1
%             str_style = '--';
%         else
%             str_style = '-.';
%         end
        line_width = 0.5;
        %line_width = rem(line_style,2)+1;
        plot(phi_vec.^2,beta_tens_data(:,rl_case,ru_case),'Color',line_color,'LineStyle',str_style,'LineWidth',line_width)
        hold on
        grid on
    end
end
lgnd = legend(legend_titles,'location','nw');
xlabel('$\hat{\phi}^2 = \frac{TL^2}{EI}$','Interpreter','latex')
ylabel('$\hat{\beta}$','Interpreter','latex')
ylim([0 30])
xlim([0 300])
temp = [lgnd; lgnd.ItemText];
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');
set(gca, 'FontName', 'Times New Roman');

%% 4. Plot beta tension - with shear (Appendix 3.2.5)
n_vec = [0 5e-3 1e-2];
rl_vec = 0.75;
ru_vec = 0.5;
beta_tens_shear_data = zeros(length(n_vec),length(phi_vec),length(rl_vec),length(ru_vec));
shear_assumption = 0; % Can change to 1 for a second plot
for n_case = 1:length(n_vec)
    neta = n_vec(n_case);
    for rl_case = 1:length(rl_vec)
        rl = rl_vec(rl_case);
        for ru_case = 1:length(ru_vec)
            ru = ru_vec(ru_case);
            for phi_case = 1:length(phi_vec)
                phi = phi_vec(phi_case);
                [beta,zeta] = Get_beta_tens_shear(phi,rl,ru,neta,shear_assumption);
                beta_tens_shear_data(n_case,phi_case,rl_case,ru_case) = beta/(1+zeta);
                if phi > 0.1/neta^0.5
                    beta_tens_shear_data(n_case,phi_case,rl_case,ru_case) = 0/0;
                end
            end
        end
    end
end
figure('Renderer', 'painters', 'Position', [40 40 500 260])
set(gca, 'FontName', 'Times New Roman');
legend_titles = {};
for n_case = 1:length(n_vec)
    for rl_case=1:length(rl_vec)
        for ru_case = rl_case:length(ru_vec)
            legend_titles(length(legend_titles)+1) = ...
                {strcat('\it r_l =\rm',num2str(rl_vec(rl_case)),...
                        ', \it r_u =\rm',num2str(ru_vec(ru_case)),...
                        ', \it \eta =\rm',num2str(n_vec(n_case)))};
            line_style = length(legend_titles);
            line_color = n_case/length(n_vec)/1.3*[1 1 1];
            str_style = '-';
%             if rem(line_style,3) == 0
%                 str_style = '-.';
%             elseif rem(line_style,3) == 1
%                 str_style = '-';
%             else
%                 str_style = '--';
%             end
            line_width = 0.5;
            %line_width = rem(line_style,2)+1;
            plot(phi_vec.^2,beta_tens_shear_data(n_case,:,rl_case,ru_case),'Color',line_color,'LineStyle',str_style,'LineWidth',line_width)
            hold on
            grid on
        end
    end
end
lgnd = legend(legend_titles,'location','se');
xlabel('$\hat{\phi}^2 = \frac{TL^2}{EI}$','Interpreter','latex')
ylabel('$\hat{\beta}/(1+\zeta)$','Interpreter','latex')
ylim([0 1.0])
xlim([0 5])
temp = [lgnd; lgnd.ItemText];
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');
set(gca, 'FontName', 'Times New Roman');

%% 5. Comparison of omega in compression (Section 3.2.2)
x_vec = 0:0.01:1;
omega_phi = zeros(1,length(x_vec));
omega_dy = zeros(1,length(x_vec));
for x_case = 1:length(x_vec)
    omega_phi(x_case) = sqrt(1+x_vec(x_case));
    omega_dy(x_case) = (1-x_vec(x_case))^(-0.5);
end
figure('Renderer', 'painters', 'Position', [40 40 500 260])
set(gca, 'FontName', 'Times New Roman');
legend_titles = {'Engesser (1891)','Haringx (1948)'};
plot(x_vec,omega_phi,'k-')
hold on
grid on
plot(x_vec,omega_dy,'k--')
lgnd = legend(legend_titles,'location','nw');
xlabel('\eta\phi^2')
ylabel('\omega')
ylim([0 5])
xlim([0 1])
temp = [lgnd; lgnd.ItemText];
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');
set(gca, 'FontName', 'Times New Roman');


%% 5. Comparison of omega in tension (Appendix 3.2.3)
x_vec = 0:0.01:1;
omega_phi = zeros(1,length(x_vec));
omega_dy = zeros(1,length(x_vec));
for x_case = 1:length(x_vec)
    omega_phi(x_case) = sqrt(1-x_vec(x_case));
    omega_dy(x_case) = (1+x_vec(x_case))^(-0.5);
end
figure('Renderer', 'painters', 'Position', [40 40 500 260])
set(gca, 'FontName', 'Times New Roman');
legend_titles = {'Engesser (1891)','Haringx (1948)'};
plot(x_vec,omega_phi,'k-')
hold on
grid on
plot(x_vec,omega_dy,'k--')
lgnd = legend(legend_titles,'location','sw');
xlabel('\eta\phi^2')
ylabel('\omega')
ylim([0 1])
xlim([0 1])
temp = [lgnd; lgnd.ItemText];
set(temp, 'FontSize', 10);
set(temp, 'FontName', 'Times New Roman');
set(gca, 'FontName', 'Times New Roman');