stream3 returns error: Sample points must be unique

Hello, I searched the community forum and tried my best to resolve the issue myself, but I could not do it.

Here is my code:

clc; close; clear all;
format compact
%% Cylindrical Coordinates
% Generate cylindrical coordinates
phi = [0 pi/2 pi 3*pi/2];
rho = [0.25 2];
z = [1 2];
% Create a meshgrid for theta, r, and z
[Phi, Rho, Z] = meshgrid(phi, rho, z);
% Convert cylindrical coordinates to Cartesian coordinates
[X, Y, Z] = pol2cart(Phi, Rho, Z);
% Combine X, Y, Z into a single matrix
points = [X(:), Y(:), Z(:)];
% Find unique rows and their indices
[unique_points, ~, ~] = unique(points, 'rows');
% Check for duplicates (optional)
num_duplicates = size(points, 1) - size(unique_points, 1);
if num_duplicates > 0
    fprintf('There are %d duplicate points.\n', num_duplicates);
else
    fprintf('All points are unique.\n');
end
% Plot the points in 3D space (optional)
figure;
% plot3(X, Y, Z, 'o','DisplayName','Sample Points');
xlabel('X');
ylabel('Y');
zlabel('Z');
title('3D Cylindrical Coordinates');
grid on; axis equal; hold on;
%% Circular magnetic field
% Current in negative z-direction (negative values mean the current
% direction is in positive z-direction)
I = 1;
% Absolute Value of the cylindrical H-field at corresponding positions of
% the radius in phi-direction (all other components are zero)
Hc_phi = 1/(2*pi*Rho);
% Calculate the x-, y- and z-components of the cylindrical field
Hc_x = Hc_phi .* sin(Phi);
Hc_y = -Hc_phi .* cos(Phi);
Hc_z = zeros(size(Z));
% Draw the vector-arrows of the cylindrical field
quiver3(X,Y,Z,Hc_x,Hc_y,Hc_z,'AutoScale','off');
%% Longitudinal magnetic field
Hl = 0.4; % Absolute field strength (+ = neg. z-direction)
Hl_x = zeros(size(X));
Hl_y = zeros(size(Y));
Hl_z = -Hl .* ones(size(Z));
% Draw the vector-arrows of the longitudinal field
quiver3(X,Y,Z,Hl_x,Hl_y,Hl_z,'AutoScale','off');
%% Helical magnetic field (sum of cylindrical and longitudinal field)
H_x = Hc_x + Hl_x;
H_y = Hc_y + Hl_y;
H_z = Hc_z + Hl_z;
% Draw the vector arrows of the helical field
quiver3(X,Y,Z,H_x,H_y,H_z,'AutoScale','off','LineWidth',1.5)
% Draw the fieldlines of the helical field
[startX,startY,startZ] = meshgrid([-1 1],[-1 1],1);
% startpoints = [startX(:), startY(:), startZ(:)]; % unique, manually checked
% samplepoints = [X(:), Y(:), Z(:)]; % unique, manually checked
verts  = stream3(X,Y,Z,H_x,H_y,H_z,startX,startY,startZ); 
lineobj = streamline(verts);
view(3)

Sorry, it might be a bit overwhelming. Ultimately, I just want to get streamlines of H_x, H_y and H_z on the X, Y, Z grid. But I get the following output of the script:

All points are unique.
Error using matlab.internal.math.interp1
Sample points must be unique.
Error in interp1 (line 188)
VqLite = matlab.internal.math.interp1(X,V,method,method,Xqcol);
Error in stream3 (line 66)
syi=interp1(yy(:),1:szu(1),sy(k));
Error in test (line 62)
verts  = stream3(X,Y,Z,H_x,H_y,H_z,startX,startY,startZ); 

What is the Problem here? Is my check for unique points wrong?

Thank you very much in advance!

[dpb]

The problem is that internally stream3 is interpolating in 1D and the yy() vector internally is not unique.

Your test shows that there are no duplicated ROWS in the array, but not that there are not duplicated y values.

Now, since meshgrid in cartesian coordinates includes more than one row with each of the two coordinates and stream3 can handle that per the examples; it appears the use of meshgrid with cylindrical points subsequently converted to cartesian broke the assumptions made internally about the structure of the input coordinates.

