How can I locate a vortex center?

58 views (last 30 days)
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:
load('vortex_velocity_data.mat')
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 :)

Accepted Answer

Riccardo Scorretti
Riccardo Scorretti on 25 Apr 2022
Edited: Riccardo Scorretti on 25 Apr 2022
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:
% Load the data
load('vortex_velocity_data.mat');
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:
% Load the data
load('vortex_velocity_data.mat');
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 Comments
Ashfaq Ahmed
Ashfaq Ahmed on 9 May 2022
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?

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by