# MagWire¶

As a first example how to use the modules and concepts that we have seen so far, we will develop two python classes that allow the calculation of magnetic fields generated by an arbitrary eletric current geometry. This is done by exploiting the law of Biot-Savart.

The source code including some examples is available here.

The frist class Wire is found in wire.py. It represents a continuous path for an electrical current (i.e. a wire). The path is specified as a sequence of 3D coordinates. The member property discretized_path returns the original path with linearly interpolated segments shorter than the specified discretization length.

  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 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 __author__ = 'wack' # part of the magwire package # calculate magnetic fields arising from electrical current through wires of arbitrary shape # with the law of Biot-Savart # written by Michael Wack # wack@geophysik.uni-muenchen.de # tested with python 3.4.3 from copy import deepcopy import numpy as np try: import visvis as vv visvis_avail = True except ImportError: visvis_avail = False print("visvis not found.") class Wire: ''' represents an arbitrary 3D wire geometry ''' def __init__(self, current=1, path=None, discretization_length=0.01): ''' :param current: electrical current in Ampere used for field calculations :param path: geometry of the wire specified as path of n 3D (x,y,z) points in a numpy array with dimension n x 3 length unit is meter :param discretization_length: lenght of dL after discretization ''' self.current = current self.path = path self.discretization_length = discretization_length @property def discretized_path(self): ''' calculate end points of segments of discretized path approximate discretization lenghth is given by self.discretization_length elements will never be combined elements longer that self.dicretization_length will be divided into pieces :return: discretized path as m x 3 numpy array ''' try: return self.dpath except AttributeError: pass self.dpath = deepcopy(self.path) for c in range(len(self.dpath)-2, -1, -1): # go backwards through all elements # length of element element = self.dpath[c+1]-self.dpath[c] el_len = np.linalg.norm(element) npts = int(np.ceil(el_len / self.discretization_length)) # number of parts that this element should be split up into if npts > 1: # element too long -> create points between # length of new sub elements sel = el_len / float(npts) for d in range(npts-1, 0, -1): self.dpath = np.insert(self.dpath, c+1, self.dpath[c] + element / el_len * sel * d, axis=0) return self.dpath @property def IdL_r1(self): ''' calculate discretized path elements dL and their center point r1 :return: numpy array with I * dL vectors, numpy array of r1 vectors (center point of element dL) ''' npts = len(self.discretized_path) if npts < 2: print("discretized path must have at least two points") return IdL = np.array([self.discretized_path[c+1]-self.discretized_path[c] for c in range(npts-1)]) * self.current r1 = np.array([(self.discretized_path[c+1]+self.discretized_path[c])*0.5 for c in range(npts-1)]) return IdL, r1 def vv_plot_path(self, discretized=True, color='r'): if not visvis_avail: print("plot path works only with visvis module") return if discretized: p = self.discretized_path else: p = self.path vv.plot(p, ms='x', mc=color, mw='2', ls='-', mew=0) def mpl3d_plot_path(self, discretized=True, show=True, ax=None, plt_style='-r'): if ax is None: fig = plt.figure(None) ax = ax3d.Axes3D(fig) if discretized: p = self.discretized_path else: p = self.path ax.plot(p[:, 0], p[:, 1], p[:, 2], plt_style) ax.set_xlabel('X') ax.set_ylabel('Y') ax.set_zlabel('Z') # make all axes the same #max_a = np.array((p[:, 0], p[:, 1], p[:, 2])).max() #ax.set_xlim3d(min(p[:, 0]), max_a) #ax.set_ylim3d(min(p[:, 1]), max_a) #ax.set_zlim3d(min(p[:, 2]), max_a) if show: plt.show() return ax def ExtendPath(self, path): ''' extends existing path by another one :param path: path to append ''' if self.path is None: self.path = path else: # check if last point is identical to avoid zero length segments if self.path[-1] == path[0]: self.path=np.append(self.path, path[1:], axis=1) else: self.path=np.append(self.path, path, axis=1) def Translate(self, xyz): ''' move the wire in space :param xyz: 3 component vector that describes translation in x,y and z direction ''' if self.path is not None: self.path += np.array(xyz) return self def Rotate(self, axis=(1,0,0), deg=0): ''' rotate wire around given axis by deg degrees :param axis: axis of rotation :param deg: angle ''' if self.path is not None: n = axis ca = np.cos(np.radians(deg)) sa = np.sin(np.radians(deg)) R = np.array([[n[0]**2*(1-ca)+ca, n[0]*n[1]*(1-ca)-n[2]*sa, n[0]*n[2]*(1-ca)+n[1]*sa], [n[1]*n[0]*(1-ca)+n[2]*sa, n[1]**2*(1-ca)+ca, n[1]*n[2]*(1-ca)-n[0]*sa], [n[2]*n[0]*(1-ca)-n[1]*sa, n[2]*n[1]*(1-ca)+n[0]*sa, n[2]**2*(1-ca)+ca]]) self.path = np.dot(self.path, R.T) return self # different standard paths @staticmethod def LinearPath(pt1=(0, 0, 0), pt2=(0, 0, 1)): return np.array([pt1, pt2]).T @staticmethod def RectangularPath(dx=0.1, dy=0.2): dx2 = dx/2.0; dy2 = dy/2.0 return np.array([[dx2, dy2, 0], [dx2, -dy2, 0], [-dx2, -dy2, 0], [-dx2, dy2, 0], [dx2, dy2, 0]]).T @staticmethod def CircularPath(radius=0.1, pts=20): return Wire.EllipticalPath(rx=radius, ry=radius, pts=pts) @staticmethod def SinusoidalCircularPath(radius=0.1, amplitude=0.01, frequency=10, pts=100): t = np.linspace(0, 2 * np.pi, pts) return np.array([radius * np.sin(t), radius * np.cos(t), amplitude * np.cos(frequency*t)]).T @staticmethod def EllipticalPath(rx=0.1, ry=0.2, pts=20): t = np.linspace(0, 2 * np.pi, pts) return np.array([rx * np.sin(t), ry * np.cos(t), 0]).T @staticmethod def SolenoidPath(radius=0.1, pitch=0.01, turns=30, pts_per_turn=20): return Wire.EllipticalSolenoidPath(rx=radius, ry=radius, pitch=pitch, turns=turns, pts_per_turn=pts_per_turn) @staticmethod def EllipticalSolenoidPath(rx=0.1, ry=0.2, pitch=0.01, turns=30, pts_per_turn=20): t = np.linspace(0, 2 * np.pi * turns, pts_per_turn * turns) return np.array([rx * np.sin(t), ry * np.cos(t), t / (2 * np.pi) * pitch]).T 

