# How can I locate a vortex center?

조회 수: 32 (최근 30일)
Ashfaq Ahmed 2022년 4월 25일
댓글: Ashfaq Ahmed 2022년 5월 9일
Hi all!
Someone in the MATLAB community once asked similar kind of question but unfortuinately, it was never answered. So, I am asking it again here.
Suppose I have a vortex field. If you please download the "vortex_velocity_data.mat" file you will be able to create the vector field which I created.
Once the code is imported to the workspace, I typed the following code, to get the quiver plot:
figure(1),
q2 = quiver(x,y,ud,vd)
box on;
set(findall(gcf,'-property','LineWidth'),'LineWidth',3);
The image is attached -
How can I locate the center of this quiver/vector plot WITHOUT using the 'Data Cursor'? What lines of code will I have to write in order to locate the center? By locating the center x and y co-ordinates I can find the relative velocity components (u,v)?
Any feedback from you will be highly appreciated :)

댓글을 달려면 로그인하십시오.

### 채택된 답변

Riccardo Scorretti 2022년 4월 25일
편집: Riccardo Scorretti 2022년 4월 25일
Nice problem! I propose the solution based on the rationale that the center of the vortex (assuming that there is only one vortex) is the point where the curl is maximum:
figure
q2 = quiver(x, y, ud, vd);
box on;
% set(findall(gcf,'-property','LineWidth'),'LineWidth',3);
hold on
% Look for the max of the curl
[curlz, cav] = curl(x, y, ud, vd);
[~, ind] = max(abs(curlz(:)));
plot(x(ind), y(ind), 'r*');
fprintf('Center position: x0 = %f , y0 = %f\n', x(ind), y(ind));
Center position: x0 = 0.076262 , y0 = 0.084869
fprintf('Speed: u = %f , v = %f\n', ud(ind), vd(ind));
Speed: u = -0.388332 , v = -1.441110
In a more general case where your vector field is given by a couple of functions Fu(x,y) and Fv(x,y) you can try to solve the problem by dicothomy, by using more or less the same argument; only the computation of the curl has to be programmed differently:
figure
q2 = quiver(x, y, ud, vd);
box on;
% set(findall(gcf,'-property','LineWidth'),'LineWidth',3);
hold on
% Define interpolators
Fu = scatteredInterpolant(x(:), y(:), ud(:));
Fv = scatteredInterpolant(x(:), y(:), vd(:));
% Try a simple 2D bisection algorithm
minx = min(x(:)) ; maxx = max(x(:));
miny = min(y(:)) ; maxy = max(y(:));
maxit = 20;
for it = 1 : maxit
mx = (minx + maxx)/2;
my = (miny + maxy)/2;
set_of_circulations = [ ...
circulation(Fu,Fv,minx,mx,miny,my) ...
circulation(Fu,Fv,minx,mx,my,maxy) ...
circulation(Fu,Fv,mx,maxx,miny,my) ...
circulation(Fu,Fv,mx,maxx,my,maxy) ...
];
[~, ind] = max(abs(set_of_circulations));
if ind == 1 || ind == 2
maxx = mx;
else
minx = mx;
end
if ind == 1 || ind == 3
maxy = my;
else
miny = my;
end
end
plot(mx, my, 'r*');
return
function val = circulation(Fx, Fy, x1, x2, y1, y2)
% Compute the circulation of the vector field (Fx,Fy) along a square path
if x1 > x2 , t = x1 ; x1 = x2 ; x2 = t ; end
if y1 > y2 , t = y1 ; y1 = y2 ; y2 = t ; end
val = ...
integral(@(t) Fx(t,y1*ones(size(t))), x1, x2) + ...
integral(@(t) Fy(x2*ones(size(t)),t), y1, y2) + ...
integral(@(t) Fx(t,y2*ones(size(t))), x2, x1) + ...
integral(@(t) Fy(x1*ones(size(t)),t), y2, y1);
end
In both cases I got more or less the same figure:
##### 댓글 수: 6이전 댓글 4개 표시이전 댓글 4개 숨기기
Riccardo Scorretti 2022년 4월 28일
편집: Riccardo Scorretti 2022년 4월 28일
If you have several vortex the problem is much more messy. Basically, it boils up in findining multiple peaks in a noisy graphic.
Perhaps it can be solved if you have an a priori information on the minimal distance between vortex. Hereafter a temptative program, of which I'm not much happy. It tries to find a number of vortex in your data. It has three parameters (which are tricky to determine) and gives many false positives:
clear all;
close all; clc;
figure
q2 = quiver(x, y, vx, vy); hold on
%normalize the field before computing the curl
nuv=sqrt(vx.^2+vy.^2);
vx=vx./nuv;
vy=vy./nuv;
% Look for the max of the curl
[curlz, cav] = curl(x, y, vx, vy);
[~, ind] = max(abs(curlz(:)));
% plot(x(ind), y(ind), 'r*','MarkerSize',18);
box on;
hold on
fprintf('Center position: x0 = %f , y0 = %f\n', x(ind), y(ind));
Center position: x0 = -2010736.133410 , y0 = 454953.143337
% These parameters are likely to be "tricky" to set up
mindist = 10000; % = minimal distance between two vortex
threshold = 0.90; % = threshold (must be in the range ]0 ; 1[ )
maxnumvortex = 10; % = max number of vortex
% Hereafter it is assume that x and y have the structure of a grid
% There are for sure smarter ways to do
[M, N] = size(x);
ii = 2:M-1;
jj = 2:N-1;
curlz = abs(curlz);
thrval = sort(curlz(:));
% figure ; plot(thrval);
thrval = thrval(round(threshold*numel(thrval)));
% figure
% pcolor(x, y, curlz); hold on ; shading interp
ind = find( ...
curlz(ii,jj) > curlz(ii+1,jj) & ...
curlz(ii,jj) > curlz(ii-1,jj) & ...
curlz(ii,jj) > curlz(ii,jj+1) & ...
curlz(ii,jj) > curlz(ii,jj-1) & ...
curlz(ii,jj) > thrval ...
);
% Mandatory tricky step
[tmp_i, tmp_j] = ind2sub([M-2 N-2], ind);
ind = sub2ind([M N], tmp_i+1, tmp_j+1);
clear tmp_i tmp_j
fprintf('Found %i candidates\n', numel(ind));
Found 825 candidates
plot(x(ind), y(ind), 'c.','MarkerSize',5);
listOf_vortex = [];
while ~isempty(ind)
[~, t] = max(curlz(ind));
i_vortex = ind(t);
listOf_vortex(end+1) = i_vortex;
% Get rid of close candidates
x0 = x(i_vortex) ; y0 = y(i_vortex);
plot(x0, y0, 'm.','MarkerSize',5);
dst2 = (x(ind) - x0).^2 + (y(ind) - y0).^2;
ind = ind(dst2 > mindist.^2);
end
fprintf('Found vortex = %i\n', numel(listOf_vortex));
Found vortex = 226
% There are still too much vortex
[~, t] = sort(curlz(listOf_vortex), 'descend');
listOf_vortex = listOf_vortex(t(1:min(maxnumvortex, numel(listOf_vortex))));
plot(x(listOf_vortex), y(listOf_vortex), 'r*','MarkerSize',18);
axis(1.0E6*[-2.0535 -1.9620 0.3917 0.5379]);
Ashfaq Ahmed 2022년 5월 9일
Hi @Riccardo Scorretti! This code is great! It's actually detecting the centers of the vortex. But I have a question: what does the cyan color points and the magenda color points indicate? Do they have something to do with the identification of curls within a distance of 10 km?

댓글을 달려면 로그인하십시오.

### 카테고리

Help CenterFile Exchange에서 Graphics Objects에 대해 자세히 알아보기

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by