Design Optimal Discrete Magnetic Coils with Wolfram Language

So, hello, my name is Noah Hardwick. I work at Wolfram Research in the. In the UX team, and I did my undergraduate physics degree at the University of Nottingham, and I still do some research with one of the groups there. And so I want to talk today about a project I've been doing with them, and that is designing optimal discrete magnetic coils with the Wolfram language. Okay, so I'll be talking about, first of all, why is this an interesting problem? Why do we want to kind of tackle this? We'll look at how to express a magnetic field mathematically, and then we'll look at how to optimize these discrete coils to minimize errors in the field. And we have packaged up all of this functionality into a packlet called Create Coil. So I will talk a bit. I'll finish by sort of going over that. Okay, so why discrete call? So, our goal is to optimally produce a desired magnetic field profile. And, you know, over 100 years ago, we were doing this by sort of directly solving the Biot Savart law for. For a pair of loops. And with that, we could generate a constant field or a linear gradient field. But optimizing beyond that was quite difficult. And then continuum coils came along, and essentially you are modeling a continuous current flow on some surface, like the surface of a cylinder. And then you are discretizing that current flow with these sort of complex, spatially varying wire patterns. And so that's what you can see here. But because of the way these are fabricated, they produce very high quality fields. But because of the way they're fabricated, they also have to be very large. And so they're good for something like an MRI machine where you have lots of space. But there are lots of applications, especially to sort of like fundamental physics experiments where you want something much, much smaller. Now, discrete coils that they're much easier to fabricate, they're less prone to breaking, and then can be far more compact. And so the question is, can we revisit the discrete coils and better optimize those? So, firstly, you know, how do we express a magnetic field? Well, the most common way to do this, to express different magnetic field profiles in free space, is to app. Apply the basis of spherical harmonics. And essentially, these are functions that are using spherical symmetry to provide the simplest possible encoding of all possible magnetic field profiles that you could generate. And so here are. These happen to be the real spherical harmonics. And yeah, an individual harmonic is described by an order, an order N and degree M, and is denoted BNM. So for example this is, you know, the B4,3 harmonic. So the lower order field harmonics, they have lower sort of spatial, lower frequency spatial oscillations. And therefore they can be expressed in Cartesian coordinates in a straightforward way. Right. So a uniform field, uniform BZ field along the Z axis is equivalent to the B10 field harmonic system. So that's order 1 and degree 0. And so we can use the harmonic field plot function from Create Coil, remember to load the package there to visualize the Cartesian components of this field along each Cartesian axis. And so we're interested in the BZ component along the z axis. And we can see that that green line there is straight as a constant field. So the sort of, the most widely employed setup for generating this constant field is a Helmholtz pair. And that's actually the call I was talking about earlier. And it's very simply just a single pair of loops which each carry the same current and then they're separated by their radius. And so we can use the loop coilplot3d function from CreateCoil to visualize this. And we'll need to give it some specific parameters for it to draw. But yeah, so these are the arguments that essentially describe the coil. And we can see what that coil looks like. Okay, so let's see what this field looks like, the field generated by this coil. And so we can use loop field plots to do that. And that's going to return the Cartesian field components of the coils field along each Cartesian axis. And so again, if we look at BZ along the z axis, we can see that it's, yeah, it's, you know, it's a fairly uniform field. We can also use Loop Field Plot 2D in Create Coil. Again taking the same parameters to describe the coil to inspect the BZ fields over the XZ plane, we take the third part of its output to get bz. And what this is actually doing is plotting. It's working out where the field is deviating. And so this inner dashed line within that the field is deviating by, by less than 1%. So we can see that, you know, there's this, there's this sort of area of high quality low deviation fields in the center of the coil. But as you move away, the field starts to deteriorate. And so if we want to optimize these fields, then we need to, you know, we need to figure out why those, why the field deteriorates away from the center of the coil. So coil primitives. And so from, for example, a cool Primitive in this case would be a pair of loops. They generate a cascading series of field harmonics, and so specifically a pair of symmetric loops where the current is flowing in the same direction in each loop that generates all of the odd ordered harmonics of 0 degree. And so the field B is composed of the sum of B10 plus B30 plus B5,0, etc. And each of these harmonics has a weighting, so the field harmonic BNM has a weighting fnm. So let's look at some of these harmonics along the z axis. Again, we can use harmonic fields plot to do that. And what we see is that as the order of these harmonics increases, the harmonics, they have steeper gradients, but they also diminish more rapidly towards the center of the coil. And so of course, what that means is that the field at the center of the coil is composed predominantly of lower order harmonics. And therefore, if we want to generate B10, which was the constant field, then we want to choose a separation of those loop pairs that nulls B30. So just a reminder, here's our field is composed. So if we want to null b30, what we want to solve is for its weighting to equal zero, so F30. Now, those weightings, they're essentially functions of the separation of the loop pairs, now expressed as derivatives with respect to that separation. It's the normalized separation chi C, where the separation has been normalized by the core radius. And so this is the expression for zero degree harmonic weighting. We have an academic paper which I'll link to at the end, which goes over like where these come from, the derivation of these. But let's just use it for now. And so we can, we can have a look at what the B30 weighting looks like. And it's this expression here. And then we want to solve that equal to zero. And this is, it gives the result that we were expecting, which is that, so that's the normalized separation above and below the axis. And remember that the coil is actually symmetric, so those two solutions are physically the same. Okay, so a Helmholtz pair only has one degree of freedom, right? That separation. And therefore we can null one unwanted variation. But to null the higher order error harmonics, we require more degrees of freedom. So in other words, we need to add more pairs of loops. So let's start with just adding one extra pair of loops. So we have two pairs in total, which will keep things sort of easy to visualize to begin with. And so we still desire that constant field B10. But now, rather than just nulling B30, we're going to null B30 and B50. So this is our objective here. You know, we want to find values of the separations chi C1 and chi C2, for which the weightings F30 and F50 are nulled. So each, each core primitive, so each pair of loops contributes to those weightings. So we now have total weightings. And so this is therefore the expression that we need to find roots of. And so we solve this by setting it equal to zero, and we are unable to find any analytic solutions. Solve, can't do this. Now, in general, an analytic solution, you know, may not exist, but even if it did, it might describe a coil that is sort of, you know, like, physically impractical. It might be very difficult to fabricate, or it might be very long or something like that. So, you know, if we want a more robust, robust method for finding optimal coils, we've kind of got to change our approach a bit. So we change to a numerical approach, and our objective now becomes we want to find. So we still have the two pairs of loops, so we still have two degrees of freedom to play with. But now we're going to purposefully under determine the system. And so we're only going to solve F3. Zero is equal to zero. So we're only going to completely null one harmonic. And so we therefore have one equation and two variables. So the solutions are going to lie on a 1D contour in 2D space. And then if we can find some solutions on that contour, we can then rank them by how Effectively they minimize F5.0. So that's the weighting of then the first leading order error harmonic. So now our system's null is just that, the, you know, the B30 weighting contribution from each of those loop pairs. And now we can just use findroot with an initial guess to find a numerical solution to this. And obviously note that, you know, because the location of that solution contour, some roots, some seeds, aren't going to converge. So if we want to kind of do this efficiently and find lots of points on the solution contour so that we can find the best coils, our approach is going to be we want to coarsely mesh our search space so that we can just kind of find a few points on the contour, get a coarse mesh of the contour, and then we can interpolate those points and then sort them. So first we want to locate the solution contour. So I'm going to Coarsely seed find root over separations over this range between 0.1 and 2 for sake of argument. And I also want to return some other information here to plot. So I can see which seed converge to which solution. But the key thing here is that I'm just making an array, a coarse array over that range and I'm applying findroot to each point. Okay, so we can plot this and so we can see we've got our coarse mesh of seeds in blue and we have our solutions, our coarse mesh of the solution contour in red. And then we can see which seeds converge to which solutions. So the first thing to notice here is that this contour is symmetric along this diagonal and that's expected given that the coil is actually symmetric. So we can actually enforce that Chi C1 is less than Chi C2 and we're not going to lose any solutions. And so let's do that. We'll call it filtered data. And then great. So now we have a coarse mesh of the solution contour and so we can now interpolate it. So in the general case we have to use like a slightly more complicated algorithm, which I'll talk about later. But in this simple 2D case, actually Mathematica's built in meshing functionality is sufficient to do this. And so I won't go through all the steps here, but you know, we're essentially using functions like reconstruction, mesh, mesh and discretized region to follow this process here. So we have our filtered points, we then connect them up into a connected mesh, we can interpolate that mesh and then finally we just apply find root to all those points to make sure that they are actually lying on the solution contour. Okay, great. So now that we have our solutions on the contour, we can rank them. And so all of those solutions are going to completely null B30, which means B50 is the leading order error harmonic. We want to minimize that and we want to maximize our desired harmonic B10. And therefore we sort by the ratio of those weightings, right? So F10 over F50, that's what this expression here is. And then we can in this case reverse sort the solutions by that value. And yeah, we have. That gives us here are the three best coils and they have a, you know, they have a desired error ratio of about A. Okay, so now we can use loop field plot to visualize these, to compare these fields. And then I'm just adding this option here to make sure that they're plotting essentially the same range and that the color bar is the same. And so here's our, here's the coil that we just optimized with two pairs of loops, and here's the Helmholtz pair. And we can see that our optimal coil, it is better because we can see, you know, this 1% deviation region on the left is larger, but it's only slightly larger. And also the coil itself is physically bigger. So we're sort of cheating a bit. And so we should, you know, it is a little bit better, but we should really be able to do better than this. So one degree of freedom that we haven't considered yet is the turn ratios of the coil primitives. And so, for example, you know, previously we kind of assumed that we had unitary current flowing in every pair of loops. But of course you could wind those, you could wind each pair of loops a different number of times. And so if you do that, each chord primitive contributes a harmonic weighting that is scaled by its turn ratio, I.e. the I chi there. Okay, so with that in mind, let's try and generate B10 again. So that's the constant Z field, but using three loop pairs. And so we want to find again, we're going to slightly under determine the problem. And so we have three separations, but we're only going to null two leading order error harmonics. So that's B30 and B50. And then rather than doing all that optimization manually ourselves, we can just use findloopcoil in createcoil to do this. And it's basically following the method that we had. So let's start by just using unitary turn ratios. And it doesn't actually find any solutions. Now if we were looking at higher order harmonics, then they have high frequency spatial oscillations. And so it could be that our coarse mesh is just a bit too coarse coarse to find any sort of like, you know, pieces of contour. But here we're looking at, you know, we're looking at sort of n equals 1. And so what this actually means is that there just isn't a piece of solution contour in the region that we're looking at between 0.1 and 1. And so we can try a different set of turn ratios, in this case 1, minus 2 and 3. And what that does is it's going to bring a solution contour into view. And so it has succeeded and it's returned the best 10 cores it's found. And now the desired to error ratio is 14,000. So this is looking much better than before. And remember as well that this is, this is now the ratio of B10 to B70 rather than B50. Okay, so now we can play the same game. We'll use loop field plot 2D to again compare this optimal 3 loop coil to the Helmholtz pair. And we see that it has a much larger, you know, 1% deviation region, is actually about twice as long along the z axis, and is relatively compact as well. So the maximum separation here is about 0.7. And so this is a far more optimal coil than the one on the right. Okay, so we can inspect the solution contour. We're expecting a 1D solution contour in 3D this time by just returning all the solutions found. And because it calculates that essentially that optimality ratio for each one, we can color the solution by how optimal the solution is by the optimality of those solutions. And so this is what we see. And we've got our optimal solution in the middle there. So let's say we couldn't find a solution. We can try a number of things, but one of them is to further under determine the problem. And so rather than nulling two leading order error harmonics, we can just null one. And so we'd only null B30. And what that means is that B50 becomes the leading order error harmonic. And our solutions are now going to lie on a 2D surface in 3D. And so we can do this by manually specifying the harmonics to null in fine loop core with this option here. Excuse me. And then we can plot these solutions. And what we see is that as expected, we now have a 2D surface of solutions. And our previous 1D contour of solutions sort of glows with optimality in this surface. Okay, so how do we do the interpolation in high dimensions? Well, we still coarsely mesh the space to try and find some points on the solution contour. And so once we've done that, let's say our solution contour, we know it's going to be two dimensional for each point on the contour. So let's consider this point in red here. I'm looking at its nearest neighbors and then picking the and drawing a vector from the point to each of its nearest neighbors. And then picking the the closest nearest neighbour so that the shortest vector, and then I'm discarding the vector that is least orthogonal to the one I just picked. Then of the ones that remain which are in black, I choose the one. So I pick the one that is most orthogonal to the one I picked last, and then I discard the one that's least orthogonal. And so in two dimensions I just have to do this once. But if you had like, I don't know, like a 4D contour or something, then you'd obviously just have to keep repeating that until you get four basis vectors. So we orthonormalize those vectors to give us a local basis at that point, and then we use that local basis to expand an array of interpolation points. And then we give those defined root to make sure that they are lying on the contour. Okay, so so far we've only looked at loops. And because loops are physically, you know, they have complete azimuthal symmetry, which means their degree is zero. They only generate field harmonics that are also complete, that have also have complete azimuthal symmetry and are of degree zero. And so for example, you can see here, so that is a slice right through the XZ plane, but you can see that going around the azimuth, there's no periodicity there. So what if we actually do want to generate tesseral fields, in other words, fields with azimuthal periodicity. Well, then we can use different call primitives. And there's a whole bunch of logic as to how these are put together. But I'm just going to show some basic examples. You can have saddle primitives and you can have ellipse primitives. And so now m is greater than zero. So here's just an example. So this saddle primitive here is going to generate so this saddle primitive and then also this ellipse primitive are going to generate this series of field harmonics. And now you can see that there is some azimuthal periodicity here. And again with an M2 example, the primitives get more complicated and the fields become more tesseract as well, as you can see on this one here. Okay, so I'll finish by looking a bit at Create Coil because it kind of packages up. Hopefully it's a nice package for if you need to kind of, you know, make a small discrete coil for a certain system, then hopefully it provides some nice functionality to do that. So it has 18 fully documented functions in it. They all have options to control, you know, the optimization process or how you're visualizing what if you're visualizing, be it a coil or a field, we've also tried to include kind of helpful error messages. So, you know, here if you're trying to specify a non physical coil, you know, it's going to try and help you and kind of explain, you know, why that coil is non physical and where you're going wrong. And yeah, so let's, let's just have a look at some of, some of the things it can do. So the scope of create coils. So one of the issues is that for ellipse and saddle primitives, it can get. As the number of primitives increase and especially as the order increases as well. So the degree of the, the harmonics, it can be quite difficult to kind of distinguish different primitives. And so we've got some nice functionality. We can say, you know, plot style goes to automatic on The Coil Plot 3D functions to color each primitive differently so we can kind of very clearly see which wire belongs to which primitive. There's a 2D version of these functions as well. So it's as if you sort of unwrapped that coil. And so here I'm looking at a saddle coil. And so this is obviously useful if you want to actually, you know, fabricate one of these. And if you hover over call primitive, if you hover over a line, it's going to show all the other wires that are part of that same call primitive. And what you might not be able to see here, because it's a bit small, is that the tooltip contains all the parameters that describe that primitive. And by the way, I should mention that also works in the 3D case as well. You can do the same thing. So, yeah, we've got some nice field inspection features as well. So we've already seen the Loop Field plot 2D and, and the Lips and Saddlefield plot 2D as well. But by default it plots deviations that are 1% and 5%. I can specify whatever deviations I want it to plot here. And by the way, I can't remember if I mentioned this before, but this is all being calculated. It's just a numerical bias of R method. And I can specify the style of those deviations as well. So rather than 1% to 5%, I'm plotting 0.11 and 10% and colors in that way. Also note that this is now a gradient field coil and Loopfield Plot 2D is kind of clever enough to know what gradient you'd be interested in. And so rather than plotting just the BX deviation, it's now plotting the deviation from the gradient, right? So it's plotting the deviation in partial BX by partial X and here it's plotting partial BZ by partial Z partial. Now, if you have a harmonic like B30, for example, where you've got a different gradient, you've got a gradient along each axis. So for example, along the Z axis here I've got a gradient and then I've also got another gradient along the X axis. When I Give this to loopfieldplot2d for the z case, it's actually going to return me two plots, right? So here I got three plots, back here I'm getting four. And the reason is it knows that I'm kind of interested in both of those gradients. And so, for example, on the left here, it's Giving me the partial 2Bz by partial X squared deviation. And here it's giving again for the BZ field, it's giving the partial 2 BZ by partial Z squared deviation. So it's kind of clever enough to sort of work out what you'd want to. What are the interesting deviations to look at. So if you want to explore more, we have our academic paper and that goes through the derivation of all the maths behind this and sort of summarizes as well what I've said here. There is also a seminal paper by Romeo and holt from the 1980s on discrete call optimization. And hopefully you'll be able to find crate call useful as well. And so you can, you can find that on. And finally, I'd just like to thank my friend and colleague from Nottingham, Peter Hobson, who is the joint author on this paper with me and worked with me on this, and also Mark Fromhold from Nottingham, who leads that research group. Thanks. Thanks very much. Okay, let's see, are there questions? There is a question. Okay, the question is for those aspects of your model where you needed to estimate phase parameters of each component. What is the algorithm you use to normalize phase values to bring them into the range of 0 to 2 PI? Also related, sorry, related, is what is your algorithm for estimating parameter fits again, where the phase values need to be expressed in the interval 0 to 2 PI? The. So. The way all those sorts of things are normalized, I probably don't want to go into that here, but I would direct you, to direct you to this paper because it sort of has the derivations for all of those sort of written that much clearly than I can kind of explain here. One thing I would say is that we derive some of our expressions slightly differently to the way Romeo and Holt derived theirs. So that's. If you're looking at both of these papers, that's just something to be aware of because it can be confusing if you're not already aware of that. Sorry, that's a slightly dodgy question there, but I feel like you'll get a clear answer from just kind of looking at this paper. Okay. There are no other questions for now, if you're welcome. If. Yeah, if anybody has any more questions, please, please send me a direct message. And, yeah, thank you very. Much.