% Step 1: Define the geometry of the rectangular beam section

% User-defined parameters
L = input('Enter the length of the beam (m): ');
w = input('Enter the width of the beam (m): ');
h = input('Enter the height of the beam (m): ');

% Select the type of beam system

beamType = input('Enter the type of beam (cantilever or simply supported): ', 's');

% Initialize support locations and types

supportLocations = [];
supportTypes = {};

if strcmpi(beamType, 'cantilever')
% Cantilever beam
supportLocations = [0];
supportTypes = {'cantilever'};
elseif strcmpi(beamType, 'simply supported')
% Simply supported beam
supportLocations = [0, L];
supportTypes = {'pin', 'roller'};
else
error('Unsupported type of beam');
end

% Plot the beam

figure;
hold on;
plot([0 L], [0 0], 'k', 'LineWidth', 2); % Beam
text(L/2, 0.05, 'Beam', 'HorizontalAlignment', 'center');

% Plot supports
for i = 1:length(supportLocations)
switch lower(supportTypes{i})
case 'cantilever'
plot(supportLocations(i), 0, 'r^', 'MarkerSize', 10, 'MarkerFaceColor',
'r'); % Cantilever support
text(supportLocations(i), -0.05, 'Cantilever', 'HorizontalAlignment',
'center');
case 'pin'
plot(supportLocations(i), 0, 'bs', 'MarkerSize', 10, 'MarkerFaceColor',
'b'); % Pin support
text(supportLocations(i), -0.05, 'Pin', 'HorizontalAlignment',
'center');
case 'roller'
plot(supportLocations(i), 0, 'go', 'MarkerSize', 10, 'MarkerFaceColor',
'g'); % Roller support
text(supportLocations(i), -0.05, 'Roller', 'HorizontalAlignment',
'center');
end
end

title('Beam with Supports');

xlabel('Length (m)');
ylabel('Height (m)');
ylim([-0.1 0.1]);
grid on;
hold off;
% Step 2: Determine the type of loading on the beam and plot the loading

% Number of loads
numLoads = input('Enter the number of loads on the beam: ');

% Initialize arrays to store load information

loadTypes = cell(numLoads, 1);
loadMagnitudes = zeros(numLoads, 2); % Stores magnitude(s)
loadLocations = zeros(numLoads, 2); % Stores location(s)

% Get load information from the user

