3 views (last 30 days)

Show older comments

I want to do a quiver3 vector plot and streamline. I am doing a test with simple data as defined in the attached file testx.m. If I define my grid with meshgrid, quiver3 gives an incorrect vector field, but streamline works expectedly on the incorrect field. If I define my grid as X(i,j,k)=x(i), etc, quiver3 works correctly but streamline returns the error below. Would appreciate some help. This may have something to do with the way grids are defined but, as it stands, it seems internally inconsistent between quiver3 and streamline.

>> h=streamline(X,Y,Z,Bx,By,Bz,0.8,0,2) ??? Error using ==> interp1 at 261 The values of X should be distinct.

Error in ==> stream3 at 67 sxi=interp1(xx(:),1:szu(2),sx(k));

Error in ==> streamline at 69 verts = stream3(x,y,z,u,v,w,sx,sy,sz,options);

Rahul Kalampattel
on 3 Mar 2017

Edited: Rahul Kalampattel
on 3 Mar 2017

Does this look like what you're expecting?

Edited code:

qpl=4;

qpd=4;

ky=2*pi/qpl;

ep=0.3;

nx=31;

ny=31;

nz=9;

dx=1/(nx-5);

dy=qpl/(ny-5);

dz=qpd/(nz-5);

% Taken out of the for loop, indices have been shifted accordingly

x=dx*(-2:nx-3);

y=dy*(-2:ny-3);

z=dz*(-2:nz-3);

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

% Taken out of the for loop, Bx calculated using the meshgrid version of Y

Bx=ep*sin(ky*Y);

By=ones(nx,ny,nz);

Bz=zeros(nx,ny,nz);

figure

quiver3(X,Y,Z,Bx,By,Bz)

streamline(X,Y,Z,Bx,By,Bz,0.8,0,2)

Rahul Kalampattel
on 4 Mar 2017

It has to do with the order in which you were allocating elements to the matrices in your for loop.

In your X matrix, in each 'layer' (e.g. z=1) the values were constant across the rows and increasing down the columns:

>> X(:,:,1)

[-0.769 -0.769 ... -0.769

-0.385 -0.385 ... -0.385

0.000 0.000 ... 0.000

. . ... .

. . ... . ]

When streamline tries to interpolate, it expects X, Y and Z to be in 'full grid' format (i.e. an output of meshgrid or ndgrid, which in this case is the transpose of the above):

>> X(:,:,1)

[-0.769 -0.385 0.000 . .

-0.769 -0.385 0.000 . .

... ... ... ... ...

-0.769 -0.385 0.000 . . ]

Looking at your original code, you could achieve this by switching a couple of indices:

qpl=4;

qpd=4;

ky=2*pi/qpl;

ep=0.3;

nx=31;

ny=31;

nz=9;

dx=1/(nx-5);

dy=qpl/(ny-5);

dz=qpd/(nz-5);

for i=1:nx;

for j=1:ny;

for k=1:nz;

x(i)=dx*(i-3);

y(j)=dy*(j-3);

z(k)=dz*(k-3);

X(i,j,k)=x(j); % Switched i for j

Y(i,j,k)=y(i); % Switched j for i

Z(i,j,k)=z(k);

Bx(i,j,k)=ep*sin(ky*y(i)); % Switched j for i

By(i,j,k)=1;

Bz(i,j,k)=0;

end;

end;

end;

(Example of full grid required for interp2: http://au.mathworks.com/help/matlab/ref/interp2.html#bt0m664-2 )

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

Start Hunting!