Setting a breakpoint in stream3(), at the point it crashes, it turns out that yy(:) is [0 0] which are, indeed, identical.

I'm not sure otomh how to create the mesh you're looking for in a form that would be compatible; maybe somebody else will have an idea.

[Florian G]

Thank you for the input! :)

This actually helped me solve my problem. In my revised script i basically just replaced this:

% Generate cylindrical coordinates

phi = [0 pi/2 pi 3*pi/2];

rho = [0.25 2];

z = [1 2];

% Create a meshgrid for theta, r, and z

[Phi, Rho, Z] = meshgrid(phi, rho, z);

% Convert cylindrical coordinates to Cartesian coordinates

[X, Y, Z] = pol2cart(Phi, Rho, Z);

With this:

% Generate cartesian coordinates

x = -5.025:0.05:5.025;

y = x;

z = 0:0.25:10;

% Create a meshgrid for x, y, and z

[X, Y, Z] = meshgrid(x, y, z);

% Convert cylindrical coordinates to Cartesian coordinates

[Phi, Rho, Z] = cart2pol(X, Y, Z);

To get the following working example:

clc; close; clear all;

format compact

%% Cartesian Coordinates

% Generate cylindrical coordinates

x = -5.025:0.05:5.025;

y = x;

z = 0:0.25:10;

% Create a meshgrid for x, y, and z

[X, Y, Z] = meshgrid(x, y, z);

% Convert cylindrical coordinates to Cartesian coordinates

[Phi, Rho, Z] = cart2pol(X, Y, Z);

% Combine X, Y, Z into a single matrix

points = [X(:), Y(:), Z(:)];

% Find unique rows and their indices (optional)

[unique_points, ~, ~] = unique(points, 'rows');

% Check for duplicates (optional)

num_duplicates = size(points, 1) - size(unique_points, 1);

if num_duplicates > 0

fprintf('There are %d duplicate points.\n', num_duplicates);

else

fprintf('All points are unique.\n');

end

% Configure 3D plot

figure;

% plot3(X, Y, Z, 'o','DisplayName','Sample Points'); % optionally plot points

xlabel('X');

ylabel('Y');

zlabel('Z');

title('3D Cylindrical Coordinates');

grid on; axis equal; hold on;

%% Circular magnetic field

% Current in negative z-direction (negative values mean the current

% direction is in positive z-direction)

I = 1;

% Absolute Value of the cylindrical H-field at corresponding positions of

% the radius in phi-direction (all other components are zero)

Hc_phi = 1/(2*pi*Rho);

% Calculate the x-, y- and z-components of the cylindrical field

Hc_x = Hc_phi .* sin(Phi);

Hc_y = -Hc_phi .* cos(Phi);

Hc_z = zeros(size(Z));

% % Draw the vector-arrows of the cylindrical field

% quiver3(X,Y,Z,Hc_x,Hc_y,Hc_z,'AutoScale','off');

%% Longitudinal magnetic field

Hl = 0.4; % Absolute field strength (+ = neg. z-direction)

Hl_x = zeros(size(X));

Hl_y = zeros(size(Y));

Hl_z = -Hl .* ones(size(Z));

% % Draw the vector-arrows of the longitudinal field

% quiver3(X,Y,Z,Hl_x,Hl_y,Hl_z,'AutoScale','off');

%% Helical magnetic field (sum of cylindrical and longitudinal field)

H_x = Hc_x + Hl_x;

H_y = Hc_y + Hl_y;

H_z = Hc_z + Hl_z;

% % Draw the vector arrows of the helical field

% quiver3(X,Y,Z,H_x,H_y,H_z,'AutoScale','off','LineWidth',1.5)

% Draw the fieldlines of the helical field

[startX,startY,startZ] = meshgrid([-0.25 0.25],[-0.25 0.25],10);

% startpoints = [startX(:), startY(:), startZ(:)]; % unique, manually checked

% samplepoints = [X(:), Y(:), Z(:)]; % unique, manually checked

verts = stream3(X,Y,Z,H_x,H_y,H_z,startX,startY,startZ,[1e-2 1e5]);

lineobj = streamline(verts);

view(3)

[dpb]

I knew the above would run; wasn't sure it would produce the grid you wished, however...