How can I locate a vortex center?

조회 수: 21(최근 30일)
Ashfaq Ahmed
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:
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 :)

채택된 답변

Riccardo Scorretti
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:
% 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
Ashfaq Ahmed
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?

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

추가 답변(0개)

범주

Find more on Vector Fields in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by