function [] = plotbeamelem(x1,y1,x2,y2,u,N)
% Plots a beam element's deformation using interpolation within it
ux1 = u(1); 
uy1 = u(2);
ux2 = u(4);
uy2 = u(5);
xp1 = x1 + ux1;
yp1 = y1 + uy1;
xp2 = x2 + ux2;
yp2 = y2 + uy2;
phi0 = atan2((y2-y1), (x2-x1));
phi = atan2((yp2-yp1), (xp2-xp1));
theta1 = u(3) + phi0 - phi;
theta2 = u(6) + phi0 - phi; 
l = sqrt( (xp2-xp1)^2 + (yp2-yp1)^2);
A = [0 0 0 1; 0 0 1 0; l^3 l^2 l 1; 3*l^2 2*l 1 0];
a = inv(A) * [0; tan(theta1); 0; tan(theta2)];
x0 = xp1;
y0 = yp1;
plot([x1 x2], [y1 y2], 'k.');
hold on
plot([x1 x2], [y1 y2], 'k--');
plot([xp1 xp2], [yp1 yp2], 'r.','LineWidth',2);
for x = l/N:l/N:l,
    y = a(1)*x^3 + a(2)*x^2 + a(3)*x + a(4);
    x1 = xp1 + cos(phi)*x - sin(phi)*y;
    y1 = yp1 + sin(phi)*x + cos(phi)*y;
    plot([x0 x1], [y0 y1], '-r','LineWidth',2);
    x0 = x1;
    y0 = y1;
end
