% This program corresponds to Section 3.3 of the thesis.
% The lateral stiffness equation for a frame containing n bays is solved,
%   and considers the axial deformations of the beams.
% Originally written by LZ ~ January 3, 2019
% Minor revision by TM ~ February 4, 2019
clc, clear, close all

% Define variables
n = 2
eqns =zeros(1, n);
syms Q
syms S1 B1 D1 Q1 eqn1
S_m = S1;                   %lateral stiffness of each column
B_m = B1;                   %stiffness of each beam
D_m = D1;                   %lateral displacement of each column
Q_m = Q1;                   %lateral force applied on each column
eqn_m = eqn1;               %define equation matrix
for i = 2:n
   syms(['S', num2str(i)]);
   S_m = [S_m, ['S', num2str(i)]];
   syms(['B', num2str(i)]);
   B_m = [B_m, ['B', num2str(i)]];
   syms(['D', num2str(i)])
   D_m = [D_m, ['D', num2str(i)]];
   syms(['Q', num2str(i)])
   Q_m = [Q_m, ['Q', num2str(i)]];
end
for i = 2:2*n
    eqn_m = [eqn_m, ['eqn', num2str(i)]];
end

% Define equation set
for i = 1:n;
    eqn_m(i) = D_m(i) == Q_m(i)/S_m(i);
    if i < n
        eqn_m(i+n) = D_m(i)-D_m(i+1) == (Q-sum(Q_m(:, 1:i), 2))/B_m(i);
    else
        eqn_m(i*2) = Q == sum(Q_m, 2);
    end
end

% Solve linear equations
[A,B] = equationsToMatrix(eqn_m, [D_m Q_m]);
X = linsolve(A,B);
X = simplify(X);
D_vec = X(1:length(X)*0.5)
Q_vec = X(length(X)*0.5+1:length(X));
S_vec = Q./D_vec;

% The results were checked against the solution obtained via the equivalent
% spring system method and found to be identical (cases n = 2 through n = 5
% were checked)