for i = 1:numLoads
loadType = input('Enter the type of load (point load, point torque, uniform
patch load, linear patch load): ', 's');
loadTypes{i} = loadType;
switch lower(loadType)
case 'point load'
loadMagnitudes(i, 1) = input('Enter the magnitude of the point load
(N): ');
loadLocations(i, 1) = input('Enter the location of the point load (m):
');
case 'point torque'
loadMagnitudes(i, 1) = input('Enter the magnitude of the point torque
(N*m): ');
loadLocations(i, 1) = input('Enter the location of the point torque
(m): ');
case 'uniform patch load'
loadMagnitudes(i, 1) = input('Enter the magnitude of the uniform patch
load (N/m): ');
loadLocations(i, 1) = input('Enter the start location of the uniform
patch load (m): ');
loadLocations(i, 2) = input('Enter the stop location of the uniform
patch load (m): ');
case 'linear patch load'
loadMagnitudes(i, 1) = input('Enter the magnitude of the linear patch
load at start (N/m): ');
loadLocations(i, 1) = input('Enter the start location of the linear
patch load (m): ');
loadMagnitudes(i, 2) = input('Enter the magnitude of the linear patch
load at stop (N/m): ');
loadLocations(i, 2) = input('Enter the stop location of the linear
patch load (m): ');
otherwise
error('Unsupported type of load');
end
end

% Plot the loads on the beam

figure;
hold on;
plot([0 L], [0 0], 'k', 'LineWidth', 2); % Beam
text(L/2, 0.05, 'Beam', 'HorizontalAlignment', 'center');

% Plot supports
for i = 1:length(supportLocations)
switch lower(supportTypes{i})
 case 'cantilever'
plot(supportLocations(i), 0, 'r^', 'MarkerSize', 10, 'MarkerFaceColor',
'r'); % Cantilever support
text(supportLocations(i), -0.05, 'Cantilever', 'HorizontalAlignment',
'center');
case 'pin'
plot(supportLocations(i), 0, 'bs', 'MarkerSize', 10, 'MarkerFaceColor',
'b'); % Pin support
text(supportLocations(i), -0.05, 'Pin', 'HorizontalAlignment',
'center');
case 'roller'
plot(supportLocations(i), 0, 'go', 'MarkerSize', 10, 'MarkerFaceColor',
'g'); % Roller support
text(supportLocations(i), -0.05, 'Roller', 'HorizontalAlignment',
'center');
end
end

% Plot loads
for i = 1:numLoads
switch lower(loadTypes{i})
case 'point load'
plot(loadLocations(i, 1), 0.1, 'rv', 'MarkerSize', 10,
'MarkerFaceColor', 'r'); % Point load
text(loadLocations(i, 1), 0.15, sprintf('P = %.2f N', loadMagnitudes(i,
1)), 'HorizontalAlignment', 'center');
case 'point torque'
plot(loadLocations(i, 1), 0.1, 'mo', 'MarkerSize', 10,
'MarkerFaceColor', 'm'); % Point torque
text(loadLocations(i, 1), 0.15, sprintf('T = %.2f N*m',
loadMagnitudes(i, 1)), 'HorizontalAlignment', 'center');
case 'uniform patch load'
plot([loadLocations(i, 1) loadLocations(i, 2)], [0.1 0.1], 'b',
'LineWidth', 2); % Uniform patch load
text(mean(loadLocations(i, :)), 0.15, sprintf('q = %.2f N/m',
loadMagnitudes(i, 1)), 'HorizontalAlignment', 'center');
case 'linear patch load'
plot([loadLocations(i, 1) loadLocations(i, 2)], [0.1 0.1], 'g',
'LineWidth', 2); % Linear patch load
text(loadLocations(i, 1), 0.15, sprintf('q1 = %.2f N/m',
loadMagnitudes(i, 1)), 'HorizontalAlignment', 'center');
text(loadLocations(i, 2), 0.15, sprintf('q2 = %.2f N/m',
loadMagnitudes(i, 2)), 'HorizontalAlignment', 'center');
otherwise
error('Unsupported type of load');
end
end

title('Beam with Supports and Loads');

xlabel('Length (m)');
ylabel('Load Magnitude');
ylim([-0.2 0.2]);
grid on;
hold off;
% Step 3: Compute forces and moments at the supports

% Initialize reaction forces and moments

R = zeros(length(supportLocations), 1); % Reaction forces
M = zeros(length(supportLocations), 1); % Reaction moments

% Sum of vertical forces (Fy)

F_total = 0;
for i = 1:numLoads
if strcmp(loadTypes{i}, 'point load')
F_total = F_total + loadMagnitudes(i, 1);
elseif strcmp(loadTypes{i}, 'uniform patch load')
F_total = F_total + loadMagnitudes(i, 1) * (loadLocations(i, 2) -
loadLocations(i, 1));
elseif strcmp(loadTypes{i}, 'linear patch load')
% Calculate the total load of linear patch load
F_total = F_total + 0.5 * (loadMagnitudes(i, 1) + loadMagnitudes(i, 2)) *
(loadLocations(i, 2) - loadLocations(i, 1));
end
end

% For a simply supported beam

if strcmpi(beamType, 'simply supported')
% Sum of moments about the first support (assuming the first support is a pin)
M_total = 0;
for i = 1:numLoads
if strcmp(loadTypes{i}, 'point load')
M_total = M_total + loadMagnitudes(i, 1) * (loadLocations(i, 1) -
supportLocations(1));
elseif strcmp(loadTypes{i}, 'point torque')
M_total = M_total + loadMagnitudes(i, 1);
elseif strcmp(loadTypes{i}, 'uniform patch load')
q_length = loadLocations(i, 2) - loadLocations(i, 1);
q_center = (loadLocations(i, 1) + loadLocations(i, 2)) / 2;
M_total = M_total + loadMagnitudes(i, 1) * q_length * (q_center -
supportLocations(1));
elseif strcmp(loadTypes{i}, 'linear patch load')
q_length = loadLocations(i, 2) - loadLocations(i, 1);
q_center = (loadLocations(i, 1) + loadLocations(i, 2)) / 2;
M_total = M_total + (1/6) * (2*loadMagnitudes(i, 1) + loadMagnitudes(i,
2)) * q_length * (loadLocations(i, 2) - supportLocations(1));
end
end
% Solve for the reaction forces
R(2) = M_total / (supportLocations(2) - supportLocations(1)); % Reaction force
at the roller support
R(1) = F_total - R(2); % Reaction force at the pin support

elseif strcmpi(beamType, 'cantilever')

% Cantilever beam has a fixed support at the start
% Calculate the reaction forces and moment
Fy = 0; % Sum of vertical forces
M_total = 0; % Sum of moments
for i = 1:numLoads
if strcmp(loadTypes{i}, 'point load')
Fy = Fy + loadMagnitudes(i, 1);
M_total = M_total + loadMagnitudes(i, 1) * (loadLocations(i, 1) -
supportLocations(1));
elseif strcmp(loadTypes{i}, 'point torque')
 M_total = M_total + loadMagnitudes(i, 1);
elseif strcmp(loadTypes{i}, 'uniform patch load')
q_length = loadLocations(i, 2) - loadLocations(i, 1);
q_center = (loadLocations(i, 1) + loadLocations(i, 2)) / 2;
Fy = Fy + loadMagnitudes(i, 1) * q_length;
M_total = M_total + loadMagnitudes(i, 1) * q_length * (q_center -
supportLocations(1));
elseif strcmp(loadTypes{i}, 'linear patch load')
q_length = loadLocations(i, 2) - loadLocations(i, 1);
Fy = Fy + 0.5 * (loadMagnitudes(i, 1) + loadMagnitudes(i, 2)) *
q_length;
M_total = M_total + (1/6) * (2*loadMagnitudes(i, 1) + loadMagnitudes(i,
2)) * q_length * (loadLocations(i, 2) - supportLocations(1));
end
end
R(1) = Fy; % Vertical reaction force at the fixed support
M(1) = M_total; % Moment at the fixed support
end

% Display the reaction forces and moments

disp('Reaction Forces and Moments at the Supports:');
for i = 1:length(supportLocations)
if strcmpi(beamType, 'cantilever') && i == 1
fprintf('Fixed support at x = %.2f m: Fy = %.2f N, M = %.2f N*m\n',
supportLocations(i), R(i), M(i));
elseif strcmpi(beamType, 'simply supported')
if strcmp(lower(supportTypes{i}), 'pin')
fprintf('Pin support at x = %.2f m: Fy = %.2f N\n',
supportLocations(i), R(i));
elseif strcmp(lower(supportTypes{i}), 'roller')
fprintf('Roller support at x = %.2f m: Fy = %.2f N\n',
supportLocations(i), R(i));
end
end
end

% Step 4: Compute Bending Moment, Shear Force, and Torsional Moment

% Number of points to evaluate along the beam

numPoints = 100;
x = linspace(0, L, numPoints); % Discrete points along the beam

% Initialize arrays to store the results

V = zeros(1, numPoints); % Shear force
M_b = zeros(1, numPoints); % Bending moment
T = zeros(1, numPoints); % Torsional moment

% Compute the shear force, bending moment, and torsional moment at each point
for j = 1:numPoints
x_j = x(j);
 % Initialize shear force and bending moment at x_j
V_j = 0;
M_j = 0;
T_j = 0;

% Consider contributions from point loads

for i = 1:numLoads
if strcmp(loadTypes{i}, 'point load')
if x_j >= loadLocations(i, 1)
V_j = V_j + loadMagnitudes(i, 1);
M_j = M_j + loadMagnitudes(i, 1) * (x_j - loadLocations(i, 1));
end
elseif strcmp(loadTypes{i}, 'point torque')
if x_j >= loadLocations(i, 1)
T_j = T_j + loadMagnitudes(i, 1);
end
elseif strcmp(loadTypes{i}, 'uniform patch load')
if x_j >= loadLocations(i, 1) && x_j <= loadLocations(i, 2)
q_length = x_j - loadLocations(i, 1);
V_j = V_j + loadMagnitudes(i, 1) * q_length;
M_j = M_j + loadMagnitudes(i, 1) * q_length * (q_length / 2);
elseif x_j > loadLocations(i, 2)
q_length = loadLocations(i, 2) - loadLocations(i, 1);
V_j = V_j + loadMagnitudes(i, 1) * q_length;
M_j = M_j + loadMagnitudes(i, 1) * q_length * (q_length / 2 + (x_j
- loadLocations(i, 2)));
end
elseif strcmp(loadTypes{i}, 'linear patch load')
if x_j >= loadLocations(i, 1) && x_j <= loadLocations(i, 2)
q1 = loadMagnitudes(i, 1);
q2 = loadMagnitudes(i, 2);
x1 = loadLocations(i, 1);
x2 = loadLocations(i, 2);
a = (q2 - q1) / (x2 - x1);
q_x = q1 + a * (x_j - x1);
V_j = V_j + (q1 * (x_j - x1) + 0.5 * a * (x_j - x1)^2);
M_j = M_j + (q1 * (x_j - x1) * (x_j - x1) / 2 + a * (x_j - x1)^3 /
6);
elseif x_j > loadLocations(i, 2)
q1 = loadMagnitudes(i, 1);
q2 = loadMagnitudes(i, 2);
x1 = loadLocations(i, 1);
x2 = loadLocations(i, 2);
a = (q2 - q1) / (x2 - x1);
V_j = V_j + (q1 * (x2 - x1) + 0.5 * a * (x2 - x1)^2);
M_j = M_j + (q1 * (x2 - x1)^2 / 2 + a * (x2 - x1)^3 / 6 + (q1 * (x2
- x1) + 0.5 * a * (x2 - x1)^2) * (x_j - x2));
end
end
end

% Consider contributions from reaction forces

if strcmpi(beamType, 'simply supported')
for k = 1:length(supportLocations)
if x_j >= supportLocations(k)
V_j = V_j - R(k);
M_j = M_j - R(k) * (x_j - supportLocations(k));
end
end
 elseif strcmpi(beamType, 'cantilever')
V_j = V_j - R(1);
M_j = M_j - M(1) - R(1) * (x_j - supportLocations(1));
end

% Store the computed values

V(j) = V_j;
M_b(j) = M_j;
T(j) = T_j;
end

% Plot the results

figure;

subplot(3, 1, 1);
plot(x, V, 'r', 'LineWidth', 2);
title('Shear Force Diagram');
xlabel('Length (m)');
ylabel('Shear Force (N)');
grid on;

subplot(3, 1, 2);
plot(x, M_b, 'b', 'LineWidth', 2);
title('Bending Moment Diagram');
xlabel('Length (m)');
ylabel('Bending Moment (N*m)');
grid on;

subplot(3, 1, 3);
plot(x, T, 'g', 'LineWidth', 2);
title('Torsional Moment Diagram');
xlabel('Length (m)');
ylabel('Torsional Moment (N*m)');
grid on;

% Step 5: Compute Maximum Stresses and Plot Them

% Initialize arrays to store maximum stresses

sigma_b_max = zeros(1, numPoints); % Maximum bending stress
tau_s_max = zeros(1, numPoints); % Maximum shear stress
tau_t_max = zeros(1, numPoints); % Maximum torsional stress

% Compute the maximum stresses at each point along the beam

for j = 1:numPoints
% Compute bending stress
sigma_b_max(j) = abs(M_b(j) * h / (w * h^3 / 12));

% Compute shear stress

tau_s_max(j) = abs(V(j) * h^2 / (2 * w * h^2 / 4));

% Compute torsional stress

r = h / 2; % Distance from the center to the outermost fiber
tau_t_max(j) = abs(T(j) * r / (w * h^3 / 3));
end
% Plot the maximum stresses along the beam
figure;

subplot(3, 1, 1);
plot(x, sigma_b_max, 'r', 'LineWidth', 2);
title('Maximum Bending Stress along the Beam');
xlabel('Length (m)');
ylabel('Bending Stress (Pa)');
grid on;

subplot(3, 1, 2);
plot(x, tau_s_max, 'b', 'LineWidth', 2);
title('Maximum Shear Stress along the Beam');
xlabel('Length (m)');
ylabel('Shear Stress (Pa)');
grid on;

subplot(3, 1, 3);
plot(x, tau_t_max, 'g', 'LineWidth', 2);
title('Maximum Torsional Stress along the Beam');
xlabel('Length (m)');
ylabel('Torsional Stress (Pa)');
grid on;