python-3.x - python中的磁偶极子

Question

我看了this code :

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()

它显示了下图:



是否有可能获得一条力线的坐标？我不明白它是如何绘制的。

Answer 1

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

线条:流线的 LineCollection

箭头:包含 FancyArrowPatch 对象的 PatchCollection(这些是三角形箭头)

c.lines.get_paths()给出所有的段。遍历这些段，可以检查它们的顶点。当一个段从前一个结束的地方开始时，两者都属于同一条曲线。注意每一段都是一条短直线；许多段一起使用以形成流线曲线。

下面的代码演示了对段的迭代。为了显示正在发生的事情，每个线段都被转换为适合于 plt.plot 的二维点数组。 .默认，plt.plot用新颜色为每条曲线着色(每 10 次重复)。圆点显示了每个短直线段的位置。

要找到一条特定的曲线，您可以将鼠标悬停在起点上，并注意该点的 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()