How to compare 2D vector fields and minimize the difference?

Question

I want to compare the field of two electrical currents and compare the resulting field to a magnetic dipol field and find magnetic momentum that minimizes the difference of the two fields.

My current approach is to calculate the field of the two currents (by superimposing the fields of the two individual currents) and then start a optimization/search algorithm that minimizes the difference of the two fields. Because I am not interested in the immediate surroundings of the dipol or wires, I mask out the area around the center (and set the values of the fields there to zero).

To compute the overall error that the search algorithm should minimize, I subtract the dipol field from the field from the two currents. I use vectorization, so I end up just substracting whole fields from each other. Feels nice, but perhaps its wrong? I assume the vecorization substracts each vector in each coordinate. I want those individual vectors to be as similar as possible and want to minimize all those individual, small errors, added up. But I fear this optimization does not take place because the fields dont look similar.

How do i compare fields effectively?

here is my code so far, so you can see all the details of how i did it.

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import minimize


# def magnetic_dipol_field_in_point(dipol_coordinates, dipol_vector, point_coordinates):
#     magnetic_permeability = 4 * np.pi * 1e-7
#     direction_vector = (point_coordinates - dipol_coordinates)
#     direction_vector_norm = np.linalg.norm(direction_vector, axis=0)
#     direction_unit_vector = np.divide(direction_vector, direction_vector_norm, out=np.zeros_like(direction_vector),
#                                       where=direction_vector_norm != 0)
#     vector_product = np.dot(dipol_vector, direction_unit_vector)
#     unit_vect_prod_ = direction_unit_vector * vector_product
#     dipol_vector = dipol_vector[:, np.newaxis]
#     dipol_diff = unit_vect_prod_ - dipol_vector
#     b_vector = (magnetic_permeability / (4 * np.pi)) * (3 * dipol_diff)\
#                / np.power(direction_vector_norm, 3)
#     return b_vector

def magnetic_dipol_field_in_point(dipol_coordinates, dipol_vector, point_coordinates):
    mu_0 = 4 * np.pi * 1e-7
    direction_vector = (point_coordinates - dipol_coordinates)
    direction_vector_norm = np.linalg.norm(direction_vector, axis=0)
    direction_unit_vector = (direction_vector / direction_vector_norm)
    b_vector = (mu_0 / (4 * np.pi)) * (3 * direction_unit_vector * np.dot(dipol_vector.T, direction_unit_vector)
                                       - dipol_vector) / np.power(direction_vector_norm, 3)

    return b_vector

def magnetic_vector_in_point(current_coordinates, current, point_coordinates):
    magnetic_permeability = 4 * np.pi * 1e-7
    direction_vector = (point_coordinates - current_coordinates)
    direction_vector_norm = np.linalg.norm(direction_vector, axis=0)
    direction_unit_vector = np.divide(direction_vector, direction_vector_norm, out=np.zeros_like(direction_vector),
                                      where=direction_vector_norm != 0)
    # orthogonal_unit_vector = np.flip(direction_unit_vector, axis=0) * np.array([[-1], [1]])
    # faster version of the above:
    # https://stackoverflow.com/questions/64795722/how-to-calculate-the-orthogonal-vector-of-a-unit-vector-with-numpy/64797015#64797015
    orthogonal_unit_vector = direction_unit_vector[::-1, :]
    np.negative(orthogonal_unit_vector[0, :], out=orthogonal_unit_vector[0, :])

    b_vector = (current * magnetic_permeability / (2 * np.pi)) * np.divide(orthogonal_unit_vector,
                                                                           direction_vector_norm,
                                                                           out=np.zeros_like(orthogonal_unit_vector),
                                                                           where=direction_vector_norm != 0)
    return b_vector


def calc_current(time, num, peak):
    return peak * np.cos(num * 2 * np.pi / 3 + time)  # time is in the interval 0..2pi


