function [K] = stiff( bar, node, area, young, free)
%
% input variables
%
% bar    integer(m,2)
%        The i-th row contains information for bar i in the following 
%        form:
%                  [ lower-left-node   upper-right-node ]
%
% node   real(n,2) 
%        The i-th row contains the coordinates of node i in the 
%        following form:
%                  [x-coordinate-of-node   y-coordinate-of-node]
%
% area   real(m)
%        area(i) is the cross sectional area of bar i.
%
% young  real(m)
%        young(i) is the Young's modulus of bar i.
%  
% free   integer(k)
%        index vector that contains the indices of displacement
%        that are not fixed. (k <= n)
%
%
% output variables
%
% K      real(k,k)
%        stiffness matrix
%


% get number of nodes n  and number of bars m
n = size(node,1);
m = size(bar,1);

% compute angles of the bar
theta = zeros(m,1);
for i = 1:m
    if( node(bar(i,2),1) - node(bar(i,1),1) == 0 )
        theta(i) = -pi/2;
    else
        theta(i) = atan( ( node(bar(i,2),2) - node(bar(i,1),2) ) / ...
                         ( node(bar(i,2),1) - node(bar(i,1),1) ) );
        if( theta(i) >= 0 )
            theta(i) = - theta(i);
        else
            theta(i) = -(pi+theta(i));
        end
    end
end

% compute lengths of the undeformed bars
length = zeros(m,1);
for i = 1:m
    length(i) = sqrt( (node(bar(i,1),1)-node(bar(i,2),1))^2 ...
                      + (node(bar(i,1),2)-node(bar(i,2),2))^2 );
end

D = diag(area.*young./length);  
        
B = zeros(2*n,m);
for i = 1:m
    % Populate row of B-TRANSPOSE, i.e., columns i of B.
    B(2*bar(i,1)-1,i) = -cos(theta(i));
    B(2*bar(i,1)  ,i) = -sin(theta(i));
    B(2*bar(i,2)-1,i) =  cos(theta(i));
    B(2*bar(i,2)  ,i) =  sin(theta(i));
end

% determine system matrix
K = B(free,:)*D*B(free,:)';

%end of stiff