How can I locate a vortex center?

31 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!