bending1(s, d, F1, F2)

function M = bending1(s, d, F1, F2)
%BENDING1 bending moments caused vehicle moving on beam.
%       See book, Section 17.1.
%    Bending moments caused by vehicle moving on
%    simply-supported beam.
%    The  input arguments and their SI units are:
%    s  - beam span, m
%    d  - car      , m
%    F1 - forward-axle load, N
%    F2 - rear axle load, N
%    The output is
%    M, m-by-n matrix of bending moments, Nm
%    Internal variables are:
%    x - position of forward axle relative to left-hand
%    support, positive forward, m
%    y - coordinate where M is calculated, measured from
%        the left-hand support, positive forward, m
%    R1 - reaction in left-hand support, N
%    R2 - reaction in right-hand support, N

s, d, F1, F2           % display the input
disp('Press ENTER to continue')
pause
step = d/10;
m = fix(s/step)            % number of y steps
n = fix((s + d)/step)      % number of x steps
x = [0: n]*step;           % equally spaced points
x = [ x (s+d) ]            % include last point
y = [0: m]*step;           % equally spaced points
y = [ y s ]                % include last point of span
M = zeros((m+2), (n+2));   % allocate space for array of moments

for k = 1:(n+1)
        % case 1           % shows where we are
        if x(k) < d
              disp('case 1')
              R1 = (1 - x(k)/s)*F1;
              for l = 1:(m+1)
                   M(l,k) = y(l)*R1;
                   if y(l) >= x(k)
                        M(l,k) = M(l,k) - (y(l) - x(k))*F1;
                   end
              end
        % case 2
        elseif x(k) < s
              disp('case 2')        % shows where we are
              R1 = (1 - x(k)/s)*F1 + (1 + (d - x(k))/s)*F2;
              for l = 1:(m+1)
                   M(l, k) = y(l)*R1;
                   if y(l) >= x(k) -d
                        M(l,k) = M(l,k) - (y(l) - (x(k)-d))*F2;
                   end
                   if y(l) >= x(k)
                        M(l,k) = M(l,k) - (y(l) - x(k))*F1;
                   end
              end
     % case 3
     else
              disp('case 3')          % shows where we are
              R1 = (1 + (d - x(k))/s)*F2;
              for l = 1:(m+1)
                        M(l,k) = y(l)*R1;
                        if y(l) > x(k) - d
                             Mdiff  = (y(l) - x(k) + d)*F2;
                             M(l,k) = M(l,k) - Mdiff;
                        end
              end
     end
end
Mmax = max(M)          % displays maximum for
                       % each vehicle position
Mmaxmax = max(Mmax)    % maximum maximorum moment

% 3-D MATLAB4.0 plot
surf(x, y, M)
xlabel('Position of forward axle, m')
ylabel('Position of section, m')
zlabel('Bending moment, Nm')
pause
print bending1 -dps
contour3(x, y, M)
xlabel('Position of forward axle, m')
ylabel('Position of section, m')
zlabel('Bending moment, Nm')
print bending3 -dps