function ell=ellipse(cx,sx,cy,sy)

%ELLIPSE  calculates tidal ellipse parameters from scalar analysis
%           results of U and V.

%	  ell=ellipse(cx,sx,cy,sy) where
%                     u=CXcos(wt)+SXsin(wt)
%	                 v=CYcos(wt)+SYsin(wt)

%	  and returns ellipse parameters:

%          ell(:,1) = umaj  = major axis.
%          ell(:,2) = umin  = minor axis, where the sign indicates the
%                             rotation of the ellipse.  A positive sign
%                             of umin indicates counterclockwise rotation.
%          ell(:,3) = inc   = counterclockwise inclination from 
%                             positive x axis (degrees). The inclination
%                             will always be in the upper half plane [0:pi].
%          ell(:,4) = phase = phase relative to t=0 (degrees).
%
%	  Uses notation of M.G.G FOREMAN'S least squares harmonic analysis (1978) 
%      Rich Signell 3-10-88
%

i=sqrt(-1);
[imax,jmax]=size(cx);
if imax>jmax
        cx=cx'; sx=sx'; cy=cy'; sy=sy';
end

template=ones(cx/cx(1))*NaN;
umaj=template;umin=template;inc=template;phase=template;
ell=zeros(length(cx),4);
%

ind=find(cx~=NaN);
cx=cx(ind);sx=sx(ind);cy=cy(ind);sy=sy(ind);
%

p=(cx+sy)/2+i*(cy-sx)/2;
m=(cx-sy)/2+i*(cy+sx)/2;
ap=abs(p);
epsp=angle(p);
am=abs(m);
epsm=angle(m);
umaj(ind)=ap+am;
umin(ind)=ap-am;
ieps=find((epsp+epsm)<0);                % choose inc to range from [0:pi]

if length(ieps)>1
         epsp(ieps)=epsp(ieps)+2*pi*ones(ieps);
else
         epsp(ieps)=epsp(ieps)+2*pi;
end

inc(ind)=.5*(epsp+epsm)*180/pi;
phase(ind)=-.5*(epsp-epsm)*180/pi;
ipha=find(phase<0);
phase(ipha)=phase(ipha)+360;

ell(:,1)=umaj';
ell(:,2)=umin';
ell(:,3)=inc';
ell(:,4)=phase';

