In [1]:
using PyPlot
using Coils
using Coils.CoilsPlot
The coil design method works on a predefined grid. The grid is defined by a set of points in space (vertices) and a directed graph g defining which are connected with one another. A connection between two vertices means, that a wire of a coil can be laid in a straight line between them.
It is absolutely possible to supply vertices (an array of arrays of x, y, z coordinates) and the graph manually. There is also fuction cuboid_system that defines a cuboid system, given it size (x, y and z) and the number of tiles along x, y and z directions.
In [2]:
g, vertex_positions = cuboid_system([1, 1, 1], [3, 3, 3])
Out[2]:
In CoilsPlot there is a number of functions for plotting the objects.
In [3]:
figure(figsize = (6, 6))
plot_vertices(vertex_positions)
plot_edges(g, vertex_positions, standalone = false)
The method optimizes the field in a finite set of Points Of Interest (poi). They can be supplied manually or a helper function cuboid_poi may be used. It has to be given the size, the position of the centre, the number of points (along each direction). The points may be created either in the whole volume (filled = true) or just on the surface.
In [4]:
poi = cuboid_poi([0.75, 0.75, 0.75], [0.0, 0.0, 0.0], [10, 10, 10], filled = false)
Out[4]:
In [5]:
figure(figsize = (6, 6))
plot_system(g, vertex_positions, poi)
Out[5]:
The goal field is defined as a function of space x (an array of size 3 - x, y and z), returning an array of size 3 - the x, y and z components of the field.
In [6]:
Bgoal(x) = [0, 100e-6, 0]
Out[6]:
Next, the graph is analysed to find cells (or tiles) - smallest loops to be found in the graph. Each will be treated as a tile coil.
In [7]:
cells = find_cells(g)
Out[7]:
In [8]:
figure(figsize = (8, 8))
plot_vertices(vertex_positions, labels = false)
plot_cells(cells, vertex_positions)
Now the matrix of the system is calculated. It is a matrix of proportionality constants between the current in each cell (columns) and the magnetic field in each of the points of interest, in each direction (rows). It is calculated using the Biot-Savart law.
In [9]:
M = system_matrix(poi, vertex_positions, cells)
Out[9]:
We evaluate the goal magnetic field in each of the points of interest and concatenate the values into one big vector.
In [10]:
Bpoi = vcat(Bgoal.(poi)...)
Out[10]:
Now we solve the system. We look for the optimal currents in each of the tile coils, that will reproduce the goal magnetic field in the points of interest best. It is a linear least-squares problem.
In [11]:
optI = M \ Bpoi
Out[11]:
In [12]:
figure(figsize = (10, 10))
plot_vertices(vertex_positions, labels = false)
current_norm = 1000 / maximum(abs, optI)
plot_cells(cells, vertex_positions, optI * current_norm)
We proceed to simplify the system. The first step is for each edge to to add the currents of the adjacent tiles. We get the net current flowing along each edge.
In [13]:
edgecurrents = find_edgecurrents(g, cells, optI)
Out[13]:
In [14]:
figure(figsize = (8, 8))
plot_vertices(vertex_positions, labels = false)
plot_edge_currents(g, edgecurrents, vertex_positions)
The result is a complicated net of current. We decompose the net into what we call simple loops.
In [15]:
simpleloops, simpleloopscurrents = find_all_simpleloops(g, edgecurrents)
Out[15]:
Now we plot the simple loops, together with the corresponding currents.
In [16]:
figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
plot_loop(simpleloops[1], vertex_positions, color = "orange")
plot_loop(simpleloops[2], vertex_positions, color = "orange")
title("Current: $(simpleloopscurrents[1] * current_norm)")
Out[16]:
In [17]:
figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
plot_loop(simpleloops[3], vertex_positions, color = "green")
plot_loop(simpleloops[4], vertex_positions, color = "green")
title("Current: $(simpleloopscurrents[3] * current_norm)")
Out[17]:
In [18]:
figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
plot_loop(simpleloops[5], vertex_positions, color = "blue")
plot_loop(simpleloops[6], vertex_positions, color = "blue")
title("Current: $(simpleloopscurrents[5] * current_norm)")
Out[18]:
In [19]:
figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
plot_loop(simpleloops[7], vertex_positions, color = "indigo")
plot_loop(simpleloops[8], vertex_positions, color = "indigo")
title("Current: $(simpleloopscurrents[7] * current_norm)")
Out[19]:
In [20]:
figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
plot_loop(simpleloops[9], vertex_positions, color = "grey")
plot_loop(simpleloops[10], vertex_positions, color = "grey")
title("Current: $(simpleloopscurrents[9] * current_norm)")
Out[20]:
We can characterise the performance of the design. For example by plotting the histogram of the difference of the goal field and the one generated with the design.
In [21]:
plot_deviation_histogram(poi, simpleloops, simpleloopscurrents, vertex_positions, Bgoal)
Out[21]:
Or we can plot the map of the difference between the goal field an the one produced by the system.
In [22]:
plot_loops_field(cells, optI, vertex_positions, [1, 2, 3], 0;
n = 50, spanA = [-0.5, 0.5], spanB = [-0.5, 0.5], Bref = Bgoal, levels = -10:10)
Out[22]:
The magnetic field can be decomposed into linearly independent cartesian harmonic polynomials. Below are the polynomials up to the first order:
n = 1: $\quad P^x_1(\mathbf{r}) = 1 \quad P^y_1(\mathbf{r}) = 0 \quad P^z_1(\mathbf{r}) = 0$
n = 2: $\quad P^x_2(\mathbf{r}) = 0 \quad P^y_2(\mathbf{r}) = 1 \quad P^z_2(\mathbf{r}) = 0$
n = 3: $\quad P^x_3(\mathbf{r}) = 0 \quad P^y_3(\mathbf{r}) = 0 \quad P^z_3(\mathbf{r}) = 1$
n = 4: $\quad P^x_4(\mathbf{r}) = x \quad P^y_4(\mathbf{r}) = 0 \quad P^z_4(\mathbf{r}) = -z$
n = 5: $\quad P^x_5(\mathbf{r}) = y \quad P^y_5(\mathbf{r}) = x \quad P^z_5(\mathbf{r}) = 0$
n = 6: $\quad P^x_6(\mathbf{r}) = 0 \quad P^y_6(\mathbf{r}) = y \quad P^z_6(\mathbf{r}) = -z$
n = 7: $\quad P^x_7(\mathbf{r}) = z \quad P^y_7(\mathbf{r}) = 0 \quad P^z_7(\mathbf{r}) = x$
n = 8: $\quad P^x_8(\mathbf{r}) = 0 \quad P^y_8(\mathbf{r}) = z \quad P^z_8(\mathbf{r}) = y$
We will design a $n=5$ gradient coil on a system that is extended in the $y$ direction and has two large openings.
We will also assume that to construct the coil we will use wires with 10A, 1A and 0.1A per 100 μT/m field.
In [23]:
g, vertex_positions = cuboid_system([1, 2, 1], [3, 6, 3],
skipfaces = [false, false, true, true, false, false])
poi = cuboid_poi([0.5, 0.5, 0.5], [0.0, 0.0, 0.0], [10, 10, 10], filled = true)
# the wires used to construct the coil
elemcurrents = [10.0, 1.0, 0.1]
figure(figsize = (6, 6))
plot_system(g, vertex_positions, poi)
Out[23]:
When there are large openings it takes a long time to find the cells automatically. It can be speeded up by providing the large cells manually.
In [24]:
# Find out the location of the openings
extrema(getindex.(vertex_positions, 2))
Out[24]:
In [25]:
# prepare the large cells
front_face_unordered = [ i for i in eachindex(vertex_positions) if vertex_positions[i][2] < -0.9 ]
front_face = order_cell(g, front_face_unordered)
back_face_unordered = [ i for i in eachindex(vertex_positions) if vertex_positions[i][2] > 0.9 ]
back_face = order_cell(g, back_face_unordered)
initialcells = [front_face, back_face]
Out[25]:
The goal field is a gradient coil, n = 5: $\quad P^x_5(\mathbf{r}) = y \quad P^y_5(\mathbf{r}) = x \quad P^z_5(\mathbf{r}) = 0$, magnitude is 100 μT/m.
In [26]:
Bgoal(x) = [√2 * 100e-6 * x[2], √2 * 100e-6 * x[1], 0]
Out[26]:
We will use the solve_system function, which performs the whole design.
In [27]:
simpleloops, simpleloopscurrents = solve_system(g, vertex_positions, poi, Bgoal,
verbose = true, initialcells = [front_face, back_face],
mincurrent = minimum(elemcurrents) / 2)
Out[27]:
We decompose the solution into the currents that we want to build. We also calculate the real currents that will be flowing in the simpleloops, when they are constructed out of 10A, 1A and 0.1A wires.
In [28]:
simpleloopscurrents_decomp = decompose_currents(simpleloopscurrents, elemcurrents)
simpleloopscurrents_real = real_decomposed_currents(simpleloopscurrents_decomp, elemcurrents)
Out[28]:
We calculate the performance of the system build using the 10A, 1A and 0.1A wires.
In [29]:
plot_deviation_histogram(poi, simpleloops, simpleloopscurrents_real, vertex_positions, Bgoal)
Out[29]:
In [30]:
plot_loops_field(simpleloops, simpleloopscurrents_real, vertex_positions, [1, 2, 3], 0;
n = 50, spanA = [-0.5, 0.5], spanB = [-1, 1], Bref = Bgoal, levels = -20:4:20)
Out[30]:
When we do not supply the reference field, we can plot the field of the coild.
In [31]:
plot_loops_field(simpleloops, simpleloopscurrents_real, vertex_positions, [1, 2, 3], 0;
n = 50, spanA = [-0.5, 0.5], spanB = [-1, 1], levels = -50:5:50)
Out[31]:
We print the decomposition of all the currents in the simple loops to the currents we want to use in the wires:
In [32]:
for i in eachindex(simpleloopscurrents)
print("#$(lpad(i,2)): ")
print("$(lpad(simpleloopscurrents[i], 20)) = ")
print("$(lpad(simpleloopscurrents_decomp[i][1], 2)) * 10A + ")
print("$(simpleloopscurrents_decomp[i][2]) * 1A + ")
println("$(simpleloopscurrents_decomp[i][3]) * 0.1A")
end
In [33]:
figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
for i in 1:8
plot_loop(simpleloops[i], vertex_positions, color = "C0")
end
title("Current: $(simpleloopscurrents[1])")
Out[33]:
In [34]:
figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
for i in 9:16
plot_loop(simpleloops[i], vertex_positions, color = "C1")
end
title("Current: $(simpleloopscurrents[9])")
Out[34]:
In [35]:
figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
for i in 17:20
plot_loop(simpleloops[i], vertex_positions, color = "C2")
end
title("Current: $(simpleloopscurrents[17])")
Out[35]:
In [36]:
figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
for i in 21:28
plot_loop(simpleloops[i], vertex_positions, color = "C3")
end
title("Current: $(simpleloopscurrents[21])")
Out[36]:
In [37]:
figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
for i in 29:30
plot_loop(simpleloops[i], vertex_positions, color = "C4")
end
title("Current: $(simpleloopscurrents[29])")
Out[37]:
In [38]:
figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
for i in 31:38
plot_loop(simpleloops[i], vertex_positions, color = "C5")
end
title("Current: $(simpleloopscurrents[31])")
Out[38]:
In [39]:
figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
for i in 39:42
plot_loop(simpleloops[i], vertex_positions, color = "C6")
end
title("Current: $(simpleloopscurrents[39])")
Out[39]:
In [40]:
figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
for i in 43:46
plot_loop(simpleloops[i], vertex_positions, color = "C7")
end
title("Current: $(simpleloopscurrents[43])")
Out[40]:
In [41]:
figure(figsize = (6, 6))
plot_vertices(vertex_positions, labels = false)
for i in 47:48
plot_loop(simpleloops[i], vertex_positions, color = "C8")
end
title("Current: $(simpleloopscurrents[47])")
Out[41]: