tworods.m

% Example: tworods
% ~~~~~~~~~~~~~~~~
% This example analyzes the effects of a load
% applied at the point of interconnect for
% two rods.  The user defines a range of
% lengths for the left rod (AB) and a range
% of angles for rod AB measured clockwise from
% the horizontal axis.
%
% Data is defined in the declaration statements
% below, where:
%
% P           - downward vertical load
% L_ac        - horizontal span for AC
% A_ab        - area of rod AB
% Gamma_ab    - density of rod AB
% Gamma_bc    - density of rod BC
% Theta_range - two element vector with the
%               starting and ending values
%               of theta to evaluate
%               (in degrees)
% No_thetas   - number of thetas to evaluate
% L_ab_range  - two element vector with the 
%               starting and ending values
%               of length for AB to evaluate
% No_L_abs    - number of lengths of AB to
%               evaluate
% A_bc_range  - two element vector with the
%               minimum and maximum acceptable
%               values for the area of BC
%
% User m functions required: 
%    genprint
%----------------------------------------------

clear;
%...Input definitions
Problem=1;
if Problem == 1
  P=5000; L_ac=100; A_ab=0.4;
  Gamma_ab=0.1; Gamma_bc=0.3;
  Theta_range=[5 85]; No_thetas=[20];
  L_ab_range=[5 95]; No_L_abs=[45];
  A_bc_range=[0.3 0.5];
end 

%...Set up some quantities
Theta=linspace(Theta_range(1), ...
               Theta_range(2),No_thetas+1);
L_ab=linspace(L_ab_range(1), ...
              L_ab_range(2),No_L_abs+1);
degrad=pi/180; thetai=Theta*degrad;
costheta=cos(thetai); sintheta=sin(thetai);
Lac2=L_ac*2; Lacsq=L_ac^2;
Lac2cos=Lac2*costheta; Labsq=L_ab.^2;

%...Loop on each theta
for i=1:No_thetas+1
  %...Loop on each length 
  for j=1:No_L_abs+1
    %...Law of cosines
    L_bc(i,j)=sqrt(Labsq(j)+Lacsq- ...
              L_ab(j)*Lac2cos(i));
    %...Rod BC slope
    Vbc=L_ab(j)*sintheta(i);
    Hbc=L_ac-L_ab(j)*costheta(i);
    %...Set up simultaneous equations and solve
    a=[-costheta(i)  Hbc/L_bc(i,j); ...
        sintheta(i)  Vbc/L_bc(i,j)];
    b=[0;P]; F=a\b;
    F_ab(i,j)=F(1); F_bc(i,j)=F(2);
    %...Calculate stress and weight
    Sigma=F(1)/A_ab; A_bc(i,j)=F(2)/Sigma;
    W(i,j)=Gamma_ab*A_ab*L_ab(j)+ ...
           Gamma_bc*A_bc(i,j)*L_bc(i,j);
  end
end

%...Plot results
clf; surf(L_ab,Theta,F_ab);
  title('Force in AB');
  xlabel('Theta (degrees)'); 
  ylabel('Length of AB');
  zlabel('Force in AB'); 
  view(135.0,45); drawnow;
% genprint('Fab');
  disp('Press key to continue'); pause;

clf; surf(L_ab,Theta,F_bc)
  title('Force in BC');
  xlabel('Theta (degrees)');
  ylabel('Length of AB');
  zlabel('Force in BC'); 
  view(135.0,45); drawnow;
% genprint('Fbc');
  disp('Press key to continue'); pause;

clf; surf(L_ab,Theta,W);
  title('Total Weight'); 
  xlabel('Theta (degrees)');
  ylabel('Length of AB');
  zlabel('Total Weight'); 
  view(135.0,45); drawnow; a1=axis;
% genprint('W');
  disp('Press key to continue'); pause;

clf; surf(L_ab,Theta,A_bc);
  title('Area of BC'); 
  xlabel('Theta (degrees)');
  ylabel('Length of AB');
  zlabel('Area of BC'); 
  view(135.0,45); drawnow; a2=axis;
% genprint('Abc');
  disp('Press key to continue'); pause;

%...Choose only prescribed range of Abc
for i=1:No_thetas+1
  for j=1:No_L_abs+1
    if A_bc(i,j) < A_bc_range(1) | ...
                      A_bc(i,j) > A_bc_range(2)
      W(i,j)=nan; A_bc(i,j)=nan;
    end
  end
end

clf; surf(L_ab,Theta,W);
  title('Weight (allowable range)');
  xlabel('Theta (degrees)');
  ylabel('Length of AB');
  zlabel('Total Weight');
  view(135.0,45); axis([a1]); drawnow;
% genprint('Wallow');
  disp('Press key to continue'); pause;

clf; surf(L_ab,Theta,A_bc);
  title('Area of BC (allowable range)');
  xlabel('Theta (degrees)');
  ylabel('Length of AB');
  zlabel('Area of BC');
  view(135.0,45); axis([a2]); drawnow;
% genprint('Abcallow');