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 205 | __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 108 109 110 | __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 89 90 91 | __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()
|