function [SEMA,  ECC, INC, PHA, w]=app2ep(Au, PHIu, Av, PHIv, plot_demo)
%
% Convert tidal amplitude and phase lag (ap-) parameters into tidal ellipse
% (e-) Parameters.
% 
% Usage:
%
% [SEMA,  ECC, INC, PHA, w]=app2ep(Au, PHIu, Av, PHIv, plot_demo)
%
% where:
%
%     Au, PHIu, Av, PHIv are the amplitudes and phase lags of 
%     u- and v- tidal current components. They can be vectors or 
%     matrices;
%            
%     plot_demo is an optional argument, when is supplied as a non-zero
%     number, the program will plot an  ellipse corresponding to
%     Au(1), PHIu(1), Av(1), and PHIv(1);
%
%     SEMA: Semi-major axes, or the maximum speed;
%     ECC:  Eccentricity, the ration of the semi-major axes over the 
%           semi-minor axis;
%     INC:  Inclination, the angles  between the semi-major axes and 
%           u-axis.                        
%     PHA:  Phase angles, the time (in angles)  when the 
%           tidal currents reach their maximum speeds,  (i.e. 
%           PHA=omega*tmax)
%          
%           These four e-parameters will have the same dimensionality 
%           (i.e., vectors, or matrices) as the input ap-parameters. 
%
%     w:    Optional. If it is requested, it will be output as matrices
%           whose rows allow for plotting ellipses and whose columns are  
%           for different ellipses corresponding columnwise to SEMA. For
%           example, plot(real(w(1,:)), imag(w(1,:))) will let you see 
%           the first ellipse.
%
% Document:   tidal_ellipse_v2.ps


  if nargin < 5
     plot_demo=0;  % by default, no plot for the ellipse
  end

% Transform the input data as column vectors
   [sizeAu_m, sizeAu_n]=size(Au); % for later transform them back.
   if sizeAu_n > 1
	   Au = Au(:);
	   Av = Av(:);
	   PHIu = PHIu(:);
	   PHIv = PHIv(:);
   end
% Assume the input phase lags are in degrees and convert them in radians.
   PHIu = PHIu/180*pi;
   PHIv = PHIv/180*pi;

% Make complex amplitudes for u and v
   i = sqrt(-1);
   u = Au.*exp(i*PHIu);
   v = Av.*exp(i*PHIv);

% Calculate complex radius of anticlockwise and clockwise circles:
   wp = (u-i*v)/2;      % for anticlockwise circles, and yes, it is (u-i*v)/2 
   wm = (u+i*v)/2;      % for clockwise circles
% and their amplitudes and angles
   Wp = abs(wp);
   Wm = abs(wm);
   THETAp = -angle(wp);
   THETAm = angle(wm);
   
% calculate e-parameters (ellipse parameters)
    SEMA = Wp+Wm;              % Semi  Major Axis, or maximum speed
    SEMI = Wp-Wm;              % Semin Minor Axis, or minimum speed
     ECC = SEMI/SEMA;          % Eccentricity

    PHA = (THETAm-THETAp)/2;   % Phase angle, the time (in angle) when 
                               % the velocity reaches the maximum
    INC = (THETAm+THETAp)/2;   % Inclination, the angle between the 
                               % semi major axis and x-axis (or u-axis).

    % convert to degrees for output
    PHA = PHA/pi*180;         
    INC = INC/pi*180;         
    THETAp = THETAp/pi*180;
    THETAm = THETAm/pi*180;
    
  % flip THETAp and THETAm, PHA, and INC in the range of 
  % [-pi, 0) to [pi, 2*pi), which at least is my convention.
    id = THETAp < 0;   THETAp(id) = THETAp(id)+360;
    id = THETAm < 0;   THETAm(id) = THETAm(id)+360;
    id = PHA < 0;      PHA(id) = PHA(id)+360;
    id = INC < 0;      INC(id) = INC(id)+360;
    
  %output e-parameter in the same matrix dimension as Au 
  if sizeAu_n > 1
       SEMA = reshape(SEMA, sizeAu_m, sizeAu_n);
       SEMI = reshape(SEMI, sizeAu_m, sizeAu_n);
       PHA  = reshape(PHA,  sizeAu_m, sizeAu_n);
       INC  = reshape(INC,  sizeAu_m, sizeAu_n);
       ECC  = reshape(ECC,  sizeAu_m, sizeAu_n);
    end

  if nargout == 5 | plot_demo
     dot=pi/18;
     ot=[0:dot:2*pi-dot];
     w=wp'*exp(i*ot)+wm*exp(-i*ot);
  end
  

