Possible bug of streamslice function or just a misuse of it?

Greetings!

The problem is simple and very straightforward: I want to plot the magnetic flow of a bar magnet (equal to a solenoid magnetic field), calculated by these equations. My MATLAB version is `R2013a`.

The function that does the job, `magnetic_field` (it returns a struct with the coordinates (r,z) and the field (Br,Bz), B(r,z) = Br(r,z) r + Bz(r,z) z is calculated), is at the end of this question.

I then call

    >> t = magnetic_field(2e-2, 16e-2, 20, -6e-2, 6e-2, 0.5e-2, -18.2e-2, 18.2e-2, 1e-2);

    >> streamslice(t.r, t.z, t.Br, t.Bz)

Notice that `t.r` and `t.z` are a `meshgrid` result in which `Br` and `Bz` are calculated (see the `for` loop in the function `magnetic_field`). Thus, all these four matrices have the same `size`.

MATLAB gives the error:

Attempted to access startgrid(NaN,13); index must be a positive integer or logical.

Error in streamslice>nicestreams (line 324) startgrid(rr,cc)=1;

Error in streamslice (line 134) [vertsout, avertsout] = nicestreams(x,y,u,v,density,arrows);

The weirdness is: if I change the grid, the streamslice gives no error:

    >> t = magnetic_field(2e-2, 16e-2, 20, -2, 2, 0.1, -2, 2, 0.1);

Although it produces a very strange result (because the solution is only defined for positive `r`, I presume).

I went to check with MATLAB sample data, like the `wind` data:

    >> w = load('wind');
    >> streamslice(w.x(:,:,5), w.y(:,:,5), w.u(:,:,5), w.v(:,:,5))

It plots correctly. However, if I change from `5` to either `2`, `3`, `10`, `12` or `15`, MATLAB displays a very similar error to the one I had:

Attempted to access startgrid(2,0); index must be a positive integer or logical.

Error in streamslice>nicestreams (line 324) startgrid(rr,cc)=1;

Error in streamslice (line 134) [vertsout, avertsout] = nicestreams(x,y,u,v,density,arrows);

the difference is that with the sample data, the error is that `streamslice` tries to access a `0` index, whilst in my case, it tries to access a `NaN` index.

Thus, I keep with either two possibilities :

1. My data (namely the variables `Br` and `Bz`) do not have the structure required by `streamslice` function
2. it is a bug in the `streamslice` function

What do you think?

Any suggestions?

    function res = magnetic_field(a, L, I, rMin, rMax, dr, zMin, zMax, dz)    
        Nr = (rMax - rMin) ./ dr + 1;
        Nz = (zMax - zMin) ./ dz + 1;

        % calculating dependent constants
        mu0div2pi = 2e-7;
        Lm1 = 1 / L; % in 1/m
        Ldiv2 = L / 2;
        constZ = -Lm1 .* I * mu0div2pi / 4;
        constR = Lm1 .* I * mu0div2pi;

        % defining functions
        kSqr = @(ar4, ar2, zeta) ar4 ./ (ar2 + zeta .^ 2);
        hSqr = @(ar4, ar2) ar4 ./ ar2;
        zetaPlus = @(z, Ldiv2) z - Ldiv2;
        zetaMinus = @(z, Ldiv2) z + Ldiv2;
        IntBr = @(kSqr, eK) ((kSqr - 2) .* eK + 2 .* ellipticE(kSqr)) ./ sqrt(kSqr);
        IntBz = @(kSqr, hSqr, r, zeta, a, eK) zeta .* sqrt(kSqr) .* (eK + ellipticPi(hSqr, kSqr) .* (a - r) ./ (a + r));

        % creating grid
        [r, z] = meshgrid(linspace(rMin, rMax, Nr), linspace(zMin, zMax, Nz));
        [m, n] = size(r);
        NN = m * n;

        Br = zeros(m,n);
        Bz = zeros(m,n);
        for i = 1:NN
            ar4 = 4 .* a .* r(i);
            ar2 = (a + r(i)) .^ 2;
            hh = hSqr(ar4, ar2);
            zP = zetaPlus(z(i), Ldiv2);
            zM = zetaMinus(z(i), Ldiv2);
            kkPlus = kSqr(ar4, ar2, zP);
            kkMinus = kSqr(ar4, ar2, zM);
            eKPlus = ellipticK(kkPlus);
            eKMinus = ellipticK(kkMinus);
            Br(i) = constR .* sqrt(a ./ r(i)) .* (IntBr(kkPlus, eKPlus) - IntBr(kkMinus, eKMinus));
            Bz(i) = constZ .* ( IntBz(kkPlus, hh, r(i), zP, a, eKPlus) - IntBz(kkMinus, hh, r(i), zM, a, eKMinus) ) ./ sqrt(a .* r(i));
        end
        res.a = a;
        res.L = L;
        res.I = I;
        res.r = r;
        res.z = z;
        res.Br = Br;
        res.Bz = Bz;
    end

EDIT

I attached the result of

    >> t = magnetic_field(2e-2, 16e-2, 20, -6e-2, 6e-2, 0.5e-2, -18.2e-2, 18.2e-2, 1e-2);

[Sean de Wolski]

Can you attach a ZIP file that contains:

• a MAT file with all of your variables
• The exact code you are using to call streamslice.

[Maurício Girardi-Schappo]

@Sean de Wolski: no other variables set in the workspace... Just start MATLAB, change directory, run the function as I explained here... The call to streamslice is also stated in the question:

    streamslice(t.r, t.z, t.Br, t.Bz)

PS: I attached the output of magnetic_field function call

Answer 1

[Maurício Girardi-Schappo]

The solution was found by me and Amro in StackOverflow website. Please refer to:

StackOverflow Answer

Answer 2

[Walter Roberson]

The NaN error can happen if your y values are all identical and (at least) one of your v values is the same as your y values. In such a case, (v - min(y)) would be 0, and that would be internally multiplied by 20/range(y) which would be 20/0 = inf if all the y are the same; the 0 times inf would generate a nan.

The 0 error can happen if one of your v values is less than your minimum y value.

[Maurício Girardi-Schappo]

I checked if any of my y values are equal to any of my v values: none of them are. I also checked if (v-min(y)) == 0, which is not true for any case. Notice that the y data has been generated by a call to meshgrid:

    [r,z] = meshgrid(linspace(...), linspace(...))

and the matrixes Br,Bz have been calculated by using the linear index (in the for loop):

    Br(i) = ...
    Bz(i) = ...

PS: I attached the output of magnetic_field function in the main question