function truss_plot(X,Y,EqMat,T,axisrange)
% truss_plot(X,Y,EqMat,T,axisrange)
% Plot a truss configuration as descibed by bar length components X,Y,
% equilibrium matrix EqMat, and, optionally truss bar forces, T.
% If truss bar forces are provided, they are displayed on the figure.
%
%    INPUT          DESCRIPTION
%   -------        -------------
%     X             The X-component dimension of each bar in a truss
%     Y             The Y-component dimension of each bar in a truss
%   EqMat           A matrix representing all the equilibrium equations.  
%                   Adjacent rows represents the equilibrium equations
%                   for forces in the X-direction and in the Y-direction.
%                   X-direction equilibrium equations are in the odd rows.
%                   Y-direction equilibrium equations are in the even rows.
%                   Each column corresponds to each bar.  
%     T             Optional vector of the truss bar forces in each bar.  
%  axisrange        Optional plot axis range [ x_min x_max y_min y_max ]
%  
%   OUTPUT          DESCRIPTION
%   -------        -------------
%   figure(9)       A plot of the truss.  red=compression, blue=tension


 [joints,bars] = size(EqMat);
 joints = joints/2;
 bar   = zeros(2,bars);                 % bar connectivity

 L = sqrt(X.^2 + Y.^2);

 % Which joints does each bar connect?

 for b = 1:bars                   % bar "b" connects joint "j1" to joint "j2"
     j1 = ceil(min(find(EqMat(:,b)))/2);  % joint "1" of bar "b"
     j2 = ceil(max(find(EqMat(:,b)))/2);  % joint "2" of bar "b"
     bar(:,b) = [ j1 ; j2 ];
 end
 % In each column of "bar", the first row should be less than the second row. 

 % What are the coordinates of each joint?

 Xj = zeros(1,joints);           % initial X- coordinates of each joint
 Yj = zeros(1,joints);           % initial Y- coordinates of each joint

 joint_ok = 1;                   % list of joints with found coordinates
                                 % assume joint 1 is located at point (0,0)
 for j=1:joints

     b1 = find(bar(1,:)==j);             % bars with "j1" == j
     j2 = bar(2,b1);                     % j2 of bars with "j1" == j
     for jj=1:length(j2)
         if ~any(j2(jj) == joint_ok)
            Xj(j2(jj)) = Xj(bar(1,b1(jj))) + L(b1(jj))*EqMat(2*j-1,b1(jj));
            Yj(j2(jj)) = Yj(bar(1,b1(jj))) + L(b1(jj))*EqMat(2*j,b1(jj));
            joint_ok = [ joint_ok bar(2,b1(jj)) ];
         end
     end

     b2 = find(bar(2,:)==j);            % bars with "j2" == j
     j1 = bar(1,b2);                    % j1 of bars with "j2" == j
     for jj=1:length(j1)
         if ~any(j1(jj) == joint_ok)
            Xj(j1(jj)) = Xj(bar(1,b2(jj))) - L(b2(jj))*EqMat(2*j-1,b2(jj));
            Yj(j1(jj)) = Yj(bar(1,b2(jj))) - L(b2(jj))*EqMat(2*j-1,b2(jj));
            joint_ok = [ joint_ok bar(1,b2(jj)) ];
         end
     end

 end
%

 figure(9)              % PLOT IT!
   clf                  % clear the prior plot

   plot(Xj,Yj,'*k');    % plot joint locations as stars


   if nargin >= 5
      axis(axisrange);
   else
      % set up default axis ranges with margins around truss figure
      Xmin = min(Xj);  Xmax = max(Xj);  Ymin = min(Yj);  Ymax = max(Yj);
      Xmrg = 0.10*(Xmax - Xmin);  Ymrg = 0.10*(Ymax - Ymin);

      axis([ Xmin-Xmrg , Xmax+Xmrg , Ymin-Ymrg , Ymax+Ymrg ]);
   end
   
   for j=1:joints	% display joint numbers
       text(Xj(j),Yj(j),sprintf('%3d',j));
   end

   hold on

   for b = 1:bars       % plot the bars one at a time and display bar forces

       j1 = bar(1,b);   j2 = bar(2,b);   % joints connected bar bar "b"
       if nargin >= 4
          if T(b) > 0         % compression = red 
             plot([Xj(j1) Xj(j2)],[Yj(j1) Yj(j2)],'-b');
          elseif T(b) < 0     % tension = blue
             plot([Xj(j1) Xj(j2)],[Yj(j1) Yj(j2)],'-r');
          else
             plot([Xj(j1) Xj(j2)],[Yj(j1) Yj(j2)],'-g');
          end
          text((Xj(j1)+Xj(j2))/2,(Yj(j1)+Yj(j2))/2,sprintf('%6.1f',T(b)));
       else
          plot([Xj(j1) Xj(j2)],[Yj(j1) Yj(j2)],'-g');
          text((Xj(j1)+Xj(j2))/2,(Yj(j1)+Yj(j2))/2,sprintf('%3d',b));
       end

   end

   axis('equal')
   drawnow;

% --------------------------------------------------------------  TRUSS_PLOT