function [SEMA,  ECC, INC, PHA, w]=app2ep(Au, PHIu, Av, PHIv, plot_demo)
% Convert u and v amp and phase lag (ap-) parameters into 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 (degrees) 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
%     integer, the program will plot an  ellipse corresponding to
%     Au(plot_demo), PHIu(plot_demo), Av(plot_demo), and PHIv(plot_demo);
%
%     SEMA: Semi-major axes, or the maximum speed.
%     ECC:  Eccentricity, the ratio of the semi-minor axis to the
%           semi-major axis.
%     INC:  Inclination, the angle (in degrees) between the x-axis and the 
%           semi-major axis.
%     PHA:  Phase angle (in degrees) 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, this matrix contains the outlines
%           of the calculated ellipses (for plotting purposes).  This matrix will
%           have one dimension greater than your input arrays.  
%           Usage example:  You input 100x50 arrays of values and you want to see
%           the ellipse at (50,4).  You type:  
%                       zi=squeeze(w(50,4,:));
%                       polar(angle(zi),abs(zi));
%
% Document:   tidal_ellipse_v2.ps

% Modified by Rich Signell (October 13, 2000):
%   1) Fixed code so it works with arbitrary array sizes on input
%   2) Changed definition of "wp" so that it correctly refers to the complex amplitude
%      of the positively rotating part

  if nargin < 5
     plot_demo=0;  % by default, no plot for the ellipse
  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 
   wm = conj(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
  % flip THETAp and THETAm, PHA, and INC in the range of 
  % [-pi, 0) to [pi, 2*pi), which at least is my convention.
    PHA = mod(PHA/pi*180,360);         
    INC = mod(INC/pi*180,360);         
    THETAp = mod(THETAp/pi*180,360);
    THETAm = mod(THETAm/pi*180,360);
    
 if nargout == 5 | plot_demo
    ndot=36; % number of points (dots) for drawing ellipses
    dot=2*pi/ndot;
    ot=[0:dot:2*pi-dot];
    w=wp(:)*exp(i*ot)+wm(:)*exp(-i*ot);
% Plot demo
    if plot_demo     
      wmax = SEMA.*exp(i*INC/180*pi);
      wmin = SEMI.*exp(i*(INC/180*pi+pi/2));
      figure(gcf)
      clf
      n=plot_demo;   %just to plot the n'th ellipse for a demo.
      plot(real(w(n,:)), imag(w(n,:)))
      axis('equal');
      hold on
      plot([0 real(wmax(n))], [0 imag(wmax(n))], 'm');
      plot([0 real(wmin(n))], [0 imag(wmin(n))], 'm');
      xlabel('u')
      ylabel('v');
      title([num2str(n) '''st ellipse'])
      THETAp=THETAp/180*pi;
      THETAm=THETAm/180*pi;
      a=Wp(n).*exp(i*(ot+THETAp(n)));
      b=Wm(n).*exp(-i*(ot-THETAm(n)));
      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
      hold off
    end
%    [m,n,p]=size(wp);
%    w=(reshape(w, [m n p ndot]));

    w=reshape(w, [size(wp) ndot]);

  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    
%_______________________________________________________________________
%

