% Solution to Problem 2.23
%
Aac = 3 % in
Abc = 5
Acd = 2
Abd = Aac
Aad = Abc
%
E = 1
L = 1
A = 1
%
ndof = 10
fdof = 2
Ksys = zeros(ndof)
%
% Set up vectors to assemble
phi = [30 60 90 180-30]*pi/180
Xinfo = [1 2 3 4; 1 2 5 6; 1 2 7 8;1 2 9 10]
pause
%%
Ksys = zeros(ndof)
for i=1:4
	k = swftrusstiff(L,E,A)
	T = swfgetT(phi(i),1)
	kg = T'*k*T
	Xmap = Xinfo(i,:)
	X = swfgetX(Xmap,ndof)
	Ksys = Ksys + X'*kg*X
	pause   
end
%%
%%
Kff = Ksys(1:fdof,1:fdof)
Kfs = Ksys(1:fdof,fdof+1:ndof)
Ksf = Ksys(fdof+1:ndof,1:fdof)
Kss = Ksys(fdof+1:ndof,fdof+1:ndof)
pause
%
Ds = zeros(ndof-fdof,1)
%%%
%%%
% Assembly load vector
Load = [1 -1]
%
% member forces
Psys = Load'
Pfm = zeros(fdof,1)
Psm = zeros(ndof-fdof,1)
%
% coordinate displacements
Pf = Psys(1:fdof)
Df = inv(Kff)*(Pf - Kfs*Ds - Pfm)
%
% reaction forces
Ps = Ksf*Df + Kss*Ds + Psm
%
D = [Df' Ds']'
% member forces
for i=1:4
	k = swftrusstiff(L,E,A)
	T = swfgetT(phi(i),1)
	Xmap = Xinfo(i,:)
	X = swfgetX(Xmap,ndof)
   Ksys = Ksys + X'*kg*X
   X*D
   d = T*X*D
   s = k*d
   bars(i) = s(3)
	pause   
end
%
['The following forces are multiplied by P.']
['positive is tension']
['negative is compression']
bars 
%