%《MATLAB及其在理工課程中的應用指南》(十一五規劃版)第六章例6－5－2程序
% 亥姆藿茨線圈的磁場驗算
% 西安電子科技大學出版社出版  陳懷琛編著 2007年5月
%
clear all; % 清工作空間及變量 (即環半徑,電流,圖形 )初始化
mu0 = 4*pi*1e-7;       % 真空導磁率 (T*m/A)
I0 = 5.0;  Rh=1;     % 環中電流(A)及環半徑,在本題中可任設,不影響結果
C0 = mu0/(4*pi) * I0;   % 歸并常數
%c=input('輸入測量點的x坐標c*Rh,(c在0.3~0.6之間選,磁場分布最勻) ');
NGx =21 ;NGy = 21;		% 設定觀測點網格數
x=linspace(-Rh,Rh, NGx);y=linspace(-Rh,Rh, 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
%plot(y,Bx),title('沿y方向的磁場Bx分布');
%xlabel('y');ylabel('Bx');grid
Bax=Bx(:,11:21)+Bx(:,1:11);		% 把x<0區域內的磁場平移疊加到x>0的區域
Bay=By(:,11:21)+By(:,1:11);		% 右邊再加一個線圈所增加的磁場
subplot(1,2,1),
mesh(x(11:21),y,Bax);xlabel('x');ylabel('y');	% 畫出其Bx分布三維圖
subplot(1,2,2),
plot(y,Bax),grid,xlabel('y');ylabel('Bx');
q=abs((Bax-Bax(11,6))/Bax(11,6))<0.02
pause,figure(2),[i,j]=find(q);plot(i,j,'*')
set(gcf,'color','w')