% This program corresponds to Seciton 7.2.1 of the thesis.
% The effective rotational stiffness contribution of a member towards the
%   columns at its ends is derived. The member has an axial force.
clc, clear, close all

% Governing differential equation
syms C1 C2 phi yi yj ti tj Pi Pj Ri Rj Yj x N EI L Rbi Rbj
y(x) = yi + C1*cos(phi/L*x) + C2*sin(phi/L*x) - Pi*Ri/N - Yj*x/N;
dy(x) = diff(y,x);

% Substitute boundary conditions and external moment equilibrium equation
eqn1 = y(0) == yi;
eqn2 = y(L) == yj;
eqn3 = dy(0) == ti-Pi;
eqn4 = dy(L) == tj-Pj;
eqn5 = Pi*Ri + Pj*Rj + N*(yj-yi) + Yj*L == 0;

% Solve system of equations, expressed as a matrix
[A,B] = equationsToMatrix([eqn1,eqn2,eqn3,eqn4,eqn5],[Pi,Pj,C1,C2,Yj]);
X = linsolve(A,B)
Pi_f = simplify(X(1))
Pj_f = simplify(X(2))

% Invert the rotational stiffness matrix
[rotate_matrix,b] = equationsToMatrix(X(1:2)==[Pi; Pj],[ti;tj;yi;yj]);
R_matrix = simplify([Ri 0; 0 Rj]*rotate_matrix)
syms ri rj
R_matrix = simplify(subs(R_matrix,[Ri Rj],3*EI/L*[ri/(1-ri) rj/(1-rj)]))

% Verify the solution for each term individually and then the final solution as well
a1 = 3*(ri*(1-rj)+rj*(1-ri));
a2 = 9*ri*rj - (1-ri)*(1-rj)*phi^2;
a3 = 18*ri*rj + a1*phi^2;
R11 = 3*EI/L*ri*phi*( (1-rj)*phi^2*sin(phi) + 3*rj*(sin(phi)-phi*cos(phi)))/...
    ( 18*ri*rj + (a1-a2)*phi*sin(phi) - a3*cos(phi) );
zero11 = simplify(subs(R11-R_matrix(1,1),phi,L*(N/EI)^0.5)) % should be zero
simplify(limit(R11,phi,0))

R12 = 9*EI/L*ri*rj*phi*(phi-sin(phi))/...
    ( 18*ri*rj + (a1-a2)*phi*sin(phi) - a3*cos(phi) );
zero12 = simplify(subs(R12-R_matrix(1,2),phi,L*(N/EI)^0.5)) % should be zero
simplify(limit(R12,phi,0))

R13 = 3*EI/L^2*ri*phi*( (1-rj)*phi^2*sin(phi) + 3*rj*phi*(1-cos(phi)))/...
    ( 18*ri*rj + (a1-a2)*phi*sin(phi) - a3*cos(phi) );
zero13 = simplify(subs(R13-R_matrix(1,3),phi,L*(N/EI)^0.5)) % should be zero
simplify(limit(R13,phi,0))

vji = tj/ti;
wji = (yj-yi)/ti/L;
R_prime = 3*EI/L*ri*phi/(18*ri*rj+(a1-a2)*phi*sin(phi)-a3*cos(phi))*...
    ( (1-rj)*phi^2*sin(phi)*(1-wji) + 3*rj*(sin(phi)*(1-vji) - phi*cos(phi)*(1-wji) + phi*(vji-wji)) ); % express the final solution
M_over_t = R11 + R12*vji - R13*wji*L;
zero_sol = simplify(subs(R_prime-M_over_t,phi,L*(N/EI)^0.5)) % should be zero
