function tunnellinercallback

% Creator: Violeta Ivanova, Office of Educational Innovation and Technology.
% TUNNELLINERCALLBACK is a function called by COMPUTE tunnelliner.m.
% Developed to support course material in 1.383 Underground Construction.
% Faculty: Prof. Herbert H. Einstein, Civil & Environmental Engineering.

fig1 = openfig('TunnellinerGUI.fig', 'reuse');
figure(fig1)

E_ground          = findobj( fig1, 'Tag', 'E_ground');      % Elastic modulus of ground
nu_ground         = findobj( fig1, 'Tag', 'nu_ground');     % Poisson ratio of ground
K_ground          = findobj( fig1, 'Tag', 'K_ground');      % K ratio of ground
gamma_ground      = findobj( fig1, 'Tag', 'gamma_ground');  % Unit weight of ground
h_ground          = findobj( fig1, 'Tag', 'h_ground');      % Depth of tunnel
E_liner           = findobj( fig1, 'Tag', 'E_liner');       % Elastic modulus of liner
nu_liner          = findobj( fig1, 'Tag', 'nu_liner');      % Poisson ratio of liner
R_liner           = findobj( fig1, 'Tag', 'R_liner');       % Radius of tunnel
theta_liner       = findobj( fig1, 'Tag', 'theta_liner');   % Angle (?) of tunnel    
tmin_liner        = findobj( fig1, 'Tag', 'tmin_liner');    % Min thickness of liner
tmax_liner        = findobj( fig1, 'Tag', 'tmax_liner');    % Max thickness of liner
t_result          = findobj( fig1, 'Tag', 't_result');      % Thickness of liner for computations
T_result          = findobj( fig1, 'Tag', 'T_result');      % Normalized Thrust
M_result          = findobj( fig1, 'Tag', 'M_result');      % Normalized Moment
P_result          = findobj( fig1, 'Tag', 'P_result');      % Ground pressure without tunnel
sT_result         = findobj( fig1, 'Tag', 'sT_result');     % Liner stress due to T 
sM_result         = findobj( fig1, 'Tag', 'sM_result');     % Liner stress due to M
slip_popup        = findobj( fig1, 'Tag', 'slip_popup');        
compute_button    = findobj( fig1, 'Tag', 'compute_button');

slip                = get( slip_popup,  'Value' );     % Gives 1 for "slip" and 2 for "no slip"
E_ground_text       = get( E_ground,    'String');     % Gets input for E of ground as string
nu_ground_text      = get( nu_ground,   'String');     % Gets input for Poisson ratio of ground as string
K_ground_text       = get( K_ground,    'String');     % Gets input for K ratio of ground as string
gamma_ground_text   = get( gamma_ground,'String');     % Gets input for gamma (unit weight) of ground as string
h_ground_text       = get( h_ground,    'String');     % Gets input for depth h of tunnel as string

E         = str2num(E_ground_text);          % Converts string to number in MPa for E of ground
nu        = str2num(nu_ground_text);         % Converts string to number for Poisson ratio of ground
K         = str2num(K_ground_text);          % Converts string to number for K ratio of ground
gamma     = str2num(gamma_ground_text);      % Converts string to number for gamma in KN/cu.m of ground
h         = str2num(h_ground_text);          % Converts string to number in m for depth h of tunnel of ground

E_liner_text       = get( E_liner,    'String');     % Gets input for E of liner as string
nu_liner_text      = get( nu_liner,   'String');     % Gets input for Poisson ratio of liner as string
R_liner_text       = get( R_liner,    'String');     % Gets input for radius R of tunnel as string
theta_liner_text   = get( theta_liner,'String');     % Gets input for theta angle of tunnel as string
tmin_liner_text    = get( tmin_liner, 'String');     % Gets input for thickness t of liner as string
tmax_liner_text    = get( tmax_liner, 'String');     % Gets input for thickness t of liner as string

E1         = str2num(E_liner_text);          % Converts string to number for E (MPa) of liner
nu1        = str2num(nu_liner_text);         % Converts string to number for Poisson ratio of liner
R          = str2num(R_liner_text);          % Converts string to number for radius R (m) of tunnel
theta      = str2num(theta_liner_text);      % Converts string to number for theta angle (deg) of tunnel
tmin       = str2num(tmin_liner_text);       % Converts string to number in cm for thickness t (cm) of liner
tmax       = str2num(tmax_liner_text);       % Converts string to number in cm for thickness t (cm) of liner

