% Linear System Example Problem - Truss Analysis
%
% This script solves a truss problem: Plesha Statics v 2 problem 6.40
% First, the matrix is explicitly defined and solved. In exercise 1, an
% altered system is used with the same geometry but a different set of
% applied loads. Finally, in exercise 2, a demonstration of linear indexing
% is given as an alternative to the unwieldy method used in the example.
%
% Joe Schoneman - 20 Jan 2013 

clc;clear; % Clear the Command Window and the Workspace

%%%%%%%% Exaxmple: Linear System %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Set up the linear system: 

c = cosd(60); % Defined for convenience
s = sind(60); 

%    Ay  Ex  Ey  Tab Tac Tbc Tbd Tcd Tce Tde
A = [0   0   0   c   1   0   0   0   0   0; ... % Fa_x (1)
     1   0   0   s   0   0   0   0   0   0; ... % Fa_y (2)
     0   0   0  -c   0   c   1   0   0   0; ... % Fb_x (3)
     0   0   0  -s   0  -s   0   0   0   0; ... % Fb_y (4)
     0   0   0   0  -1  -c   0   c   1   0; ... % Fc_x (5)
     0   0   0   0   0   s   0   s   0   0; ... % Fc_y (6) [= P]
     0   0   0   0   0   0  -1  -c   0   c; ... % Fd_x (7)
     0   0   0   0   0   0   0  -s   0  -s; ... % Fd_y (8)
     0   1   0   0   0   0   0   0  -1  -c; ... % Fe_x (9)
     0   0   1   0   0   0   0   0   0   s];    % Fe_y (10)
     
P = 1; % Solve the system in terms of 'P' for now
p_vec = [0; 0; 0; 0; 0; 1; 0; 0; 0; 0]; % Applied force vector, 'P' in row 6

t1 = inv(A)*p_vec; % Literal syntax to solve (MATLAB throws a warning)
t2 = A\p_vec; % Matlab operator to solve

% The next few lines simply label the answers for convenience, they're not
% relevant to solving the system.

Labels = [' Ay'; ' Ex'; ' Ey'; 'Tab'; 'Tac'; 'Tbc'; 'Tbd'; 'Tcd'; 'Tce'; 'Tde'];
disp('Example Answers:');
for i = 1:length(t2);
    fprintf('%s: %.4f\n',Labels(i,:),t2(i)); % Google 'fprintf MATLAB' for more info
end
disp(' ');

%%%%%%%% Exercise 1: Altered System %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The only changes that must be made are to the 'p' vector:
% Row 3 = P/2; Row 6 = P; Row 7 = -P/3*cosd(theta); Row 8 = -P/3*sind(theta)

P = 5000;   % [lb]
theta = 30; % [deg]

% Define the new vector of applied forces
p_ex1 = [0; 0; P/2; 0; 0; P; -P/3*cosd(theta); -P/3*sind(theta); 0; 0];

% Solve the linear system
t_ex1 = A\p_ex1;

% Print Answers to Command Window:
disp('Exercise 1 Answers:');
for i = 1:length(t_ex1);
    fprintf('%s: %.4f\n',Labels(i,:),t_ex1(i));
end
disp(' ');

%%%%%%%% Exercise 2: Matrix Indexing %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 A1 = zeros(10); % Make a 10-by-10 matrix of 'zero' entries
 A1(1,[4 5]) = [c 1];               % Fa_x (1)
 A1(2,[1 4]) = [1 s];               % Fa_y (2)
 A1(3,[4 6 7]) = [-c c 1];          % Fb_x (3)
 A1(4,[4 6]) = [-s -s];             % Fb_y (4)
 A1(5,[5 6 8 9]) = [-1 -c c 1];     % Fc_x (5)
 A1(6,[6 8]) = [s s];               % Fc_y (6) [= P]
 
 % Using the 'end' syntax for demonstration
 A1(end-3,[end-3 end-2 end]) = [-1 -c c];  % Fd_x (7)
 A1(end-2,[end-2 end]) = [-s -s];          % Fd_y (8)
 A1(end-1,[2 end-1 end]) = [1 -1 -c];      % Fe_x (9)
 A1(end,[3 end]) = [1 s];                  % Fe_y (10);
 
 P = 1; % P was changed in exercise 1, change it back for comparison
 t_ex2 = A1\p_vec;
 
 % Print answers to Command Window
 disp('Exercise 2 Answers:');
for i = 1:length(t_ex2);
    fprintf('%s: %.4f\n',Labels(i,:),t_ex2(i));
end
disp(' ');

commandwindow % Bump focus to the command window

 %%%%%%% End of script %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%