passive_shimming/simulation_data_test.py at main · imr-framework/passive_shimming · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
import numpy as np
import magpylib as magpy
from make_shim_rings import make_shim_ring_template
from utils import get_field_pos, display_scatter_3D, get_magnetic_field, load_magnets_in_rings
from colorama import Style, Fore
from target_B0_2_shim_locations_v2 import shimming_problem
from pymoo.core.mixed import MixedVariableMating, MixedVariableGA, MixedVariableSampling, MixedVariableDuplicateElimination
from pymoo.algorithms.moo.nsga2 import NSGA2, RankAndCrowdingSurvival
from pymoo.optimize import minimize
from pymoo.core.population import Population
from pymoo.core.evaluator import Evaluator
import matplotlib.pyplot as plt

# Read magnetic field and positions
fname = './fmr_data/calibration_jig/test_data.npy'
data = np.load(fname)
resolution = 2 #mm
x, y, z, B = get_field_pos(data)
print(Fore.GREEN + 'Done reading data')
x = (np.float64(x).transpose() - 0.5 * np.max(x))  * 1e-3 #conversion to m
y = (np.float64(y).transpose() - 0.5 * np.max(y)) * 1e-3 #conversion to m
z = (np.float64(z).transpose() - 0.5 * np.max(z)) * 1e-3 #conversion to m
B = B * 1e-3 # mT to T

# Map robot space to magpy space
x_magpy = z # length
y_magpy = x # depth
z_magpy = -y # height


# Display measured field as scattered data - plot3
display_scatter_3D(x_magpy, y_magpy, z_magpy, B, center=False)
print(Fore.CYAN + 'Mean B0 before shimming is:' + str(np.round(np.mean(B),2)) + ' mT') # What decimal should we round off to? 1mT - 85kHz
pos = np.zeros((x.shape[0], 3))
pos[:, 0] = x_magpy
pos[:, 1] = y_magpy
pos[:, 2] = z_magpy

dsv_sensors = magpy.Collection(style_label='sensors')
sensor1 = magpy.Sensor(position=pos,style_size=2)
dsv_sensors.add(sensor1)
print(Fore.GREEN + 'Done creating position sensors')

# Specify geometry of the magnet calibration jig
magnet_dims_x = 6.35 *1e-3 # m
magnet_dims_y = 6.35 *1e-3 # m
magnet_dims_z = 6.35 *1e-3 * 0.5 # m
magnet_dims = (magnet_dims_x, magnet_dims_y, magnet_dims_z)
diameter = 60 * 1e-3 # m
magnetization = [0, 0 ,8* 1e5] # 1.34, 0.7957
magnetization2 = [0, 0 , 8 * 1e5] # 1.34, 0.7957
position1 = np.multiply([0, 0, 8 + (0.5 * magnet_dims_z)], 1e-3)
position2 = np.multiply([0, 0, -8 - (0.5 * magnet_dims_z)], 1e-3)

# Make the jig
magnet_calib_jig = magpy.Collection()
cube1 = magpy.magnet.Cuboid(
                        dimension=magnet_dims,
                        # polarization=self.polarization,
                        position=position1,
                        # polarization=polarization,
                        magnetization=magnetization,
                        style_color='red',
                    )
magnet_calib_jig.add(cube1)
cube2 = magpy.magnet.Cuboid(
                        dimension=magnet_dims,
                        # polarization=self.polarization,
                        position=position2,
                        # polarization=polarization,
                        magnetization=magnetization2,
                        style_color='red',
                    )

magnet_calib_jig.add(cube2)
magnet_calib_jig.show()
B_sim = get_magnetic_field(magnet_calib_jig, dsv_sensors, axis = 2)
display_scatter_3D(x_magpy, y_magpy, z_magpy, B_sim, center=False)

# Figure how to export this to CAD
plt.plot(B, label = 'Measured')
plt.plot(B_sim, label = 'Simulated')
plt.legend()
plt.show()

# Visualize differences
display_scatter_3D(x_magpy, y_magpy, z_magpy, (B - B_sim), center=False)