% This program corresponds to Appendix 4.2.5 of the thesis.
% The deformation equation for a Timoshenko column with initial imperfections
%    is derived.
clc, clear, close all

% Governing differential equations
syms EI Rl tl Q N kAG D0 d0 PI L
syms y(x) psi(x)
ode1 = EI*diff(psi) == Rl*tl - Q*x - N*(y+d0*sin(PI/L*x)+D0*x/L)
ode2 = Q + N*(psi + d0*PI/L*cos(PI/L*x) + D0/L) == kAG*(-psi + diff(y))
odes = [ode1; ode2];
S = dsolve(odes);
y = simplify(S.y)
psi = simplify(S.psi)
pretty(y)
pretty(psi)

% Simplify y(x) and psi(x) by hand, and verify the solution by hand
syms D tu Ru
syms C1 C2 phi w n
y(x) = C1*cos(w*phi/L*x) + C2*sin(w*phi/L*x) + tl*Rl/N + ...
    d0*sin(PI*x/L)*(phi^2*(n*PI^2+w^2)/(PI^2-w^2*phi^2)) - D0/L*x - Q/N*x;
psi(x) = 1/w^2*(C2*w*phi/L*cos(w*phi/L*x) - C1*w*phi/L*sin(w*phi/L*x))...
    -D0/L - Q/N + d0*PI/L*cos(PI*x/L)*w^2*phi^2/(PI^2-phi^2*w^2); 

% Boundary conditions
%y(x) = subs(y,[n w],[0 1])
%psi(x) = subs(psi,[n w],[0 1])
eqn1 = y(0) == 0;
eqn2 = y(L) == D;
eqn3 = psi(0) == tl;
eqn4 = psi(L) == tu;
eqn5 = tu*Ru + tl*Rl == N*(D+D0) + Q*L;

% Solve the for the deflection, integration constants and end rotations
[A,B] = equationsToMatrix([eqn1,eqn2,eqn3,eqn4,eqn5],[D,C1,C2,tu,tl]);
X = linsolve(A,B)

%% Verify all of the solutions
syms rl ru
% Do these before the expressions get redefined
Delta_sol = subs(X(1),PI,pi);
Delta_sol = subs(Delta_sol,w,sqrt(1+phi^2*n));
Delta_sol = subs(Delta_sol,n,EI/L^2 / kAG);
Delta_sol = subs(Delta_sol,phi,sqrt(N/EI)*L);
Delta_sol = subs(Delta_sol,[rl ru], [1/(1+3*EI/L/Rl) 1/(1+3*EI/L/Ru)]);
C1_sol = subs(X(2),PI,pi);
C1_sol = subs(C1_sol,w,sqrt(1+phi^2*n));
C1_sol = subs(C1_sol,n,EI/L^2 / kAG);
C1_sol = subs(C1_sol,phi,sqrt(N/EI)*L);
C1_sol = subs(C1_sol,[rl ru], [1/(1+3*EI/L/Rl) 1/(1+3*EI/L/Ru)]);
C2_sol = subs(X(3),PI,pi);
C2_sol = subs(C2_sol,w,sqrt(1+phi^2*n));
C2_sol = subs(C2_sol,n,EI/L^2 / kAG);
C2_sol = subs(C2_sol,phi,sqrt(N/EI)*L);
C2_sol = subs(C2_sol,[rl ru], [1/(1+3*EI/L/Rl) 1/(1+3*EI/L/Ru)]);
tu_sol = subs(X(4),PI,pi);
tu_sol = subs(tu_sol,w,sqrt(1+phi^2*n));
tu_sol = subs(tu_sol,n,EI/L^2 / kAG);
tu_sol = subs(tu_sol,phi,sqrt(N/EI)*L);
tu_sol = subs(tu_sol,[rl ru], [1/(1+3*EI/L/Rl) 1/(1+3*EI/L/Ru)]);
tl_sol = subs(X(5),PI,pi);
tl_sol = subs(tl_sol,w,sqrt(1+phi^2*n));
tl_sol = subs(tl_sol,n,EI/L^2 / kAG);
tl_sol = subs(tl_sol,phi,sqrt(N/EI)*L);
tl_sol = subs(tl_sol,[rl ru], [1/(1+3*EI/L/Rl) 1/(1+3*EI/L/Ru)]);

% Check hand solution for Q
phi = sqrt(N/EI)*L;
n = EI/L^2 / kAG;
w = sqrt(1+phi^2*n);
a1 = 3*rl+3*ru-6*rl*ru;
a2 = 9*rl*ru - (w*phi)^2*(1-ru)*(1-rl);
a3 = 18*rl*ru+a1*(w*phi)^2;
zeta = (w^2-1)*(18*rl*ru*(1-cos(w*phi)) + a1*w*phi*sin(w*phi))/...
        (18*rl*ru-a3*cos(w*phi)+(a1-a2)*w*phi*sin(w*phi));
beta = (w*phi)*phi^2/12 * (a1*w*phi*cos(w*phi) + a2*sin(w*phi) ) /...
        (18*rl*ru-a3*cos(w*phi)+(a1-a2)*w*phi*sin(w*phi));
S = 12*EI/L^3*beta/(1+zeta);
chi = 3*(w*phi)^3*w^2*PI*sin(w*phi)*(ru-rl) /...
    (PI^2-(w*phi)^2) / ...
    (1+zeta) / ...
    (18*rl*ru-a3*cos(w*phi)+(a1-a2)*w*phi*sin(w*phi));
