function [def] = simpbeam(P,a,L,E,I,incr)
% [def] = SIMPBEAM(P,a,L,E,I,incr) returns an array with
% the deflection [def].  [def] is also plotted.
% 
% ALL CALCULATIONS ARE FOR A SIMPLE BEAM SUPPORTED AT EITHER END
%
%	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 arrary 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

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

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


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

	def = def + tempdef;
	%slp = slp + tempslp;
end %for j


x = [0:stepp:L];
% theta = atan(slp);
figure(1);
clf;
plot(x,def.*1000);
hold on;
grid on;
title('SIMPLE BEAM 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


% figure(2);
% clf;

% plot(x,theta.*(180/pi));
% hold on;
% grid on;
% title('SIMPLE BEAM slope vs position');
% xlabel('position [meters]');
% ylabel('angle from horizontal [degrees]');

% q=0.2;
% stdlength = (max(theta.*(180/pi))-min(theta.*(180/pi)))*q;
% for i = 1:size(a,2)
% 	k = find(x<a(i));
% 	g=max(k);
% 	if P(i)>0
% 		h = line([a(i) a(i)],[theta(g).*(180/pi)+stdlength/2 theta(g).*(180/pi)-stdlength/2]);
% 	else
% 		h = line([a(i) a(i)],[theta(g).*(180/pi)-stdlength/2 theta(g).*(180/pi)+stdlength/2]);
% 	end
% 	set(h,'color','red');
% end %for