% This program corresponds to Appendix 7.2.3 of the thesis.
% The shape parameters corresponding to the effective rotational stiffness
%   contribution of a member towards the columns at its ends are derived
% Note that only the derivation for Euler-Bernoulli members is presented as
%   as on this solution the Timoshenko solution is trivial and was obtained
%   via inspection (compare with the solution in Section 3)
clc, clear, close all

% Goverening differential equation (Euler-Bernoulli)
syms C1 C2 phi P D L Q tl Rl tu Ru x EI phi
y(x) = C1*cos(phi*x/L) + C2*sin(phi*x/L) + tl*Rl/P - Q*x/P;
dy(x) = diff(y,x);

% Substitute boundary conditions and external moment equilibrium
eqn1 = y(0) == 0;
eqn2 = y(L) == D;
eqn3 = dy(0) == tl;
eqn4 = dy(L) == tu;
eqn5 = tu*Ru + tl*Rl == Q*L + P*D;

% Solve system of equations
[A,B] = equationsToMatrix([eqn1,eqn2,eqn3,eqn4,eqn5],[D,tu,tl,C1,C2]);
X = linsolve(A,B)
syms ru rl
X = subs(X,[Rl Ru],3*EI/L*[rl/(1-rl) ru/(1-ru)]);
D_f = simplify(X(1))
tu_f = simplify(X(2))
tl_f = simplify(X(3))

% Verify the hand-simplified solutions
a1 = 3*(ru*(1-rl)+rl*(1-ru));
a2 = 9*ru*rl-(phi)^2*(1-ru)*(1-rl);
a3 = 18*rl*ru + a1*(phi)^2;
tu_sol = Q*L^2/(phi*EI)*((1-ru)*(1-rl)*phi*sin(phi) + 3*rl*(1-ru)*(1-cos(phi)))/...
    ( a1*phi*cos(phi) + a2*sin(phi)  );
tl_sol = Q*L^2/(phi*EI)*((1-ru)*(1-rl)*phi*sin(phi) + 3*ru*(1-rl)*(1-cos(phi)))/...
    ( a1*phi*cos(phi) + a2*sin(phi)  );
delta_sol = Q*L^3/(EI*phi^3)*(18*rl*ru - a3*cos(phi) + (a1-a2)*phi*sin(phi))/...
    ( a1*phi*cos(phi) + a2*sin(phi)  );
zero = simplify(subs(tu_sol-tu_f,phi,sqrt(P/EI)*L)) % should be zero
zero = simplify(subs(tl_sol-tl_f,phi,sqrt(P/EI)*L)) % should be zero
zero = simplify(subs(delta_sol-D_f,phi,sqrt(P/EI)*L)) % should be zero

% Solve the limits for first order analysis
vul = tu_sol/tl_sol;
vul_0 = simplify(limit(vul,phi,0))
vlu = tl_sol/tu_sol;
vlu_0 = simplify(limit(vlu,phi,0))
wul = delta_sol/L/tl_sol;
wul_0 = simplify(limit(wul,phi,0))
wlu = delta_sol/L/tu_sol;
wlu_0 = simplify(limit(wlu,phi,0))