Delta = (Q + N*D0/L + N*d0/L*chi)/S;
Delta = subs(Delta,PI,pi);
Delta = subs(Delta,[rl ru], [1/(1+3*EI/L/Rl) 1/(1+3*EI/L/Ru)]);
Delta_zero = simplify(Delta_sol-Delta) % Delta 
subs(Delta_zero,d0,0)
subs(Delta_zero,D0,0) 
subs(Delta_zero,Q,0) % should be zero

% Check hand solution for C1
C1_mine = (N*D0/L+Q)*(w*L^3/(EI*phi^3)) *...
    (9*rl*ru*(cos(w*phi)-1)-3*w*phi*(1-ru)*rl*sin(w*phi)) / ...
    (a1*w*phi*cos(w*phi) + a2*sin(w*phi) ) + ...
 d0*PI*phi*w^3/(PI^2-phi^2*w^2)*...
    (3*rl*(1-ru)*w*phi*sin(w*phi) - 9*ru*rl*(cos(w*phi)+1))/...
    (a1*w*phi*cos(w*phi) + a2*sin(w*phi) );
C1_mine = subs(C1_mine,PI,pi);
C1_mine = subs(C1_mine,[rl ru], [1/(1+3*EI/L/Rl) 1/(1+3*EI/L/Ru)]);
C1_zero = simplify(C1_sol-C1_mine) % C1
simplify(subs(C1_zero,d0,0))
simplify(subs(C1_zero,D0,0))
simplify(subs(C1_zero,Q,0)) % should be zero
double(subs(C1_zero,[Q D0 d0 N EI kAG Ru Rl L],[2e5 0.002 0.003 8e5 2e11*152e-6 0.44*4510e-6*76e9 1e6 2e6 4.285]))

% Check hand solution for C2
C2_mine = (N*D0/L+Q)*(w*L^3/(EI*phi^3)) *...
    (9*rl*ru*sin(w*phi) + 3*w*phi*ru*(1-rl) + 3*w*phi*rl*(1-ru)*cos(w*phi))/...
    (a1*w*phi*cos(w*phi) + a2*sin(w*phi) ) + ...
 d0*PI*phi*w^3/(PI^2-phi^2*w^2)*...
    (3*w*phi*ru*(1-rl)-3*w*phi*rl*(1-ru)*cos(w*phi)-9*ru*rl*sin(w*phi))/...
    (a1*w*phi*cos(w*phi) + a2*sin(w*phi) );
C2_mine = subs(C2_mine,PI,pi);
C2_mine = subs(C2_mine,[rl ru], [1/(1+3*EI/L/Rl) 1/(1+3*EI/L/Ru)]);
C2_zero = simplify(C2_sol-C2_mine) % C2
simplify(subs(C2_zero,d0,0))
simplify(subs(C2_zero,D0,0))
simplify(subs(C2_zero,Q,0)) % should be zero
double(subs(C2_zero,[Q D0 d0 N EI kAG Ru Rl L],[2e5 0.002 0.003 8e5 2e11*152e-6 0.44*4510e-6*76e9 1e6 2e6 4.285]))

% Check hand solution for tu
tu_mine = (w*L^2/(EI*phi))*(...
    (N*D0/L+Q)*(3*rl*(1-ru)*(1-cos(w*phi)) + (1-ru)*(1-rl)*w*phi*sin(w*phi))/...
                (a1*w*phi*cos(w*phi) + a2*sin(w*phi) ) - ...
    N*d0/L*(PI*w^2*phi^2)/(PI^2-w^2*phi^2) * ...
               (3*rl*(1-ru)*(1+cos(w*phi)) - (1-ru)*(1-rl)*w*phi*sin(w*phi))/...
                (a1*w*phi*cos(w*phi) + a2*sin(w*phi) )...
        );
tu_mine = subs(tu_mine,PI,pi);
tu_mine = subs(tu_mine,[rl ru], [1/(1+3*EI/L/Rl) 1/(1+3*EI/L/Ru)]);
tu_zero = simplify(tu_sol-tu_mine) % tu
simplify(subs(tu_zero,d0,0))
simplify(subs(tu_zero,D0,0))
simplify(subs(tu_zero,Q,0)) % should be zero
double(subs(tu_zero,[Q D0 d0 N EI kAG Ru Rl L],[2e5 0.002 0.003 8e5 2e11*152e-6 0.44*4510e-6*76e9 1e6 2e6 4.285]))
    
% Check hand solution for tl
tl_mine = (w*L^2/(EI*phi))*(...
    (N*D0/L+Q)*(3*ru*(1-rl)*(1-cos(w*phi)) + (1-ru)*(1-rl)*w*phi*sin(w*phi))/...
                (a1*w*phi*cos(w*phi) + a2*sin(w*phi) ) + ...
    N*d0/L*(PI*w^2*phi^2)/(PI^2-w^2*phi^2) * ...
               (3*ru*(1-rl)*(1+cos(w*phi)) - (1-ru)*(1-rl)*w*phi*sin(w*phi))/...
                (a1*w*phi*cos(w*phi) + a2*sin(w*phi) )...
        );
tl_mine = subs(tl_mine,PI,pi);
tl_mine = subs(tl_mine,[rl ru], [1/(1+3*EI/L/Rl) 1/(1+3*EI/L/Ru)]);
tl_zero = simplify(tl_sol-tl_mine) % tl
simplify(subs(tl_zero,d0,0))
simplify(subs(tl_zero,D0,0))
simplify(subs(tl_zero,Q,0)) % should be zero
double(subs(tl_zero,[Q D0 d0 N EI kAG Ru Rl L],[2e5 0.002 0.003 8e5 2e11*152e-6 0.44*4510e-6*76e9 1e6 2e6 4.285]))
    
    
    
