Finding axial, radial velocity and core radius of a 3D vortex from PTV data

I have obtained PTV data and would like to calculate various components of the vortex, is there a better code than the one i have used below.Not certain about the accuracy of the code.
% select a CSV file to import
[file, path] = uigetfile(‘*.csv’, ‘Select the CSV file to import’);
if file == 0
return; % user canceled file selection
end
%Load the data from the CSV file
data = readmatrix(fullfile(path, file), ‘HeaderLines’, 1);
% x,y coordinate vectors (mm) and Vx and Vy fields (m/s)
x = data(:, 1);
y = data(:, 2);
z = data(:, 3);
u = data(:, 4);
v = data(:, 5);
w = data(:, 6);
vorticity =data(:,20);
% Example data (replace these with your actual data)
% Example data (replace these with your actual data)
y = linspace(-5, 5, 100) / 1000; % y-coordinates (converted to meters)
z = linspace(-5, 5, 100) / 1000; % z-coordinates (converted to meters)
x = linspace(-5, 5, 100) / 1000; % x-coordinates (converted to meters)
[Y, Z, X] = meshgrid(y, z, x); % Create meshgrid for y, z, x

% Step 1: Calculate the vorticity magnitude
vorticity_magnitude=max(abs(vorticity(:)));

% Step 2: Define the threshold for vortex core (e.g., 10% of max vorticity)
max_vorticity = max(vorticity_magnitude(:));
threshold = 0.1 * max_vorticity; % Threshold for vortex core (10% of max vo

% Step 3: Find the vortex core radius (where vorticity exceeds threshold)
% Compute the radial distance from the center (0, 0, 0)
radius = sqrt(X.^2 + Y.^2 + Z.^2); % Radial distance in 3D

% Identify the region where vorticity exceeds threshold
core_radius = max(radius(vorticity_magnitude >= threshold)); % Vortex core radius

% Step 4: Find the location of maximum vorticity (core of the vortex)
[~, core_idx] = max(vorticity_magnitude(:)); % Find the index of max vorticity
[core_x, core_y, core_z] = ind2sub(size(vorticity_magnitude), core_idx); % Get the coordinates of the core

% Step 5: Radial and Axial Velocities at the vortex core
% The radial velocity is the component of velocity in the direction of the radial vector
% The axial velocity is the component of velocity along the axis of the vortex (e.g., the x-axis)

% Calculate the radial velocity at the core
r_vec = [X(core_x, core_y, core_z), Y(core_x, core_y, core_z), Z(core_x, core_y, core_z)]; % Position vector at core
r_mag = norm(r_vec); % Magnitude of position vector

% Calculate the radial velocity at core
radial_velocity = (u(core_x, core_y, core_z) * r_vec(1) + …
v(core_x, core_y, core_z) * r_vec(2) + …
w(core_x, core_y, core_z) * r_vec(3)) / r_mag; % Dot product of velocity and radial direction

% The axial velocity is the velocity component along the x-axis (assuming the vortex is aligned with the x-axis)
axial_velocity = u(core_x, core_y, core_z); % Axial velocity along x-axis at the core

% Step 6: Find the velocity at the outer boundary (where radius is maximum)
% Find the index at the outermost boundary (max radial distance)
[~, boundary_idx] = min(abs(radius(:) – max(radius(:)))); % Find the outer boundary index
[boundary_x, boundary_y, boundary_z] = ind2sub(size(radius), boundary_idx); % Get boundary coordinates

boundary_velocity = sqrt(u(boundary_x, boundary_y, boundary_z)^2 + v(boundary_x, boundary_y, boundary_z)^2 + w(boundary_x, boundary_y, boundary_z)^2);

% Display Results
fprintf(‘Vortex core radius: %.4f metersn’, core_radius);
fprintf(‘Radial velocity at the core: %.4f m/sn’, radial_velocity);
fprintf(‘Axial velocity at the core: %.4f m/sn’, axial_velocity);

%I have obtained PTV data and would like to calculate various components of the vortex, is there a better code than the one i have used below.Not certain about the accuracy of the code.
% select a CSV file to import
[file, path] = uigetfile(‘*.csv’, ‘Select the CSV file to import’);
if file == 0
return; % user canceled file selection
end
%Load the data from the CSV file
data = readmatrix(fullfile(path, file), ‘HeaderLines’, 1);
% x,y coordinate vectors (mm) and Vx and Vy fields (m/s)
x = data(:, 1);
y = data(:, 2);
z = data(:, 3);
u = data(:, 4);
v = data(:, 5);
w = data(:, 6);
vorticity =data(:,20);
% Example data (replace these with your actual data)
% Example data (replace these with your actual data)
y = linspace(-5, 5, 100) / 1000; % y-coordinates (converted to meters)
z = linspace(-5, 5, 100) / 1000; % z-coordinates (converted to meters)
x = linspace(-5, 5, 100) / 1000; % x-coordinates (converted to meters)
[Y, Z, X] = meshgrid(y, z, x); % Create meshgrid for y, z, x

% Step 1: Calculate the vorticity magnitude
vorticity_magnitude=max(abs(vorticity(:)));

% Step 2: Define the threshold for vortex core (e.g., 10% of max vorticity)
max_vorticity = max(vorticity_magnitude(:));
threshold = 0.1 * max_vorticity; % Threshold for vortex core (10% of max vo

% Step 3: Find the vortex core radius (where vorticity exceeds threshold)
% Compute the radial distance from the center (0, 0, 0)
radius = sqrt(X.^2 + Y.^2 + Z.^2); % Radial distance in 3D

% Identify the region where vorticity exceeds threshold
core_radius = max(radius(vorticity_magnitude >= threshold)); % Vortex core radius

% The axial velocity is the velocity component along the x-axis (assuming the vortex is aligned with the x-axis)
axial_velocity = u(core_x, core_y, core_z); % Axial velocity along x-axis at the core

boundary_velocity = sqrt(u(boundary_x, boundary_y, boundary_z)^2 + v(boundary_x, boundary_y, boundary_z)^2 + w(boundary_x, boundary_y, boundary_z)^2);

% I have obtained PTV data and would like to calculate various components of the vortex, is there a better code than the one i have used below.Not certain about the accuracy of the code.
% select a CSV file to import
[file, path] = uigetfile(‘*.csv’, ‘Select the CSV file to import’);
if file == 0
return; % user canceled file selection
end
%Load the data from the CSV file
data = readmatrix(fullfile(path, file), ‘HeaderLines’, 1);
% x,y coordinate vectors (mm) and Vx and Vy fields (m/s)
x = data(:, 1);
y = data(:, 2);
z = data(:, 3);
u = data(:, 4);
v = data(:, 5);
w = data(:, 6);
vorticity =data(:,20);
% Example data (replace these with your actual data)
% Example data (replace these with your actual data)
y = linspace(-5, 5, 100) / 1000; % y-coordinates (converted to meters)
z = linspace(-5, 5, 100) / 1000; % z-coordinates (converted to meters)
x = linspace(-5, 5, 100) / 1000; % x-coordinates (converted to meters)
[Y, Z, X] = meshgrid(y, z, x); % Create meshgrid for y, z, x

% Step 1: Calculate the vorticity magnitude
vorticity_magnitude=max(abs(vorticity(:)));

% Step 2: Define the threshold for vortex core (e.g., 10% of max vorticity)
max_vorticity = max(vorticity_magnitude(:));
threshold = 0.1 * max_vorticity; % Threshold for vortex core (10% of max vo

% Step 3: Find the vortex core radius (where vorticity exceeds threshold)
% Compute the radial distance from the center (0, 0, 0)
radius = sqrt(X.^2 + Y.^2 + Z.^2); % Radial distance in 3D

% Identify the region where vorticity exceeds threshold
core_radius = max(radius(vorticity_magnitude >= threshold)); % Vortex core radius

% The axial velocity is the velocity component along the x-axis (assuming the vortex is aligned with the x-axis)
axial_velocity = u(core_x, core_y, core_z); % Axial velocity along x-axis at the core

boundary_velocity = sqrt(u(boundary_x, boundary_y, boundary_z)^2 + v(boundary_x, boundary_y, boundary_z)^2 + w(boundary_x, boundary_y, boundary_z)^2);

% 3d vortex, ptv, matlab MATLAB Answers — New Questions

Cart

Cart