Commit ddc4b0da authored by Jean-Luc Shaw's avatar Jean-Luc Shaw

Added the ctrspd.m routine with calculates the speed at which a line in 2D is moving.

parent fa2630c0
function m_circle(lon,lat,r,label)
if isempty(label)
rng = [0:360 0] ;
else
rng = [290:360 0:250] ;
m_text(lon,lat-r,label,'horizontalalignment','center','fontsize',6) ;
end
x = lon + r*cosd(rng)/cosd(lat) ;
y = lat + r*sind(rng) ;
m_plot(x,y,'k-') ;
end
function [u,v] = ctrspd(x1,y1,x2,y2,t)
% Synthax : [u,v] = ctrspd(x1,y1,x2,y2,t)
%
% Calculates the speed at which a contour travels perpendicularly to itself.
% The instantaneous speed of a vertex is calculated as the distance between
% itself at time t_1, and the intersection of the contour at time t_2 with a
% line normal to contour 1 and crossing the vertex.
%
% Inputs are coordinates at t_1 : x1,y1
% coordinates at t_2 : x2,y2
% times 1 & 2 : t
%
% Where t is size two vector. Contours may be periodic or non periodic. The
% beginning and ends of the contour will be handled differently in these
% cases.
% PARAMETERS
NA = numel(x1) ;
NB = numel(x2) ;
per1 = false;
per2 = false;
warning('off','MATLAB:singularMatrix') ;
warning('off','MATLAB:polyfit:RepeatedPointsOrRescale') ;
warning('off','MATLAB:nearlySingularMatrix') ;
% MAKE CONTOURS PERIODIC according to prior periodicity
if x1(1) == x1(end) & y1(1) == y1(end) & x2(1) == x2(end) & y2(1) == y2(end)
x1 = [x1(end-1) x1(1:end-1) x1(1)] ;
y1 = [y1(end-1) y1(1:end-1) y1(1)] ;
x2 = [x2(end-1) x2(1:end-1) x2(1)] ;
y2 = [y2(end-1) y2(1:end-1) y2(1)] ;
per1 = true;
per2 = true;
elseif x1(1) == x1(end) & y1(1) == y1(end) & x2(1) ~= x2(end) & y2(1) ~= y2(end)
x1 = [x1(end-1) x1(1:end-1) x1(1)] ;
y1 = [y1(end-1) y1(1:end-1) y1(1)] ;
x2 = [x2(end) x2 x2(1)] ;
y2 = [y2(end) y2 y2(1)] ;
per1 = true;
elseif x1(1) ~= x1(end) & y1(1) ~= y1(end) & x2(1) ~= x2(end) & y2(1) ~= y2(end)
x1 = [x1(end) x1 x1(1)] ;
y1 = [y1(end) y1 y1(1)] ;
x2 = [x2(end-1) x2(1:end-1) x2(1)] ;
y2 = [y2(end-1) y2(1:end-1) y2(1)] ;
per2 = true;
else
x1 = [x1(end) x1 x1(1)] ;
y1 = [y1(end) y1 y1(1)] ;
x2 = [x2(end) x2 x2(1)] ;
y2 = [y2(end) y2 y2(1)] ;
end
% OUTPUT VARIABLES
u = nan(size(x1)) ;
v = nan(size(x1)) ;
% LOOP OVER CONTOUR VERTICES
for ii = 2:numel(x1) - 1
gam = range_pm180_2_360( atan2d(y1(ii+1)-y1(ii-1),x1(ii+1)-x1(ii-1))) ;
if ~per1 & ii == 2
gam = range_pm180_2_360( atan2d(y1(ii+1)-y1(ii),x1(ii+1)-x1(ii))) ;
elseif ~per1 & ii == numel(x1) - 1
gam = range_pm180_2_360( atan2d(y1(ii)-y1(ii-1),x1(ii)-x1(ii-1))) ;
end
bet = gam - 90 ;
if bet < 0 ; bet+360 ; end
xp = x1(ii) + cosd(bet) ;
yp = y1(ii) + sind(bet) ;
% Make the parametric equation for this line
m1 = polyfit([x1(ii) xp],[y1(ii) yp],1) ;
vec = @(x)m1(1)*x + m1(2) ;
% Draw for the next contours x values
yvec = vec(x2) ;
% Find intersections
I = false(NB,1) ;
for jj = 2:numel(x2) - 1
if y2(jj) > yvec(jj) & y2(jj-1) < yvec(jj-1) | y2(jj) < yvec(jj) & y2(jj-1) > yvec(jj-1)
I(jj) = true ;
end
end
I = find(I) ;
% Narrow it down to the nearest intersection
d = nan(size(I)) ;
for jj = 1:numel(I)
d(jj) = sqrt( (x1(ii) - x2(I(jj))).^2 + (y1(ii) - y2(I(jj))).^2 ) ;
end
[~,Imin] = min(d) ;
Imin = I(Imin) ;
% Fit this line segment
m2 = polyfit([x2(Imin) x2(Imin-1)],[y2(Imin) y2(Imin-1)],1) ;
% Solve the linear system of these two lines
A = [ m1(1) -1; m2(1) -1];
B = [-m1(2); -m2(2)];
X = A\B ;
% Calculate angle and magnitude of velocity
dist = sqrt( (x1(ii) - X(1)).^2 + (y1(ii) - X(2)).^2 ) ;
the = range_pm180_2_360( atan2d(X(2)-y1(ii),X(1)-x1(ii))) ;
% Decompose into components
spd = dist/(t(2) - t(1)) ;
u(ii) = spd*cosd(the) ;
v(ii) = spd*sind(the) ;
end
% DROP PERIODICITY
u([1 end]) = [];
v([1 end]) = [];
x1([1 end]) = [];
y1([1 end]) = [];
% MIRROR INPUT PERIODICITY
if per1 ; u = app(u,u(1)) ; v = app(v,v(1)) ; end
function y = app(x,v)
if iscolumn(x)
y = [x;v] ;
else
y = [x v] ;
end
end
end
function [A,HS] = genopt(Q,C,dom,NP,NI,pc,pm)
% Synthax: [A,HS] = genopt(@Q,@X,dom,NP,NI,pc,pm)
function [A,S] = genopt(Q,C,dom,NP,NI,pc,pm)
% Synthax: [A,S] = genopt(@Q,@X,dom,NP,NI,pc,pm)
%
% Optimisation routine following a genetic algorithm. Returns the best parameter
% set which solves the problem expressed by the anonymous function 'Q', when
......@@ -61,9 +61,9 @@ for kk = 1:NI
end
end
end
disp(['Last iteration: ' num2str(kk)]) ;
%disp(['Last iteration: ' num2str(kk)]) ;
% OUTPUT
A = MOM ;
HS = HS(1) ;
S = HS(1) ;
end
......@@ -56,9 +56,9 @@ for ii = 1:2:nargin-5
case 'section'
switch varargin{ii+1}
case 'top'
W = zc/z;
W = zc./z;
case 'bottom'
W = z/zc;
W = z./zc;
case 'all'
W = ones(size(z));
end
......
function points2grid(x,y,m_map,varargin)
function [varargout] = points2grid(x,y,m_map,varargin)
% Synthax : points2grid(x,y,m_map,varargin)
%
% Takes as input a grid of coordinates contained in 'x' and 'y' with the
......@@ -16,6 +16,10 @@ function points2grid(x,y,m_map,varargin)
% plotted grid itself.
%
% Options
% default
%
% input size
[M,N] = size(x) ;
......@@ -36,6 +40,10 @@ INTERPy = scatteredInterpolant(2*I,2*J,yvec(:), ...
X = INTERPx(indexx',indexy') ;
Y = INTERPy(indexx',indexy') ;
% Save full grid for output
FX = X ;
FY = Y ;
% Replace initial nan values
[I,J] = find(isnan(x) | isnan(y)) ;
for ii = 1:numel(I)
......@@ -184,4 +192,18 @@ for jj = 1:2:2*N+1
plot(X(:,jj),Y(:,jj),varargin{:})
end
end
% Manage output
if nargout==2
varargout{1} = X ;
varargout{2} = Y ;
elseif nargout==4
varargout{1} = X ;
varargout{2} = Y ;
varargout{3} = FX ;
varargout{4} = FY ;
elseif nargout==0
else
varargout = [] ;
end
end
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment