1、理论基础（待续）

2、算法实现

import numpy as np
import matplotlib.pyplot as plt                
R=2
I0=12
mu0=4*np.pi*1e-7
C0=mu0/(4*np.pi)
N=50
x, y = np.meshgrid(np.linspace(-5,5,N), np.linspace(-5,5,N))
#x = np.linspace(-3,3,N)
#x=linspace(-3,3,N)
#y=x

theta0 = np.linspace(0,2*np.pi,N+1)
#theta0=linspace(0,2*pi,N+1)
theta1=theta0[0:N]
y1=R*np.cos(theta1)
z1=R*np.sin(theta1)
theta2=theta0[1:N+1]
y2=R*np.cos(theta2)
z2=R*np.sin(theta2)
xc=0
yc=(y2+y1)/2
zc=(z2+z1)/2
dlx=0
dly=y2-y1
dlz=z2-z1

NGx=N
NGy=NGx
Bx =[[0] * N for _ in range(N)]
By =[[0] * N for _ in range(N)]


for i in range(1,NGy):
    for j in range(1,NGx):
        rx=x[0][j]-xc
        ry=y[i][0]-yc
        rz=0-zc
        r3=np.sqrt(rx**2+ry**2+rz**2)**3
        dlXr_x=dly*rz-dlz*ry
        dlXr_y=dlz*rx-dlx*rz
        
        Bx[i][j] = sum(C0*I0*dlXr_x/r3)
        By[i][j] = sum(C0*I0*dlXr_y/r3)
        #Bx[i].append(sum(C0*I0*dlXr_x/r3))
        #By[i].append(sum(C0*I0*dlXr_y/r3))
       

Bx = np.array(Bx)
By = np.array(By)
B=(Bx**2+By**2)**0.5
plt.figure(1)
temp = plt.quiver(x,y,Bx,By,angles="xy",color="#666666")
plt.show()

plt.figure(2)
ax = plt.subplot(111, projection='3d') 
ax.plot_surface(x, y, B, rstride=1, cstride=1,cmap='rainbow', alpha=0.9)#绘面
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()

3、结果