%%% COMPUTE SIGMA_T AND SIGMA_M V. THICKNESS %%% 
t = linspace(tmin, tmax, 50);

% Compute Compressibility Ratio C 
% Reference: C* in Einstein & Schwartz 1979
A = t; % Cross-sectional area per unit length of tunnel. A (sq.m.)
C = E*R*(1-nu1^2)./(E1.*A*(1-nu^2));  % dimensionless/m

% Compute Flexibility Ratio F 
% Reference: F* in Einstein & Schwartz 1979
I = t.^3/12; % Moment of intertia per unit length of tunnel. I (m^4)
F = E*(R^3)*(1-nu1^2)./(E1.*I*(1-nu^2));   % dimensionless/m

% Compute coefficients
% Reference: Eqs. 20 in Einstein & Schwartz 1979.
a0      = C.*F*(1-nu) ./ (C + F + C.*F*(1-nu));

% Compute ground stress P
P = gamma * h /1000; % ground stress P in MPa

% Compute normalized thrust T/(PR) and moment M/(PR^2) 
% Reference: Eqs. 22 and 23 in Einstein & Schwartz 1979.

if slip == 1 % For full slip conditions   
    a2 = (F+6)*(1-nu) ./ (2*F*(1-nu) + 6*(5-6*nu));
    T_norm = 0.5*(1+K)*(1-a0) + 0.5*(1-K)*(1-2*a2)*cos(2*theta*pi/180);
    M_norm = 0.5*(1-K)*(1-2*a2)*cos(2*theta*pi/180);
    
elseif slip == 2 % For no-slip conditions 
    beta    = ((6+F).*C*(1-nu) + 2.*F*nu) ./ (3.*F + 3.*C + 2.*C.*F*(1-nu));
    b2      = C*(1-nu) ./ (2.*(C*(1-nu) + 4*nu - 6*beta - 3*beta.*C*(1-nu)));
    a2      = beta.*b2;
    T_norm = 0.5*(1+K)*(1-a0) + 0.5*(1-K)*(1+2*a2)*cos(2*theta*pi/180);
    M_norm = 0.25*(1-K)*(1-2*a2+2*b2)*cos(2*theta*pi/180);
end

% Compute Thrust T and Moment M 
T = T_norm*P*R*1000;    % T (kN)
M = M_norm*P*R^2*1000;  % M (kN.m)

% Compute stresses due to Thrust T and Moment M 
sT = (T./t)/1000;        % sT (MPa)
sM = (6.*M./t.^2)/1000;  % abs(sM) (MPa)

% Compute stresses in inner and outer fiber of liner sIn, sOut
% sIn  = sT - sM;      % sIn (MPa)
% sOut = sT + sM;      % sOut (MPa)

% Convert P, T and M to strings to display in the GUI
set( t_result,   'String', [tmin_liner_text, ' ']);           % set t for computations to tmin 
set( P_result,   'String', [num2str(P), ' ']);                % num2str(P) converts rounded P to a string
set( T_result,   'String', [num2str(round(T(1))), ' ']);      % num2str(T) converts T to a string
set( M_result,   'String', [num2str(round(M(1))), ' ']);      % num2str(M) converts M to a string
set( sT_result,  'String', [num2str(sT(1), 3), ' ']);            % num2str(sT) converts sT to a string
set( sM_result,  'String', [num2str(sM(1), 3), ' ']);            % num2str(sM) converts sM to a string
% set( sIn_result, 'String', [num2str(round(sIn(1))), ' ']);           % num2str(sIn) converts sIn to a string
% set( sOut_result,'String', [num2str(round(sOut(1))), ' ']);          % num2str(sOut) converts sOut to a string

% Plot sT v. t amd abs(sM) v. t

figure(2);
set(2, 'Position', [50 50 1000 500]); 
set(2, 'Name', 'TUNNELLINER: Liner Stresses Due to Thrust and Moment', ...
       'Color', [0.95 0.87 0.73], ... 
       'NumberTitle', 'off');


