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

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

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()
• Sounds like the goal here is to find a magnetic dipole field approximating the magnetic field from a given system of stationary electric currents. Instead of a numerical optimization procedure, the right way to do it is to apply multipole expansion to your coil current, that would produce the dipole field as the lowest-order term. – Maxim Umansky Dec 19 '20 at 1:52
• @MaximUmansky reading about multipole expansion, it sounds as if it was the standard tool for the job. Its new to me, though. Can you please help with the approximation at hand? – Andreas Schuldei Dec 19 '20 at 7:44
• For localized current it will boil down to taking an integral along the path of the current, as explained, e.g., in farside.ph.utexas.edu/teaching/jk1/lectures/node53.html. What is the current distribution? Is it a closed loop? – Maxim Umansky Dec 19 '20 at 16:33
• @MaximUmansky the two line currents (very long, parallel, just in one point) are distance 2a apart and are I1 = -I2. it can be considered a closed loop, but an infietly long one. – Andreas Schuldei Dec 19 '20 at 16:48
• If the currents are long but of finite length L then the dipole moment is $I \cdot 2a \cdot L$. If the length is infinite that is more subtle but probably some approximation can be worked out for the magnetic field using the symmetry of the system. – Maxim Umansky Dec 19 '20 at 18:32

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$$

• You write that the total current is the sum of all currents. Your vector potential is that of one current. With the two currents in opposit directions, the sum of the currents is zero and hence the dipol moment would be zero, too. Do I understand that correctly? – Andreas Schuldei Dec 25 '20 at 10:06
• No, the monopole moment (the one that gives $\ln(\rho)$ dependence for $A_z$) vanishes. Dipole moment is the second one that gives $1/\rho$ dependence. – Maxim Umansky Dec 25 '20 at 15:27
• working my way through this. Shouldn't it be $A_z = \frac{\mu I}{4\pi}ln(|r|)$ ? – Andreas Schuldei Dec 26 '20 at 8:57
• In SI units there will be $\mu_0$ there, in CGS there will be speed of light $c$. I'll fix the units asap. – Maxim Umansky Dec 26 '20 at 16:30