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));
fprintf('Speed: u = %f , v = %f\n', ud(ind), vd(ind));
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:
