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:

  • MultipoleInversion from the multipole_inversion.multipole_inversion library

  • MagneticSample from the multipole_inversion.magnetic_sample library

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#

[2]:
%matplotlib inline
[3]:
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
[4]:
# 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
[5]:
# 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?

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)
../_images/tutorial_tutorial_basics_15_1.png

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 0x7efd16bdb1f0>
../_images/tutorial_tutorial_basics_17_1.png

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 0x7efd10ac1d60>
../_images/tutorial_tutorial_basics_22_1.png
[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.

[15]:
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')
Scanning array size = 200 x 201
Generating forward matrix
Generation of Q matrix took: 0.1300 s
Using scipy.linalg.pinv for inversion
../_images/tutorial_tutorial_basics_25_1.png

We can also check the residual is small compared to the inverted field:

[16]:
f, ax = plt.subplots()
cf, c1 = minv.plot_tools.plot_difference_Bz(ax, qinv)
colorbar(cf)
ax.set_aspect('equal')
../_images/tutorial_tutorial_basics_27_0.png

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

[17]:
qinv.inv_multipole_moments
[17]:
array([[ 4.80000000e-13,  3.00753220e-30,  1.05818295e-29,
         1.42164522e-31, -2.74340099e-29, -3.90631829e-29,
         4.08979310e-30, -1.42587698e-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.

[6]:
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()
[7]:
# 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

[8]:
f, ax = plt.subplots()
cf, c1, c2 = sample.plot_sample(ax)
colorbar(cf)
ax.set_aspect('equal')
../_images/tutorial_tutorial_basics_36_0.png

Now we redefine the dipole_positions in the sample instance in order to make a single quadrupole source rather than two dipoles:

[9]:
# 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
[10]:
sample.save_data(filename='quadrupole_y-orientation')
[11]:
!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": "20220929-160148", "Scan origin x": 0.0, "Scan 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:

[22]:
qinv = MultipoleInversion('./MetaDict_quadrupole_y-orientation.json',
                          './MagneticSample_quadrupole_y-orientation.npz',
                          expansion_limit='quadrupole',
                          sus_functions_module='maxwell_cartesian_polynomials')
Scanning array size = 200 x 200
[23]:
qinv.compute_inversion(method='sp_pinv2')
Generating forward matrix
Generation of Q matrix took: 0.2340 s
Using scipy.linalg.pinv for 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:

[24]:
f, ax = plt.subplots()
cf, c1, c2 = minv.plot_tools.plot_inversion_Bz(ax, qinv)
colorbar(cf)
ax.set_aspect('equal')
../_images/tutorial_tutorial_basics_45_0.png

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.

[25]:
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:

[16]:
f, ax = plt.subplots()
cf, c1 = minv.plot_tools.plot_difference_Bz(ax, qinv)
colorbar(cf)
ax.set_aspect('equal')
../_images/tutorial_tutorial_basics_49_0.png

We can compare these results if we used the orthogonal basis given by spherical harmonics polynomials defined by Burnham and English:

[26]:
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')
Scanning array size = 200 x 200
Generating forward matrix
Generation of Q matrix took: 0.2705 s
Using scipy.linalg.pinv for inversion
[27]:
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]]
[28]:
print('Magnetization: ', np.linalg.norm(qinv.inv_multipole_moments[:3]) / 1e-18)
Magnetization:  91.23148267085703

Quadrupole x-direction#

We repeat the same calculations here, but setting the two dipoles in the \(x\)-direction

[34]:
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()
[35]:
f, ax = plt.subplots()
p, *_ = sample.plot_sample(ax)
colorbar(p)

ax.set_aspect('equal')
../_images/tutorial_tutorial_basics_57_0.png
[36]:
# 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')
[38]:
qinv = minv.MultipoleInversion('./MetaDict_quadrupole_x-orientation.json',
                               './MagneticSample_quadrupole_x-orientation.npz',
                               sus_functions_module='maxwell_cartesian_polynomials')
qinv.compute_inversion(method='sp_pinv2')
Scanning array size = 200 x 200
Generating forward matrix
Generation of Q matrix took: 0.1622 s
Using scipy.linalg.pinv for inversion

We again obtain the highest quadrupole for the \(Q_{xy}\) component

[39]:
qinv.inv_multipole_moments
[39]:
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)

[41]:
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()
[42]:
f, ax = plt.subplots()
p, *_ = sample.plot_sample(ax)
colorbar(p)

ax.set_aspect('equal')
../_images/tutorial_tutorial_basics_65_0.png
[43]:
# 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')
[44]:
qinv = minv.MultipoleInversion('./MetaDict_quadrupole_xy-orientation.json',
                               './MagneticSample_quadrupole_xy-orientation.npz',
                               expansion_limit='quadrupole')
qinv.compute_inversion(method='sp_pinv2')
Scanning array size = 200 x 200
Generating forward matrix
Generation of Q matrix took: 0.2268 s
Using scipy.linalg.pinv for inversion

Now the highest moments are the \(Q_{1}=Q_{11}=Q_{xx}\) and the \(Q_{4}=Q_{22}=Q_{yy}\) components

[45]:
qinv.inv_multipole_moments
[45]:
array([[ 5.40141200e-16, -5.40141200e-16,  2.79421928e-27,
        -1.13264993e-32, -2.09919164e-21,  2.09919164e-21,
         2.64190432e-18,  2.87650521e-32]])
[46]:
f, ax = plt.subplots()
cf, c1, c2 = minv.plot_tools.plot_inversion_Bz(ax, qinv)
colorbar(cf)
ax.set_aspect('equal')
../_images/tutorial_tutorial_basics_70_0.png

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.

[48]:
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()
[49]:
f, ax = plt.subplots()
p, *_ = sample.plot_sample(ax)
colorbar(p)

ax.set_aspect('equal')
../_images/tutorial_tutorial_basics_74_0.png
[50]:
# 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')
[51]:
# 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')
[55]:
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)
Scanning array size = 200 x 200
Generating forward matrix
Generation of Q matrix took: 0.2239 s
Using scipy.linalg.pinv for inversion
[56]:
f, ax = plt.subplots()
cf, c1, c2 = minv.plot_tools.plot_inversion_Bz(ax, oinv)
ax.set_aspect('equal')
colorbar(cf)
[56]:
<matplotlib.colorbar.Colorbar at 0x7efd13fbe0a0>
../_images/tutorial_tutorial_basics_78_1.png
[58]:
f, ax = plt.subplots()
cf, c1 = minv.plot_tools.plot_difference_Bz(ax, oinv)
ax.set_aspect('equal')
colorbar(cf)
[58]:
<matplotlib.colorbar.Colorbar at 0x7efd13e68bb0>
../_images/tutorial_tutorial_basics_79_1.png
[59]:
qinv.inv_multipole_moments
[59]:
array([[ 3.53852928e-28, -3.53854975e-28, -1.42195925e-35,
        -2.38432549e-33,  4.61478886e-23, -4.61478886e-23,
         4.21851577e-24, -3.03183503e-34,  1.63992684e-38,
        -3.03166486e-28,  3.03166486e-28, -1.54723019e-29,
         3.18554377e-39, -2.48650698e-27, -2.48650698e-27]])

Multiple particles sample#

[61]:
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')
[62]:
sample.Bz_array.shape
[62]:
(100, 100)
[63]:
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')
../_images/tutorial_tutorial_basics_84_0.png
[64]:
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()
[64]:
<matplotlib.legend.Legend at 0x7efd13ccfd00>
../_images/tutorial_tutorial_basics_85_1.png

Inversion#

[65]:
qinv = minv.MultipoleInversion('./MetaDict_seed42.json',
                               './MagneticSample_seed42.npz')
Scanning array size = 100 x 100
[68]:
qinv.compute_inversion(rcond=1e-10, method='np_pinv')
Using numpy.pinv for inversion
[69]:
qinv.inv_multipole_moments[:10]
[69]:
array([[-1.00275987e-11, -1.09300850e-11,  4.00512755e-12,
        -4.48719898e-17,  1.30387106e-17,  5.43854803e-17,
         2.94773220e-17, -3.13438666e-18],
       [-6.96637235e-12,  1.06269116e-12, -3.33727440e-12,
         1.75710238e-20, -1.00402229e-19, -1.10267772e-19,
         4.81786070e-21, -8.00986367e-20],
       [-7.18369963e-12, -1.19192257e-11, -6.45404249e-12,
         1.77954991e-18, -4.55280174e-20,  6.92800314e-19,
         2.48954902e-18, -2.32297886e-19],
       [ 1.00804803e-11, -4.10555104e-12, -6.68273378e-12,
         3.81543339e-19, -2.92837371e-18, -7.16780422e-19,
        -8.76046453e-19,  9.35422321e-19],
       [-6.43159997e-13, -5.85191859e-12, -1.34308431e-11,
        -1.20492365e-20, -7.22968400e-22, -1.35846767e-20,
        -3.86755605e-21,  3.73394350e-21],
       [ 8.87989770e-13,  1.84399215e-12, -6.35868735e-12,
        -1.13144600e-18, -3.89993965e-18, -7.34107560e-19,
         7.48552073e-20,  1.67301350e-19],
       [ 2.46072703e-12,  1.22550665e-12,  1.58046839e-11,
        -6.45820268e-18, -2.57599415e-19, -2.59381088e-18,
         3.87006865e-18, -1.70070822e-18],
       [-1.16912666e-12,  1.40821923e-11, -2.10319835e-12,
         1.66776910e-19,  1.79100282e-19,  1.77720034e-19,
        -1.84216176e-19, -1.50965374e-19],
       [-1.83151421e-11, -1.24620339e-11, -5.26831161e-12,
        -2.11857666e-21, -1.67303816e-20, -1.03491845e-21,
        -1.26636974e-21, -7.18513872e-23],
       [ 7.30050202e-12,  2.44423307e-12,  2.96386084e-12,
         2.90652343e-22, -1.68420131e-21,  6.48879963e-22,
        -8.82403559e-22, -3.49515465e-22]])
[71]:
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')
../_images/tutorial_tutorial_basics_90_0.png
[72]:
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()
[72]:
<matplotlib.legend.Legend at 0x7efd138c28b0>
../_images/tutorial_tutorial_basics_91_1.png
[74]:
f, ax = plt.subplots(figsize=(6, 6))
cf, c1 = minv.plot_tools.plot_difference_Bz(ax, qinv)
colorbar(cf)
ax.set_aspect('equal')
../_images/tutorial_tutorial_basics_92_0.png
[75]:
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)
[75]:
[<matplotlib.lines.Line2D at 0x7efd025e5d90>]
../_images/tutorial_tutorial_basics_93_1.png
[76]:
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
[77]:
# 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.8623405853795806e-09
[78]:
NRRMS = RRMS / (Bzinv_minus_Bz.max() - Bzinv_minus_Bz.min())
print(f'Normalised RRMS: {NRRMS:.4f} %')
Normalised RRMS: 0.0679 %
[79]:
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}}$')
[79]:
Text(0, 0.5, '$m_{i}^{\\mathrm{predicted}} - m_{i}^{\\mathrm{data}}$')
../_images/tutorial_tutorial_basics_97_1.png