% This program corresponds to Section 3.2.2 of the thesis.
% The lateral stiffness of a Timoshenko column is derived under the
%   assumption of Haringx (1948) 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*diff(y) == kAG*(-psi + diff(y))
odes = [ode1; ode2];
S = dsolve(odes);

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

% Verify the solution for y(x), psi(x)
syms C2 C1
phiL = sqrt(P/EI);
w = 1/sqrt(1-P/kAG);
y(x) = Rl*tl/P - Q*x/P + C1*cos(x*w*phiL) + C2*sin(x*w*phiL);
psi(x) = -Q/P - 1/w^2*(C1*phiL*w*sin(phiL*w*x) - C2*phiL*w*cos(phiL*w*x));
simplify(EI*diff(psi) - Rl*tl + Q*x + P*y)
simplify(Q+P*diff(y) - kAG*(-psi+diff(y)))

syms D L tu Ru
%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 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 == -Q/P + 1/w^2*(C2*phi*w/L)
eqn4 = tu == -Q/P - 1/w^2*(C1*phi*w/L*sin(phi*w) - C2*phi*w/L*cos(phi*w));
eqn5 = tu*Ru + tl*Rl == P*D + Q*L;

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

% Manually simplify Q from above, and write the solution below to verify that it matches the above expression
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 =  (phi^2)/12 * ( a1*wphi*cos(wphi) + a2*sin(wphi) ) / ( 18*rl*ru*w/phi - (18*rl*ru*w/phi + a1*wphi)*cos(wphi) + (a1*w^2-a2)*sin(wphi) );
%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 * D;
Q = subs(Q,phi,sqrt(P/EI)*L);
Q_result = subs(Q_result,phi,sqrt(P/EI)*L);

% Compare results
error = simplify(Q_result-Q) % Should be zero
%error_val = subs(error,w, 1/sqrt(1-P/kAG));
%error_val = double(subs(error_val,[P,EI,L,Rl,Ru,kAG,D],[5E5,2E7,5,4E6,3E6,77E9*2E-3,1]))