Plotting multiple current loops for magnetic field

5 ビュー (過去 30 日間)
Eric Wilckens
Eric Wilckens 2022 年 11 月 26 日
編集済み: Voss 2022 年 11 月 28 日
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

採用された回答

Voss
Voss 2022 年 11 月 26 日
編集済み: Voss 2022 年 11 月 27 日
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
  2 件のコメント
Voss
Voss 2022 年 11 月 28 日
You're welcome!

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Find more on Vector Fields in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by