% Plot demo
  if plot_demo     
    wmax = SEMA.*exp(i*INC/180*pi);
    wmin = SEMI.*exp(i*(INC/180*pi+pi/2));
  %  wmax1 = wp'.*exp(i*(PHA/180*pi))     +  wm.*exp(-i*(PHA/180*pi));
  %  wmin1 = wp'.*exp(i*(PHA/180*pi+pi/2))+  wm.*exp(-i*(PHA/180*pi+pi/2));
  %  wmax1 and wmax should be equal to each other, so are wmin and wmin1.  

  figure(gcf)
  clf
     n=1;   %just to plot the first ellipse for a demo.
     plot(real(w(n,:)), imag(w(n,:)))
     axis('equal');
     hold on
     plot([0 real(wmax)], [0 imag(wmax)], 'm');
     plot([0 real(wmin)], [0 imag(wmin)], 'm');
     xlabel('u')
     ylabel('v');
     title([num2str(n) '''st ellipse'])
  
     THETAp=THETAp/180*pi;
     THETAm=THETAm/180*pi;
     a=Wp.*exp(i*(ot+THETAp));
     b=Wm.*exp(-i*(ot-THETAm));
     c=a+b;

     plot(real(a), imag(a), 'r');
     plot(real(b), imag(b), 'g');
     hnd_a=line([0 real(a(1))], [0 imag(a(1))], 'color', 'r', 'marker','o');
     hnd_b=line([0 real(b(1))], [0 imag(b(1))], 'color', 'g', 'marker','o');
     hnd_c=line([0 real(c(1))], [0 imag(c(1))], 'color', 'b', 'marker','o');
     plot(real(a(1)), imag(a(1)), 'ro');
     plot(real(b(1)), imag(b(1)), 'go');
     plot(real(c(1)), imag(c(1)), 'bo');
     hnd_ab=line(real([a(1) a(1)+b(1)]), imag([a(1) a(1)+b(1)]), ... 
              'linestyle', '--', 'color', 'g');
     hnd_ba=line(real([b(1) a(1)+b(1)]), imag([b(1) a(1)+b(1)]), ... 
              'linestyle', '--', 'color', 'r');
  
     for n=1:length(ot);
         set(hnd_a, 'xdata', [0 real(a(n))], 'ydata', [0 imag(a(n))]);
         set(hnd_b, 'xdata', [0 real(b(n))], 'ydata', [0 imag(b(n))]);
         set(hnd_c, 'xdata', [0 real(c(n))], 'ydata', [0 imag(c(n))]);     
              hold on
         plot(real(a(n)), imag(a(n)), 'ro');
         plot(real(b(n)), imag(b(n)), 'go');
         plot(real(c(n)), imag(c(n)), 'bo');
         set(hnd_ab, 'xdata',real([a(n) a(n)+b(n)]), 'ydata', ... 
         imag([a(n) a(n)+b(n)]))
         set(hnd_ba, 'xdata',real([b(n) a(n)+b(n)]), 'ydata', ... 
         imag([b(n) a(n)+b(n)]))
     end


  end


%Authorship Copyright:
%
%    The author of this program retains the copyright of this program, while
% you are welcome to use and distribute this program as long as you credit 
% the author properly. Particularly, you are expected to retain the
% origianl author's name in this original version of the program 
% or any of its modifcations that you might make. Enjoying my program!
%
%Author Info:
%_______________________________________________________________________
%  Zhigang Xu, Ph.D.                            
%  (pronounced as Tsi Gahng Hsu)
%  Research Scientist
%  Coastal Circulation                   
%  Bedford Institute of Oceanography     
%  1 Challenge Dr.
%  P.O. Box 1006                    Phone  (902) 426-2307 (o)       
%  Dartmouth, Nova Scotia           Fax    (902) 426-7827            
%  CANADA B2Y 4A2                   email zhigangx@emerald.bio.dfo.ca    
%_______________________________________________________________________
%



