% This function computes the maximal extent of a trajectory (MET) a new
% Lagrangian diagnostic family.
%
% USAGE:
%
% INPUT:
% X: x position
% Y: y position
% T: time span
%
% OUTPUT:
% AD - Absolute Dispersion   :    ||x(t)-x(0)||
% D - current location       :    ||x(t)||
% MAD - Maximal distance     :     max(||x(t)||)
% XMaxD  Maximal extent in x :     Max(x_1(t)) (in direction X)
% XMinD  Minimal extent in x :     Min(x_1(t))
% XMAD   Maximal Shift in x  :     XMaxD-XMinD
% XMean  Mean extent in x    :    (XMaxD+XMinD)/2
% YMaxD,YMinD,YMAD,YMean - same as for X in direction Y
% XYMaxD,XYMinD,XYMAD,XYMean - same as for X,Y, in the direction "o"=[1 1]
%
% <LICENSE
%   Copyright (c) 2013, 2014 Erick Fredj
%
%   This file is part of Metlab.
%
%   Metlab is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   Metlab is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with Metlab.  If not, see <http://www.gnu.org/licenses/>.
% LICENSE>

function [AD,D,MAD,XMaxD,XMinD,XMAD,XMean,YMaxD,YMinD,YMAD,YMean,XYMaxD,XYMinD,XYMAD,XYMean] = compute_met(X,Y,T, varargin)

% AD - Absolute Dispersion   :    ||x(t)-x(0)||
% D - current location       :    ||x(t)||
% MAD - Maximal distance     :     max(||x(t)||)  
% XMaxD  Maximal extent in x :     Max(x_1(t)) (in direction X)
% XMinD  Minimal extent in x :     Min(x_1(t))
% XMAD   Maximal Shift in x  :     XMaxD-XMinD
% XMean  Mean extent in x    :    (XMaxD+XMinD)/2
% YMaxD,YMinD,YMAD,YMean - same as for X in direction Y
% XYMaxD,XYMinD,XYMAD,XYMean - same as for X,Y, in the direction "o"=[1 1]

% default local definition
numFrames = length(T);
numPositions = size(X,1);

% default options
startFrame = 1;
finalFrame = numFrames;
o = [1,1]; % projection [1,1]

% look at any remaining options
k=1;
while k<length(varargin),
    optn=lower(varargin{k});
    switch optn,
        case 'startframe',
            startFrame = varargin{k+1};
        case 'finalframe',
            finalFrame = varargin{k+1};
        case 'direction'
            o = varargin{k+1};
    end;
    k=k+2;
end

if finalFrame > numFrames
    error(['finalFrame should be little to ',num2str(numFrames)]);
end

if startFrame > numFrames
    error(['startFrame should be little to ',num2str(numFrames)]);
end

if startFrame <= 0
    error('startFrame should be bigger to 0');
end

if finalFrame <= 0
    error('finalFrame should be bigger to 0');
end

% local definition
localNumFrames = finalFrame-startFrame+1;

% memory allocation
AD  =NaN(numPositions,localNumFrames);             % unknown value
D   =NaN(numPositions,localNumFrames);             % unknown value
MAD =NaN(numPositions,localNumFrames);             % unknown value
XMAD=NaN(numPositions,localNumFrames);             % unknown value
YMAD=NaN(numPositions,localNumFrames);             % unknown value
XYMAD=NaN(numPositions,localNumFrames);            % unknown value
XD =NaN(numPositions,localNumFrames);              % unknown value
YD = NaN(numPositions,localNumFrames);             % unknown value

XMaxD = NaN(numPositions,localNumFrames);          % unknown value
YMaxD = NaN(numPositions,localNumFrames);          % unknown value
XYMaxD = NaN(numPositions,localNumFrames);         % unknown value
XMinD = NaN(numPositions,localNumFrames);          % unknown value
YMinD = NaN(numPositions,localNumFrames);          % unknown value
XYMinD = NaN(numPositions,localNumFrames);         % unknown value

XYMean = NaN(numPositions,localNumFrames);         % unknown value
XMean = NaN(numPositions,localNumFrames);          % unknown value
YMean = NaN(numPositions,localNumFrames);          % unknown value

% first Frame
AD(:,1)=sqrt((X(:,startFrame)-X(:,startFrame)).^2+(Y(:,startFrame)-Y(:,startFrame)).^2);  %%% this is just zero!
D(:,1)=sqrt(X(:,startFrame).^2+Y(:,startFrame).^2);
XD(:,1) =X(:,startFrame);
YD(:,1) =Y(:,startFrame);
XMaxD(:,1) = XD(:,1);
XMinD(:,1) = XD(:,1);
YMaxD(:,1) = YD(:,1);
YMinD(:,1) = YD(:,1);
XMean(:,1) = 0.5*(XMaxD(:,1)+XMinD(:,1));
YMean(:,1) = 0.5*(YMaxD(:,1)+YMinD(:,1));
XYMaxD(:,1) = o(1)*XD(:,1)+o(2)*YD(:,1);
XYMinD(:,1) = o(1)*XD(:,1)+o(2)*YD(:,1);
XYMAD(:,1) = XYMaxD(:,1)-XYMinD(:,1);
XYMean(:,1) = 0.5*(XYMaxD(:,1)+XYMinD(:,1));

for n = 2:localNumFrames
globaln=n+startFrame-1;
fprintf('%s: Processing frame%02d\n',mfilename,n);
AD(:,n)=sqrt((X(:,globaln)-X(:,startFrame)).^2+(Y(:,globaln)-Y(:,startFrame)).^2);
D(:,n)=sqrt(X(:,globaln).^2+Y(:,globaln).^2);
MAD(:,n)=max(D(:,n),MAD(:,n-1));

XD(:,n)=X(:,globaln);
XMaxD(:,n)=max(XD(:,n),XMaxD(:,n-1));
XMinD(:,n)=min(XD(:,n),XMinD(:,n-1));
XMAD(:,n)=XMaxD(:,n)-XMinD(:,n);
XMean(:,n)= 0.5*(XMaxD(:,n)+XMinD(:,n));

YD(:,n)=Y(:,globaln);
YMaxD(:,n)=max(YD(:,n),YMaxD(:,n-1));
YMinD(:,n)=min(YD(:,n),YMinD(:,n-1));
YMAD(:,n)=YMaxD(:,n)-YMinD(:,n);
YMean(:,n)= 0.5*(YMaxD(:,n)+YMinD(:,n));

XYMaxD(:,n) = max(o(1)*XD(:,n)+o(2)*YD(:,n),XYMaxD(:,n-1));
XYMinD(:,n) = min(o(1)*XD(:,n)+o(2)*YD(:,n),XYMinD(:,n-1));
XYMAD=XYMaxD(:,n)-XYMinD(:,n);
XYMean(:,n) = 0.5*(XYMaxD(:,n)+XYMinD(:,n));

end
end