% problem 2.26
%
Av = .76
Ah = .43
Ad = 1.01
%
E = 29*10^3
Lv = 70
Lh = 48
Ld = sqrt(Lv^2 + Lh^2)
%
Ksys = zeros(10)
ndof = 10
fdof = 7
%
% Set up vectors to assemble
Area = [Av Av Ah Ah Ad Ad Av Ad]
Length = [Lv Lv Lh Lh Ld Ld Lv Ld]
E = E*ones(1,8)
phi = [pi/2 pi/2 0 0 -atan(Lv/Lh) -atan(Lv/Lh) pi/2 atan(Lv/Lh)] 
Xinfo = [8 1 2 3; 2 3 9 10; 8 1 4 5; 2 3 6 7; 2 3 4 5; 9 10 6 7; 4 5 6 7; 8 1 6 7]
%%
Ksys = zeros(ndof)
for i=1:8
	k = swftrusstiff(Length(i),E(i),Area(i))
	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)
%
Ds = zeros(ndof-fdof,1)
%%%
%%%
% Assembly load vector
Lphi = [0 350 0 290 120]*pi/180
Load = [0 1 0 1 1]
XLinfo = [1 8; 2 3; 9 10; 4 5; 6 7]
Psys = zeros(ndof,1)
for i=1:5
	T = swfgetT(Lphi(i),1)
   r = T(1:2,1:2)
   Pg = r'*[Load(i) 0]'
	Xmap = XLinfo(i,:)
	X = swfgetX(Xmap,ndof)
	Psys = Psys + X'*Pg
	pause   
end
%
% member forces
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 + Psys(fdof+1:ndof)
%
D = [Df' Ds']'
['member forces']
for i=1:8
	k = swftrusstiff(Length(i),E(i),Area(i))
	T = swfgetT(phi(i),1)
	Xmap = Xinfo(i,:)
	X = swfgetX(Xmap,ndof)
   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