%《MATLAB及其在理工課程中的應用指南》(十一五規劃版)第六章例6－5－1程序
% 由畢奧－薩瓦定律計算磁場
% 西安電子科技大學出版社出版  陳懷琛編著 2007年7月
%
clear all; % 清工作空間及變量 (即環半徑,電流,圖形 )初始化
mu0 = 4*pi*1e-7;       % 真空導磁率 (T*m/A)
I0 = 5.0;       % 環中電流(A)
Rh = input('輸入環半徑Rh(m):(書上取2) ');
C0 = mu0/(4*pi) * I0;   % 歸并常數
xMax = 3; yMax = 3;              % 規定圖的范圍
NGx = 21; NGy = 21;          		% 規定觀測點網格線數
x=linspace(-xMax, xMax, NGx);		% 確定觀測點的x,y坐標數組
y=linspace(-yMax, yMax, NGy);
Nh = 20;		% 電流環分段數 
% 計算每段的端點,環在x=0平面上,其坐標x1,x2均為零
theta0 = linspace(0,2*pi, Nh+1);	% 環的圓周角分段
theta1 = theta0(1:Nh);
y1 = Rh*cos(theta1);		% 環各段向量的起點坐標y1,z1
z1 = Rh*sin(theta1);
theta2 = theta0(2:Nh+1);
y2 = Rh*cos(theta2);		% 環各段向量的終點坐標y2,z2
z2 = Rh*sin(theta2);  
dlx = 0;   		% 計算環各段向量dl的三個分量
dly = y2-y1;
dlz = z2-z1;  
xc = 0; 	  	% 計算環各段向量中點的三個坐標分量
yc = (y2+y1)/2;
zc = (z2+z1)/2;
% 循環計算各網格點上的B(x,y) 值
for i=1:NGy
   for j=1:NGx 
      % 對yz平面內的電流環分段作元素群運算,先算環上某段與觀測點之間的向量r
      rx = x(j) - xc;
      ry = y(i) - yc;
      rz = -zc;             % 觀測點在z=0平面上	 
      r3 = sqrt(rx.^2 + ry.^2 + rz.^2).^3;	 	 % 計算r^3
      dlXr_x = dly.*rz - dlz.*ry; 	 % 計算叉乘積dl X r的 x 和 y 分量 
      dlXr_y = dlz.*rx - dlx.*rz;	 
      Bx(i,j) = sum(C0*dlXr_x./r3); 	 % 把磁場各段的x 和 y 分量累加
      By(i,j) = sum(C0*dlXr_y./r3);
   end
end
clf;   quiver(x,y,Bx,By);    	%  用quiver 畫磁場向量圖
hold on;						
plot(0,Rh,'bo');plot(0,-Rh,'rx');  	% 在圖上畫出電流環
xlabel('x');  ylabel('y');
hold off;