% problem 4.27
%
I = 200*10^-6 %m^4
A = 1000000 %m^2
E = 200*10^9 %Pa
Lab = sqrt(4^2+3^2) %m
Lbc = 4 %m
L = [Lab Lbc]'
type = 2
%
ndof = 9
fdof = ndof-6
nelem = 2
%
% Set up vectors to assemble
phi = [atan(4/3) 0] 
Xinfo = [4 5 6 1 2 3;1 2 3 7 8 9]
%
Ksys = zeros(ndof)
for i=1:nelem
	k = swffrmstiff(L(i),E,I,A)
	T = swfgetT(phi(i),type)
	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
Pf = zeros(fdof,1)
%
% member forces
% Assembly member load vector
Pmember = zeros(ndof,1)
Loads = [6 9/2 9*5/8 6 9/2 -9*5/8; 0 12*4/2 12*4^2/12 0 12*4/2 -12*4^2/12]*1000
for i=1:2
	T = swfgetT(phi(i),type)
   Pm = T'*[Loads(i,:)]'
	Xmap = Xinfo(i,:)
	X = swfgetX(Xmap,ndof)
	Pmember = Pmember + X'*Pm
	pause   
end
%
%
Pfm = Pmember(1:fdof)
Psm = Pmember(fdof+1:ndof)
%
% coordinate displacements
Df = inv(Kff)*(Pf - Kfs*Ds - Pfm)
%
% reaction forces
Ps = Ksf*Df + Kss*Ds + Psm
%
D = [Df' Ds']'
% member forces
for i=1:nelem
	k = swffrmstiff(L(i),E,I,A)
	T = swfgetT(phi(i),type)
	Xmap = Xinfo(i,:)
	X = swfgetX(Xmap,ndof)
   X*D
   d = T*X*D
   s = k*d + Loads(i,:)'
   beams(i,:) = s'
	pause   
end
%
%

['The following displacements are in global coordinates.']
['Note that translations are close to zero.']
Df

['The following forces are in member coordinates. Column i = Member i.']
beams'./1000

['The reaction forces are in global coordinates']
Ps./1000