python中的磁偶极子

import numpy as np from matplotlib import pyplot as plt def dipole(m, r, r0): """ Calculation of field B in point r. B is created by a dipole moment m located in r0. """ # R = r - r0 - subtraction of elements of vectors r and r0, transposition of array R = np.subtract(np.transpose(r), r0).T # Spatial components of r are the outermost axis norm_R = np.sqrt(np.einsum("i...,i...", R, R)) # einsum - Einsteinova sumace # Dot product of R and m m_dot_R = np.tensordot(m, R, axes=1) # Computation of B B = 3 * m_dot_R * R / norm_R**5 - np.tensordot(m, 1 / norm_R**3, axes=0) B *= 1e-7 # abbreviation for B = B * 1e-7, multiplication B of 1e-7, permeability of vacuum: 4\pi * 10^(-7) # The result is the magnetic field B return B X = np.linspace(-1, 1) Y = np.linspace(-1, 1) Bx, By = dipole(m=[0, 1], r=np.meshgrid(X, Y), r0=[-0.2,0.8]) plt.figure(figsize=(8, 8)) plt.streamplot(X, Y, Bx, By) plt.margins(0, 0) plt.show()

这个 streamplot 返回包含两部分的容器对象“StreamplotSet”:

线条:流线的线条集合
箭头:包含FancyArrowPatch对象的PatchCollection(这些是三角形箭头)

c.lines.get_paths() 给出所有片段。通过迭代这些段,可以检查它们的顶点。当一段从上一段结束的地方开始时,两者都属于同一条曲线。请注意,每个线段都是一条短直线;许多线段一起用于形成流线型曲线。

下面的代码演示如何遍历这些段。为了显示正在发生的事情,每个段都被转换成一个二维点数组,适合于 plt.plot . 违约, plt.绘图

要查找一条特定曲线,可以将鼠标悬停在起点上,并注意该点的x坐标。然后测试代码中的坐标。例如,开始于 x=0.48 是以一种特殊的方式画出来的。

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import patches

def dipole(m, r, r0):
    R = np.subtract(np.transpose(r), r0).T
    norm_R = np.sqrt(np.einsum("i...,i...", R, R))
    m_dot_R = np.tensordot(m, R, axes=1)
    B = 3 * m_dot_R * R / norm_R**5 - np.tensordot(m, 1 / norm_R**3, axes=0)
    B *= 1e-7
    return B

X = np.linspace(-1, 1)
Y = np.linspace(-1, 1)

Bx, By = dipole(m=[0, 1], r=np.meshgrid(X, Y), r0=[-0.2,0.8])

plt.figure(figsize=(8, 8))
c = plt.streamplot(X, Y, Bx, By)
c.lines.set_visible(False)
paths = c.lines.get_paths()
prev_end = None
start_indices = []
for index, segment in enumerate(paths):
    if not np.array_equal(prev_end, segment.vertices[0]):  # new segment
        start_indices.append(index)
    prev_end = segment.vertices[-1]
for i0, i1 in zip(start_indices, start_indices[1:] + [len(paths)]):
    # get all the points of the curve that starts at index i0
    curve = np.array([paths[i].vertices[0] for i in range(i0, i1)] + [paths[i1 - 1].vertices[-1]])
special_x_coord = 0.48
for i0, i1 in zip(start_indices, start_indices[1:] + [len(paths)]):
    # get all the points of the curve that starts at index i0
    curve = np.array([paths[i].vertices[0] for i in range(i0, i1)] + [paths[i1 - 1].vertices[-1]])
    if abs(curve[0,0] - special_x_coord) < 0.01:  # draw one curve in a special way
        plt.plot(curve[:, 0], curve[:, 1], '-', lw=10, alpha=0.3)
    else:
        plt.plot(curve[:, 0], curve[:, 1], '.', ls='-')

plt.margins(0, 0)
plt.show()