Plotting multiple current loops for magnetic field

I'm attempting to plot mutiple loops of current to develop a magnetic field. The overall shape of the multiple loops should appear to be sections of a larger torus. The code that I am using ends up plotting only 4 of the expected 9 loops, but I can't determine why the other loops are not being plotted. I have attempted to plot the loop then rotate it around the Z-axis, but the quiver function only plots for the initial loop, not the additional rotated ones. Full code listed below.

% Set Parameters
I = 2.0;
u0 = 1.577e-6;
R0=1;
r_out = 1.5;
r_in = 0.5;
nR = 20;
alpha=360;
nloops=10;
% Graph Window
xmin = -2;
xmax = 2;
ymin = -2;
ymax = 2;
zmin = -2;
zmax = 2;
dx = .5;
x=xmin:dx:xmax;
y=ymin:dx:ymax;
z=zmin:dx:zmax;
plotOrigin([-2 2],[-2 2],[-2 2], [0,0,0]);
theta=linspace(0,2*pi,nR+1);
phi=linspace(0,2*pi,nR+1);
sigma=linspace(0,alpha, nloops);
r=.5*(r_out-r_in)*ones(1,nR+1);
Rz = r.*cos(theta);
Ry=r.*sin(theta)+r_in+.5*(r_out-r_in);
Rx=zeros(1,nR+1);
[X,Y,Z]=meshgrid(x,y,z);
Bx=zeros(size(X));
By=zeros(size(X));
Bz=zeros(size(X));
Rx1=[];
Ry1=[];
Rz1=[];
for nt=1:nloops-1
    temp_CS = rotationMatrix('z', sigma(nt))*[Rx;Ry;Rz];
    Rx=[Rx1;temp_CS(1,:)];
    Ry=[Ry1;temp_CS(2,:)];
    Rz=[Rz1;temp_CS(3,:)];
    for it=1:nR
        Points=transpose([Rx;Ry;Rz]);
        rx=X-Points(it,1);
        ry=Y-Points(it,2);
        rz=Z-Points(it,3);
        R=sqrt(rx.^2+ry.^2+rz.^2);
        R=R.*(R>=dx)+dx*(R<dx);
        dLx = Points(it+1,1)-Points(it,1);
        dLy=Points(it+1,2)-Points(it,2);
        dLz=Points(it+1,3)-Points(it,3);
        Bx=Bx+u0*I/(4*pi)*(dLy.*rz-dLz*ry).*R.^(-3);
        By=By+u0*I/(4*pi)*(dLz.*rx-dLx*rz).*R.^(-3);
        Bz=Bz+u0*I/(4*pi)*(dLx.*ry-dLy*rx).*R.^(-3);  
    end
    plot3(Points(:,1),Points(:,2),Points(:,3), 'black')
    axis([xmin xmax ymin ymax zmin zmax])
    hold on
end
quiver3(X,Y,Z, Bx,By,Bz,'b')
hold off
rotate3d on
axis square
xlabel('x')
ylabel('y')
zlabel('z')
view(-30,30)
title('test')
function []=plotOrigin(rangeX, rangeY, rangeZ, offset)
x0 = offset(1);
y0 = offset(1);
z0 = offset(1);
plot3(rangeX+x0,[y0 y0],[z0 z0], '--r')
axis equal
hold on
plot3([x0 x0], rangeY+y0,[z0 z0], '--b')
plot3([x0 x0], [y0 y0],rangeZ+z0, '--g')

end
function rotation_matrix=rotationMatrix(type, angle)
angle= angle*(2*pi)/360;
if type == 'x'
    rotation_matrix=[1   0   0
        0 cos(angle) -sin(angle)
        0 sin(angle) cos(angle)];
elseif type == 'y' 
    rotation_matrix=[cos(angle) 0 sin(angle)
        0          1          0
        -sin(angle) 0      cos(angle)];
    
elseif type == 'z'
    rotation_matrix=[cos(angle) -sin(angle) 0
        sin(angle)  cos(angle) 0
        0           0          1];
else
    fprintf("invalid input for type");
end
end

Answer 1

[Voss]

The problem is that Rx, Ry, and Rz are overwritten in the outer for loop:

% before the loop, calculate Rx, Ry, Rz:
Rz = r.*cos(theta);
Ry=r.*sin(theta)+r_in+.5*(r_out-r_in);
Rx=zeros(1,nR+1);

% ...

for nt=1:nloops-1
    
    % use Rx, Ry, Rz to calculate tempCS (OK the first time):
    temp_CS = rotationMatrix('z', sigma(nt))*[Rx;Ry;Rz];
    
    % then overwriting Rx, Ry, Rz (these values will be used 
    % to calculate temp_CS on the next iteration of this loop):
    Rx=[Rx1;temp_CS(1,:)];
    Ry=[Ry1;temp_CS(2,:)];
    Rz=[Rz1;temp_CS(3,:)];
    
    % ...
    
end

(If you inspect the value of Points each time through the loop, you'll see Points is the same on iteration #1 as it is on iteration #9, Points is the same in iterations #2, #5, and #8, Points is the same on iterations #3 and #7, and Points is the same on iterations #4 and #6. All 9 loops are plotted, but in fact there are only four unique values of Points, so you see only 4 distinct loops. This is happening because the Rx, Ry, and Rz are not correct after the first iteration.)

To fix it, don't change Rx, Ry, or Rz inside the loop, e.g.:

% ... 

Rz = r.*cos(theta);
Ry=r.*sin(theta)+r_in+.5*(r_out-r_in);
Rx=zeros(1,nR+1);

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

Bx=zeros(size(X));
By=zeros(size(X));
Bz=zeros(size(X));

% (no need for Rx1, Ry1, Rz1)

for nt=1:nloops-1
    
    temp_CS = rotationMatrix('z', sigma(nt))*[Rx;Ry;Rz];
    
    % note that Points can be calculated here, outside the inner loop
    Points = transpose(temp_CS);
    
    for it=1:nR
        % the inner loop in the same as before, except that Points is not
        % calculated (since it is already calculated before the loop)
        rx=X-Points(it,1);
        ry=Y-Points(it,2);
        rz=Z-Points(it,3);
        % ...
    end
    % ...
end

[Eric Wilckens]

That did it! Thanks for the assistance!

[Voss]

You're welcome!