 function [ D, nodes, links] = truss_load(N,a,b,theta)
% Construct matrix for internal loads on o triangular truss structure
%
% D      load matrix (D * loads + V_ext = 0)
% nodes(i).x   coordinates of ith node
% links.nodes(k) = [i j] end nodes for link k
% links.prop(k)  property of link k (batten, etc)
%
% N      number of bays
% theta  interior angle in y-z triangles
%  a     length of longerons
%  b     length of battens
%
if nargin < 4; theta = pi/3;
    if nargin < 3; b = 1;
        if nargin < 2; a = 1;
            if nargin < 1; N = 2;
            end
        end
    end
end
YZ = b*[0 0; 1 0; cos(theta) sin(theta)]; % basic triangle
%
nodes = struct([]);  % initialize empty data structure
in = 0; xa = -a;
for ii=1:N+1
    xa = xa + a; 
    for kk=1:3
        in = in+1;
        nodes(in).x = [xa  YZ(kk,:)];
    end
end
%
links = assem_links(N); % build the links data structure
%
% Each bay adds 3 nodes and 9 links
% There are 9*N unknowns - the tension in each link
% There are 3*(3*N) equations - force equilibrium at each node
%
D = zeros(9*N,9*N);
% rows 1:3 characterize (x,y,z) equilib for the 1st unknown node (4)
id = [-2:1:0]';
for ii=4:3*(N+1) % loop through the nodes
    id = id + 3;
    for jj=1:9*N % loop through the links
            D(id,jj) = dir_cosine(ii,links(jj).nodes,nodes);
    end % jj
end % ii
%
%
% internal sub-functions
function links = assem_links(N)
% Build the  links data structure
%
% Nodes 1 -3 are built in (anchor points)
% First build an adjaceny matrix
adj = reshape(blanks((3*(N+1))^2),[3*(N+1),3*(N+1)]);
B = ['   ';'b  ';'bb '];
M = ['ld ';' l ';'ddl'];
ia = [1:3];
adj(ia,ia) = B;
%
for ii=1:N
    ia = ia + 3;
    adj(ia,ia) = B;
    adj(ia,ia-3) = M;
    M = M';
end
% Now scan the sub-diagonals of rows 4 --> end
%
links = struct([]);
il = 0;
for ii=4:size(adj,1)
    for jj=1:ii-1
        if ~strncmp(adj(ii,jj),' ',1)
            il = il + 1;
            links(il).nodes = [jj ii];
            links(il).prop = adj(ii,jj);
        end
    end
end
function alpha = dir_cosine(ii,ends,nodes) 
% Compute the 3 \times 1 direction cosine vector
% ii is the node at which we sum forces
% ends is an array with the end nodes of a link
%
if ends(1) == ii
   alpha = nodes(ii).x - nodes(ends(2)).x;
   alpha = alpha(:) / norm(alpha);
elseif ends(2) == ii
   alpha = nodes(ii).x - nodes(ends(1)).x;
   alpha = alpha(:) / norm(alpha);
else
    alpha = zeros(3,1);
end
