Магнитное поле квадратной рамки с током

Вычисляем магнитное поле простым интегрированием закона Био-Саварра-Лапласа по прямолинейному участку провода. \(d\vec{B}=\frac{\mu_0 I [\vec{dl} \times \vec{r}]}{4 \pi r^3}\)

def mf(I, x0, y0, z0, a1, a2, d):
 Bx=np.zeros_like(x0)
 By=np.zeros_like(y0)
 Bz=np.zeros_like(z0)
 lngth=int((((a2[0]-a1[0])**2+(a2[1]-a1[1])**2+(a2[2]-a1[2])**2)**0.5)/d)
 dr=(a2-a1)/(lngth)
 for i in range(lngth-1):
  q1=a1+i*dr 
  rx=q1[0]-x0
  ry=q1[1]-y0
  rz=q1[2]-z0    

  r3=(rx**2+rz**2+ry**2)**1.5
  K=I*(1.0e-3)
  Bx+=(ry*dr[2]-rz*dr[1])*K/r3
  By+=(rz*dr[0]-rx*dr[2])*K/r3
  Bz+=(rx*dr[1]-ry*dr[0])*K/r3
 return Bx, By, Bz

Выполнить описанный пример

#!/usr/bin/python 
# -*- coding: utf-8
import numpy as np
import math

# начальные данные
I=15.0 # ток, втекающий в квадрат, Ампер
d=0.05 # интервал разбиения
a1=np.array([0.0, 0.0, 0.0]) # координаты начала отрезка с током
a2=np.array([2.0, 0.0, 0.0]) # координаты конца отрезка с током
a3=np.array([2.0, 2.0, 0.0])
a4=np.array([0.0, 2.0, 0.0])


#функция вычисляющая магнитное поле в точках x0, y0, z0
def mf(I, x0, y0, z0, a1, a2, d):
 Bx=np.zeros_like(x0)
 By=np.zeros_like(y0)
 Bz=np.zeros_like(z0)
 lngth=int((((a2[0]-a1[0])**2+(a2[1]-a1[1])**2+(a2[2]-a1[2])**2)**0.5)/d)
 dr=(a2-a1)/(lngth)
 for i in range(lngth-1):
  q1=a1+i*dr 
  rx=q1[0]-x0
  ry=q1[1]-y0
  rz=q1[2]-z0    

  r3=(rx**2+rz**2+ry**2)**1.5
  K=I*(1.0e-3)
  Bx+=(ry*dr[2]-rz*dr[1])*K/r3
  By+=(rz*dr[0]-rx*dr[2])*K/r3
  Bz+=(rx*dr[1]-ry*dr[0])*K/r3
 return Bx, By, Bz

from mayavi import mlab
x0,y0,z0 = np.mgrid[-1.1:3.1:0.2, -1.1:3.1:0.2, -1.1:1.1:0.2]  #[1.0:3.0:0.2, 1.0:3.0:0.2, 1.0:3.0:0.2]
print '0'
Vx1, Vy1, Vz1 = mf(I, x0, y0, z0, a1, a2, d)
print '1'
Vx2, Vy2, Vz2 = mf(I, x0, y0, z0, a2, a3, d)
print '2'
Vx3, Vy3, Vz3 = mf(I, x0, y0, z0, a3, a4, d)
print '3'
Vx4, Vy4, Vz4 = mf(I, x0, y0, z0, a4, a1, d)
print '4'
Vx=Vx1+Vx3+Vx4
Vy=Vy1+Vy2+Vy3+Vy4
Vz=Vz1+Vz2+Vz3+Vz4
mlab.clf()
#mlab.quiver3d(x0, y0, z0, Vx, Vy, Vz, scale_factor=0.1, colormap='YlOrRd')
field = mlab.pipeline.vector_field(x0, y0, z0, Vx, Vy, Vz, name='B field')
vectors = mlab.pipeline.vectors(field)
vcp = mlab.pipeline.vector_cut_plane(field)

mlab.plot3d(np.array([0.0, 2.0]),np.array([0.0, 0.0]),np.array([0.0, 0.0]), tube_radius=0.025, colormap='Spectral')
mlab.plot3d(np.array([2.0, 2.0]),np.array([0.0, 2.0]),np.array([0.0, 0.0]), tube_radius=0.025, colormap='Spectral')
mlab.plot3d(np.array([2.0, 0.0]),np.array([2.0, 2.0]),np.array([0.0, 0.0]), tube_radius=0.025, colormap='Spectral')
mlab.plot3d(np.array([0.0, 0.0]),np.array([2.0, 0.0]),np.array([0.0, 0.0]), tube_radius=0.025, colormap='Spectral')
mlab.show()


# python -i magn5.py