This is a simple approximation to the Biot-Savart law.
I've implemented the integral(sum) in the function calc(),
If the number of spatial points is big, say 10^7 or 10^8 -ish, can calc be written to use NumPy arrays more efficiently? Thanks for your suggestions!
def calc(points, x_seg, idl_seg):
r = points[:, None, :] - x_seg[None, :, :] # START CALCULATION
bottom = ((r**2).sum(axis=-1)**1.5)[...,None] # 1/|r|**3 add axis for vector
top = np.cross(idl_seg[None,:,:], r) # np.cross defaults to last axis
db = (mu0 / four_pi) * top / bottom
b = db.sum(axis=-2) # sum over the segments of the current loop
return b
EDIT: So for example, I can do this. Now there are just two arrays (r and hold) of size nx * ny * nz * nseg * 3. Maybe I should pass smaller chunks of points at a time, so it can all fit in cache at once?
def calc_alt(points, x_seg, idl_seg):
r = points[:, None, :] - x_seg[None, :, :]
hold = np.ones_like(r)*((r**2).sum(axis=-1)**-1.5)[...,None] # note **-1.5 neg
b = (hold * np.cross(idl_seg[None,:,:], r)).sum(axis=-2)
return b * (mu0 / four_pi)
The rest of the code is posted to show how calc is used.
import numpy as np
import matplotlib.pyplot as plt
pi, four_pi = np.pi, 4. * np.pi
mu0 = four_pi * 1E-07 # Tesla m/A exact, defined
r0 = 0.05 # meters
I0 = 100.0 # amps
nx, ny, nz = 48, 49, 50
x,y,z = np.linspace(0,2*r0,nx), np.linspace(0,2*r0,ny), np.linspace(0,2*r0,nz)
xg = np.zeros((nx, ny, nz, 3)) # 3D grid of position vectors
xg[...,0] = x[:, None, None] # fill up the positions
xg[...,1] = y[None, :, None]
xg[...,2] = z[None, None, :]
xgv = xg.reshape(nx*ny*nz, 3) # flattened view of spatial points
nseg = 32 # approximate the current loop as a set of discrete points I*dl
theta = np.linspace(0, 2.*pi, nseg+1)[:-1] # get rid of the repeat
xdl = np.zeros((nseg, 3)) # these are the position vectors
idl = np.zeros((nseg, 3)) # these are the current vectors
xdl[:,0], xdl[:,1] = r0 * np.cos(theta), r0 * np.sin(theta)
idl[:,0], idl[:,1] = I0 * -np.sin(theta), I0 * np.cos(theta)
b = calc(xgv, xdl, idl) # HERE IS THE CALCULATION
bv = b.reshape(nx, ny, nz, 3) # make a "3D view" again to use for plotting
bx, by, bz = bv[...,0], bv[...,1], bv[...,2] # make component views
bperp = np.sqrt(bx**2 + by**2) # new array for perp field
zround = np.round(z, 4)
iz = 5 # choose a transverse plane for a plot
fields = [ bz, bperp, bx, by]
names = ['Bz', 'Bperp', 'Bx', 'By']
titles = ["approx " + name + " at z = " + str(zround[iz])
for name in names]
plt.figure()
for i, field in enumerate(fields):
print i
plt.subplot(2, 2, i+1)
plt.imshow(field[..., iz], origin='lower') # fields at iz don't use Jet !!!
plt.title(titles[i])
plt.colorbar()
plt.show()
The plotting at the end is just to see that it appears to be working. In reality, never use the default colormap. Bad, awful, naughty Jet! In this case, a divergent cmap with symmetric vmin = -vmax might be good. (see Jake VanderPlas' post, and the matplotlib documentation, and there's some lovely demos down here.
I have just run across np.numexpr which does (among other things) what I suggested in the edit - breaks the arrays into "chunks" so that they can fit into cache, including all temporary arrays needed to evaluate expressions.
There are nice explanations here and especially in this wiki.