function scatter(x, y, R, r) 
% function scatter(x, y, R, r)  Generates a scatter plot of y against x.  R 
%    is an optional matrix parameter representing the correlation matrix 
%    between x and y.  R must be 2 by 2.  If a valid R is given, a 
%    concentration ellipse will be plotted.  r is the radius of the ellipse. 
%    If no r is specified, 1 is used. 
if nargin>=3 
   if min((size(R))==[2 2])==0 
      error('correlation matrix not 2 by 2') 
   end 
   R=inv(R); 
   if nargin==4 
      rsqd=r^2; 
   else 
      rsqd=1; 
   end 
   CB=R(1,2)+R(2,1); 
   A=R(1,1); 
   D=R(2,2); 
   if -4*D*rsqd/(CB.^2-4*A*D)<=0 
      error('improper correlation matrix') 
   end 
   x1range=sqrt(-4*D*rsqd/(CB.^2-4*A*D)); 
   x1=[-x1range:x1range/50:x1range]; 
   x2p=(-CB.*x1+sqrt((CB^2).*(x1.^2)-4*D.*(A.*(x1.^2)-rsqd)))./(2*D); 
   x2m=(-CB.*x1-sqrt((CB^2).*(x1.^2)-4*D.*(A.*(x1.^2)-rsqd)))./(2*D); 
   plot(x1,x2p); 
   hold on 
   plot(x1,x2m); 
end 
plot(x, y, '.') 
hold off