% This program corresponds to Section 4.2.1 of the thesis.
% The deformation equation for a Euler-Bernoulli column with initial imperfections
%    is derived.
clc, clear, close all

% Governing differential equation
syms C1 C2 phi 
syms L d D F S 
syms P PI tl tu 
syms Ru Rl EI ru rl
syms y Ml Mu
eqn(1) = 0 == C1 + Ml/P;
eqn(2) = 0 == C1*cos(phi) + C2*sin(phi) - Mu/P;
eqn(3) = tl == C2*phi/L - (d*phi^2*PI/L)/(phi^2-PI^2) - D/L - S/P;
eqn(4) = tu == C2*phi/L*cos(phi) - C1*phi/L*sin(phi) + (d*phi^2*PI/L)/(phi^2-PI^2) - D/L - S/P;
eqn(5) = Mu + Ml == P*(D+F)+S*L;
eqn = subs(eqn, [Ml Mu], [Rl*tl Ru*tu]);
eqn = subs(eqn, [Rl Ru], 3*EI/L*[rl/(1-rl) ru/(1-ru)]);
vec = [S tl tu C1 C2];

% Solve for S and lower end rotation
X = solve(eqn,vec)
X.S
X.tl

% Simplify the solution for X.tl and verify the solution
tl_sol =( F*phi^2*(PI^2-phi^2)*(1-rl)*((1-ru)*phi*sin(phi) + 3*ru*(1-cos(phi))) +...
         d*PI*phi^3*(1-rl)*(phi^2*sin(phi)*(1-ru) + 3*ru*(2*sin(phi) - phi*(1+cos(phi)))) ) /...
         ( L*(PI^2-phi^2)*(phi^3*sin(phi)*(1-ru)*(1-rl) + 3*phi*(2*rl*ru-ru-rl)*(phi*cos(phi)-sin(phi)) + 9*rl*ru*(2-2*cos(phi)-phi*sin(phi))) );
tl_sol = subs(tl_sol,phi,sqrt(P/EI)*L);
zero = subs(X.tl,phi,sqrt(P/EI)*L)-tl_sol;
simplify(zero)
test = subs(zero, [F EI L P ru rl D d PI], [1 2E11*1.29E-4 4.877 5E5 0.716912189334201 1 4.877/500 0.1 3.1415926536]);
test = double(test)
     
% Simplify the solution for X.S and verify the solution
test = X.S;
test = subs(test, phi, sqrt(P/EI)*L);
test = subs(test, [F EI L P ru rl D d PI], [1 2E11*1.29E-4 4.877 5E5 0.716912189334201 1 4.877/500 0.1 3.1415926536]);
test = double(test)
syms beta a1 a2 a3 Xi
S_eq = 12*EI/L^3*beta - P/L*D/F - P/L*d/F*Xi;
%Y = S_eq*F;
subs(S_eq,beta, 1/12*phi^3*(a1*phi*cos(phi)+a2*sin(phi))/(18*rl*ru-a3*cos(phi)+(a1-a2)*phi*sin(phi)));
subs(S_eq,Xi, 3*phi^3*PI*sin(phi)*(ru-rl)/(PI^2-phi^2)/(18*rl*ru-a3*cos(phi)+(a1-a2)*phi*sin(phi)));
subs(S_eq,a3, 18*rl*ru + a1*phi^2);
subs(S_eq,a1, 3*rl +3*ru - 6*ru*rl);
subs(S_eq,a2, 9*rl*ru-(1-ru)*(1-rl)*phi^2);
S_eq = subs(test, [F EI L P ru rl D d PI], [1 2E11*1.29E-4 4.877 5E5 0.716912189334201 1 4.877/500 0.1 3.1415926536]);
S_eq = double(S_eq) % this matched the results of a numerical example

% Verify the solutions for C1 and C2
X.C1
X.C2
format long
a1 = 3*rl+3*ru-6*rl*ru;
a2 = 9*rl*ru - phi^2*(1-ru)*(1-rl);
a3 = 18*rl*ru+a1*phi^2;
C_denom = 18*rl*ru - a3*cos(phi)+(a1-a2)*phi*sin(phi);
C1_eq = 3*rl/C_denom*(...
    d*PI/(PI^2-phi^2)*(3*phi^2*ru*(1+cos(phi)) - 6*phi*ru*sin(phi) - phi^3*sin(phi)*(1-ru)) -...
    F*( 3*ru*(1-cos(phi)) + phi*sin(phi)*(1-ru) ) );
C2_eq = 3/C_denom*(...
    d*PI/(PI^2-phi^2)*(3*phi^2*rl*ru*sin(phi) - 6*phi*rl*ru*(1-cos(phi)) - phi^3*(ru*(1-rl) - rl*cos(phi)*(1-ru))) + ...
    F*(phi*(ru*(1-rl)+rl*cos(phi)*(1-ru))+3*rl*ru*sin(phi))...
    );

zero = simplify( subs(X.C1,phi,sqrt(P/EI)*L)-subs(C1_eq,phi,sqrt(P/EI)*L) )
zero = simplify( subs(X.C2,phi,sqrt(P/EI)*L)-subs(C2_eq,phi,sqrt(P/EI)*L) )