function [def] = cantbeam(P,a,L,E,I,incr)
% [def] = CANTBEAM(P,a,L,E,I,incr) returns arrays with
% the deflection [def].  [def] is also plotted.
% 
% ALL CALCULATIONS ARE FOR A CANTILEVERED BEAM
%
%	P is the force applied to the beam in NEWTONS
%	a is the distance from the left end that the force is applied in METERS
%	L is the length of the beam in METERS
%	E is the young's modulus of the material in Pa
%	I is the cross sectional inertia in METER^4
%	icnr is the number of increments to sample along the beam
%
% Multiple forces can be entered in P, however, a must be the same length
% to give a position for each force.  P can be positive or negative.
%
% all values of a may not be greater then L
% E, I, a, incr, and L must be greater then 0
%
% def is in meters (plotted in mm)
%
%By Roger Cortesi (Course 2 Class of 1999) for 2.007

if size(a)==size(P)

else
	error('P and a must be the same length')
end

if (E < 0 | E ==0)
	error('E must be greater then 0');
end

if (I < 0 | I ==0)
	error('I must be greater then 0');
end

if (incr < 0 | incr ==0)
	error('incr must be greater then 0');
end
if (L < 0 | L ==0)
	error('L must be greater then 0');
end

for j = 1:size(a,2)
	if a(j) > L
		error('a must be less then or equal then L');
	end
	if a(j) < 0
		error('a may not be negative');
	end
	
end %for

stepp = L/incr;
x = [0:stepp:L];
def=zeros(size(x));


for j = 1:size(a,2)
	
	c1 = P(j)/(6*E*I);
	c2 = P(j)/(2*E*I);
	c3 = P(j)/(6*E*I);

	i=1;
	for x= 0:stepp:L
		if x < a(j)
			tempdef(i) = c1*(x^2)*(3*a(j)-x);
		elseif (x==a(j) | x>a(j))
			tempdef(i) =c3*(a(j)^2)*(3*x-a(j));
		end %if
			i=i+1;
	end %for i

	def = def + tempdef;
end %for j


x = [0:stepp:L];

figure(1);
clf;
plot(x,def.*1000);
hold on;
grid on;
title('CANTILEVER deflection vs position');
xlabel('position [meters] Forces shown at location, in Newtons');
ylabel('deflection [mm]');

q=0.2;
stdlength = (max(def)-min(def))*1000*q;

for i = 1:size(a,2)
	k = find(x<a(i));
	g = max(k);
	if P(i)>0
		ty=def(g)*1000-stdlength;
	else
		ty=def(g)*1000+stdlength;
	end
	
	h = line([a(i) a(i)],[def(g)*1000 ty]);
	text(a(i),ty,num2str(P(i),5));
	set(h,'color','red');
end %for