function varargout = showVarianceEllipseWithAxes(varargin)
% function varargout = showVarianceEllipseWithAxes(varargin)
%
% A plotting routine that displays an ellipsoid for a particular variance
% together with its axes.
%
% Default params:
% variances = [1 2 3];
% xaxiscol = 'm';
% yaxiscol = 'c';
% zaxiscol = 'y';
% axiswidth = 1;
% ellipcenter = [0 0 0];
% color = [132 132 255]/255;
% axish = [];
% std_lim =2;
%
% Example:
%
% h = showVarianceEllipseWithAxes;
%
% Example:
%
% h = showVarianceEllipseWithAxes('variances', [1 2 3]);
%
% Dr. A. I. Hanna (2007)

error(nargchk(0,inf,nargin));

variances = [1 2 3];
xaxiscol = 'm';
yaxiscol = 'c';
zaxiscol = 'y';
axiswidth = 1;
ellipcenter = [0 0 0];
color = [132 132 255]/255;
axish = [];
std_lim =2;
if mod(length(varargin),2) ~= 0
    % input args have not com in pairs, woe is me
    error(['Arguments to ', mfilename, ' must come param/value in pairs.'])
end
for i=1:2:length(varargin)
    switch lower(varargin{i})
        case 'variances'
            variances = varargin{i+1};
        case 'color'
            color = varargin{i+1};
        case 'parent'
            axish = varargin{i+1};
        case 'std_lim'
            std_lim = varargin{i+1};
        case 'offset'
            ellipcenter = varargin{i+1};
        case 'xaxiscol'
            xaxiscol = varargin{i+1};
        case 'yaxiscol'
            yaxiscol = varargin{i+1};
        case 'zaxiscol'
            zaxiscol = varargin{i+1};
        otherwise
            error(['Unknown parameter name passed to ', mfilename, '.  Name was ' varargin{i}])
    end
end
if isempty(axish)
    axish = gca;
end

[x,y,z] = sphere(20);
[mm,nn] = size(x);
xsd = sqrt(variances(1));
ysd = sqrt(variances(2));
zsd = sqrt(variances(3));
C = [std_lim*x(:)*sqrt(variances(1)) std_lim*y(:)*sqrt(variances(2)) std_lim*z(:)*sqrt(variances(3))];
x = C(:,1); y = C(:,2); z = C(:,3);
x = reshape(x,mm,nn);
y = reshape(y,mm,nn);
z = reshape(z,mm,nn);
x = x + ellipcenter(1);
y = y + ellipcenter(2);
z = z + ellipcenter(3);
ellip_h = surf(axish, x, y, z);
set(ellip_h, 'FaceAlpha', .2, 'FaceColor', color, 'EdgeColor', color, 'Visible', 'on');
lighting(axish, 'phong');
axis(axish, 'equal');
hold(axish, 'on');
xaxis_h = plot3(axish, [-std_lim*xsd, std_lim*xsd], [0 0], [0 0], 'Color', xaxiscol, 'LineWidth', axiswidth);
yaxis_h = plot3(axish, [0 0], [-std_lim*ysd, std_lim*ysd], [0 0], 'Color', yaxiscol, 'LineWidth', axiswidth);
zaxis_h = plot3(axish, [0 0], [0 0], [-std_lim*zsd, std_lim*zsd], 'Color', zaxiscol, 'LineWidth', axiswidth);
ellip_handles.ellip_h =ellip_h;
ellip_handles.xaxis_h =xaxis_h;
ellip_handles.yaxis_h =yaxis_h;
ellip_handles.zaxis_h =zaxis_h;
varargout{1} = ellip_handles;
return