% *******************************************************************************************************
%                               						CAPSLD.M
% *******************************************************************************************************

% Distribution of lateral forces between planar substructures in  Multistory systems

% **************
%      Input
% **************

% jc = vector cell array with the joints at the ends of each bar for each one of the different
% types of frames involved. The joint numbers at the foundation are taken as zero. The nodes can be
% selected in any order but the first 'ns' nodes (where ns is the number of stories) must be along a
% column line that streches the full height of the structure. For example, in a 3 story frame 
% nodes 1 2 and 3 are in the 1st 2nd and 3rd floors - after that the nodes can be placed anywhere.

% eai = cell array where cell j is a matrix where each row lists E,A,I the length (l) and the angle
% fi (measured in degress counterclockwise from the local bar axis) in the bars of frame type j. 

% kc = cell array, each row has the three coefficients of the local flexural stiffness (typically 4 4 2)
% for the bars of frame type j.

% D = matrix that positions each frame in the floor plan. Each row corresponds to one frame and the 
% entries are: x y fi and type - where x and y are points in the line of action of the frame, fi is angle
% of the frame measured counterclockwise with respect to the global x-x axis and type is the type of 
% frame (1 2 etc). 

% L global load vector applied at the origin of the global coordinates.

% *****************************
%           Output
% *****************************

% F = cell array with cell j listing the lateral forces in frame j.
% def = cell array with cell j listing the lateral deformations of frame j.

function [F,delta,K,lm]=capsld(jc,eai,kc,D,L);

% Compute the number of stories (ns), the number of frames (nf), and the number of frame types (nt).
ns=length(L)/3;
[dum,nt]=size(jc);
[nf,dum]=size(D);

% Compute the lateral stiffness of the various frame types
for j=1:nt  
   
   % Define the values associated with the frame in question
   jcl=jc{j};
   lm=lmg1(jcl);
   eail=eai{j};
   kcl=kc{j};

% Compute the stiffness matrix for all the coordinates
[stiff]=sma1(eail,kcl,lm);

% Condense to the lateral DOF
flex=inv(stiff);
flex=flex(1:ns,1:ns);
kl=inv(flex);

% Cell array with the local stiffness matrices
kf{j}=kl;
end;

% Assemble the various frames into the global stiffness

% Initialize the global stiffness
K=zeros(3*ns,3*ns);

for j=1:nf;
Dl=D(j,:);   
% Compute the transformation matrix.
c=cos(Dl(3)*pi/180);
s=sin(Dl(3)*pi/180);
x=Dl(1);
y=Dl(2);
type=Dl(4);
d=[c s x*s-y*c]; 
dz=d*0;
T=zeros(ns,3*ns);

for jj=1:ns;
   ind1=3*(jj-1)+1;
   ind2=ind1+2;
   T(jj,ind1:ind2)=d;
end;
TT{j}=T;

K=K+T'*kf{type}*T;
end;

% Solve for the motion at the global DOF
delta=inv(K)*L;

% Calculate the lateral forces and the lateral deformations in each one of the frames
for j=1:nf;
   type=D(j,4);
   def{j}=TT{j}*delta;
   F{j}=kf{type}*def{j};
end;

% ***** SUBPROGRAMS ********* 

% Subprogram to generate lm
function [lm]=lmg1(jc);
[nb dum]=size(jc);
a=max(jc);
nc=a(2);
lm=[];

for j=1:nb;
   n1=jc(j,1);
   n2=jc(j,2);
   if n1==0;
      e1=[0 0 0];
   e2=[n2 n2+nc n2+2*nc];
else;
   e1=[n1 n1+nc n1+2*nc];
   e2=[n2 n2+nc n2+2*nc];
   end;
   y=[e1 e2];
   lm=[lm;y];
end;

% Subprogram to compute the stiffness of the substuctures
function [stif]=sma1(eai,kc,lm);

[nb,q]=size(lm);
dof=max(lm);
dof=max(dof);
stif=zeros(size(dof));

% Assembly
for i=1:nb
% At the level of distortions
kl=[kc(i,1) kc(i,3);kc(i,3) kc(i,2)];
c=eai(i,1)*eai(i,3)/eai(i,4);
kl=c*kl;
kl=[kl [0 0]'];
comp=[0 0 eai(i,1)*eai(i,2)/eai(i,4)];
kl=[kl;comp];

% Transfer to 6*6 in local coordinates
B=[0 1/eai(i,4) 1 0 -1/eai(i,4) 0;0 1/eai(i,4) 0 0 -1/eai(i,4) 1;-1 0 0 1 0 0];
k1=B'*kl*B;

% Rotate
D=[cos(eai(i,5)*pi/180) sin(eai(i,5)*pi/180) 0; -sin(eai(i,5)*pi/180) cos(eai(i,5)*pi/180) 0; 0 0 1];
dc=zeros(3);
D=[D dc;dc D];
k2=D'*k1*D;

% Assemble the locator matrix
G=zeros(6,dof);
for j=1:6
a=lm(i,j);
if a==0;
G=G;
else;
G(j,a)=1;
end;
end;
% Contribution to global stiffness
kgs=G'*k2*G;
stif=stif+kgs;
end;