Plotting multiple current loops for magnetic field

조회 수: 13 (최근 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
Eric Wilckens
Eric Wilckens 2022년 11월 28일
That did it! Thanks for the assistance!
Voss
Voss 2022년 11월 28일
You're welcome!

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Line Plots에 대해 자세히 알아보기

제품

Community Treasure Hunt

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

Start Hunting!

Translated by