% This program corresponds to Appendix 7.2.1 and Section 3.2.7 of the thesis.
% The effective rotational stiffness contribution of a Timoshenko member to
%  a connecting column at its end is derived. The member has an axial force.
% The limit as the axial force approaches zero is solved in order to verify 
%  the result in Section 3.2.7
clc, clear, close all

% Governing differential equation
syms EI Rl tl Yj N kAG yi yj 
syms y(x) psi(x)
ode1 = -EI*diff(psi) == Rl*tl + Yj*x + N*(y-yi)
ode2 = Yj + N*psi == kAG*(-psi + diff(y))
odes = [ode1; ode2];
S = dsolve(odes);

% Simplify these expressions by hand and verify the solution by hand
y = simplify(S.y)
psi = simplify(S.psi)

% Boundary conditions
syms D L tu Ru
syms C1 C2 phi w Tl Tu
y(x) = yi - Rl*tl/N - Yj*x/N + C1*cos(w*phi/L*x) + C2*sin(w*phi/L*x);
psi(x) = -Yj/N + 1/w^2 * ( C2*w*phi/L*cos(w*phi/L*x) - C1*w*phi/L*sin(w*phi/L*x) );
eqn1 = y(0) == yi;
eqn2 = y(L) == yj;
eqn3 = psi(0) == Tl-tl;
eqn4 = psi(L) == Tu-tu;
eqn5 = tu*Ru + tl*Rl + N*(yj-yi) + Yj*L == 0;

% Solve for key variables [Yj,tu,tl,C1,C2] <- these won't show up in the expressions
[A,B] = equationsToMatrix([eqn1,eqn2,eqn3,eqn4,eqn5],[Yj,tu,tl,C1,C2]);
X = linsolve(A,B)
ti_f = simplify(X(3));
tj_f = simplify(X(2));

% Invert the rotational stiffness matrix
[rotate_matrix,b] = equationsToMatrix([X(3);X(2)]==[tl; tu],[Tl;Tu;yi;yj]);
R_matrix = simplify([Rl 0; 0 Ru]*rotate_matrix)
syms ri rj n
R_matrix = simplify(subs(R_matrix,[Rl Ru],3*EI/L*[ri/(1-ri) rj/(1-rj)]))

% Verify the solution
R_EB = subs(R_matrix,w,1) % Should match p. 4 of research notes pt. 3 (Rbeam with N using M&Wu method) -- seems to match
a1 = 3*(ri*(1-rj)+rj*(1-ri));
a2 = 9*ri*rj - (1-ri)*(1-rj)*(w*phi)^2;
a3 = 18*ri*rj + a1*(w*phi)^2;
DEN = 18*ri*rj + (a1-a2)*w*phi*sin(w*phi) - a3*cos(w*phi);
zeta = (w^2-1)*(a1*w*phi*sin(w*phi)+18*ri*rj*(1-cos(w*phi))) / DEN;
G11 = (w^2-1)*3*rj*sin(w*phi) / ...
    ( (1-rj)*(w*phi)^2*sin(w*phi) + 3*rj*(sin(w*phi)-w*phi*cos(w*phi)));
R11 = 3*EI/L*ri*w*phi*( (1-rj)*(w*phi)^2*sin(w*phi) + 3*rj*(sin(w*phi)-w*phi*cos(w*phi)))/ DEN * (1+G11)/(1+zeta);
zero11 = simplify(subs(R11-R_matrix(1,1),phi,L*(N/EI)^0.5)) % should be zero
simplify(limit(subs(R11,w,sqrt(1+n*phi^2)),phi,0)) % Matches Section 3.2.7

G12 = (w^2-1)*sin(w*phi)/(sin(w*phi)-w*phi);
R12 = 9*EI/L*ri*rj*w*phi*(w*phi-sin(w*phi))/DEN* (1+G12)/(1+zeta);
zero12 = simplify(subs(R12-R_matrix(1,2),phi,L*(N/EI)^0.5))% should be zero
simplify(limit(subs(R12,w,sqrt(1+n*phi^2)),phi,0))  % Matches Section 3.2.7

R13 = 3*EI/L^2*ri*w*phi*( (1-rj)*(w*phi)^2*sin(w*phi) + 3*rj*w*phi*(1-cos(w*phi)))/DEN * 1/(1+zeta);
zero13 = simplify(subs(R13-R_matrix(1,3),phi,L*(N/EI)^0.5))% should be zero
simplify(limit(subs(R13,w,sqrt(1+n*phi^2)),phi,0)) % Matches Section 3.2.7

% Get the limit corresponding to Euler-Bernoulli (no shear deformations)
syms vji wji
Rij = R11+R12*vji-R13*wji*L;
simplify(diff(limit(subs(Rij,w,sqrt(1+n*phi^2)),phi,0),n)) % this matches the Euler-Bernoulli derivation from Monforton & Wu