function [VEL,DIR,VLIM,DLIM] = velocity_triangulation(X,Y,T,TOL,varargin)
%VELOCITY_TRIANGULATION Triangulates 2-dimensional planar waves velocity.
%
% SYNTAX:
% VEL = velocity_triangulation(X,Y,T,TOL);
% VEL = velocity_triangulation(X,Y,T,TOL,...PN,PV,...);
% [VEL,DIR] = velocity_triangulation(...);
% [VEL,DIR,VLIM,DLIM] = velocity_triangulation(...);
%
% INPUT:
% X - Horizontal position [uL units] or longitude (if 'coord').
% Y - Vertical position [uL units] or latitude (if 'coord').
% T - Arrival time [uT units] at (X,Y) position.
% TOL - Tolerance time [uT units].
% PN/PV - Property name and value. One of:
% -----------------------------------------------------------
% NAME VALUES DEFINITION
% -----------------------------------------------------------
% 'Type' 'positions', Treates X,Y as LON,LAT
% or 'type' 'coordinates' coordinates. See NOTE
% DEFAULT: 'pos' below for details.
%
% 'Velocity' A vector or an Estimations of velocity
% or 'vel' integer. in uL/uT units. See NOTE
% DEFAULT: 100 below for details.
%
% 'Direction' A vector or an Estimations of direction
% or 'dir' integer. in degrees. See NOTE
% DEFAULT: 100 below for details.
%
% 'VLimits' true or false Estimates velocity
% or 'vlim' DEFAULT: true limits.
%
% 'DLimits' true or false Estimates direction
% or 'dlim' DEFAULT: true limits. True if 'vlim'
% is true.
%
% 'Figure' true or false Outputs a polar figure.
% or 'fig' DEFAULT: false See NOTE below for
% details.
% -----------------------------------------------------------
%
% OUTPUT:
% VEL - Wave velocity estimation in units of uL/uT, or meter/uT if
% 'coord' was used.
% DIR - Wave propagation angle or tracking if 'coord' was used, in
% degrees.
% VLIM - Estimated velocity limits.
% DLIM - Estimated direction limits.
%
% DESCRIPTION:
% By trial and error estimates the velocity magnitude and direction
% from a planar wave by its arrival time and positions of at least
% 2 sensors.
%
% Because it deals with planar waves, its phase is a line that should
% touch the first point before the others. This is used to estimate
% direction limits of this phase. Besides, the wave speed is faster
% when its direction its right from one sensor to another, and decrease
% down to zero when not. This is also used to estimate velocity limits.
%
% NOTE:
% * Optional outputs may or not be called.
% * In the output figure, a polar plot is used with the first arriving
% point at its center. The direction a velocity limits are drawn by
% green and red arrows, besides of its respective perpendicular wave
% phase at first point. If any velocity were found, it is plotted in
% blue. The region where those blue arrows should be is the one swept
% by the red arrow with increasing direction up to the green one.
% Remember that tracking increase clockwise. Velocities are
% normalized to fit the area.
% * A tolerance of 5% is added to the direction limits. This is
% observable in the figure.
% * By default, a vector of possible directions and velocities are
% calculated from:
% eDir = linspace(D1,D2,Nd);
% eVel = linspace(V1,V2,Nv);
% where:
% D1,D2 - are estimated from the extrema directions a planar wave
% could have after arriving the first point.
% V1,V2 - are estimated from the D1,D2 directions limits and the
% direction of the line joining the first arrival point to
% the other ones.
% Nd,Nv - By default are set as 100, but the user may give its own
% to increase the precision of the estimation.
% * All input directions must be between 0 and 360. The outputs will
% have this range also.
% * VEL and DIR outputs would be one of the given or estimated velocity
% and direction vectors.
% * When using 'coordinates', the program uses the M-MAP Mapping
% Toolbox by Rich Pawlowicz, which is found at
% http://www.eos.ubc.ca/~rich/
% Or, it uses the Mapping Toolbox. If any one is found, it results an
% error. Note, the user may define its own function. Look for
% lonlat2distrk handle function inside this file.
% * If more than one velocity was estimated, the mean velocity speed
% and direction may be better estimated by using my PMEAN function
% from the MATLAB FileExchange.
%
% EXAMPLE:
% % a) Wave traveling at 45:
% X = [0 1 1]; % m
% Y = [0 0 1]; % m
% T = [0 1 2]; % s
% TOL = 0.01; % s
% [VEL,DIR] = velocity_triangulation(X,Y,T,TOL,'fig',true);
% mean(VEL),mean(DIR)
% % RESULT: 0.71 m/s and 44.6
% % b) Wave traveling at 90:
% T = [0 0 1];
% TOL = 0.01;
% [VEL,DIR] = velocity_triangulation(X,Y,T,TOL,'Dir',500);
% mean(VEL),mean(DIR)
% % RESULT: 1.00 m/s and 90.0
%
% SEE ALSO:
% M_IDIST on M_MAP Toolbox by Rich Pawlowicz
% at http://www.eos.ubc.ca/~rich/
% and
% PMEAN by Carlos Vargas
% at http://www.mathworks.com/matlabcentral/fileexchange
%
%
% ---
% MFILE: velocity_triangulation.m
% VERSION: 1.0 (Sep 30, 2009) (<a href="matlab:web('http://www.mathworks.com/matlabcentral/fileexchange/authors/11258')">download</a>)
% MATLAB: 7.7.0.471 (R2008b)
% AUTHOR: Carlos Adrian Vargas Aguilera (MEXICO)
% CONTACT: nubeobscura@hotmail.com
% REVISIONS:
% 1.0 Released. (Sep 30, 2009)
% DISCLAIMER:
% velocity_triangulation.m is provided "as is" without warranty of any
% kind, under the revised BSD license.
% Copyright (c) 2009 Carlos Adrian Vargas Aguilera
% INPUTS CHECK-IN
% -------------------------------------------------------------------------
% Default.
eVel = [];
eDir = [];
Nv = 100;
Nd = 100;
vlim = true;
dlim = true;
coord = false;
draws = false;
VEL = [];
DIR = [];
VLIM = [NaN NaN];
DLIM = [NaN NaN];
% Parameter.
dirTol = 5; % Direction tolerance in percentaje of the
% difference between direction limits: diff(DLIM).
ellipsoid = 'wgs84'; % Datum system for the coordinates.
% Checks number of inputs and outputs.
if nargin<4
error('CVARGAS:velocity_triangulation:notEnoughInputs',...
'At least 4 inputs are required.')
elseif nargin>16
error('CVARGAS:velocity_triangulation:tooManyInputs',...
'At most 16 inputs are allowed.')
elseif nargout>4
error('CVARGAS:velocity_triangulation:tooManyOutputs',...
'At most 4 outputs are allowed.')
end
% Checks size.
Np = length(X);
if (Np<2) || (Np~=numel(X))
error('CVARGAS:velocity_triangulation:incorrectVectorInput',...
'Inputs must be vectors.')
end
if (Np~=numel(Y)) || (Np~=numel(T))
error('CVARGAS:velocity_triangulation:incorrectVectorInputSize',...
'Vector length must be equal.')
end
if numel(TOL)~=1
error('CVARGAS:velocity_triangulation:incorrectTolInput',...
'TOL must be an scalar.')
end
% Parses properties.
[eVel,eDir,Nv,Nd,vlim,dlim,coord,draws] = ...
parseProperties(eVel,eDir,Nv,Nd,vlim,dlim,coord,draws,varargin{:});
% Sets function to convert coordinates to distance.
if coord
if exist('m_idist','file')==2
lonlat2distrk = ...
@(lon1,lat1,lon2,lat2) m_idist(lon1,lat1,lon2,lat2,ellipsoid);
elseif exist('distance','file')==2
lonlat2distrk = ...
@(lon1,lat1,lon2,lat2) distance(lat1,lon1,lat2,lon2,...
almanac('earth',ellipsoid,'meters'));
else
error('CVARGAS:velocity_triangulation:notFoundMmapOrMappingToolbox',...
'M_MAP or Mapping Toolbox were not found.')
end
end
% -------------------------------------------------------------------------
% MAIN
% -------------------------------------------------------------------------
% Forces column vectors.
X = X(:);
Y = Y(:);
T = T(:);
% Forces positive TOL.
TOL = abs(TOL);
% Forces 'dlim'
if vlim
dlim = true;
end
% Sorts time.
[T,ind] = sort(T);
X = X(ind);
Y = Y(ind);
T = T-T(1);
% =========================================================================
% Distance, time lapses and angles between points.
Npairs = Np*(Np-1)/2;
dis = zeros(Npairs,1);
ang = dis;
tim = dis;
cont = 0;
if ~coord
for i = 1:Np
for j = i+1:Np
cont = cont+1;
tim(cont) = T(j)-T(i);
dis(cont) = sqrt((X(j)-X(i)).^2 + (Y(j)-Y(i)).^2);
ang(cont) = atan2((Y(j)-Y(i)),(X(j)-X(i)))*180/pi;
end
end
else
for i = 1:Np
for j = i+1:Np
cont = cont+1;
tim(cont) = T(j)-T(i);
[dis(cont),ang(cont)] = lonlat2distrk(X(i),Y(i),X(j),Y(j));
end
end
end
% =========================================================================
% EXTREME DIRECTIONS
if dlim
% Gets reference angle to convert all angles between points in the range
% (-90,90], using as an starting guess the angle between the first and
% last point. Then, corrects reference angle if any point is outside the
% required range.
dirRef = ang(Np-1);
range = 1:Np-1;
temp = ang(range)-dirRef;
ind2 = (temp)>90;
ind1 = (temp)<-90;
if any(ind1) && any(ind2)
error('CVARGAS:velocity_triangulation:incorrectTimeArrival',...
'Incorrect arrival times for planar waves.')
end
if any(ind1)
dirRef = mod(dirRef - min(temp(ind1)+90),360);
elseif any(ind2)
dirRef = mod(dirRef + max(temp(ind2)-90),360);
end
% Find direction limits.
D1 = DLIM(1);
D2 = DLIM(2);
for k = 1:Np-1
r = range(1:Np-k) + (k-1)*Np - sum(0:k-1);
% Converts all angles between points to (-90,90] range.
dir180 = ang(r)-dirRef;
if any(abs(dir180)>180)
error('CVARGAS:velocity_triangulation:incorrectTimeArrival',...
'Incorrect arrival times for planar waves.')
end
% Gets limits.
D1 = min(D1,min(dir180));
D2 = max(D2,max(dir180));
end
% Now converts those phase wave limits to velocity direction limits, wich
% are penpendicular.
D1 = D1+90;
D2 = D2-90;
% Gets extreme directions of the wave.
temp = max([D1 D2]);
D1 = min([D1 D2]);
D2 = temp;
% Adds some error to this extremes.
D1 = D1-(D2-D1)*(dirTol/100);
D2 = D2+(D2-D1)*(dirTol/100);
% Remember that both D1,D2 are in (-90,90] range, referenced to dirRef.
DLIM = mod([D1 D2]+dirRef,360);
end
% Checks direction vector.
if isempty(eDir)
% Generates directions vector.
if ~dlim
error('CVARGAS:velocity_triangulation:notGivenDir',...
'When ''DLim'' is false, ''Dir'' vector must be provided.')
end
% a) Linear space between limits directions.
eDir = linspace(D1,D2,Nd)';
% b) Changes eDir reference from dirRef to 0.
eDir = eDir+dirRef;
% c) Sets eDir range to (0,360].
eDir = mod(eDir,360);
% d) Saves direction limits referenced to dirRef.
elseif dlim
% Discriminates directions.
% a) Changes eDir reference from 0 to dirRef.
temp = eDir-dirRef;
% b) Sets eDir range to (-180,180].
temp = temp-360*(eDir>=180);
% c) Eliminates directions outside direction limits.
bad = (temp<=D1) | (temp>=D2);
eDir(bad) = [];
end
% =========================================================================
%EXTREME VELOCITIES
if vlim
V1 = VLIM(1);
V2 = VLIM(2);
for k = 1:Np-1
% a) Gets angles in (-90,90] range.
r = range(1:Np-k) + (k-1)*Np - sum(0:k-1);
dir180 = ang(r)-dirRef;
% b) Gets factor of distance traveled by the phase wave from first point to
% the others.
f1 = abs(cos((D1-dir180)*pi/180));
f2 = abs(cos((D2-dir180)*pi/180));
% c) Estimates velocity from first point to the others.
w1 = dis(r)./(tim(r)+TOL); % Velocity direct to the points.
w2 = dis(r)./(tim(r)-TOL); % Velocity direct to the points.
bad = (tim(r)<=TOL); % Gets points in same phase.
w1(bad) = Inf;
w2(bad) = Inf;
u11 = w1.*f1; % Velocity to the direction D1
u12 = w1.*f2; % Velocity to the direction D2
u21 = w2.*f1; % Velocity to the direction D1
u22 = w2.*f2; % Velocity to the direction D2
% d) Eliminates direct directions to other points when the point does not
% fall between direction limits.
bad = (dir180<D1) | (dir180>D2);
w1(bad) = NaN;
w2(bad) = NaN;
% e) Eliminates points in the same phase as the first one.
w1 (~isfinite(w1 )) = NaN;
w2 (~isfinite(w2 )) = NaN;
u11(~isfinite(u11)) = NaN;
u12(~isfinite(u12)) = NaN;
u21(~isfinite(u21)) = NaN;
u22(~isfinite(u22)) = NaN;
% f) Gets limits of this velocities, from 0 to Inf.
v1 = min([w1(:) u11(:) u12(:)],[],2);
v2 = max([w2(:) u21(:) u22(:)],[],2);
v1(~isfinite(v1)) = 0;
v2(~isfinite(v2)) = +Inf;
% g) Gets velocity magnitude limits.
v1 = max(v1);
v2 = min(v2) ;
% Saves new limits.
V1 = max(V1,v1);
V2 = min(V2,v2);
end
VLIM = [V1 V2];
end
% Velocity vector.
if isempty(eVel)
% Generates velocity vector.
if ~vlim
error('CVARGAS:velocity_triangulation:notGivenVel',...
'When ''VLim'' is false, ''Vel'' vector must be provided.')
end
if any(~isfinite([V1 V2]))
error('CVARGAS:velocity_triangulation:unableToGetVelLimits',...
['Unable to estimate correct velocity limits. '...
'Please provide velocity estimation.'])
end
eVel = linspace(V1,V2,Nv)';
elseif vlim
% Discriminates velocities.
eVel(eVel<V1) = [];
eVel(eVel>V2) = [];
end
% =========================================================================
% Starts estimations by comparing the time of arrival.
kcontinue = false;
for k = 1:length(eVel) % MAIN LOOP
% Checks arrival time.
ind = abs((dis(1)*cos((ang(1)-eDir)*pi/180)/eVel(k))-tim(1))<=TOL;
if ~any(ind)
continue
end
ind = find(ind);
% Following estimations.
for l = 2:Npairs
% Checks arrival time.
ind(~(abs((dis(l)*cos(abs(ang(l)-eDir(ind))*pi/180)/eVel(k))-tim(l))<=TOL)) = [];
if isempty(ind)
kcontinue = true;
break
end
end
if kcontinue
kcontinue = false;
continue
end
% Saves good estimations.
M = length(ind);
VEL(end+(1:M),1) = repmat(eVel(k),M,1);
DIR(end+(1:M),1) = eDir(ind);
end
% OUTPUTS CHECK-OUT
% -------------------------------------------------------------------------
% Draws.
if draws
if ~dlim
error('CVARGAS:velocity_triangulation:unableToDraw',...
'When ''figure'' is true, ''DLim'' cannot be false.')
end
disRef = max(dis(1:Np-1));
velMax = disRef/max(VEL);
color = flipud(jet(Np));
if ~coord
% POSITIONS.
polar(NaN,disRef)
hold on
h1 = polar((D1+dirRef + 90*[-1 1])*pi/180,disRef*[1 1],'r');
h2 = polar((D2+dirRef + 90*[-1 1])*pi/180,disRef*[1 1],'g');
if vlim
h3 = quiver([0 0],[0 0],VLIM*velMax*cos((D1+dirRef)*pi/180),...
VLIM*velMax*sin((D1+dirRef)*pi/180),0,'r');
h4 = quiver([0 0],[0 0],VLIM*velMax*cos((D2+dirRef)*pi/180),...
VLIM*velMax*sin((D2+dirRef)*pi/180),0,'g');
else
h3 = quiver(0,0,disRef*cos((D1+dirRef)*pi/180),...
disRef*sin((D1+dirRef)*pi/180),0,'r');
h4 = quiver(0,0,disRef*cos((D2+dirRef)*pi/180),...
disRef*sin((D2+dirRef)*pi/180),0,'g');
end
h5 = quiver(zeros(size(DIR)),zeros(size(DIR)),...
(VEL*velMax).*cos(DIR*pi/180),...
(VEL*velMax).*sin(DIR*pi/180),0,'b');
h6 = polar(0,0,'ko');
set(h6,'MarkerFaceColor',color(1,:))
for k = 1:Np-1
h6 = polar(ang(k)*pi/180,dis(k),'ko');
set(h6,'MarkerFaceColor',color(k+1,:))
end
set([h1;h2;h3;h4;h5],'LineWidth',2)
hold off
if ~isempty(DIR)
legend([h1(1);h2(1);h5(1);h6(1)],...
'Extrema wave','Extrema wave','Velocity','Points',...
'Location','EastOutside')
else
legend([h1(1);h2(1);h6(1)],...
'Extrema wave','Extrema wave','Points',...
'Location','EastOutside')
end
else
% COORDINATES.
polar(NaN,disRef)
hold on
h1 = polar((90-(D1+dirRef + 90*[-1 1]))*pi/180,disRef*[1 1],'r');
h2 = polar((90-(D2+dirRef + 90*[-1 1]))*pi/180,disRef*[1 1],'g');
if vlim
h3 = quiver([0 0],[0 0],VLIM*velMax*cos((90-(D1+dirRef))*pi/180),...
VLIM*velMax*sin((90-(D1+dirRef))*pi/180),0,'r');
h4 = quiver([0 0],[0 0],VLIM*velMax*cos((90-(D2+dirRef))*pi/180),...
VLIM*velMax*sin((90-(D2+dirRef))*pi/180),0,'g');
else
h3 = quiver(0,0,disRef*cos((90-(D1+dirRef))*pi/180),...
disRef*sin((90-(D1+dirRef))*pi/180),0,'r');
h4 = quiver(0,0,disRef*cos((90-(D2+dirRef))*pi/180),...
disRef*sin((90-(D2+dirRef))*pi/180),0,'g');
end
h5 = quiver(zeros(size(DIR)),zeros(size(DIR)),...
(VEL*velMax).*cos((90-DIR)*pi/180),...
(VEL*velMax).*sin((90-DIR)*pi/180),0,'b');
h6 = polar(0,0,'ko');
set(h6,'MarkerFaceColor',color(1,:))
for k = 1:Np-1
h6 = polar((90-ang(k))*pi/180,dis(k),'ko');
set(h6,'MarkerFaceColor',color(k+1,:))
end
set([h1;h2;h3;h4;h5],'LineWidth',2)
hold off
if ~isempty(DIR)
legend([h1(1);h2(1);h5(1);h6(1)],...
'Extrema wave','Extrema wave','Velocity','Coordinates',...
'Location','EastOutside')
else
legend([h1(1);h2(1);h6(1)],...
'Extrema wave','Extrema wave','Coordinates',...
'Location','EastOutside')
end
title North
end
colormap(color)
caxis([0 1])
colorbar('YTick',[0 1],'YTickLabel',['First';'Last ']);
end
% =========================================================================
% SUBFUNCTIONS
% -------------------------------------------------------------------------
function [eVel,eDir,Nv,Nd,vlim,dlim,coord,draws] = ...
parseProperties(eVel,eDir,Nv,Nd,vlim,dlim,coord,draws,varargin)
% Parses properties.
if rem(nargin-8,2)~=0
error('CVARGAS:velocity_triangulation:incorrectPairedInputLength',...
'Extra inputs must be paired.')
end
while ~isempty(varargin)
if isempty(varargin{1}) || ~ischar(varargin{1}) || isempty(varargin{2})
error('CVARGAS:velocity_triangulation:incorrectPairedInput',...
'Incorrect paired input.')
end
switch lower(varargin{1})
case {'type'}
coord = varargin{2};
if strncmpi(coord,'positions' ,min(length(coord),9 ))
coord = false;
elseif strncmpi(coord,'coordinates',min(length(coord),11))
coord = true;
else
error('CVARGAS:velocity_triangulation:incorrectTypeInput',...
'"Type" must be one of ''positions'' or ''coordinates''.')
end
case {'velocity','velocit','veloci','veloc','velo','vel'}
if numel(varargin{2})==1
Nv = round(varargin{2}); % Forces integer.
if Nv<=0
error('CVARGAS:velocity_triangulation:incorrectVelIntegerInput',...
'Single "Velocity" value must be a positive integer.')
end
else
eVel = varargin{2}(:);
if any(eVel)<0
error('CVARGAS:velocity_triangulation:incorrectVelVectorInput',...
'Vector "Velocity" must have positive elements.')
end
end
case {'direction','directio','directi','direct','direc','dire','dir'}
if numel(varargin{2})==1
Nd = round(varargin{2}); % Forces integer.
if Nd<=0
error('CVARGAS:velocity_triangulation:incorrectDirIntegerInput',...
'Single "Direction" value must be a positive integer.')
end
else
eDir = varargin{2}(:);
if any(eDir)<0 || any(eDir)>360
warning('CVARGAS:velocity_triangulation:incorrectDirVectorInput',...
['Vector "Direction" has elements outside (0,360] degrees.'...
' Used mod(x,360).'])
end
eDir = mod(eDir,360);
end
case {'dlimits','dlimit','dlimi','dlim'}
dlim = logical(varargin{2}); % Forces logical.
case {'vlimits','vlimit','vlimi','vlim'}
vlim = logical(varargin{2}); % Forces logical.
case {'figure','figur','figu','fig'}
draws = logical(varargin{2}); % Forces logical.
otherwise
error('CVARGAS:velocity_triangulation:unrecognizedProperty',...
['Unrecognized ' varargin{1} ' property.'])
end
varargin(1:2) = [];
end
% [EOF] velocity_triangulation.m
