diff --git a/as5013-test/sim/CMD b/as5013-test/sim/CMD
new file mode 100644
index 0000000000000000000000000000000000000000..7390fc96192caceb79f522e9b05aa4bfed6f686f
--- /dev/null
+++ b/as5013-test/sim/CMD
@@ -0,0 +1,2 @@
+python bfield_fd.py --mode single_magnet -nt 40000 -mu 100 -nx 200 -ny 200
+python bfield_fd.py --mode double_magnet -nt 40000 -mu 100 -nx 200 -ny 200
diff --git a/as5013-test/sim/bfield_fd.py b/as5013-test/sim/bfield_fd.py
new file mode 100644
index 0000000000000000000000000000000000000000..9086061fe1d52d53c521e5c9e0881c48484bc47f
--- /dev/null
+++ b/as5013-test/sim/bfield_fd.py
@@ -0,0 +1,199 @@
+#! /usr/bin/env python
+from __future__ import division
+from numpy import *
+import argparse
+from matplotlib import pyplot as plt
+from matplotlib import animation
+plt.style.use('bmh')
+from matplotlib.patches import Rectangle
+import matplotlib.colors as colors
+
+#using sor to solve for electric fields
+eps_0 = 8.85e-12 #permittivity of free space
+mu_0 = 1.25663706e-6 #permittivity of free space (m*kg/s^2/A^2)
+
+fn_ind = 0
+contours = None
+
+def compute_a(mu):
+ #compute square contour integrals for permeability mu
+ a0 = mu[1:,1:] + mu[:-1,1:] + mu[1:,:-1] + mu[:-1,:-1]
+ a1 = .5*(mu[1:,1:] + mu[1:,:-1])
+ a2 = .5*(mu[1:,1:] + mu[:-1,1:])
+ a3 = .5*(mu[:-1,:-1] + mu[:-1,1:])
+ a4 = .5*(mu[:-1,:-1] + mu[1:,:-1])
+ return a0,a1,a2,a3,a4
+def solve_poisson(args,V,rho,a,neuman_bcs=[None,None,None,None]):
+ for ti in range(0,args.nt):
+ V[0] = V[1]
+ if neuman_bcs[0] is not None:
+ V[1,0,:] = V[1,1,:] - neuman_bcs[0]
+ if neuman_bcs[1] is not None:
+ V[1,-1,:] = V[1,-2,:] + neuman_bcs[1]
+ if neuman_bcs[2] is not None:
+ V[1,:,0] = V[1,:,1] - neuman_bcs[2]
+ if neuman_bcs[3] is not None:
+ V[1,:,-1] = V[1,:,-2] + neuman_bcs[3]
+ V[1,1:-1,1:-1] = (1-args.alpha)*V[0,1:-1,1:-1] + \
+ args.alpha*( a[1]*V[0,2:,1:-1] + a[2]*V[0,1:-1,2:] + a[3]*V[0,:-2,1:-1] + a[4]*V[0,1:-1,:-2] )/a[0] + \
+ args.alpha*rho[1:-1,1:-1]/mu_0/a[0]
+
+ R = absolute(V[0,1:-1,1:-1] - V[1,1:-1,1:-1])
+ if amax(R) < 1e-6:
+ print "broke after %d iterations"%ti
+ break
+ return V
+def compute_H(V,dx,dy):
+ #take derivatives of V and shift indices
+ Hx = -.5*(V[-1,1:,1:] - V[-1,1:,:-1])/dx - .5*(V[-1,:-1,1:] - V[-1,:-1,:-1])/dx
+ Hy = -.5*(V[-1,1:,1:] - V[-1,:-1,1:])/dy - .5*(V[-1,1:,:-1] - V[-1,:-1,:-1])/dy
+ Hmag = sqrt(Hx*Hx + Hy*Hy)
+ return Hx,Hy,Hmag
+
+
+def run_single_magnet(args):
+ nx = args.nx; dx = args.xs/nx
+ ny = args.ny; dy = args.ys/ny
+ py1 = int(args.ny/2 - int(args.l/2/dy)) #index for plate location
+ py2 = int(args.ny/2 + int(args.l/2/dy)) #index for plate location
+ py_sep1 = int(args.ny/2 - int(args.l/2/dy) - int(args.sep/dy))
+ px1 = int(args.nx/2 - int(args.w/2/dx)) #index for plate edge
+ px2 = int(args.nx/2 + int(args.w/2/dx))
+ V = zeros((2,ny,nx)) #use twice time steps for half-stepping
+ V[:,:,0] = zeros(ny)
+ V[:,0,:] = zeros(nx)
+ V[:,:,-1] = zeros(ny)
+ V[:,-1,:] = zeros(nx)
+
+ mu = ones((ny-1,nx-1))
+ mu[ py1:py2, px1:px2+1 ] = args.mu
+ #mu[ 0:py1, px1:px2+1 ] = args.mu
+
+ divM = zeros_like(V[0]) #divergence of magnetic dipole moment density M = B_r / mu_0
+ vol = args.w*args.w*args.l
+ divM[py1,px1:px2+1] = -args.Br/mu_0/(px2-px1)*vol
+ divM[py2,px1:px2+1] = args.Br/mu_0/(px2-px1)*vol
+
+ V = solve_poisson(args,V,divM,compute_a(mu),neuman_bcs=[0,0,0,0])
+ Hx,Hy,Hmag = compute_H(V,dx,dy)
+
+
+ fig, axs = plt.subplots(2, 1, gridspec_kw = {'height_ratios':[2.5, 1]})
+ axs[0].set_aspect(1.)
+ x = linspace(-.5*args.xs, .5*args.xs, nx); y = linspace(-.5*args.ys, .5*args.ys, ny);
+ CS = axs[0].contourf(x[1:], y[1:], Hmag, cmag = plt.cm.autumn)
+ #cbar = fig.colorbar(CS, ax=axs[0])
+ #cbar.ax.set_ylabel('field strength A/m')
+ #plt.contourf(x, y, V[-1])
+ X,Y = meshgrid(x[1:],y[1:])
+ axs[0].set_xlim([x[0],x[-1]])
+ axs[0].set_ylim([args.pf*y[0],args.pf*y[-1]])
+ strm = axs[0].streamplot(X,Y, Hx,Hy,color=Hmag, minlength=args.w/2,linewidth=1, cmap=plt.cm.summer)
+ axs[0].add_patch(Rectangle( (-args.w/2, -args.l/2), args.w, args.l, facecolor="grey", alpha=.4 ))
+ axs[0].plot([-args.xs/2,args.xs/2],[-args.l/2-args.sep,-args.l/2-args.sep],c='k',lw=2,ls='--')
+
+ axs[1].plot(x[:-1],Hy[py_sep1])
+ axs[1].set_xlabel('x (m)')
+ axs[1].set_ylabel('B_z')
+ axs[1].set_yticks([0.])
+ axs[1].set_xlim([-2*args.w,2*args.w]) #set x range based on magnet width
+
+ plt.savefig('run_single_magnet.png', bbox_inches='tight', pad_inches=0)
+ plt.show()
+
+
+def run_double_magnet(args):
+ nx = args.nx; dx = args.xs/nx
+ ny = args.ny; dy = args.ys/ny
+ py1 = int(args.ny/2 - int(args.l/2/dy)) #index for plate location
+ py2 = int(args.ny/2 + int(args.l/2/dy)) #index for plate location
+ py_sep1 = int(args.ny/2 - int(args.l/2/dy) - int(args.sep/dy))
+ px1 = int(args.nx/2 - int(args.w/dx)) #index for plate edge
+ pxm = int(args.nx/2)
+ px2 = int(args.nx/2 + int(args.w/dx))
+ V = zeros((2,ny,nx)) #use twice time steps for half-stepping
+ V[:,:,0] = zeros(ny)
+ V[:,0,:] = zeros(nx)
+ V[:,:,-1] = zeros(ny)
+ V[:,-1,:] = zeros(nx)
+
+ mu = ones((ny-1,nx-1))
+ mu[ py1:py2, px1:px2+1 ] = args.mu
+ #mu[ 0:py1, px1:px2+1 ] = args.mu
+
+ divM = zeros_like(V[0]) #divergence of magnetic dipole moment density M = B_r / mu_0
+ vol = args.w*args.w*args.l
+ divM[py1,px1:pxm] = -args.Br/mu_0/(px2-px1)*vol
+ divM[py1,pxm:px2] = args.Br/mu_0/(px2-px1)*vol
+ divM[py2,px1:pxm] = args.Br/mu_0/(px2-px1)*vol
+ divM[py2,pxm:px2] = -args.Br/mu_0/(px2-px1)*vol
+
+ V = solve_poisson(args,V,divM,compute_a(mu),neuman_bcs=[0,0,0,0])
+ Hx,Hy,Hmag = compute_H(V,dx,dy)
+
+
+ fig, axs = plt.subplots(2, 1, gridspec_kw = {'height_ratios':[2.5, 1]})
+ axs[0].set_aspect(1.)
+ x = linspace(-.5*args.xs, .5*args.xs, nx); y = linspace(-.5*args.ys, .5*args.ys, ny);
+ CS = axs[0].contourf(x[1:], y[1:], Hmag, cmag = plt.cm.autumn)
+ #cbar = fig.colorbar(CS, ax=axs[0])
+ #cbar.ax.set_ylabel('field strength A/m')
+ #plt.contourf(x, y, V[-1])
+ X,Y = meshgrid(x[1:],y[1:])
+ axs[0].set_xlim([x[0],x[-1]])
+ axs[0].set_ylim([args.pf*y[0],args.pf*y[-1]])
+ strm = axs[0].streamplot(X,Y, Hx,Hy,color=Hmag, minlength=args.w/2,linewidth=1, cmap=plt.cm.summer)
+ axs[0].add_patch(Rectangle( (-args.w, -args.l/2), .99*args.w, args.l, facecolor="grey", alpha=.4 ))
+ axs[0].add_patch(Rectangle( (0.01*args.w, -args.l/2), args.w, args.l, facecolor="grey", alpha=.4 ))
+ axs[0].plot([-args.xs/2,args.xs/2],[-args.l/2-args.sep,-args.l/2-args.sep],c='k',lw=2,ls='--')
+
+ axs[1].plot(x[:-1],Hy[py_sep1])
+ #axs[1].plot(x[:-1],Hy[py_sep2])
+ #axs[1].plot(x[:-1],Hy[py_sep3])
+ axs[1].set_xlabel('x (m)')
+ axs[1].set_ylabel('B_z')
+ axs[1].set_yticks([0.])
+ axs[1].set_xlim([-2*args.w,2*args.w]) #set x range based on magnet width
+
+ plt.savefig('run_double_magnet.png', bbox_inches='tight', pad_inches=0)
+ plt.show()
+
+
+if __name__ == '__main__':
+ parser = argparse.ArgumentParser()
+ parser.add_argument('-M','--mode',choices=('single_magnet','double_magnet'))
+ parser.add_argument("-nx","--nx", type=int, default=200, help="grid size")
+ parser.add_argument("-ny","--ny", type=int, default=200, help="grid size")
+ parser.add_argument("-xs","--xs", type=float, default=.05, help="space extent x")
+ parser.add_argument("-ys","--ys", type=float, default=.04, help="space extent y")
+ parser.add_argument("-pf","--pf", type=float, default=.75, help="plot fraction (how much of simulated space to plot)")
+ parser.add_argument("-l","--l", type=float, default=.01, help="magnet length")
+ parser.add_argument("-w","--w", type=float, default=.002, help="magnet width")
+ parser.add_argument("-sep","--sep", type=float, default=.001, help="separation between magnet bottom and graphing plane")
+ parser.add_argument("-mu","--mu", type=float, default=1., help="relative permeability of magnet")
+ parser.add_argument("-Br","--Br", type=float, default=1., help="residual flux density of magnet (Teslas)")
+ parser.add_argument("-nt","--nt", type=int, default=50, help="num time steps")
+ parser.add_argument("-alpha","--alpha", type=float, default=1., help="alpha, relaxation parameter")
+ args = parser.parse_args()
+ if args.mode == 'single_magnet':
+ run_single_magnet(args)
+ if args.mode == 'double_magnet':
+ run_double_magnet(args)
+ '''elif args.mode == 'moving_dialectric_sweep':
+ Cs = []
+ seps = [(args.ys/args.ny)*i for i in range(3,15)]
+ for sep in seps:
+ args.sep = sep
+ Cs.append(run_moving_dialectric(args))
+ fn_ind += 1
+ print "capacitances: ",list(Cs)
+ print "seps",list(seps)
+ plt.figure()
+ plt.plot(1e3*asarray(seps),1e12*asarray(Cs))
+ plt.title('Adjacent trace Capacitance vs. separation')
+ plt.xlabel('separation distance (mm)')
+ plt.ylabel('capacitance (pF)')
+ plt.ylim([0,1.1*amax(1e12*asarray(Cs))])
+ plt.show()
+ '''
+
diff --git a/as5013-test/sim/run_double_magnet.png b/as5013-test/sim/run_double_magnet.png
new file mode 100644
index 0000000000000000000000000000000000000000..afc422fd0d22af85a64aba1fc2f19cd68eb98a85
Binary files /dev/null and b/as5013-test/sim/run_double_magnet.png differ
diff --git a/as5013-test/sim/run_single_magnet.png b/as5013-test/sim/run_single_magnet.png
new file mode 100644
index 0000000000000000000000000000000000000000..f6d3e85b73ebd2e1b02624f1b6382c665e417007
Binary files /dev/null and b/as5013-test/sim/run_single_magnet.png differ