%
% function ellipse(ox,oy,maj,mnr,ang,nptspq,filcol,scal)
%
% routine to plot an ellipse given
%
%  ox, oy  = the x&y origin
%  maj     = major axis
%  mnr     = minor axis
%  ang     = angle counterclockwise from due east
%  nptspq  = number of points per quarter ellipse (default=20)
%  filcol  = fill color;
%  scal    = 0 no scale
%          = 1 add a scale
%
%
% assume a 1 degree latitude plot and scale the magnitudes accordingly
% so that say a 100 cm/s current ellipse is plotted over
% 0.2 degree of geographic space
%
%
function ellipse(ox,oy,maj,mnr,ang,nptspq,filcol,scal)
  scale=.6/100;
  maj=maj*scale;
  mnr=mnr*scale;
if nargin<6;
  nptspq=20;
end;
if nargin<7;
  filcol=2;
end;
if nargin<8;
  scal=0;
end;
%
% find the minumum latitude given info input
% first assume the major axis is n/s
% increment through longitude
%
  inc=mnr/nptspq;
%
for i=-nptspq:1:nptspq; 
  x(i+nptspq+1)=inc*i;
  y(i+nptspq+1)=(maj/mnr)*sqrt(mnr^2-x(i+nptspq+1)^2);
end;
%
% so far, we only have the upper have of the ellipse
%
for i=nptspq*2+2:nptspq*4+2;
  x(i)=x(nptspq*4+2-i+1);
  y(i)=y(nptspq*4+2-i+1)*(-1);
end;

plot(x,y);
%pause(3);
%
% then rotate this vector by and angle "ang"
%
  ang=(ang-90)*(-1);
for i=1:length(x);
  [nx(i),ny(i)]=rot(ang,x(i),y(i));
end;
%
%plot(nx,ny);
%fill(nx,ny,1)
%
% now add the origin
  nx2=nx+ox;
  ny2=ny+oy;
  plot(nx2,ny2);
  if length(filcol)==1;
    fill(nx2,ny2,['c',num2str(filcol)])
  else;
    fill(nx2,ny2,filcol)
  end;

if scal==1; % add a scal to the ellipse
  [minx,miny]=rot(ang,mean(x),min(y));
  [maxx,may]=rot(ang,mean(x),max(y));
  plot([minx+ox maxx+ox],[miny+oy  may+oy],...
   'linewidth',3,'color','w');
  text(ox,oy,[num2str(maj/scale,2),' cm/s'],'color','w',...
    'horizontalalignment','left',...
    'verticalalignment','bottom');
end;
