[Key,Imax,Imin,Iavg, ...

function [Key,Imax,Imin,Iavg, ...
  Ixymxmn,Imax_ang,Imin_ang, ...
  Ixymax_ang,Ixymin_ang]= ...
  prinert(Ix,Iy,Ixy,epsilon)
%
% [Key,Imax,Imin,Iavg,Ixymxmn, ...
%   Imax_ang,Imin_ang,Ixymax_ang, ...
%   Ixymin_ang]=prinert(Ix,Iy,Ixy,epsilon)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function calculates the principal
% inertia axes and values.
% 
% Ix         - inertia about y axis
% Iy         - inertia about x axis
% Ixy        - product of inertia
% epsilon    - tolerance for determining if two 
%              numbers are equal (optional)
%
% Key        - =0, normal exit
%              =1, all axes are principal axes
% Imax       - Maximum principal inertia
% Imin       - Minimum principal inertia
% Iavg       - I for max/min Ixy plane
% Ixymxmn    - Max/min for product of inertia
% Imax_ang   - Angle for Imax 
% Imin_ang   - Angle for Imin
% Ixymax_ang - Angle for max Ixymxmn
% Ixymin_ang - Angle for min Ixymxmn
%
% NOTES:  1)  Positive angles measured CCW from
%             positive x-axis in degrees.
%
% User m functions called:  none
%----------------------------------------------

%...Initialize
raddeg=180/pi; Key=0;

%...Max/min values
Idiff=0.5*(Ix-Iy); Iavg=0.5*(Ix+Iy);
Ixymxmn=sqrt(Idiff^2+Ixy^2);
Imax=Iavg+Ixymxmn; Imin=Iavg-Ixymxmn;

if nargin ~=4
  if Imax == Imin
    epsilon=1e-3;
  else
    epsilon=abs(Imax-Imin)*1e-3;
  end
end

if abs(Imax-Imin) < epsilon
  %...All axes are principal axes
  Key=1; Imax_ang=0; Imin_ang=0; 
  Ixymax_ang=0; Ixymin_ang=0;
  Ixymxmn=0; Imax=Ix; Imin=Ix;
else
  Imax_ang=raddeg*1/2*atan2(-2*Ixy,Ix-Iy);
  Imin_ang=Imax_ang+90; Ixymax_ang=Imax_ang+45;
  Ixymin_ang=Imax_ang-45;
end