function [M,R,M_Diag,V_Diag,d_Diag,D_max] = three_moment(L,I,w,E)
% [M,R,M_Diag,V_Diag,d_Diag,D_max] = three_moment(L,I,w,E)
% solve the three-moment equations for a continuous beam of N spans.
%
% input: 
%   L is a vector of length N containing the lengths of each span.
%   I is a vector of length N containing the moments of inertia of each span.
%   w is a vector of length N containing the uniform loads on each span.
%   E is a scalar constant for the modulus of elasticity.  
%
% output:
%   M is a vector of lenght N+1 containing the moments at each support
%   R is a vector of lenght N+1 containing the reactions at each support
%   M_Diag is a vector of the moment diagram
%   V_Diag is a vector of the shear diagram
%   d_Diag is a vector of the displacement diagram
%   D_max  is a vector of the max absolute displacement of each span
%
% This program assumes that none of the supports are moment resisting, 
% that there are no hinges in the beam, that all the spans are made of 
% the same material, and that each spans is prismatic.  

% H.P. Gavin, Civil and Environmental Engineering, Duke University, 1/24/00

 N = length(L);				% number  of spans

 F = zeros(N-1,N-1);			% initialize the matrix to be zero

 for j=1:N-1				% create the flexibility matrix (8)-(10)
	if ( j > 1 ) 
		F(j,j-1) = L(j) / I(j);
	end

	F(j,j) = 2 * (  L(j) / I(j) + L(j+1) / I(j+1)  );

	if ( j < N-1 )
		F(j,j+1) = L(j+1) / I(j+1);
	end
 end
	
 for j=1:N-1				% create the right-hand-side vector (11)

	d(j) = - w(j)*L(j)^3 / ( 4 * I(j) ) - w(j+1)*L(j+1)^3 / (4 * I(j+1) );
 end

 m = inv(F) * d(:);			% compute the internal moments  (7)

 for j=1:N-1				% compute the reactions  (12)

	r(j) = w(j)*L(j)/2 + w(j+1)*L(j+1)/2 - m(j)/L(j) - m(j)/L(j+1);

	if ( j > 1 )
		r(j) = r(j) + m(j-1)/L(j);
	end

	if ( j < N-1 )
		r(j) = r(j) + m(j+1)/L(j+1);
	end
 end

 R = [	w(1)*L(1)/2 + m(1)/L(1)   r(:)'   w(N)*L(N)/2 + m(N-1)/L(N) ];

 M = [ 0  m' 0 ];			% end moments are zero.

 for j=1:N				% compute the slopes  (15)
	slope(j) = -( w(j)*L(j)^3 / 24 + ... 
			M(j+1)*L(j) / 6 + M(j)*L(j) / 3 ) / (E*I(j));
 end

 if (  abs ( sum(R) - sum (w .* L) ) < 1e-9  )
	disp (' yes! ')		% equilibrium check ... should be close to zero
 end

% ----------  shear, moment, slope, and deflection data and plots  -----------

 for j=1:N	% x-axis data for shear, moment, slope, and deflection diagrams
	x(:,j) = [ 0:L(j)/55:L(j) ]' ;
 end

 [row,col] = size(x);		

 for j=1:N
	Vo = ( M(j) - M(j+1) ) / L(j) - w(j)*L(j)/2;		% shear at left
	V_Diag(:,j) = Vo + w(j)*x(:,j);
	M_Diag(:,j) = M(j) - Vo*x(:,j) - (1/2)*w(j)*x(:,j).^2 ;
	s_Diag(:,j) = cumtrapz( M_Diag(:,j) ) * x(2,j) / (E*I(j)) + slope(j) ;
	d_Diag(:,j) = cumtrapz( s_Diag(:,j) ) * x(2,j) ;
 end

% ------------- print key results to screen --------------
 fprintf('------------------------------------------------------------------\n');
 fprintf('                 Moment              Shear             Deflection\n');
 fprintf('   Maximum      %12.5e        %12.5e      %12.5e\n', ...
             max(max(M_Diag)), max(max(-V_Diag)), max(max(d_Diag)) );
 fprintf('   Minimum      %12.5e        %12.5e      %12.5e\n', ...
             min(min(M_Diag)), min(min(-V_Diag)), min(min(d_Diag)) );
 fprintf('------------------------------------------------------------------\n');
 
 for j=1:N		% x-axis data for shear and moment diagram plots
	x(:,j) = x(:,j) + sum( L(1:j-1) ) ;
 end

 V_Diag(row+1,:) = zeros(1,col);	% close the polygon for shaded plots 
 M_Diag(row+1,:) = zeros(1,col);
 V_Diag(row+2,:) = zeros(1,col);
 M_Diag(row+2,:) = zeros(1,col);
 x(row+1,:)      = x(row,:);
 x(row+2,:)      = x(1,:);

% Plotting

figure(1);
 subplot(4,1,1)	 	
% fill ( x , -V_Diag , -V_Diag )		% shaded plot
% grid on
  plot ( x , -V_Diag, '-b' )			% line plot
  ylabel('Internal Shear (lbf)')

 subplot(4,1,2)  		
% fill ( x , M_Diag , M_Diag )			% shaded plot
% grid on
  plot ( x , M_Diag , '-b' )			% line plot
  ylabel('Internal Moment (lbf.in)')

 subplot(4,1,3)	 
  plot ( x(1:row,:) , s_Diag , '-b' )		% line plot
% plot ( x(1:row,:) , s_Diag , '-b', 'LineWidth', 2 )		% line plot
% grid on
  ylabel('Slope')

 subplot(4,1,4)  	
% plot ( x(1:row,:) , d_Diag , '-b', 'LineWidth', 2 )		% line plot
% grid on
  plot ( x(1:row,:) , d_Diag , '-b' )		% line plot
  xlabel('span (in)')
  ylabel('Deflection (in)')

 D_max  = max(abs(d_Diag)) 
 M_Diag = M_Diag(:);
 V_Diag = V_Diag(:);
 d_Diag = d_Diag(:);

% ------------------------------------------------------------- THREE_MOMENT