Магнети c диполь в python - PullRequest
       25

Магнети c диполь в python

4 голосов
/ 17 февраля 2020

Я посмотрел на этот код :

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

Показывает следующий рисунок:

enter image description here

Можно ли получить координаты одной силовой линии? Я не понимаю, как это строится.

1 Ответ

3 голосов
/ 18 февраля 2020

streamplot возвращает контейнерный объект 'StreamplotSet' с двумя частями:

  • lines: LineCollection линий потока
  • стрелки: PatchCollection, содержащая FancyArrowPatch объекты (это стрелки tri angular)

c.lines.get_paths() дает все сегменты. Итерируя по этим сегментам, можно проверить их вершины. Когда сегмент начинается там, где закончился предыдущий, оба относятся к одной и той же кривой. Обратите внимание, что каждый сегмент представляет собой короткую прямую линию; многие сегменты используются вместе для формирования кривой тока.

Код ниже демонстрирует итерацию по сегментам. Чтобы показать, что происходит, каждый сегмент преобразуется в массив точек 2D, подходящих для 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()

resulting plot

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...