% This Matlab code can be used for a simply supported beam with uniformly distributed load to find the
% * Support reactions 
% * Shear force and the bending moment at a specific point
% * Maximum Bending Moment 
% * Shear force diagram
% * Bending Moment diagram
% Also, It verifies that the derivative of M is equal to the shear force V
% and the value of M at a specific point through the polynomial is the same with the previous one.
% Be careful measurements units must be kN,m.  
clc; clear; close all;
disp('Simply Supported Beam'); 


% Data input section 
disp(' ');
L  = input('Length of the beam in meters= ');
disp(' ');
q = input('Uniformly distributed load in kN/m = ');
disp(' ');
x0 = input('give the specific point that you want to find the intensity measures x0='); 
disp(' ');
while x0<0 ||(L<2||L>5)||q>20
    disp('error inputs')
    disp(' ');
    L  = input('Length of the beam in meters= ');
    disp(' ');
    q = input('Uniformly distributed load in kN/m = ');
    disp(' ');
    x0 = input('give the specific point that you want to find the intensity measures x0='); 
    disp(' ');
end
R1 = q*L/2; % Left Support Reaction. 
R2 = -q*L/2; % Right Support Reaction.

% Discretization of x axis.
n = 1000;           % Number of discretization of x axis.
dx= L/n;            % Increment for discretization of x axis.
x = 0:dx:L;         % Generate column array for x-axis.  

V = zeros(size(x, 1), 1);   % Shear force function of x. 
M = zeros(size(x, 1), 1);   % Bending moment function of x. 


% Data processing section;
number1=x(1);
Mmax=M(1);
for i = 1:n+1
   V(i) = R1-q*x(i);
   M(i) = R1*x(i)-q*x(i)^2/2;
   if Mmax<abs(M(i))
       Mmax=abs(M(i));
       number1=i;
   end
   if x(i)==x0
       V_value=V(i);
       M_value=M(i);  
   end
end
disp('Results:')
disp(' ');disp (['Left support Reaction'  '  = ' num2str(R1) ' ' 'kN'])
disp (['Right support Reaction'  '  = ' num2str(R2) ' ' 'kN'])

disp(' ');disp (['Maximum bending moment'  ' = ' num2str(Mmax) ' ' 'kNm' ' at the point ', num2str(x(number1)),'m '])

disp(' '); disp(['the value of V at the x=',num2str(x0),'m is ',num2str(V_value),' kN '])
disp(['the value of M at the x=',num2str(x0),'m is ', num2str(M_value),' kN '])


figure
subplot(2,1,1);
plot(x, V, 'r','linewidth',1.5);
refline(0,0)
% grid
% line([x(1) x(end)],[0 0],'Color','k');
% line([0 0],[0 V(1)],'Color','r','linewidth',1.5);
% line([x(end) x(end)],[0 V(end)],'Color','r','linewidth',1.5);
% title('Shear Force Diagram','fontsize',16)
% text(0,V(1)+V(1)/6,num2str(V(1)),'HorizontalAlignment','center','FontWeight','bold','fontsize',16)
% text(L,V(end)+V(end)/6,num2str(V(end)),'HorizontalAlignment','center','FontWeight','bold','fontsize',16)
% axis off

subplot(2,1,2);
plot(x, M, 'r','linewidth',1.5);
% grid
% line([x(1) x(end)],[0 0],'Color','k');
% line([x(number1) x(number1)],[0 Mmax],'LineStyle','--','Color','b');
% title('Bending Moment Diagram','fontsize',16);
% text(x(number1),Mmax+Mmax/10,num2str(roundn(Mmax,-2)),'HorizontalAlignment','center','FontWeight','bold','fontsize',16);
% text(x(number1),-Mmax/10,[num2str(roundn(x(number1),-2)) ' m'],'HorizontalAlignment','center','FontWeight','bold','fontsize',16);
% axis off


%Verification
disp(' ');disp('Verification:')
 V(i) = R1-q*x(i);
   M(i) = R1*x(i)-q*x(i)^2/2;
V = [-q R1];
M = [-q/2 R1 0];
diffM=polyder(M);
ver=0;
if diffM==V
    disp(' ');disp('Indeed, the derivate of M is equal to V')
else
    disp(' ');disp('The derivative of M is not equal to V. Error input')
end
Mvalue=polyval(M,x0);
if Mvalue==M_value
    disp(['Indeed, the value of M at x=',num2str(x0),'m is the same with the previous one']);disp(' ')
else
    disp('The values are not the same. Error input'); disp(' ')
end