% This program corresponds to Section 3.2.2 of the thesis.
% The lateral stiffness of a Timoshenko column is derived under the
%   assumption of Engesser (1891) where the shear angle is dy/dx
clc, clear, close all

% Governing system of differential equations
syms EI Rl tl Q P kAG
syms y(x) psi(x)
ode1 = EI*diff(psi) == Rl*tl - Q*x - P*y
ode2 = Q + P*psi == kAG*(-psi + diff(y))
odes = [ode1; ode2];
S = dsolve(odes);

% Simplify these manually by hand
y = simplify(S.y)
psi = simplify(S.psi)

syms D L tu Ru
% Normally this would work but MATLAB didn't simplify to cos and sin
% Also C1 and C2 are defined differently if left in imaginary
%eqn1 = subs(y,x,0) == 0;
%eqn2 = subs(y,x,L) == D;
%eqn3 = subs(psi,x,0) == tl;
%eqn4 = subs(psi,x,L) == tu;
%eqn5 = tu*Ru + tl*Rl == P*D + Q*L;

% Substitute the boundary conditions
syms C1 C2 phi w
eqn1 = 0 == Rl*tl/P + C1;
eqn2 = D == Rl*tl/P - Q*L/P + C1*cos(w*phi) + C2*sin(w*phi);
eqn3 = tl == 1/w^2*( C2*phi*w/L - Q*(1/kAG + 1/P) );
eqn4 = tu == 1/w^2*( C2*phi*w/L*cos(w*phi) - C1*phi*w/L*sin(w*phi) - Q*(1/kAG+1/P) );
eqn5 = tu*Ru + tl*Rl == P*D + Q*L;

% Solve the system for Q
[A,B] = equationsToMatrix([eqn1,eqn2,eqn3,eqn4,eqn5],[Q,tu,tl,C1,C2]);
X = linsolve(A,B)
Q = simplify(X(1))

% Simplify the above solution for Q and substitute the result to verify its correctness
ru = 1/(1+3*EI/(Ru*L));
rl = 1/(1+3*EI/(Rl*L));
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;
beta = (w*phi^3)/12 * ( a1*wphi*cos(wphi) + a2*sin(wphi) ) / ( 18*rl*ru - a3*cos(wphi) + (a1-a2)*wphi*sin(wphi) );
zeta = P/kAG*(18*rl*ru*(1-cos(wphi)) + a1*wphi*sin(wphi))/   ( 18*rl*ru - a3*cos(wphi) + (a1-a2)*wphi*sin(wphi) );
Q_result = 12*EI/L^3*beta/(1+zeta) * D;
Q = subs(Q,phi,sqrt(P/EI)*L);
Q_result = subs(Q_result,phi,sqrt(P/EI)*L);

% Compare results for any example:
error = simplify(Q_result-Q) % Should be zero