function K = beam2d(A,E,Iz,P1,P2)
%   BEAM2D
%
%   This command generates the 6 x 6 stiffness matrix for a two
%   dimensional beam element in global coordinates. Slender beams
%   are assumed such that shear deformations can be neglected.
%   The syntax is:
%              K = beam2d(A,E,Iz,P1,P2)
%   where: A is the cross-sectional area;
%          E is the Young's modulus;
%          Iz is the cross-section moment of inertia about the
%              z axis;
%          P1 and P2 are vectors of the {x,y} coordinates of the
%              endpoints.

L = norm(P2-P1);               % Find length of beam.

lox = (P2(1)-P1(1))/L;         % Find direction cosines.
mox = (P2(2)-P1(2))/L;
loy = -mox;
moy = lox;
lam = [lox mox 0; loy moy 0; 0 0 1];
zerom = zeros(3);
lambda = [lam zerom; zerom lam];

kl = [zerom zerom; zerom zerom];     % Create local matrix.
kl(1,1) = A*(L^2)/Iz;
kl(2,2) = 12;
kl(3,2) = 6*L;
kl(3,3) = 4*(L^2);
kl(4,1) = -kl(1,1);
kl(4,4) = kl(1,1);
kl(5,2) = -12;
kl(5,3) = -kl(3,2);
kl(5,5) = 12;
kl(6,2) = kl(3,2);
kl(6,3) = 2*(L^2);
kl(6,5) = -kl(6,2);
kl(6,6) = kl(3,3);
for i = 1:5
   for j = i+1:6
      kl(i,j) = kl(j,i);
   end
end

K = E*Iz/(L^3)*lambda'*kl*lambda;      % Create global matrix.