def mask_area(field, grid, coord, radius):
    x_center = coord[0]
    y_center = coord[1]
    field[0] = np.where(np.sqrt((grid[0] - x_center) ** 2 + (grid[1] - y_center) ** 2) > radius, field[0], 0)
    field[1] = np.where(np.sqrt((grid[0] - x_center) ** 2 + (grid[1] - y_center) ** 2) > radius, field[1], 0)
    return field


def main():
    size = 4
    a = 2.67e-2

    X, Y = np.mgrid[-size:size:.1, -size:size:.1]
    grid = np.vstack([X.ravel(), Y.ravel()])

    current_coordinates = np.array([[a, -a / np.sqrt(3), -a / np.sqrt(3)],
                                    [0, a, -a]])
    magnetic_fields = np.zeros([2 + 1, grid.shape[0], grid.shape[1]])

    currents = np.zeros(3)

    for i in range(0, 2):
        currents[i] = calc_current(1 + np.pi / 2, i, 1)
        currents = [5, -5, 0]
        magnetic_fields[i] = magnetic_vector_in_point(current_coordinates[:, np.newaxis, i], currents[i],
                                                      grid)
        magnetic_fields[i] = mask_area(magnetic_fields[i], grid, [0.00564237, 0.01335], 1)

    def find_dipol_moment(x, fields, slot_cnt, end_points, koordinaten_feld):
        # calculate middle between currents, where dipol should be located.
        dipol_coord = np.mean(end_points[:, 0:2], axis=1)
        # calculate direction of the dipol moment, it is orthogonal to the line between the currents
        vect = end_points[:, 0] - end_points[:, 1]
        dipol_vektor = np.flip(vect, axis=0) * np.array([-1, 1])
        #
        dipol_vektor = dipol_vektor[:, np.newaxis]
        dipol_unit_vektor = - dipol_vektor / np.linalg.norm(dipol_vektor)
        dipol_coord = dipol_coord[:, np.newaxis]
        fields[slot_cnt] = magnetic_dipol_field_in_point(dipol_coord, x[0] * dipol_unit_vektor, koordinaten_feld)
        fields[slot_cnt] = -mask_area(fields[slot_cnt], grid, [0.00564237, 0.01335], 1)
        error_field = np.sum(magnetic_fields, axis=0)
        error = np.nansum(np.linalg.norm(error_field, axis=1))
        print("Dipolmoment:", x,"\tFehler: ", error)
        return error

    res = minimize(find_dipol_moment, -5,  # 1 ist mein anfangswert bei der suche
                   args=(magnetic_fields, 2, current_coordinates, grid))
    x = res.x
    print("Ergebnis: äquivalentes Dipolmoment", x)
    magnetic_field = np.sum(magnetic_fields[0:2], axis=0)
    dipol_field = magnetic_fields[2]
    # overall_field = np.sum(magnetic_fields, axis=0)
    #    magnetic_field = magnetic_fields[2]

    fig, ax = plt.subplots(figsize=(10, 10), dpi=180, facecolor='w', edgecolor='k')
    ax.scatter(0.00564237, 0.01335, marker="o")
    ax.scatter(a, 0, marker="+")
    ax.scatter(-a / np.sqrt(3), a, marker="+")
    q = ax.quiver(grid[0, :], grid[1, :], magnetic_field[0, :], magnetic_field[1, :], color='b')
    q = ax.quiver(grid[0, :], grid[1, :], dipol_field[0, :], dipol_field[1, :], color = 'k')
    # q = ax.quiver(coordiante_field[0, :], coordiante_field[1, :], overall_field.T[0, :], overall_field.T[1, :],
    #              color='r')
    # print(str(np.nansum(np.linalg.norm(overall_field, axis=1 )) ))
    ax.quiverkey(q, X=0.3, Y=1.1, U=10,
                 label='Quiver key, length = 10', labelpos='E')
    ax.set_aspect('equal', 'box')
    plt.show()


