Either quiver3 or streamline work for me, but not both, for the same test data

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);

Answer 1

[Rahul K]

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)

[adil hassam]

Yes, this works, thanks! Not quite sure why what I was doing before is not self-consistent.

[adil hassam]

Would appreciate any suggestions as to what I was doing inconsistently earlier.

[Rahul K]

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 )