The second class BiotSavart performs the actual calculations. An arbitrary number Wire class instances can be added to define the current distribution in space. The member function CalculateB accepts a list of 3D coordinates and returns a list of 3D vectors corresponding to the B field at those points

  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 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 __author__ = 'wack' # part of the magwire package # calculate magnetic fields arising from electrical current through wires of arbitrary shape # with the law of Biot-Savart # written by Michael Wack # wack@geophysik.uni-muenchen.de # tested with python 3.4.3 import numpy as np import time #import multiprocessing as mp from matplotlib import pyplot as plt import mpl_toolkits.mplot3d.axes3d as ax3d class BiotSavart: ''' calculates the magnetic field generated by currents flowing through wires ''' def __init__(self, wire=None): self.wires = [] if wire is not None: self.wires.append(wire) def AddWire(self, wire): self.wires.append(wire) def CalculateB(self, points): """ calculate magnetic field B at given points :param points: numpy array of n points (xyz) :return: numpy array of n vectors representing the B field at given points """ print("found {} wire(s).".format(len(self.wires))) c = 0 # generate list of IdL and r1 vectors from all wires for w in self.wires: c += 1 _IdL, _r1 = w.IdL_r1 print("wire {} has {} segments".format(c, len(_IdL))) if c == 1: IdL = _IdL r1 = _r1 else: IdL = np.vstack((IdL, _IdL)) r1 = np.vstack((r1, _r1)) print("total number of segments: {}".format(len(IdL))) print("number of field points: {}".format(len(points))) print("total number of calculations: {}".format(len(points)*len(IdL))) # now we have # all segment vectors multiplied by the flowing current in IdL # and all vectors to the central points of the segments in r1 # calculate vector B*1e7 for each point in space t1 = time.process_time() # simple list comprehension to calculate B at each point r B = np.array([BiotSavart._CalculateB1(r, IdL, r1) * 1e-7 for r in points]) # multi processing # slower than single processing? #pool = mp.Pool(processes=16) #B = np.array([pool.apply(_CalculateB1, args=(r, IdL, r1)) for r in points]) t2 = time.process_time() print("time needed for calculation: {} s".format(t2-t1)) return B def vv_PlotWires(self): for w in self.wires: w.vv_plot_path() def mpl3d_PlotWires(self, ax): for w in self.wires: w.mpl3d_plot_path(show=False, ax=ax) @staticmethod def _CalculateB1(r, IdL, r1): ''' calculate magnetic field B for one point r in space :param r: 3 component numpy array representing the location where B will be calculated :param IdL: all segment vectors multiplied by the flowing current :param r1: all vectors to the central points of the segments :return: numpy array of 3 component vector of B multiplied by 1e7 ''' # calculate law of biot savart for all current elements at given point r r2 = r - r1 r25 = np.linalg.norm(r2, axis=1)**3 r3 = r2 / r25[:, np.newaxis] cr = np.cross(IdL, r3) # claculate sum of contributions from all current elements s = np.sum(cr, axis=0) return s 

