prplot(Sigma_max,Sigma_min, ...

function prplot(Sigma_max,Sigma_min, ...
  Sigma_avg,Tau_mxmn,Sigma_max_ang, ...
  Tau_max_ang,Sigma_x,Sigma_y,Tau_xy)
%
% prplot(Sigma_max,Sigma_min, ...
%   Sigma_avg,Tau_mxmn,Sigma_max_ang, ...
%   Tau_max_ang,Sigma_x,Sigma_y,Tau_xy)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function plots the element, the
% principal normal stresses and the shear 
% stresses for an element.
%
% Sigma_max     - Sigma maximum
% Sigma_min     - Sigma minimum
% Sigma_avg     - Sigma on planes of max/min
%                 shear stress
% Tau_mxmn      - Tau maximum, Tau minimum is
%                 negative of Tau maximum
% Sigma_max_ang - angle to Sigma_max (degs)
% Tau_max_ang   - angle to Tau_max (degs)
% Sigma_x       - sigma in the x direction
% Sigma_y       - sigma in the y direction
% Tau_xy        - tau in the xy plane
%
% NOTE: To make the plotting perform 
%       as desired an undocumented MATLAB
%       capability has been used (Renderlimits).
%       MathWorks does not promise that it
%       will be supported in the future.
%       Renderlimits allows the two side-by-
%       side plots to maintain the same 
%       y-axes. (no longer applicable!)
%
% User m functions called:  arrows
%----------------------------------------------

degrad=pi/180; pi2=pi/2;
%...Define a unit square for element
Element=[-1 -1; 1 -1; 1 1; -1 1; -1 -1];
Raxis=[0 0;1 0];

%...................
%...Original element
%...................
%...Arrowhead direction
dirSx=sign(Sigma_x); dirSy=sign(Sigma_y);
dirTxy=sign(Tau_xy);

clf; ax=subplot(2,2,1);
  plot(Element(:,1),Element(:,2),'-');
    hold on;
    plot(Raxis(:,1),Raxis(:,2),':');
    if Sigma_x ~= 0
      arrows(dirSx,0,0,Sigma_x,0,0,1,0);
      arrows(dirSx,0,pi,0,0,0,0,0);
    end
    if Sigma_y ~= 0
      arrows(dirSy,0,pi2,0,0,0,0,0);
      arrows(dirSy,0,-pi2,0,Sigma_y,0,2,0);
    end
    if Tau_xy ~= 0
      arrows(-dirTxy,-1,0,0,0,Tau_xy,3,Element);
      arrows(-dirTxy,-1,pi,0,0,0,0,0);
      arrows(dirTxy,-1,pi2,0,0,0,0,0);
      arrows(dirTxy,-1,-pi2,0,0,0,0,0);
    end
    title('Original Element');
    axis('equal'); axis('off');
    hold off; drawnow;

%..............
%...Text angles
%..............
str1=['Sigma-max angle = ', ...
      num2str(Sigma_max_ang),' degs'];
str2=['Tau-max angle = ', ...
      num2str(Tau_max_ang),' degs'];
subplot(2,2,2);
  axis([0 1 0 1]); axis('off');
  text(0.0,0.6,str1);
  text(0.0,0.4,str2);
  drawnow;

%.....................
%...Principal stresses
%.....................
angle=Sigma_max_ang*degrad;
S1=Sigma_max; S2=Sigma_min;

%...Arrowhead direction
dirSx=sign(Sigma_max); dirSy=sign(Sigma_min);

%...Rotate by angle
trans=[cos(angle) -sin(angle);
       sin(angle)  cos(angle)];
temp=trans*Element'; Elementt=temp';
temp=trans*Raxis'; Raxist=temp';

subplot(2,2,3);
  plot(Elementt(:,1),Elementt(:,2),'-');
    hold on;
    plot(Raxist(:,1),Raxist(:,2),':');
    if S1 ~= 0
      arrows(dirSx,0,angle,S1,0,0,1,Elementt);
      arrows(dirSx,0,pi+angle,0,0,0,0,0);
    end
    if S2 ~= 0
      arrows(dirSy,0,pi2+angle,0,0,0,0,0);
      arrows(dirSy,0,-pi2+angle,0,S2,0,2,0);
    end
    title('Principal Normal Stresses');
    set(gca,'DataAspectRatio',[1 1 1]);
    axis('off'); hold off;
    drawnow;

%.........................
%...Max/min shear stresses
%.........................
angle=Tau_max_ang*degrad;
S1=Sigma_avg; T=Tau_mxmn;

%...Arrowhead direction
dirSx=sign(Sigma_avg); dirSy=sign(Sigma_avg);
dirTxy=sign(Tau_mxmn);

%...Rotate by angle
trans=[cos(angle) -sin(angle);
       sin(angle)  cos(angle)];
temp=trans*Element'; Elementt=temp';
temp=trans*Raxis'; Raxist=temp';

subplot(2,2,4);
  plot(Elementt(:,1),Elementt(:,2),'-');
    hold on;
    plot(Raxist(:,1),Raxist(:,2),':');
    arrows(dirSx,0,angle,S1,0,0,1,0);
    arrows(dirSx,0,pi+angle,0,0,0,0,0);
    arrows(dirSy,0,pi2+angle,0,0,0,0,0);
    arrows(dirSy,0,-pi2+angle,0,S1,0,2,0);
    if T ~= 0
      arrows(-dirTxy,-1,angle,0,0,0,0,0);
      arrows(-dirTxy,-1,pi+angle,0,0,T,3, ...
             Elementt);
      arrows(dirTxy,-1,pi2+angle,0,0,0,0,0);
      arrows(dirTxy,-1,-pi2+angle,0,0,0,0,0);
    end
    title('Maximum Shear Stresses');
    axis('off'); hold off;
    set(gca,'DataAspectRatio',[1 1 1]);
    drawnow;