求截面为矩形的圆线圈周围产生的磁场
一、数值方法
（一）数学模型：所研究的电流圆线圈产生磁场的问题在柱坐标系下研究， 根据磁场强度跟矢势之间的关系，得到磁场；
磁场为B ，矢势为A
B A =∇⨯
r r z z A A e A e A e θθ=++ A e θθ= (,)A r z e θθ=
（由A 具有轴对称得到）
所以B A =∇⨯ A e θθ=∇⨯
在柱坐标系中，由公式1()()11()()r r z z
z r
r z r
z f f e f e f e f f f r z
f f f z r
f
f rf r r r θθθ
θθθθ
⎧=++⎪∂∂⎪∇⨯=-⎪∂∂⎪⎨∂∂∇⨯=-⎪∂∂⎪∂∂
⎪∇⨯=-⎪∂∂⎩
-得
B A =∇⨯ 1()r z f e rf e z r r θθ∂∂
=-+∂∂
即r A B z θ
∂=-∂，1()z B rA r r θ∂
=∂
（1）先求矢势A
4L Idl A r
μπ=⎰
一个电流为I ，半径为a 的线圆环周围空间产生的磁场，其矢势表示为
202220cos (,)42cos Ia A r z d r z a ar πθμϕ
ϕπϕ=++-⎰
推广到截面为矩形的圆环线圈中
22
11202220
cos (,)4()2cos R z R z I r A r z d dz dr s r z z r r r πθμϕ
ϕπϕ'''='''+-+-⎰⎰⎰
其中S 为矩形截面的面积，12,R R 为矩形截面的两边距圆环中心的距离，12,z z 为矩形截面的上下面的z 轴坐标。
（二）数值模型离散化（均匀网格有限差分）
（1）高斯方法计算三重积分（参考书:徐士良常用算法程序集第二版）
（2）根据一阶三点公式来求解磁场的分布：
磁场分量r B ,z B
,1,1()()2i j i j r A A A B z z
θθθ-+-∂=-=∂∆ 1,1,,,()()()1()2i j i j i j z i j A A A B rA r r r r θθθθ+--∂=
=+∂∆ 对于0r →点处
1,1,0
00()()1lim lim ()lim()2i j i j z r r r A A A A A B rA r r r r r r θθθθθθ+-→→→-∂∂∂==+==∂∂∂∆ 22z r B B B +=
二 数值计算程序
SUBROUTINE FGAUS(N,JS,X,FS,F,S,Z,R)
DIMENSION JS(N),X(N)
DIMENSION T(5),C(5),D(2,11),CC(11),IS(2,11)
DATA T/-0.90617,-0.538469,0.0,0.538469,0.90617/
DATA C/0.2369,0.47862,0.568889,0.4786,0.2369/
M=1
D(1,N+1)=1.0
D(2,N+1)=1.0
10 DO 20 J=M,N
CALL FS (J,N,X,DN,UP)
D(1,J)=0.5*(UP-DN)/JS(J)
CC(J)=D(1,J)+DN
X(J)=D(1,J)*T(1)+CC(J)
D(2,J)=0.0
IS(1,J)=1
IS(2,J)=1
20 CONTINUE
J=N
30 K=IS(1,J)
IF(J.EQ.N)THEN
P=F(N,X,Z,R)
ELSE
P=1.0
ENDIF
D(2,J)=D(2,J+1)*D(1,J+1)*P*C(K)+D(2,J) IS(1,J)=IS(1,J)+1
IF(IS(1,J).GT.5)THEN
IF(IS(2,J).GE.JS(J)) THEN
J=J-1
IF(J.EQ.0) THEN
S=D(2,1)*D(1,1)
RETURN
ENDIF
GOTO 30
ENDIF
IS(2,J)=IS(2,J)+1
CC(J)=CC(J)+D(1,J)*2.0
IS(1,J)=1
ENDIF
K=IS(1,J)
X(J)=D(1,J)*T(K)+CC(J)
IF(J.EQ.N) GOTO 30
M=J+1
GOTO 10
END
EXTERNAL FS,F
DIMENSION
JS(3),X(3),Z(30),R(30),a0(30,30),BB(30,30),BR(30,30),BZ(30,30) DATA JS/4,4,4/
N=3
c=1
H=0.5
Z(1)=1.5
DO I=1,29
Z(I+1)=Z(I)+H
ENDDO
R(1)=1.5
do j=1,29
R(J+1)=R(J)+H
ENDDO
DO I=1,30
DO J=1,30
CALL FGAUS(N,JS,X,FS,F,S,Z(I),R(J))
a0(i,j)=c*s
WRITE(1,*)Z(I),R(J),S
OPEN(1,FILE='DUHAI.DAT')
ENDDO
ENDDO
!求解磁场的Br
do i=1,30
Br(i,1)=(-3*A0(i,1)+4*A0(i,2)-A0(i,3))/2*h
Br(i,30)=(A0(i,28)-4*A0(i,29)+3*A0(i,30))/2*h
do j=2,29
Br(i,j)=(A0(i,j-1)-A0(i,j+1))/2*h
end do
end do
