% This program corresponds to Section 6.2.6 of the thesis.
% The rotational buckling load of a Euler-Bernoulli column with three segments is derived in following the Hoblit (1951) method.
clc, clear, close all

% Write out Hoblit (1951) equations for the column
syms zu z1 z2 zL
syms tu t1 t2 tL
syms Mu M1 M2 ML
syms Lu Lm LL L
syms pu pm pL
syms lbu lbm lbL
syms P S
syms Ru RL
syms E I
syms F
eqn(1) = z1 == zL+ ML/P*(1-cos(pL*LL/L)) + (tL-S/P)/(pL/L)*sin(pL*LL/L) + S*LL/P;
eqn(2) = z2 == z1+ M1/P*(1-cos(pm*Lm/L)) + (t1-S/P)/(pm/L)*sin(pm*Lm/L) + S*Lm/P;
eqn(3) = zu == z2+ M2/P*(1-cos(pu*Lu/L)) + (t2-S/P)/(pu/L)*sin(pu*Lu/L) + S*Lu/P;
eqn(4) = t1 == ML*pL/(P*L)*sin(pL*LL/L) + (tL-S/P)*cos(pL*LL/L) + S/P;
eqn(5) = t2 == M1*pm/(P*L)*sin(pm*Lm/L) + (t1-S/P)*cos(pm*Lm/L) + S/P;
eqn(6) = tu == M2*pu/(P*L)*sin(pu*Lu/L) + (t2-S/P)*cos(pu*Lu/L) + S/P;
eqn(7) = M1 == ML*cos(pL*LL/L) - (tL-S/P)*(P*L/pL)*sin(pL*LL/L);
eqn(8) = M2 == M1*cos(pm*Lm/L) - (t1-S/P)*(P*L/pm)*sin(pm*Lm/L);
eqn(9) = Mu == M2*cos(pu*Lu/L) - (t2-S/P)*(P*L/pu)*sin(pu*Lu/L);
eqn(10) = zL == F;
eqn(11) = zu == 0;
eqn(12) = ML == +RL*tL;
%eqn(13) = Mu == -Ru*tu; % this is the final equation, but the thirteen equations can't be solved unless the first twelve are solved for the shear force

% Solve for the shear force (X.S is designated as 'V' in the thesis)
vec = [S zu z2 z1 zL tu t2 t1 Mu M2 M1 ML];
X = solve(eqn,vec)
X.S


