super.m

% Example: super 
% ~~~~~~~~~~~~~~
% This example calculates the shear, moment,
% rotation, and deflection for a simply 
% supported beam using an approximate technique
% based on superposition.  Various loads are 
% replaced with a system of concentrated loads 
% which approximate the original loading 
% configuation.  The results from these 
% calculations are compared with the true
% results.  The beam equations used in this
% example were developed using discontinuity 
% functions.
%
% All the input to this program is defined
% in function setup.
%
% User m functions required: 
%    setup, truebeam, genprint
%----------------------------------------------

clear; 
Problem=1;      % Multi-load example
%Problem=2;      % Uniformly distributed load 
%Problem=3;      % Linearly distributed load 

%...Parameter definitions for specific problem
[L,E,I,no_defl_pts,start_x,end_x, ...
  start_val,end_val,no_sub_loads,no_loads]= ...
  setup(Problem);

%...Initialize
right_react=0; left_react=0;
EI=E*I; dx=L/(no_defl_pts-1); x=[0:dx:L];
V=zeros(1,no_defl_pts); M=zeros(1,no_defl_pts);
EItheta=zeros(1,no_defl_pts);
EIdelta=zeros(1,no_defl_pts);

%...Loop on each loading
%.......................
for i=1:no_loads
  %...x spacing and slope of substitute loads
  if no_sub_loads(i) == 1
    %...not a distributed load
    dx=1; slope=0; loop_limit=1; shift_x=0;
  else
    dx=(end_x(i)-start_x(i))/no_sub_loads(i);
    slope=(end_val(i)-start_val(i))/ ...
          (end_x(i)-start_x(i));
    loop_limit=no_sub_loads(i); shift_x=dx/2;
  end
  %...Instantaneous load value
  if start_val(i) == end_val(i)
    %...Point loads & uniformly distributed
    load_val=start_val(i)*dx;
  end

  %...Substitute load loop
  %.......................
  for j = 1:loop_limit
    % fprintf('.');    % Let you know work is 
                       %   being done
    %...Instantaneous location
    load_x=start_x(i)+dx*(j-1)+shift_x;
    if start_val(i) ~= end_val(i)
      %...linearly distributed load
      load_val=0.5*slope*dx^2*(2*j-1);
    end
    %...Superpose reaction contributions
    left_react=left_react-load_val*(L-load_x)/L;
    right_react=right_react-load_val*load_x/L;

    %...Loop on each x along entire beam
    %...................................
    bl6=L*6; b=L-load_x; Cv=-load_val*b/L;
    Ctheta=(load_val*b/6)*(L-b^2/L);
    V=V+Cv; M=M+Cv*x;
    EItheta=EItheta+Cv/2*x.^2+Ctheta;
    EIdelta=EIdelta+Cv/6*x.^3+Ctheta*x;
    for k = 1:no_defl_pts
      % fprintf('.');
      if (x(k)-load_x) > 0 
        xd=x(k)-load_x; V(k)=V(k)+load_val;
        M(k)=M(k)+load_val*xd;
        EItheta(k)=EItheta(k)+load_val/2*xd^2;
        EIdelta(k)=EIdelta(k)+load_val/6*xd^3;
      end
    end
  end
end

%...Calculate the true values
[V_true,M_true,EItheta_true,EIdelta_true, ...
  left_react_true,right_react_true]= ...
  truebeam(Problem,L,x);

%...Calculate the difference between the
%...true and approximate solutions
V_diff=V_true-V; M_diff=M_true-M;
theta_diff=EItheta_true-EItheta;
delta_diff=EIdelta_true-EIdelta;
V_diff=abs(V_diff); M_diff=abs(M_diff);
theta_diff=abs(theta_diff);
delta_diff=abs(delta_diff);

%...Divide by EI
theta=EItheta/EI; theta_true=EItheta_true/EI;
delta=EIdelta/EI; delta_true=EIdelta_true/EI;
theta_diff=theta_diff/EI;
delta_diff=delta_diff/EI;

%...Output 
fprintf('\n\nApproximations for Simply ');
fprintf('Supported Beam\n');
fprintf(    '--------------------------');
fprintf('--------------\n');
fprintf( ...
  '\nBeam length:                  %g', L);
fprintf( ...
  '\nModulus of elasticity:        %g', E);
fprintf( ...
  '\nMoment of inertia:            %g', I);
fprintf( ...
  '\nNumber of deflection points:  %g', ...
  no_defl_pts);
fprintf( ...
  '\nNumber of idealized loadings: %g', ...
  no_loads);
for i=1:no_loads
  fprintf('\n\nLoad number %g', i);
  fprintf(  '\n  Starting position');
  fprintf(  '\n    x:              %g', ...
    start_x(i));
  fprintf(  '\n    load magnitude: %g', ...
    start_val(i));
  fprintf(  '\n  Ending position');
  fprintf(  '\n    x:              %g', ...
    end_x(i));
  fprintf(  '\n    load magnitude: %g', ...
    end_val(i));
  fprintf(  '\n  # subdivisions:   %g', ...
    no_sub_loads(i));
end
fprintf('\n\nLeft pin reaction:');
fprintf(  '\n  Approximate: %g',left_react);
fprintf(  '\n  True:        %g', ...
  left_react_true);
fprintf('\n\nRight pin reaction:');
fprintf(  '\n  Approximate: %g',right_react);
fprintf(  '\n  True:        %g', ...
  right_react_true);
fprintf('\n\n')

%...Plot the shear results
clf;
subplot(2,1,1);
  plot(x,V,'o',x,V_true,'-');
  xlabel('x'); ylabel('Shear Force');
  title('Shear Diagram'); 
subplot(2,1,2);
  plot(x,V_diff)
  xlabel('x'); ylabel('Difference');
  title('Difference'); drawnow;
% genprint('shear');
disp('Press a key to continue'); pause;

%...Plot the moment results
clf;
subplot(2,1,1);
  plot(x,M,'o',x,M_true,'-');
  xlabel('x'); ylabel('Moment');
  title('Moment Diagram'); 
subplot(2,1,2);
  plot(x,M_diff)
  xlabel('x'); ylabel('Difference');
  title('Difference'); drawnow;
% genprint('moment');
disp('Press a key to continue'); pause;

%...Plot the rotation results
clf;
subplot(2,1,1);
  plot(x,theta,'o',x,theta_true,'-');
  xlabel('x'); ylabel('Rotation (radians)');
  title('Rotation Diagram'); 
subplot(2,1,2);
  plot(x,theta_diff)
  xlabel('x'); ylabel('Error Measure');
  title('Difference'); drawnow;
% genprint('theta');
disp('Press a key to continue'); pause;

%...Plot the displacement results
clf;
subplot(2,1,1);
  plot(x,delta,'o',x,delta_true,'-');
  xlabel('x'); ylabel('Displacement');
  title('Deflection Diagram'); 
subplot(2,1,2);
  plot(x,delta_diff)
  xlabel('x'); ylabel('Error Measure');
  title('Difference'); drawnow;
% genprint('displ');