% This program corresponds to Section 3.2.2 of the thesis
% The limits for zeta and beta in the Timoshenko solution of the lateral stiffness
%   are solved in the absence of lateral load (phi -> 0 or kAG -> inf)
clc, clear, close all
syms P kAG rl ru phi EI L
w = sqrt(1+(phi^2*EI/L^2/kAG));
a1 = 3*(ru*(1-rl)+rl*(1-ru));
a2 = 9*ru*rl-(w*phi)^2*(1-ru)*(1-rl);
a3 = 18*rl*ru + a1*(w*phi)^2;
wphi = phi*w;

% Solve the limit for zeta0
zeta = (w^2-1)*(18*rl*ru*(1-cos(wphi)) + a1*wphi*sin(wphi))/   ( 18*rl*ru - a3*cos(wphi) + (a1-a2)*wphi*sin(wphi) );
zeta_0 = simplify(limit(zeta,phi,0)) % Get limit of zeta
zeta_0_noshear = simplify(limit(zeta_0,kAG,inf))

% Solve the limit for beta0
beta = (w*phi^3)/12 * ( a1*wphi*cos(wphi) + a2*sin(wphi) ) / ( 18*rl*ru - a3*cos(wphi) + (a1-a2)*wphi*sin(wphi) );
beta_0 = simplify(limit(beta,phi,0)) % Get limit of beta
beta_0_noshear = simplify(limit(beta_0,kAG,inf))