function testellipsef (d1,d2,mu1,mu2,ci)
%%% The goal is to plot the 95% confidence ellipse
%%% for n points of bivariate (2-parameter) data.  The mean and covariance is
%%% determined, as well as the inverse covariance matrix.  Equations using the eigenvalues and
%%% eigenvectors are used to determine the confidence ellipse axes half-
%%% lengths.  The standard ellipse equation is used (furnished by
%%% Danielle Rabina of the MathWorks company), and then rotated by the
%%% "rotate" command.  The angle of rotation is determined by using the
%%% equation in "Data Analysis for Scientists and Engineers," by Stuart
%%% L. Meyer, 1975, by John Wiley & Sons, page 289, right side, second
%%% equatioon for theta as defined in the textbook.  The size of the
%%% ellipse is determined by setting the desired probability, 1 - alpha,
%%% and using the F-Distribution Tables.
%%%
%%% ********************************************************************



cov_x = cov(d1);
cov_x_inv = inv(cov_x);

cov_y = cov(d2);
cov_y_inv = inv(cov_y);

XY = [d1'; d2'];


M = mean(XY')';
S = cov(d1, d2);
% S=[20 100;100 1000];
S_inv = inv(S);

d = eig(S')'
[V,D] = eig(S')

%%  --------------------------------------------------------

%%  n = 42;  p = 2;
%%  The 95% confidence ellipse for mu consist of all values
%%  (mu1, mu2) satisfying (5.7) pp. 119 in Rencher:
%%  n * [mean1 - mu1; mean2 - mu2] * inv(S) *
%%  [mean1 - mu1; mean2 - mu2]'
%%     <=  (p(n - 1)/(n - p)) * F_p,n-p(alpha)
%%
%%  ----------------------------------------------------------

%%  Confidence Ellipsoid:

n = length(d1)-1;
p = 2;   
mean1 = M(1,1)   
mean2 = M(2,1)
%%  F_p_n-p is from can be determined from any set of F-tables but the
%%  inverse F is not as readily accessible. 
%%  MATLAB has an inverse F-Distribution statement to look
%%  into the F-Distribution Tables.

% F_2_40 = 3.23; % from F table

F_2_40 = finv(ci,p,n-p);%95% confidence,nu1=p; nu2=n-p

x_bar = M';

term1 = (p*(n-1))/(n*(n-p));

axis1 = sqrt(d(1,1))*sqrt(term1*F_2_40); %half-axis
axis2 = sqrt(d(1,2))*sqrt(term1*F_2_40); %half-axis

%%  ----------------------------------------------------------
%%  Ellipsoid Continued . . .

h = mean1;     k = mean2;
h= mu1; k=mu2;
a = 2*axis1;   b = 2*axis2;    ellipse_area = pi*a*b;

th=0:2*pi/100:2*pi;
x=a*cos(th)+h;
y=b*sin(th)+k;
figure(2)
clf
plot(x,y)
%   axis equal
grid on
hold on

sig_x = sqrt(S(1,1));
sig_y = sqrt(S(2,2));
sig_x_2 = S(1,1);
sig_y_2 = S(2,2);

rho1 = corrcoef(XY(1,:),XY(2,:));
rho1 = rho1(1,2);
term2 = (2 * rho1 * sig_x * sig_y);
term3 = sig_x_2 - sig_y_2;

angle = 0.5 * atan2(term2,term3)*180/pi

%%  These next 2 lines grab the ellipse and rotate it:

h1 = findobj('type','line');
rotate(h1,[0 0 1],-angle)

%%  -----------------------------------------------------------------

x1 = mu1;
y1 = mu2;
% plot(d1,d2,'o',x1,y1,'r*')
plot(x1,y1,'r*')
title('The 95% Confidence Ellipse for Mean Vectors')
xlabel('Parameter Number 1')
ylabel('Parameter Number 2')
gtext('Mean Vector under test (Asterisk):')
gtext('Ellipse Center:  Diamond')

x2 = mean1;
y2 = mean2;
plot(x2,y2,'bd')

hold off