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 :
- My data (namely the variables `Br` and `Bz`) do not have the structure required by `streamslice` function
- 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);
Best Answer