% sT v. t
     subplot(121)
     p1 = plot(t, sT, 'ro', 'MarkerFaceColor','r');
     hold on
     p1 = plot(t, sM, 'bs');
     hold off
     legend ('Stress due to Thrust [MPa]', 'Stress due to Moment [MPa]');
     titlestr = ['STRESSES FOR THETA = ', num2str(theta), ' DEG'];
     title(titlestr, ...
                'FontWeight', 'bold', ...
                'FontSize', 12);
     xlabel('LINER THICKNESS [m]', 'FontSize', 12); 
     ylabel('STRESS [MPa]', 'FontSize', 12);


%%% COMPUTE SIGMA_IN AND SIGMA_OUT V. THETA %%%

% Coefficients for t = tmin
Ft  = F(1);
Ct  = C(1);
a0t = a0(1);

% theta ranges from 0 to 2*pi
Theta = [0 : pi/48 : 2*pi];

% Compute normalized thrust T/(PR) and moment M/(PR^2)
% Reference: Eqs. 22 and 23 in Einstein & Schwartz 1979.

if slip == 1 % For full slip conditions   
    a2 = (Ft+6)*(1-nu) / (2*Ft*(1-nu) + 6*(5-6*nu));
    T_norm = 0.5*(1+K)*(1-a0t) + 0.5*(1-K)*(1-2*a2)*cos(2*Theta);
    M_norm = 0.5*(1-K)*(1-2*a2)*cos(2*Theta);
    
elseif slip == 2 % For no-slip conditions 
    beta    = ((6+Ft)*Ct*(1-nu) + 2*Ft*nu) / (3*Ft + 3*Ct + 2*Ct*Ft*(1-nu));
    b2      = Ct*(1-nu) / (2*(Ct*(1-nu) + 4*nu - 6*beta - 3*beta*Ct*(1-nu)));
    a2      = beta*b2;
    T_norm  = 0.5*(1+K)*(1-a0t) + 0.5*(1-K)*(1+2*a2)*cos(2*Theta);
    M_norm  = 0.25*(1-K)*(1-2*a2+2*b2)*cos(2*Theta);
end

% Compute Thrust T and Moment M 
T = T_norm*P*R*1000;    % T (kN)
M = M_norm*P*R^2*1000;  % M (kN.m)

% Compute stresses due to Thrust T and Moment M 
sT = (T/tmin)/1000;      % sT (MPa)
sM = (6*M/tmin^2)/1000;  % sM (MPa)

% Compute stresses in inner and outer fiber of liner sIn, sOut
% sIn  = sT - sM;      % sIn (MPa)
% sOut = sT + sM;     % sOut (MPa)

% Plot sT and sM v. Theta

% Define axes for the polar plot
sTmax = max(sT);
sMmax = max(sM);
     
% Plot sT and sM v. theta in polar coordinates
     subplot(122);
     
     if (sTmax > sMmax)
        p1 = polar(Theta(1:25), sT(1:25), 'r.');   % sT v. Theta=0-pi/4
        hold on
        p2 = polar(Theta(1:25), sM(1:25), 'bs');  % sM v. Theta=pi/4
        hold on
        p3 = polar(Theta(25:97), sT(25:97), 'r-');  % sT v. Theta=pi/4-2*pi
        hold on
        p4 = polar(Theta(25:97), sM(25:97), 'b--'); % sM v. Theta=pi/4-2*pi
        hold off
        legend ('Stress due to Thrust for theta from 0 to 90 deg [MPa]', ...
                'Stress due to Moment for theta from 0 to 90 deg [MPa]', ...
                'Location', 'SouthOutside');
     else
        p2 = polar(Theta(1:25), sM(1:25), 'bs');  % sM v. Theta=pi/4
        hold on
        p1 = polar(Theta(1:25), sT(1:25), 'r.');   % sT v. Theta=0-pi/4
        hold on
        p3 = polar(Theta(25:97), sT(25:97), 'r-');  % sT v. Theta=pi/4-2*pi
        hold on
        p4 = polar(Theta(25:97), sM(25:97), 'b--'); % sM v. Theta=pi/4-2*pi
        hold off
        legend ('Stress due to Moment for theta from 0 to 90 deg', ...
                'Stress due to Thrust for theta from 0 to 90 deg', ...
                'Location', 'SouthOutside');
     end
    
     titlestr = ['STRESSES FOR LINER THICKNESS ', num2str(tmin), ' M'];
     title(titlestr, ...
                'FontWeight', 'bold', ...
                'FontSize', 12);

return