Tutorial: Basics¶
This tutorial shows how to perform a multipole inversion from a magnetic scan surface into one or more magnetic sources that are represented as physical point sources. The multipole_inversion library contains two main classes to perform the numerical inversions:
MultipoleInversionfrom themultipole_inversion.multipole_inversionlibraryMagneticSamplefrom themultipole_inversion.magnetic_samplelibrary
With the MagneticSample module it is possible to define the scan grid dimensions and the location of the magnetic sources. From this class we save this information in json and npz files that can be inputted into the MultipoleInversion class. Both of these classes have extensive docstrings that can be read in this notebook for more information on the input parameters/arguments.
Additional tools include plotting functions defined in the class libraries or in the multipole_inversion.plot_tools module.
Import and definitions¶
[1]:
%matplotlib inline
[2]:
import numpy as np
import matplotlib.pyplot as plt
# from palettable.cartocolors.diverging import Geyser_5
from mpl_toolkits.axes_grid1 import make_axes_locatable
import importlib as imp # to reload libraries if implementing new features
[3]:
# Load the libraries for the calculation of dipole fields
import mmt_multipole_inversion as minv
from mmt_multipole_inversion import MultipoleInversion
from mmt_multipole_inversion import MagneticSample
[4]:
# Define a colorbar for the plots
def colorbar(mappable, ax=None, location='right', size='5%', pad=0.05,
orientation='vertical', ticks_pos='right', **kwargs):
"""
Note: append_axes() reduces the size of ax to make room for the colormap
ticks_pos :: if orientation is vertical -> 'right' or 'left'
if orientation is horizontal -> 'top' or 'bottom'
"""
if not ax:
ax = plt.gca()
# fig = ax.figure
divider = make_axes_locatable(ax)
cax = divider.append_axes(location, size=size, pad=pad)
# cbar = fig.colorbar(mappable, cax=cax)
cbar = plt.colorbar(mappable, cax=cax, orientation=orientation,
**kwargs)
if orientation == 'vertical':
cax.yaxis.set_ticks_position(ticks_pos)
elif orientation == 'horizontal':
cax.xaxis.set_ticks_position(ticks_pos)
# plt.sca(ax)
return cbar
An example of docstring for the MultipoleInversion class:
[5]:
MultipoleInversion?
Init signature:
MultipoleInversion(
sample_config_file: Union[str, pathlib.Path],
sample_arrays: Union[str, pathlib.Path, NoneType],
expansion_limit: Literal['dipole', 'quadrupole', 'octupole'] = 'quadrupole',
sus_functions_module: Literal['spherical_harmonics_basis', 'maxwell_cartesian_polynomials', 'cartesian_spherical_harmonics', 'spherical_harmonics_basis_area', 'spherical_harmonics_basis_volume'] = 'spherical_harmonics_basis',
) -> None
Docstring:
Class to perform multipole inversions
Class to perform multipole inversions of a magnetic scan surface into
multiple magnetic sources located within a sample. Specifications of
the scan grid and the magnetic particles in the sample can be generated
using the `MagneticSample` class. The sensors of the magnetic scan surface
are modelled either as point sensors or sensor with cuboid shape (volume).
Init docstring:
Parameters
----------
sample_config_file
Path to a `json` file with the specifications of the scan grid and
the magnetic particles. The following keys are mandatory::
Scan height Hz
Scan area x-dimension Sx
Scan area y-dimension Sy
Scan x-step Sdx
Scan y-step Sdy
Number of particles
Time stamp
The following are optional::
Sensor origin x
Sensor origin y
Sensor dimensions
Sensor dimensions are required if a 2D or 3D sensor is used. The
sensor origins are the coordinates of the lower left sensor center,
in the scanning surface. By default it is `(0.0, 0.0)`. The sensor
grid in each dimension (`x` or `y`) is computed as multiples of
`Scan area / Scan step`
sample_arrays
An `npz` file containing the scan signal Bz and the particle
positions (magnetic sources). The file can contain other
information as well but it is not read here. If the two arrays are
not specified they can be set manually using the `Bz_array` and
`particle_positions` class variables.
expansion_limit
Higher order multipole term to compute the field contribution from
the potential of the magnetic particles. Options:
`dipole`, `quadrupole`, `octupole`
sus_functions_module
Spherical harmonic basis for the susceptibility matrix used for the
multipole inversion. The fully orthogonal and linearly independent
basis is the `spherical_harmonics_basis`. Other options are not
orthogonal but might be necessary for comparison. These modules
populate the forward matrix with the assumption that sensors are a
point, which can be seen as sensors with an area or volume where
the magnetic flux from the sources is treated as constant.
Alternatively, susceptibility modules ending with `_area` or
`_volume` model the sensors using a geometry with a higher
dimension, e.g. a rectangle, where the magnetic flux is integrated
within it. For details see the comments in the libraries in the
`sus_functions_module/` directory and the Notes.
Notes
-----
Mathematical details for the multipole inversion can be found in::
D. Cortés-Ortuño, K. Fabian, L. V. de Groot
Single Particle Multipole Expansions From Micromagnetic Tomography
G^3, 22(4), e2021GC009663 (2021)
https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2021GC009663
File: ~/git/mmt_multipole_inversion/mmt_multipole_inversion/multipole_inversion.py
Type: type
Subclasses:
Testing the dipole field function¶
Here we are testing the dipole field from the dipole.py library. We set the dipole as close to the origin as possible, and oriented in the \(\hat{x}\) direction.
[6]:
dip_r = np.array([[0., 0., -1e-10]])
dip_m = np.array([[1, 0., 0.]])
[7]:
# we set the space grid in the -1,1 range in both x and y directions
x = np.linspace(-1, 1, 150)
X, Y = np.meshgrid(x, x)
positions = np.column_stack([X.ravel(), Y.ravel(), np.zeros_like(X.ravel())])
[8]:
# Bz = msp.dipole_Bz(dip_r, dip_m, positions)
B = minv.magnetic_sample.dipole_field(dip_r, dip_m, positions)
We can now plot the dipole field around the origin, isolines are plotted for the \(B_z\) component
[9]:
plt.contour(X, Y, B[:, 2].reshape(-1, len(x)),
levels=[-1e-11, -1e-12, 1e-12, 1e-11],
colors='k', linestyles='-')
# plt.contourf(X, Y, B[:, 1].reshape(-1, len(x)).T)
# Normalise arrows
U, V = B[:, 0], B[:, 1]
norm = np.sqrt(U ** 2 + V ** 2)
# print(norm)
U, V = U / norm, V / norm
p = plt.quiver(positions[:, 0], positions[:, 1], U, V, norm, scale=25,
cmap='magma', width=.005, edgecolor='k', linewidth=.5)
plt.colorbar(p)
plt.scatter(dip_r[:, 0], dip_r[:, 1], c='C2', s=50)
plt.xlim(-0.085, 0.085)
plt.ylim(-0.085, 0.085)
[9]:
(-0.085, 0.085)
The same using streamlines:
[10]:
# Normalise arrows
U, V = B[:, 0], B[:, 1]
plt.streamplot(X, Y, B[:, 0].reshape(-1, len(x)), B[:, 1].reshape(-1, len(x)),
density=2, linewidth=1, color='k')
plt.scatter(dip_r[:, 0], dip_r[:, 1], c='C2', s=50)
# plt.xlim(-0.085, 0.085)
# plt.ylim(-0.085, 0.085)
[10]:
<matplotlib.collections.PathCollection at 0x7f75509b5640>
Inversion of a single dipole source¶
In this Section we test the inversion of the field flux of a single dipole measured at a surface located at \(H_z\) above a \(Lx\times Ly\times Lz\) rectangular sample region, which contains the sipole at its centre.
[11]:
Hz = 2e-6 # Scan height in m
Sx = 20e-6 # Scan area x - dimension in m
Sy = 20.1e-6 # Scan area y - dimension in m
Sdx = 0.1e-6 # Scan x - step in m
Sdy = 0.1e-6 # Scan y - step in m
Lx = Sx * 0.9 # Sample x - dimension in m
Ly = Sy * 0.9 # Sample y - dimension in m
Lz = 5e-6 # Sample thickness in m
# Initialise the dipole class
sample = MagneticSample(Hz, Sx, Sy, Sdx, Sdy, Lx, Ly, Lz)
# Generate two random particles in the sample region, which we are going to
# redefine (this might not be necessary, we need to add more methods to the class)
# sample.generate_particles(N_particles=1)
# Manually set the positions and magnetization of the two dipoles
Ms = 4.8e5
dipole_positions = np.array([[sample.Lx * 0.5, sample.Ly * 0.5, -sample.Lz * 0.5]])
magnetization = Ms * (1 * 1e-18) * np.array([[1., 0., 0.]])
volumes = np.array([1e-18])
sample.generate_particles_from_array(dipole_positions, magnetization, volumes)
print('Magnetization:', sample.dipole_moments)
# Generate the dipole field measured as the Bz field flux through the
# measurement surface
sample.generate_measurement_mesh()
Magnetization: [[4.8e-13 0.0e+00 0.0e+00]]
We can visualise the field generated by the dipole at the measurement surface. The colormap refers to the \(B_z\) flux, and we add a streamplot to observe the dipole field direction \((B_x, B_y)\) at the surface.
[12]:
f, ax = plt.subplots()
cf, c1, c2 = sample.plot_sample(ax)
colorbar(cf)
ax.set_aspect('equal')
c2.set_color('C3')
# Streamplot: take the measurement surface range and generate a regular
# rectangular mesh grid. We take these mesh points to compute the field in them
x, y = sample.Sx_range, sample.Sy_range
X, Y = np.meshgrid(x, y)
positions = np.column_stack([X.ravel(), Y.ravel(), sample.Hz * np.ones_like(X.ravel())])
B = minv.magnetic_sample.dipole_field(sample.dipole_positions, sample.dipole_moments, positions)
# Generate random seed points from where streamlines emerge (density not
# necessary if random seeds are used)
seed_points_x = sample.Lx * np.random.random(150)
seed_points_y = sample.Ly * np.random.random(150)
ax.streamplot(X, Y,
B[:, 0].reshape(-1, len(x)),
B[:, 1].reshape(-1, len(x)),
density=1.5, linewidth=0.5, color='k',
start_points=np.column_stack((seed_points_x, seed_points_y)),
# start_points=[[0, 0]],
)
# plt.xlim(5e-6, 15e-6)
# plt.ylim(5e-6, 15e-6)
[12]:
<matplotlib.streamplot.StreamplotSet at 0x7f754e26b9e0>
[13]:
sample.save_data(filename='dipole_y-orientation')
Now we use the inv_quadrupole.py library to load the dipole field data and inverse the signal into the particle position, which gives us the dipole and quadrupole moments. The latter should be close to zero.
In the plot we observe that the inverte dsignal reproduces the original field accurately.
[14]:
qinv = MultipoleInversion('./MetaDict_dipole_y-orientation.json',
'./MagneticSample_dipole_y-orientation.npz',
expansion_limit='quadrupole')
qinv.compute_inversion(method='sp_pinv2', atol=1e-25)
f, ax = plt.subplots()
cf, c1, c2 = minv.plot_tools.plot_inversion_Bz(ax, qinv)
colorbar(cf)
ax.set_aspect('equal')
2025-09-08 16:51:29,044 | INFO | multipole_inversion :: 🧲 Initialized new MultipoleInversion simulation
2025-09-08 16:51:29,045 | WARNING | multipole_inversion :: Parameter "Sensor dimensions" not found in json file. Setting class attribute sensor_dims to default ()
2025-09-08 16:51:29,045 | INFO | multipole_inversion :: Scanning array sizes (row x col)
2025-09-08 16:51:29,046 | INFO | multipole_inversion :: Computed Sx x Sy sizes : 201 x 200
2025-09-08 16:51:29,047 | INFO | multipole_inversion :: Scan positions array memory: 0.9201 Mb
2025-09-08 16:51:29,049 | INFO | multipole_inversion :: Bz data array size : 201 x 200
2025-09-08 16:51:29,049 | INFO | multipole_inversion :: Bz data memory : 0.3067 Mb
2025-09-08 16:51:29,049 | WARNING | multipole_inversion :: Bz array set. Setting raw fieldMask array
2025-09-08 16:51:29,051 | INFO | multipole_inversion :: Generating forward matrix
2025-09-08 16:51:29,051 | INFO | multipole_inversion :: Green matrix memory: 2.4536 Mb
2025-09-08 16:51:31,045 | INFO | multipole_inversion :: Generation of Q matrix took: 1.9942 s
2025-09-08 16:51:31,047 | INFO | multipole_inversion :: Using scipy.linalg.pinv for inversion
2025-09-08 16:51:31,058 | INFO | multipole_inversion :: Finished inversion
We can also check the residual is small compared to the inverted field:
[15]:
f, ax = plt.subplots()
cf, c1 = minv.plot_tools.plot_difference_Bz(ax, qinv)
colorbar(cf)
ax.set_aspect('equal')
The inverted magnetic moments (3 dipole moments and 5 quadrupole moments) shows us the original dipole moment in the \(x\)-direction. The other moments are significantly small with respecto to \(m_x\)
[16]:
qinv.inv_multipole_moments
[16]:
array([[ 4.80000000e-13, 1.20724993e-29, 2.43596434e-29,
1.42226394e-31, -2.74334360e-29, -3.90631751e-29,
4.08983467e-30, -1.42569894e-29]])
Quadrupole¶
Here we define a quadrupole by specifying two magnetic dipoles oriented in opposite directions and located close to the center of the sample. Before saving the data, we redefine the two dipoles as a single particle at the center of the sample. The purpose of this idea is to analyse the strength of a magnetic quadrupole when solving the inversion problem. Accordingly, the dipole moments should be close to zero and one or more quadrupole moments should be stronger.
Quadrupole y-direction¶
The first example is a quadrupole oriented in the \(+\hat{y}\) and \(-\hat{y}\) directions, located at 4.5 micrometers from the measurement surface of the sample.
[17]:
Hz = 2e-6 # Scan height in m
Sx = 20e-6 # Scan area x - dimension in m
Sy = 20e-6 # Scan area y - dimension in m
Sdx = 0.1e-6 # Scan x - step in m
Sdy = 0.1e-6 # Scan y - step in m
Lx = Sx * 0.9 # Sample x - dimension in m
Ly = Sy * 0.9 # Sample y - dimension in m
Lz = 5e-6 # Sample thickness in m
# Initialise the dipole class
sample = MagneticSample(Hz, Sx, Sy, Sdx, Sdy, Lx, Ly, Lz)
# Manually set the positions and magnetization of the two dipoles
Ms = 4.8e5
dipole_positions = np.array([[sample.Lx * 0.5 - 1e-6, sample.Ly * 0.5, -sample.Lz * 0.5],
[sample.Lx * 0.5 + 1e-6, sample.Ly * 0.5, -sample.Lz * 0.5]])
magnetization = Ms * (1 * 1e-18) * np.array([[0., 1., 0], [0., -1, 0]])
volumes = np.array([1e-18, 1e-18])
sample.generate_particles_from_array(dipole_positions, magnetization, volumes)
# Generate the dipole field measured as the Bz field flux through the
# measurement surface
sample.generate_measurement_mesh()
[18]:
# DEBUG:
# pos_r = np.array([sample.Sx_range[len(sample.Sx_range) // 2],
# sample.Sy_range[len(sample.Sy_range) // 2],
# sample.Hz])
# r = pos_r - sample.dipole_positions
# print(r)
Here we show the dipole field at the measurement surface, generated from the two dipoles in the sample region
[19]:
f, ax = plt.subplots()
cf, c1, c2 = sample.plot_sample(ax)
colorbar(cf)
ax.set_aspect('equal')
Now we redefine the dipole_positions in the sample instance in order to make a single quadrupole source rather than two dipoles:
[20]:
# Hack the positions array making a single particle at the centre
# (ideal quadrupole)
sample.dipole_positions = np.array([[sample.Lx * 0.5, sample.Ly * 0.5, -sample.Lz * 0.5]])
# This magnetisation direction should not matter (?)
sample.magnetization = Ms * (1 * 1e-18) * np.array([[0., 1., 0]])
# Update the N of particles to update the internal dict
sample.N_particles = 1
[21]:
sample.save_data(filename='quadrupole_y-orientation')
[22]:
!cat MetaDict_dipole_y-orientation.json
{"Scan height Hz": 2e-06, "Scan area x-dimension Sx": 2e-05, "Scan area y-dimension Sy": 2.01e-05, "Scan x-step Sdx": 1e-07, "Scan y-step Sdy": 1e-07, "Time stamp": "20250908-165127", "Sensor origin x": 0.0, "Sensor origin y": 0.0, "Number of particles": 1}
At this point we can load the data for the inversion of the measurement generated in the previous steps, in the inversion code/class. Notice we are going to use the maxwell_cartesian_polynomials as a basis for the multipole expansion in order to physically interepret the results. Strictly, this basis is not orthogonal so it is not the most robust basis for the expansion. Nevertheless, for a single magnetic source it should solve the problem:
[23]:
qinv = MultipoleInversion('./MetaDict_quadrupole_y-orientation.json',
'./MagneticSample_quadrupole_y-orientation.npz',
expansion_limit='quadrupole',
sus_functions_module='maxwell_cartesian_polynomials')
2025-09-08 16:51:31,799 | INFO | multipole_inversion :: 🧲 Initialized new MultipoleInversion simulation
2025-09-08 16:51:31,800 | WARNING | multipole_inversion :: Parameter "Sensor dimensions" not found in json file. Setting class attribute sensor_dims to default ()
2025-09-08 16:51:31,801 | INFO | multipole_inversion :: Scanning array sizes (row x col)
2025-09-08 16:51:31,801 | INFO | multipole_inversion :: Computed Sx x Sy sizes : 200 x 200
2025-09-08 16:51:31,802 | INFO | multipole_inversion :: Scan positions array memory: 0.9155 Mb
2025-09-08 16:51:31,804 | INFO | multipole_inversion :: Bz data array size : 200 x 200
2025-09-08 16:51:31,804 | INFO | multipole_inversion :: Bz data memory : 0.3052 Mb
2025-09-08 16:51:31,805 | WARNING | multipole_inversion :: Bz array set. Setting raw fieldMask array
[24]:
qinv.compute_inversion(method='sp_pinv2')
2025-09-08 16:51:31,808 | INFO | multipole_inversion :: Generating forward matrix
2025-09-08 16:51:31,809 | INFO | multipole_inversion :: Green matrix memory: 2.4414 Mb
2025-09-08 16:51:34,057 | INFO | multipole_inversion :: Generation of Q matrix took: 2.2477 s
2025-09-08 16:51:34,058 | INFO | multipole_inversion :: Using scipy.linalg.pinv for inversion
2025-09-08 16:51:34,069 | INFO | multipole_inversion :: Finished inversion
We can compute the inverted measurement grid to compare it with the original measurement grid. We also notice we have now a single particle at the centre of the sample:
[25]:
f, ax = plt.subplots()
cf, c1, c2 = minv.plot_tools.plot_inversion_Bz(ax, qinv)
colorbar(cf)
ax.set_aspect('equal')
From the inversion we can now check the magnitude of the inverted multipole moments. In this case, the inverted magnetization array is
\[\texttt{inv magnetization} = [m_x, m_y, m_z, Q_1, Q_2, \ldots, Q_5]\]
where \(m_i\) are the dipole moments and \(Q_i\) are the quadrupole moments.
In this example, we can see that \(Q_2=Q_{xy}\) has the highest magnitude among the quadrupoles, and the magnitude is around \(10^{-18}\) (check units). The dipole moments should be around \(\approx 10^{-12}\) if we had only a dipolar field.
[26]:
with np.printoptions(precision=2, suppress=False):
print(qinv.inv_multipole_moments)
[[-4.00e-17 -8.20e-17 6.61e-19 3.46e-25 -1.38e-18 1.11e-22 2.09e-24
2.27e-22]]
And finally the difference between the measured field \(B_z\) and the field from the inversion. The residual has an octupole character:
[27]:
f, ax = plt.subplots()
cf, c1 = minv.plot_tools.plot_difference_Bz(ax, qinv)
colorbar(cf)
ax.set_aspect('equal')
We can compare these results if we used the orthogonal basis given by spherical harmonics polynomials defined by Burnham and English:
[28]:
qinv = MultipoleInversion('./MetaDict_quadrupole_y-orientation.json',
'./MagneticSample_quadrupole_y-orientation.npz',
expansion_limit='quadrupole',
sus_functions_module='spherical_harmonics_basis')
qinv.compute_inversion(method='sp_pinv2')
2025-09-08 16:51:34,458 | INFO | multipole_inversion :: 🧲 Initialized new MultipoleInversion simulation
2025-09-08 16:51:34,459 | WARNING | multipole_inversion :: Parameter "Sensor dimensions" not found in json file. Setting class attribute sensor_dims to default ()
2025-09-08 16:51:34,459 | INFO | multipole_inversion :: Scanning array sizes (row x col)
2025-09-08 16:51:34,459 | INFO | multipole_inversion :: Computed Sx x Sy sizes : 200 x 200
2025-09-08 16:51:34,461 | INFO | multipole_inversion :: Scan positions array memory: 0.9155 Mb
2025-09-08 16:51:34,462 | INFO | multipole_inversion :: Bz data array size : 200 x 200
2025-09-08 16:51:34,462 | INFO | multipole_inversion :: Bz data memory : 0.3052 Mb
2025-09-08 16:51:34,462 | WARNING | multipole_inversion :: Bz array set. Setting raw fieldMask array
2025-09-08 16:51:34,464 | INFO | multipole_inversion :: Generating forward matrix
2025-09-08 16:51:34,464 | INFO | multipole_inversion :: Green matrix memory: 2.4414 Mb
2025-09-08 16:51:34,502 | INFO | multipole_inversion :: Generation of Q matrix took: 0.0374 s
2025-09-08 16:51:34,502 | INFO | multipole_inversion :: Using scipy.linalg.pinv for inversion
2025-09-08 16:51:34,520 | INFO | multipole_inversion :: Finished inversion
[29]:
with np.printoptions(precision=2, suppress=False):
print(qinv.inv_multipole_moments)
[[-4.00e-17 -8.20e-17 6.61e-19 -2.98e-24 1.57e-22 3.21e-22 -1.23e-24
-1.95e-18]]
[30]:
print('Magnetization: ', np.linalg.norm(qinv.inv_multipole_moments[:3]) / 1e-18)
Magnetization: 91.23148267088258
Quadrupole x-direction¶
We repeat the same calculations here, but setting the two dipoles in the \(x\)-direction
[31]:
Hz = 2e-6 # Scan height in m
Sx = 20e-6 # Scan area x - dimension in m
Sy = 20e-6 # Scan area y - dimension in m
Sdx = 0.1e-6 # Scan x - step in m
Sdy = 0.1e-6 # Scan y - step in m
Lx = Sx * 0.9 # Sample x - dimension in m
Ly = Sy * 0.9 # Sample y - dimension in m
Lz = 5e-6 # Sample thickness in m
sample = MagneticSample(Hz, Sx, Sy, Sdx, Sdy, Lx, Ly, Lz)
# Manually set positions
Ms = 4.8e5
dipole_positions = np.array([[sample.Lx * 0.5, sample.Ly * 0.5 - 1e-6, -sample.Lz * 0.5],
[sample.Lx * 0.5, sample.Ly * 0.5 + 1e-6, -sample.Lz * 0.5]])
magnetization = Ms * (1 * 1e-18) * np.array([[-1., 0., 0], [1., 0, 0]])
volumes = np.array([1e-18, 1e-18])
sample.generate_particles_from_array(dipole_positions, magnetization, volumes)
sample.generate_measurement_mesh()
[32]:
f, ax = plt.subplots()
p, *_ = sample.plot_sample(ax)
colorbar(p)
ax.set_aspect('equal')
[33]:
# Now hack the positions array making a single particle at the centre
# (ideal quadrupole)
sample.dipole_positions = np.array([[sample.Lx * 0.5, sample.Ly * 0.5, -sample.Lz * 0.5]])
# This magnetisation direction should not matter (?)
sample.magnetization = Ms * (1 * 1e-18) * np.array([[1., 0., 0]])
sample.N_particles = 1 # Need to modify the JSON file!
sample.save_data(filename='quadrupole_x-orientation')
[34]:
qinv = minv.MultipoleInversion('./MetaDict_quadrupole_x-orientation.json',
'./MagneticSample_quadrupole_x-orientation.npz',
expansion_limit='quadrupole',
sus_functions_module='maxwell_cartesian_polynomials')
qinv.compute_inversion(method='sp_pinv2')
2025-09-08 16:51:34,768 | INFO | multipole_inversion :: 🧲 Initialized new MultipoleInversion simulation
2025-09-08 16:51:34,768 | WARNING | multipole_inversion :: Parameter "Sensor dimensions" not found in json file. Setting class attribute sensor_dims to default ()
2025-09-08 16:51:34,769 | INFO | multipole_inversion :: Scanning array sizes (row x col)
2025-09-08 16:51:34,769 | INFO | multipole_inversion :: Computed Sx x Sy sizes : 200 x 200
2025-09-08 16:51:34,770 | INFO | multipole_inversion :: Scan positions array memory: 0.9155 Mb
2025-09-08 16:51:34,771 | INFO | multipole_inversion :: Bz data array size : 200 x 200
2025-09-08 16:51:34,771 | INFO | multipole_inversion :: Bz data memory : 0.3052 Mb
2025-09-08 16:51:34,772 | WARNING | multipole_inversion :: Bz array set. Setting raw fieldMask array
2025-09-08 16:51:34,773 | INFO | multipole_inversion :: Generating forward matrix
2025-09-08 16:51:34,773 | INFO | multipole_inversion :: Green matrix memory: 2.4414 Mb
2025-09-08 16:51:34,812 | INFO | multipole_inversion :: Generation of Q matrix took: 0.0388 s
2025-09-08 16:51:34,813 | INFO | multipole_inversion :: Using scipy.linalg.pinv for inversion
2025-09-08 16:51:34,834 | INFO | multipole_inversion :: Finished inversion
We again obtain the highest quadrupole for the \(Q_{xy}\) component
[35]:
qinv.inv_multipole_moments
[35]:
array([[ 8.19900901e-17, 3.99570646e-17, -6.61299707e-19,
-2.08881684e-24, 1.37917152e-18, -2.27255428e-22,
-3.45930322e-25, -1.11141388e-22]])
Quadrupole xy-direction¶
And the calculation for the dipoles in the \(xy\) direction (\(\phi = \pi/4\) in polar)
[36]:
Hz = 2e-6 # Scan height in m
Sx = 20e-6 # Scan area x - dimension in m
Sy = 20e-6 # Scan area y - dimension in m
Sdx = 0.1e-6 # Scan x - step in m
Sdy = 0.1e-6 # Scan y - step in m
Lx = Sx * 0.9 # Sample x - dimension in m
Ly = Sy * 0.9 # Sample y - dimension in m
Lz = 5e-6 # Sample thickness in m
sample = MagneticSample(Hz, Sx, Sy, Sdx, Sdy, Lx, Ly, Lz)
# Manually set positions
Ms = 4.8e5
dipole_positions = np.array([[sample.Lx * 0.5 - 1e-6, sample.Ly * 0.5 - 1e-6, -sample.Lz * 0.5],
[sample.Lx * 0.5 + 1e-6, sample.Ly * 0.5 + 1e-6, -sample.Lz * 0.5]])
n = np.sqrt(2)
magnetization = Ms * (1 * 1e-18) * np.array([[-1 / n, 1 / n, 0],
[1 / n, -1 / n, 0]])
volumes = np.array([1e-18, 1e-18])
sample.generate_particles_from_array(dipole_positions, magnetization, volumes)
sample.generate_measurement_mesh()
[37]:
f, ax = plt.subplots()
p, *_ = sample.plot_sample(ax)
colorbar(p)
ax.set_aspect('equal')
[38]:
# Now hack the positions array making a single particle at the centre
# (ideal quadrupole)
sample.dipole_positions = np.array([[sample.Lx * 0.5, sample.Ly * 0.5, -sample.Lz * 0.5]])
# This magnetisation direction should not matter (?)
sample.magnetization = Ms * (1 * 1e-18) * np.array([[-1 / n, 1 / n, 0]])
sample.N_particles = 1
sample.save_data(filename='quadrupole_xy-orientation')
[39]:
qinv = minv.MultipoleInversion('./MetaDict_quadrupole_xy-orientation.json',
'./MagneticSample_quadrupole_xy-orientation.npz',
sus_functions_module='maxwell_cartesian_polynomials',
expansion_limit='quadrupole')
qinv.compute_inversion(method='sp_pinv2')
2025-09-08 16:51:35,084 | INFO | multipole_inversion :: 🧲 Initialized new MultipoleInversion simulation
2025-09-08 16:51:35,084 | WARNING | multipole_inversion :: Parameter "Sensor dimensions" not found in json file. Setting class attribute sensor_dims to default ()
2025-09-08 16:51:35,085 | INFO | multipole_inversion :: Scanning array sizes (row x col)
2025-09-08 16:51:35,085 | INFO | multipole_inversion :: Computed Sx x Sy sizes : 200 x 200
2025-09-08 16:51:35,086 | INFO | multipole_inversion :: Scan positions array memory: 0.9155 Mb
2025-09-08 16:51:35,087 | INFO | multipole_inversion :: Bz data array size : 200 x 200
2025-09-08 16:51:35,087 | INFO | multipole_inversion :: Bz data memory : 0.3052 Mb
2025-09-08 16:51:35,088 | WARNING | multipole_inversion :: Bz array set. Setting raw fieldMask array
2025-09-08 16:51:35,089 | INFO | multipole_inversion :: Generating forward matrix
2025-09-08 16:51:35,090 | INFO | multipole_inversion :: Green matrix memory: 2.4414 Mb
2025-09-08 16:51:35,128 | INFO | multipole_inversion :: Generation of Q matrix took: 0.0380 s
2025-09-08 16:51:35,129 | INFO | multipole_inversion :: Using scipy.linalg.pinv for inversion
2025-09-08 16:51:35,139 | INFO | multipole_inversion :: Finished inversion
Now the highest moments are the \(Q_{1}=Q_{11}=Q_{xx}\) and the \(Q_{4}=Q_{22}=Q_{yy}\) components
[40]:
qinv.inv_multipole_moments
[40]:
array([[ 5.40141200e-16, -5.40141200e-16, -7.35806220e-28,
1.86810846e-18, 1.59916575e-32, -1.48435264e-21,
-1.86810846e-18, 1.48435264e-21]])
[41]:
f, ax = plt.subplots()
cf, c1, c2 = minv.plot_tools.plot_inversion_Bz(ax, qinv)
colorbar(cf)
ax.set_aspect('equal')
Octupole¶
In this case we generate an artificial octupole using four dipoles. Below it can be seen that the inversion fails to produce an optimal solution for the octupole. It is possible to obtain a better solution by defining the octupole using two quadrupoles with two point sources, rather than a single point source.
[42]:
Hz = 2e-6 # Scan height in m
Sx = 20e-6 # Scan area x - dimension in m
Sy = 20e-6 # Scan area y - dimension in m
Sdx = 0.1e-6 # Scan x - step in m
Sdy = 0.1e-6 # Scan y - step in m
Lx = Sx * 0.9 # Sample x - dimension in m
Ly = Sy * 0.9 # Sample y - dimension in m
Lz = 5e-6 # Sample thickness in m
sample = MagneticSample(Hz, Sx, Sy, Sdx, Sdy, Lx, Ly, Lz)
# Manually set positions
Ms = 4.8e5
dipole_positions = np.array([[sample.Lx * 0.5 + 1e-6, sample.Ly * 0.5 + 1e-6, -sample.Lz * 0.5],
[sample.Lx * 0.5 - 1e-6, sample.Ly * 0.5 + 1e-6, -sample.Lz * 0.5],
[sample.Lx * 0.5 - 1e-6, sample.Ly * 0.5 - 1e-6, -sample.Lz * 0.5],
[sample.Lx * 0.5 + 1e-6, sample.Ly * 0.5 - 1e-6, -sample.Lz * 0.5]])
n = np.sqrt(2)
magnetization = Ms * (1 * 1e-18) * np.array([[-1 / n, 1 / n, 0],
[-1 / n, -1 / n, 0],
[1 / n, -1 / n, 0],
[1 / n, 1 / n, 0]])
volumes = np.array([1e-18, 1e-18, 1e-18, 1e-18])
sample.generate_particles_from_array(dipole_positions, magnetization, volumes)
sample.generate_measurement_mesh()
[43]:
f, ax = plt.subplots()
p, *_ = sample.plot_sample(ax)
colorbar(p)
ax.set_aspect('equal')
[44]:
# Now hack the positions array making a single particle at the centre
# (ideal octupole)
sample.dipole_positions = np.array([[sample.Lx * 0.5, sample.Ly * 0.5, -sample.Lz * 0.5]
])
# This magnetisation direction should not matter (?)
sample.magnetization = Ms * (1 * 1e-18) * np.array([[-1 / n, 1 / n, 0]])
sample.N_particles = 1 # Need to modify the JSON file!
sample.save_data(filename='octupole')
[45]:
# We could try to use 2 quadrupoles as well:
# # Now hack the positions array making a single particle at the centre
# # (ideal octupole)
# sample.dipole_positions = np.array([[sample.Lx * 0.5 + 1e-6, sample.Ly * 0.5 + 1e-6, -sample.Lz * 0.5],
# [sample.Lx * 0.5 - 1e-6, sample.Ly * 0.5 - 1e-6, -sample.Lz * 0.5]
# ])
# # This magnetisation direction should not matter (?)
# sample.magnetization = Ms * (1 * 1e-18) * np.array([[-1 / n, 1 / n, 0],
# [1 / n, -1 / n, 0]
# ])
# sample.N_particles = 2
# sample.save_data(filename='octupole')
[46]:
oinv = minv.MultipoleInversion('./MetaDict_octupole.json',
'./MagneticSample_octupole.npz',
expansion_limit='octupole',
sus_functions_module='spherical_harmonics_basis'
)
oinv.compute_inversion(method='sp_pinv2', atol=1e-20)
2025-09-08 16:51:35,528 | INFO | multipole_inversion :: 🧲 Initialized new MultipoleInversion simulation
2025-09-08 16:51:35,529 | WARNING | multipole_inversion :: Parameter "Sensor dimensions" not found in json file. Setting class attribute sensor_dims to default ()
2025-09-08 16:51:35,529 | INFO | multipole_inversion :: Scanning array sizes (row x col)
2025-09-08 16:51:35,530 | INFO | multipole_inversion :: Computed Sx x Sy sizes : 200 x 200
2025-09-08 16:51:35,531 | INFO | multipole_inversion :: Scan positions array memory: 0.9155 Mb
2025-09-08 16:51:35,532 | INFO | multipole_inversion :: Bz data array size : 200 x 200
2025-09-08 16:51:35,532 | INFO | multipole_inversion :: Bz data memory : 0.3052 Mb
2025-09-08 16:51:35,532 | WARNING | multipole_inversion :: Bz array set. Setting raw fieldMask array
2025-09-08 16:51:35,534 | INFO | multipole_inversion :: Generating forward matrix
2025-09-08 16:51:35,534 | INFO | multipole_inversion :: Green matrix memory: 4.5776 Mb
2025-09-08 16:51:37,814 | INFO | multipole_inversion :: Generation of Q matrix took: 2.2795 s
2025-09-08 16:51:37,816 | INFO | multipole_inversion :: Using scipy.linalg.pinv for inversion
2025-09-08 16:51:37,840 | INFO | multipole_inversion :: Finished inversion
[47]:
f, ax = plt.subplots()
cf, c1, c2 = minv.plot_tools.plot_inversion_Bz(ax, oinv)
ax.set_aspect('equal')
colorbar(cf)
[47]:
<matplotlib.colorbar.Colorbar at 0x7f754cd705f0>
[48]:
f, ax = plt.subplots()
cf, c1 = minv.plot_tools.plot_difference_Bz(ax, oinv)
ax.set_aspect('equal')
colorbar(cf)
[48]:
<matplotlib.colorbar.Colorbar at 0x7f7550e6c110>
[49]:
qinv.inv_multipole_moments
[49]:
array([[ 5.40141200e-16, -5.40141200e-16, -7.35806220e-28,
1.86810846e-18, 1.59916575e-32, -1.48435264e-21,
-1.86810846e-18, 1.48435264e-21]])
Multiple particles sample¶
[50]:
Hz = 5e-6 # Scan height in m
Sx = 200e-6 # Scan area x - dimension in m
Sy = 300e-6 # Scan area y - dimension in m
Sdx = 2e-6 # Scan x - step in m
Sdy = 3e-6 # Scan y - step in m
Lx = Sx * 0.9 # Sample x - dimension in m
Ly = Sy * 0.9 # Sample y - dimension in m
Lz = 30e-6 # Sample thickness in m
sample = MagneticSample(Hz, Sx, Sy, Sdx, Sdy, Lx, Ly, Lz)
sample.generate_random_particles(seed=42)
# print(sample.dipole_positions)
sample.generate_measurement_mesh()
sample.save_data(filename='seed42')
[51]:
sample.Bz_array.shape
[51]:
(100, 100)
[52]:
f, ax = plt.subplots(figsize=(6, 6))
cf, c1, c2 = sample.plot_sample(ax)
c2.set_sizes(c2.get_sizes() / 4)
colorbar(cf)
ax.set_aspect('equal')
[53]:
plt.plot(sample.dipole_moments[:, 0], 'o', label=r'$m_x$')
plt.plot(sample.dipole_moments[:, 1], 'v', label=r'$m_y$')
plt.plot(sample.dipole_moments[:, 2], 's', label=r'$m_z$')
plt.legend()
[53]:
<matplotlib.legend.Legend at 0x7f7550e6fbf0>
Inversion¶
[54]:
qinv = minv.MultipoleInversion('./MetaDict_seed42.json',
'./MagneticSample_seed42.npz')
2025-09-08 16:51:38,874 | INFO | multipole_inversion :: 🧲 Initialized new MultipoleInversion simulation
2025-09-08 16:51:38,875 | WARNING | multipole_inversion :: Parameter "Sensor dimensions" not found in json file. Setting class attribute sensor_dims to default ()
2025-09-08 16:51:38,875 | INFO | multipole_inversion :: Scanning array sizes (row x col)
2025-09-08 16:51:38,876 | INFO | multipole_inversion :: Computed Sx x Sy sizes : 100 x 100
2025-09-08 16:51:38,876 | INFO | multipole_inversion :: Scan positions array memory: 0.2289 Mb
2025-09-08 16:51:38,878 | INFO | multipole_inversion :: Bz data array size : 100 x 100
2025-09-08 16:51:38,878 | INFO | multipole_inversion :: Bz data memory : 0.0763 Mb
2025-09-08 16:51:38,878 | WARNING | multipole_inversion :: Bz array set. Setting raw fieldMask array
[55]:
qinv.compute_inversion(rcond=1e-10, method='np_pinv')
2025-09-08 16:51:38,882 | INFO | multipole_inversion :: Generating forward matrix
2025-09-08 16:51:38,883 | INFO | multipole_inversion :: Green matrix memory: 61.0352 Mb
2025-09-08 16:51:38,942 | INFO | multipole_inversion :: Generation of Q matrix took: 0.0584 s
2025-09-08 16:51:38,942 | INFO | multipole_inversion :: Using numpy.pinv for inversion
2025-09-08 16:51:40,055 | INFO | multipole_inversion :: Finished inversion
[56]:
qinv.inv_multipole_moments[:10]
[56]:
array([[-1.00275988e-11, -1.09300849e-11, 4.00512749e-12,
-4.48719900e-17, 1.30387095e-17, 5.43854790e-17,
2.94773226e-17, -3.13438694e-18],
[-6.96637233e-12, 1.06269115e-12, -3.33727443e-12,
1.75710002e-20, -1.00402300e-19, -1.10267875e-19,
4.81796358e-21, -8.00986253e-20],
[-7.18369943e-12, -1.19192258e-11, -6.45404262e-12,
1.77955087e-18, -4.55294651e-20, 6.92800738e-19,
2.48954897e-18, -2.32297867e-19],
[ 1.00804804e-11, -4.10555107e-12, -6.68273379e-12,
3.81543973e-19, -2.92837420e-18, -7.16779564e-19,
-8.76046515e-19, 9.35422636e-19],
[-6.43159998e-13, -5.85191859e-12, -1.34308431e-11,
-1.20492735e-20, -7.22946120e-22, -1.35846748e-20,
-3.86754647e-21, 3.73394883e-21],
[ 8.87989998e-13, 1.84399235e-12, -6.35868718e-12,
-1.13144737e-18, -3.89994261e-18, -7.34107707e-19,
7.48578428e-20, 1.67300667e-19],
[ 2.46072815e-12, 1.22550665e-12, 1.58046833e-11,
-6.45819817e-18, -2.57603206e-19, -2.59381034e-18,
3.87006504e-18, -1.70070991e-18],
[-1.16912658e-12, 1.40821922e-11, -2.10319842e-12,
1.66777301e-19, 1.79099948e-19, 1.77720275e-19,
-1.84215598e-19, -1.50965932e-19],
[-1.83151421e-11, -1.24620339e-11, -5.26831162e-12,
-2.11852558e-21, -1.67303589e-20, -1.03494297e-21,
-1.26636608e-21, -7.18386313e-23],
[ 7.30050202e-12, 2.44423307e-12, 2.96386084e-12,
2.90640436e-22, -1.68418980e-21, 6.48873385e-22,
-8.82406587e-22, -3.49512049e-22]])
[57]:
f, ax = plt.subplots(figsize=(6, 6))
cf, c1, c2 = minv.plot_tools.plot_inversion_Bz(ax, qinv)
c2.set_sizes(c2.get_sizes() / 4)
colorbar(cf)
ax.set_aspect('equal')
[58]:
plt.plot(qinv.inv_multipole_moments[:, 0], 'o', label=r'$m_x$')
plt.plot(qinv.inv_multipole_moments[:, 1], 'v', label=r'$m_y$')
plt.plot(qinv.inv_multipole_moments[:, 2], 's', label=r'$m_z$')
plt.legend()
[58]:
<matplotlib.legend.Legend at 0x7f759c777680>
[59]:
f, ax = plt.subplots(figsize=(6, 6))
cf, c1 = minv.plot_tools.plot_difference_Bz(ax, qinv)
colorbar(cf)
ax.set_aspect('equal')
[60]:
f, ax = plt.subplots(figsize=(6, 6))
ax.plot(qinv.Bz_array.flatten() - qinv.inv_Bz_array.flatten(),
'o', label=r'$m_x$', ms=1)
[60]:
[<matplotlib.lines.Line2D at 0x7f754e000320>]
[61]:
print(f'B_z (max) = {np.max(qinv.Bz_array)} | B_z (min) = {np.min(qinv.Bz_array)}')
B_z (max) = 0.0033845869641065352 | B_z (min) = -0.004771340039837247
[62]:
# Residual root mean square
L = len(qinv.Bz_array.flatten())
Bzinv_minus_Bz = (qinv.Bz_array.flatten() - qinv.inv_Bz_array.flatten())
RRMS = np.sqrt(np.sum(Bzinv_minus_Bz ** 2)) / np.sqrt(L)
print(RRMS)
1.8623405813152847e-09
[63]:
NRRMS = RRMS / (Bzinv_minus_Bz.max() - Bzinv_minus_Bz.min())
print(f'Normalised RRMS: {NRRMS:.4f} %')
Normalised RRMS: 0.0679 %
[64]:
f, ax = plt.subplots(figsize=(6, 6))
ax.plot(qinv.inv_multipole_moments[:, 0] - qinv.dipole_moments[:, 0], 'o-', label=r'$m_x$')
ax.plot(qinv.inv_multipole_moments[:, 1] - qinv.dipole_moments[:, 1], 'v-', label=r'$m_y$')
ax.plot(qinv.inv_multipole_moments[:, 2] - qinv.dipole_moments[:, 2], 's-', label=r'$m_z$')
ax.legend()
ax.set_ylabel(r'$m_{i}^{\mathrm{predicted}} - m_{i}^{\mathrm{data}}$')
[64]:
Text(0, 0.5, '$m_{i}^{\\mathrm{predicted}} - m_{i}^{\\mathrm{data}}$')