if __name__ == '__main__':
    main()

Answer 1

We consider two thin infinite straight parallel wires at distance $2 a$ apart and carrying equal currents $I$ in opposite directions. We need to find an approximate description of the system in terms of multipole expansion moments.

Let's first make a few preliminary comments. The magnetic field is described by the magnetic vector potential $\vec{A}$ that is related to the magnetic field $\vec{B}$ according to $\nabla \times \vec{A} = \vec{B}$, and $\nabla^2 \vec{A} = - \mu_0 \vec{j}$, where $\vec{j}$ is the electric current (using SI units).

Let's use the coordinate $z$ in the direction of the wires. For an infinite straight wire in the $\hat{e}_z$ direction the vector potential is also in the $\hat{e}_z$ direction. Due to the symmetry, our problem is 2D and the magnetic field is fully described by a scalar quantity $A_z(x,y)$.

For a single wire located in the 2D plane at $\vec{r}_1 =(x_1,y_1)$ and carrying current $I_1$ the vector potential is found from $\nabla^2 A_z = -\mu_0 I_1 \delta(\vec{r}-\vec{r}_1)$, where on the right hand side is the delta-function, and the solution is well known to be $A_z = -(\mu_0/2\pi) I_1 \ln(|\vec{r}-\vec{r}_1|)$, where the distance between the observation point and the wire is $|\vec{r}-\vec{r}_1|$.

Further it is convenient to use the polar coordinates $(\rho,\theta)$, thus a wire is described by the current $I_1$ and its location ($\rho_1,\theta_1$). Then at the observation point $(\rho,\theta)$, using the cosine theorem, we have

$ A_z = -\frac{\mu_0}{4\pi} I_1 \ln(\rho^2 + \rho_1^2 - 2 \rho \rho_1 \cos(\theta-\theta_1)) $

For the far field, assuming $\rho_1 \ll \rho$, we expand the expression as

$ A_z = -\frac{\mu_0}{4\pi} I_1 \left[ \ln(\rho^2) - \ln(1 - \frac{2 \rho_1}{\rho} \cos(\theta-\theta_1)) \right] \approx -\frac{\mu_0}{4\pi}I_1 \ln(\rho^2) -\frac{\mu_0}{4\pi} I_1 \frac{2 \rho_1}{\rho} \cos(\theta-\theta_1) $

The second term on the right-hand side is the dipole term which can be written as $ -(\mu_0/2\pi) \vec{p}_1 \cdot \vec{r} / r^2$, where the dipole moment is $\vec{p}_1 = I_1 \vec{r}_1$

Now, using superposition, we can immediately generalize for N parallel wires at locations described in the 2D plane by vectors $\vec{r}_k$ and carrying currents $I_k$. The far field at location $\vec{r}$, assuming $|r| \gg |r_k|$, is given by

$ A_z = -\frac{\mu_0}{2\pi} I_{tot} \ln(|r|) -\frac{\mu_0}{2\pi} (\vec{p}_{tot} \cdot \vec{r})/r^2 $

where the total current is $I_{tot} = \sum_{k=1}^N I_k$ and the total dipole moment is $\vec{p}_{tot} = \sum_{k=1}^N I_k \vec{r}_k$.

Of course the expansion procedure could be continued further and the exact answer would be an infinite sum of multipole moments; but for the far-away field only the lowest-order terms matter.

Back to the situation with two wires at distance $2a$ apart, the total current is zero so the first (monopole) term vanishes. The total dipole moment is $p_{tot} = 2 a I$ and is in the direction from one wire to the other one, in the plane perpendicular to the wires. The result for $A_z$ is very similar to the electric potential from a regular electric dipole; only the dependence on the distance is different.

Finally, the magnetic field is found from $A_z$ as

$ \vec{B} = \nabla \times \vec{A} = \hat{e}_z \times \nabla A_z $