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