solenoid_demo.py demonstrates the calculation and plotting of B fields generated by a simple solenoid.

  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 __author__ = 'wack' # part of the magwire package # calculate magnetic fields arising from electrical current through wires of arbitrary shape # with the law of Biot-Savart # written by Michael Wack 2015 # wack@geophysik.uni-muenchen.de # tested with python 3.4.3 # some basic calculations for testing import numpy as np import matplotlib.pyplot as plt import wire import biotsavart # simple solenoid # approximated analytical solution: B = mu0 * I * n / l = 4*pi*1e-7[T*m/A] * 100[A] * 10 / 0.5[m] = 2.5mT w = wire.Wire(path=wire.Wire.SolenoidPath(pitch=0.05, turns=10), discretization_length=0.01, current=100).Rotate(axis=(1, 0, 0), deg=90) #.Translate((0.1, 0.1, 0)). sol = biotsavart.BiotSavart(wire=w) resolution = 0.04 volume_corner1 = (-.2, -.8, -.2) volume_corner2 = (.2+1e-10, .3, .2) # matplotlib plot 2D # create list of xy coordinates grid = np.mgrid[volume_corner1[0]:volume_corner2[0]:resolution, volume_corner1[1]:volume_corner2[1]:resolution] # create list of grid points points = np.vstack(map(np.ravel, grid)).T points = np.hstack([points, np.zeros([len(points),1])]) # calculate B field at given points B = sol.CalculateB(points=points) Babs = np.linalg.norm(B, axis=1) # remove big values close to the wire cutoff = 0.005 B[Babs > cutoff] = [np.nan,np.nan,np.nan] #Babs[Babs > cutoff] = np.nan for ba in B: print(ba) #2d quiver # get 2D values from one plane with Z = 0 fig = plt.figure() ax = fig.gca() ax.quiver(points[:, 0], points[:, 1], B[:, 0], B[:, 1], scale=.15) X = np.unique(points[:, 0]) Y = np.unique(points[:, 1]) cs = ax.contour(X, Y, Babs.reshape([len(X), len(Y)]).T, 10) ax.clabel(cs) plt.xlabel('x') plt.ylabel('y') plt.axis('equal') plt.show() # matplotlib plot 3D grid = np.mgrid[volume_corner1[0]:volume_corner2[0]:resolution*2, volume_corner1[1]:volume_corner2[1]:resolution*2, volume_corner1[2]:volume_corner2[2]:resolution*2] # create list of grid points points = np.vstack(map(np.ravel, grid)).T # calculate B field at given points B = sol.CalculateB(points=points) Babs = np.linalg.norm(B, axis=1) fig = plt.figure() # 3d quiver ax = fig.gca(projection='3d') sol.mpl3d_PlotWires(ax) ax.quiver(points[:, 0], points[:, 1], points[:, 2], B[:, 0], B[:, 1], B[:, 2], length=0.04) plt.show()