#! /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
fn_ind = 0
contours = None

def compute_a(e):
	#compute square contour integrals for permittivity e
	a0 = e[1:,1:] + e[:-1,1:] + e[1:,:-1] + e[:-1,:-1]
	a1 = .5*(e[1:,1:] + e[1:,:-1])
	a2 = .5*(e[1:,1:] + e[:-1,1:])
	a3 = .5*(e[:-1,:-1] + e[:-1,1:])
	a4 = .5*(e[:-1,:-1] + e[1:,:-1])
	return a0,a1,a2,a3,a4
def solve_poisson(args,V,rho,a):
	#TODO: implement a sparse matrix formulation of the solver -- should work better for higher permittivities
	for ti in range(0,args.nt):
		V[0] = V[1]
		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]/eps_0/a[0]
		R = absolute(V[0,1:-1,1:-1] - V[1,1:-1,1:-1])
		#V[ti+1,1:-1,1:-1] = new
		if amax(R) < 1e-6:
			print "broke after %d iterations"%ti
			break

	return V
def compute_E(V,dx,dy):
	#take derivatives of V and shift indices
	Ex = -.5*(V[-1,1:,1:] - V[-1,1:,:-1])/dx - .5*(V[-1,:-1,1:] - V[-1,:-1,:-1])/dx
	Ey = -.5*(V[-1,1:,1:] - V[-1,:-1,1:])/dy - .5*(V[-1,1:,:-1] - V[-1,:-1,:-1])/dy
	Emag = sqrt(Ex*Ex + Ey*Ey)
	return Ex,Ey,Emag

def run_moving_dialectric(args):
	nx = args.nx; dx = args.xs/nx
	ny = args.ny; dy = args.ys/ny
	py = int(args.ny/3) #index for plate location
	py2 = int( py + int(args.sep/dy) )
	py3 = int( py + int(args.sep/dy) + int(3*args.dialectric_thickness/dy) )
	py4 = int(py - int(args.dialectric_thickness/dy))

	px1 = int(args.nx/2 - int(args.w/dx) - int(args.d/2/dx)) #index for plate edge
	px2 = int(args.nx/2 - int(args.d/2/dx)) #index for plate edge
	px3 = int(args.nx/2 + int(args.d/2/dx)) #index for plate edge
	px4 = int(args.nx/2 + int(args.w/dx) + int(args.d/2/dx)) #index for plate edge
	
	px5 = int(args.nx/2 + int(args.w/dx) + int(3*args.d/2/dx))
	px6 = int(args.nx/2 + 2*int(args.w/dx) + int(3*args.d/2/dx))
	px7 = int(args.nx/2 - 2*int(args.w/dx) - int(3*args.d/2/dx)) #index for plate edge
	px8 = int(args.nx/2 - int(args.w/dx) - int(3*args.d/2/dx)) #index for plate edge

	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)

	e = ones((ny-1,nx-1))
	e[ py2:py3+1, : ] = 1e6*args.eps
	e[ py4:py+1, : ] = args.eps #base layer

	rho = zeros_like(V[0]) #total charge at each grid point (enclosed by square contour)
	rho[py,px1:px2+1] = args.Q/(px2-px1) #charge per point in plate
	rho[py,px3:px4+1] = -args.Q/(px4-px3)
	rho[py,px5:px6+1] = args.Q/(px4-px3)
	rho[py,px7:px8+1] = -args.Q/(px4-px3)

	V = solve_poisson(args,V,rho,compute_a(e))
	Ex,Ey,Emag = compute_E(V,dx,dy)

	#fig = plt.figure()
	fig = plt.figure(figsize=(20,10),dpi=600)
	axs = [plt.gca()]
	axs[0].set_aspect(1.)
	x = linspace(-.5*args.xs, .5*args.xs, nx); y = linspace(-.5*args.ys, .5*args.ys, ny);
	#if contours is None:
	#CS = axs[0].contourf(x[1:], y[1:], Emag,  20, norm=colors.LogNorm(vmin=Emag.min(), vmax=Emag.max()/10.), cmap = plt.cm.jet)
	#print Emag.min(), Emag.max()

	#CS = axs[0].contourf(x[1:], y[1:], Emag, levels=logspace(-3.5,5,num=15), norm=colors.LogNorm(vmin=Emag.min(), vmax=Emag.max()), cmap = plt.cm.jet)
	CS = axs[0].contourf(x[1:], y[1:], Emag, levels=logspace(-3.5,5,num=15), 
		norm=colors.LogNorm(vmin=0.0273988722319, vmax=4204.10388092), cmap = plt.cm.jet)
	#else:
	#	CS = axs[0].contourf(x[1:], y[1:], Emag, contours, norm=colors.LogNorm(vmin=Emag.min(), vmax=Emag.max()), cmap = plt.cm.jet)
	print "levels=", list(CS.levels)
	#cbar = fig.colorbar(CS, ax=axs[0])
	#cbar.ax.set_ylabel('field strength V/m')	
	
	X,Y = meshgrid(x[1:],y[1:])
	
	strm = axs[0].streamplot(X,Y, Ex,Ey,color=Emag/amax(Emag), minlength=args.d/2,linewidth=1, cmap=plt.cm.summer)
	
	axs[0].plot([x[px1],x[px2]],[y[py],y[py]],c='k',lw=2)
	axs[0].plot([x[px3],x[px4]],[y[py],y[py]],c='k',lw=2)
	axs[0].plot([x[px5],x[px6]],[y[py],y[py]],c='k',lw=2)
	axs[0].plot([x[px7],x[px8]],[y[py],y[py]],c='k',lw=2)

	def add_rect(x0,x1,y0,y1):
		axs[0].add_patch(Rectangle( (x0, y0), x1-x0, y1-y0, facecolor="grey", alpha=.4 ))

	add_rect(x[0],x[-1],y[py2], y[py3] )
	add_rect(x[0],x[-1],y[py4], y[py] )

	axs[0].set_xlim([x[0],x[-1]])
	axs[0].set_ylim([y[0],y[-1]])
	
	plt.grid(False)

	q = sum(rho[py,px1:px2+1]) + sum(rho[py,px5:px6+1])
	V_average = average(V[-1,py,px1:px2])
	#print "Total charge: %.2e Coloumbs"%q
	#print "Average Voltage: %.2f Volts"%V_average
	C = q/V_average
	#print "Approximate Capacitance: %.2e F"%C
	axs[0].set_title(r'Moving dialectric ($\epsilon_r=%.1f$), interleaved traces, $C=%.2e F$'%(args.eps,C),fontsize=20)
	plt.savefig('capacitance-%03d'%fn_ind, bbox_inches='tight', pad_inches=0)
	plt.close()
	#plt.show()
	return C


def run_moving_conductor(args):
	nx = args.nx; dx = args.xs/nx
	ny = args.ny; dy = args.ys/ny
	py = int(args.ny/3) #index for plate location
	py2 = int( py + int(args.sep/dy) )
	py3 = int( py + int(args.sep/dy) + int(2*args.dialectric_thickness/dy) )
	py4 = int(py - int(args.dialectric_thickness/dy))

	print py,py2,py3,py4

	px1 = int(args.nx/2 - int(args.w/dx) - int(args.d/2/dx)) #index for plate edge
	px2 = int(args.nx/2 - int(args.d/2/dx)) #index for plate edge
	px3 = int(args.nx/2 + int(args.d/2/dx)) #index for plate edge
	px4 = int(args.nx/2 + int(args.w/dx) + int(args.d/2/dx)) #index for plate edge
	
	#px5 = int(args.nx/2 + int(args.w/dx) + int(3*args.d/2/dx))
	#px6 = int(args.nx/2 + 2*int(args.w/dx) + int(3*args.d/2/dx))
	#px7 = int(args.nx/2 - 2*int(args.w/dx) - int(3*args.d/2/dx)) #index for plate edge
	#px8 = int(args.nx/2 - int(args.w/dx) - int(3*args.d/2/dx)) #index for plate edge

	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)

	e = ones((ny-1,nx-1))
	#e[ py2:py3+1, : ] = 1e6*args.eps
	e[ py4:py+1, : ] = args.eps #base layer

	rho = zeros_like(V[0]) #total charge at each grid point (enclosed by square contour)
	rho[py,px1:px2+1] = args.Q/(px2-px1) #charge per point in plate
	rho[py,px3:px4+1] = -args.Q/(px4-px3)
	rho[py2,px1:px2+1] = -args.Q/(px4-px3)
	rho[py2,px3:px4+1] = args.Q/(px4-px3)

	#V = solve_poisson(args,V,rho,compute_a(e))
	a = compute_a(e)
	for ti in range(0,args.nt):
		V[0] = V[1]
		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]/eps_0/a[0]
		#impose von neumann condition on conductor
		V[1,py2:py3+1,:] = average(V[1,py2:py3+1,:])
		R = absolute(V[0,1:-1,1:-1] - V[1,1:-1,1:-1])
		#V[ti+1,1:-1,1:-1] = new
		if amax(R) < 1e-6:
			print "broke after %d iterations"%ti
			break
	print "Residual max size: ",amax(R)
	Ex,Ey,Emag = compute_E(V,dx,dy)

	#fig = plt.figure()
	fig = plt.figure(figsize=(20,10),dpi=600)
	axs = [plt.gca()]
	axs[0].set_aspect(1.)
	x = linspace(-.5*args.xs, .5*args.xs, nx); y = linspace(-.5*args.ys, .5*args.ys, ny);
	#if contours is None:
	#CS = axs[0].contourf(x[1:], y[1:], Emag,  20, norm=colors.LogNorm(vmin=Emag.min(), vmax=Emag.max()/10.), cmap = plt.cm.jet)
	print Emag.min(), Emag.max()
	Emag += 1e-6

	#CS = axs[0].contourf(x[1:], y[1:], Emag, levels=logspace(-3.5,5,num=15), norm=colors.LogNorm(vmin=Emag.min(), vmax=Emag.max()), cmap = plt.cm.jet)
	CS = axs[0].contourf(x[1:], y[1:], Emag, levels=logspace(-3.5,5,num=15), 
		norm=colors.LogNorm(vmin=1e-6, vmax=18372), cmap = plt.cm.jet)
	#else:
	#	CS = axs[0].contourf(x[1:], y[1:], Emag, contours, norm=colors.LogNorm(vmin=Emag.min(), vmax=Emag.max()), cmap = plt.cm.jet)
	print "levels=", list(CS.levels)
	#cbar = fig.colorbar(CS, ax=axs[0])
	#cbar.ax.set_ylabel('field strength V/m')	
	
	X,Y = meshgrid(x[1:],y[1:])
	
	strm = axs[0].streamplot(X,Y, Ex,Ey,color=Emag/amax(Emag), minlength=args.d/2,linewidth=1, cmap=plt.cm.summer)
	
	axs[0].plot([x[px1],x[px2]],[y[py],y[py]],c='k',lw=2)
	axs[0].plot([x[px3],x[px4]],[y[py],y[py]],c='k',lw=2)
	axs[0].plot([x[px1],x[px2]],[y[py2],y[py2]],c='k',lw=2)
	axs[0].plot([x[px3],x[px4]],[y[py2],y[py2]],c='k',lw=2)

	def add_rect(x0,x1,y0,y1):
		axs[0].add_patch(Rectangle( (x0, y0), x1-x0, y1-y0, facecolor="grey", alpha=.4 ))

	add_rect(x[0],x[-1],y[py2], y[py3] )
	add_rect(x[0],x[-1],y[py4], y[py] )

	axs[0].set_xlim([x[0],x[-1]])
	axs[0].set_ylim([y[0],y[-1]])
	
	plt.grid(False)

	q = sum(rho[py,px1:px2+1]) #+ sum(rho[py,px5:px6+1])
	V_average = average(V[-1,py,px1:px2])
	#print "Total charge: %.2e Coloumbs"%q
	#print "Average Voltage: %.2f Volts"%V_average
	C = q/V_average
	#print "Approximate Capacitance: %.2e F"%C
	axs[0].set_title(r'Moving conductor, adjacent traces, $C=%.2e F$'%(C),fontsize=20)
	plt.savefig('capacitance-%03d'%fn_ind, bbox_inches='tight', pad_inches=0)
	plt.close()
	#plt.show()
	return C





def run_moving_plate(args):
	nx = args.nx; dx = args.xs/nx
	ny = args.ny; dy = args.ys/ny
	py = int(args.ny/3) #index for plate location
	py2 = int( py + int(args.sep/dy) )
	py3 = int( py + int(args.sep/dy) + int(args.dialectric_thickness/dy) )
	py4 = int(py - int(args.dialectric_thickness/dy))

	px1 = int(args.nx/2 - int(args.w/dx))  #index for left plate edge
	px2 = int(args.nx/2 + int(args.w/dx))  #index for right plate edge

	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)

	e = ones((ny-1,nx-1))
	e[ py2:py3+1, px1: ] = args.eps
	e[ py4:py+1, px1: ] = args.eps #base layer

	rho = zeros_like(V[0]) #total charge at each grid point (enclosed by square contour)
	rho[py,px1:px2+1] = args.Q/(px2-px1) #charge per point in plate
	rho[py2,px1:px2+1] = -args.Q/(px2-px1)

	V = solve_poisson(args,V,rho,compute_a(e))
	Ex,Ey,Emag = compute_E(V,dx,dy)

	#fig = plt.figure()
	fig = plt.figure(figsize=(20,10),dpi=600)
	axs = [plt.gca()]
	axs[0].set_aspect(1.)
	x = linspace(-.5*args.xs, .5*args.xs, nx); y = linspace(-.5*args.ys, .5*args.ys, ny);
	#if contours is None:
	#CS = axs[0].contourf(x[1:], y[1:], Emag,  20, norm=colors.LogNorm(vmin=Emag.min(), vmax=Emag.max()/10.), cmap = plt.cm.jet)
	#print Emag.min(), Emag.max()
	#CS = axs[0].contourf(x[1:], y[1:], Emag, levels=logspace(-3.5,5,num=15), norm=colors.LogNorm(vmin=Emag.min(), vmax=Emag.max()), cmap = plt.cm.jet)
	CS = axs[0].contourf(x[1:], y[1:], Emag, levels=logspace(-3.5,5,num=15), 
		norm=colors.LogNorm(vmin=2.12416889361e-05, vmax=1456.04532414), cmap = plt.cm.jet)
	#else:
	#	CS = axs[0].contourf(x[1:], y[1:], Emag, contours, norm=colors.LogNorm(vmin=Emag.min(), vmax=Emag.max()), cmap = plt.cm.jet)
	print "levels=", list(CS.levels)
	#cbar = fig.colorbar(CS, ax=axs[0])
	#cbar.ax.set_ylabel('field strength V/m')	
	
	X,Y = meshgrid(x[1:],y[1:])
	
	strm = axs[0].streamplot(X,Y, Ex,Ey,color=Emag/amax(Emag), minlength=args.d/2,linewidth=1, cmap=plt.cm.summer)
	
	axs[0].plot([x[px1],x[px2]],[y[py],y[py]],c='k',lw=2)
	axs[0].plot([x[px1],x[px2]],[y[py2],y[py2]],c='k',lw=2)

	def add_rect(x0,x1,y0,y1):
		axs[0].add_patch(Rectangle( (x0, y0), x1-x0, y1-y0, facecolor="grey", alpha=.4 ))

	add_rect(x[px1],x[-1],y[py2], y[py3] )
	add_rect(x[px1],x[-1],y[py4], y[py] )

	axs[0].set_xlim([x[0],x[-1]])
	axs[0].set_ylim([y[0],y[-1]])
	
	plt.grid(False)

	q = sum(rho[py,px1:px2+1])
	V_average = average(V[-1,py,px1:px2])
	print "Total charge: %.2e Coloumbs"%q
	print "Average Voltage: %.2f Volts"%V_average
	C = q/V_average
	print "Approximate Capacitance: %.2e F"%C
	axs[0].set_title(r'Moving plate ($\epsilon_r=%.1f$), $C=%.2e F$'%(args.eps,C),fontsize=20)
	plt.savefig('capacitance-%03d'%fn_ind, bbox_inches='tight', pad_inches=0)
	plt.close()
	#plt.show()
	return C


def run_parallel_plate(args):
	nx = args.nx; dx = args.xs/nx
	ny = args.ny; dy = args.ys/ny
	py1 = int(args.ny/2 - int(args.d/2/dy)) #index for plate location
	py2 = int(args.ny/2 + int(args.d/2/dy)) #index for plate offset
	px1 = int(args.nx/2 - int(args.w/2/dx)) #index for plate edge
	px2 = int(args.nx/2 + int(args.w/2/dx))
	print py1, py2, px1, px2
	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)
	
	e = ones((ny-1,nx-1))
	#define permittivity distribution
	e[ py1:py2+1, px1:px2+1 ] = args.eps

	rho = zeros_like(V[0]) #total charge at each grid point (enclosed by square contour)
	rho[py1,px1:px2+1] = args.Q/(px2-px1) #charge per point in plate
	rho[py2,px1:px2+1] = -args.Q/(px2-px1)

	V = solve_poisson(args,V,rho,compute_a(e))
	Ex,Ey,Emag = compute_E(V,dx,dy)

	fig, axs = plt.subplots(2, 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:], Emag, cmag = plt.cm.autumn)
	cbar = fig.colorbar(CS, ax=axs[0])
	cbar.ax.set_ylabel('field strength V/m')	
	#plt.contourf(x, y, V[-1])
	X,Y = meshgrid(x[1:],y[1:])

	q = sum(rho[py1,px1:px2+1])
	V_average = average(V[-1,py1,px1:px2])
	print "Total charge: %.2e Coloumbs"%q
	print "Average Voltage: %.2f Volts"%V_average
	print "Approximate Capacitance: %.2e F"%(q/V_average)
	eps = average(e[ py1:py2+1, px1:px2+1 ])
	print "Parallel plate approximation: %.2e F"%(eps*eps_0*args.w/args.d)


	axs[0].set_title(r'Parallel plate capacitor, $\epsilon_r=%.1f$'%(args.eps))
	axs[1].set_title('Parallel plate capacitance: %.2e F, Calculated Capacitance: %.2e F'%((eps*eps_0*args.w/args.d),(q/V_average)),fontsize=12)


	strm = axs[0].streamplot(X,Y, Ex,Ey,color=Emag, minlength=args.d/2,linewidth=1, cmap=plt.cm.summer)
	
	axs[0].plot([-args.w/2,args.w/2],[-args.d/2,-args.d/2],c='k',lw=2)
	axs[0].plot([-args.w/2,args.w/2],[args.d/2,args.d/2],c='k',lw=2)
	axs[0].set_xlim([x[0],x[-1]])
	axs[0].set_ylim([y[0],y[-1]])
	
	axs[0].add_patch(Rectangle( (-args.w/2, -args.d/2), args.w, args.d, facecolor="grey", alpha=.4 ))

	axs[1].plot(x[px1:px2],V[-1,py1,px1:px2])
	axs[1].plot(x[px1:px2],V[-1,py2,px1:px2])
	axs[1].set_xlabel('distance along plate (m)')
	axs[1].set_ylabel('potential (volts)')


	plt.savefig('parallel-plate.png', bbox_inches='tight', pad_inches=0)

	plt.show()


if __name__ == '__main__':
	parser = argparse.ArgumentParser()
	parser.add_argument('-M','--mode',choices=('parallelplate','moving_dialectric','moving_dialectric_sweep', 'moving_plate', 'moving_plate_sweep','moving_conductor','moving_conductor_sweep'))
	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=1., help="space extent x")
	parser.add_argument("-ys","--ys", type=float, default=.5, help="space extent y")
	parser.add_argument("-d","--d", type=float, default=.2, help="plate separation")
	parser.add_argument("-w","--w", type=float, default=.5, help="plate width")
	parser.add_argument("-sep","--sep", type=float, default=.001, help="separation between load cell dialectric and pads")
	parser.add_argument("-dialectric_thickness","--dialectric_thickness", type=float, default=.001, help="thickness of dialectric layer")
	parser.add_argument("-eps","--eps", type=float, default=1., help="relative permittivity of dialectric")
	parser.add_argument("-Q","--Q", type=float, default=1.e-5, help="plate charge")
	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 == 'parallelplate':
		run_parallel_plate(args)
	elif args.mode == 'moving_dialectric':
		#contours = [0.0, 35., 75., 150., 300., 600.0, 1200.0, 2400.0, 4800.0]
		#hstack((array([0]),logspace(1,4,num=7)))#[200*i for i in range(24)]#[0.0, 60000000.0, 120000000.0, 180000000.0, 240000000.0, 300000000.0, 360000000.0, 420000000.0, 480000000.0]
		run_moving_dialectric(args)
	elif args.mode == 'moving_conductor':
		run_moving_conductor(args)
	elif args.mode == 'moving_plate':
		run_moving_plate(args)
	elif args.mode == 'moving_dialectric_sweep':
		Cs = []
		#contours = [0.0, 300000.0, 600000.0, 1200000.0, 1800000.0, 2400000.0, 3000000.0, 3600000.0, 4200000.0, 4800000.0]
		#seps = linspace(.00005,.0005,25)
		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()
	elif args.mode == 'moving_conductor_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_conductor(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 (nF)')
		plt.ylim([0,1.1*amax(1e12*asarray(Cs))])
		plt.show()
	elif args.mode == 'moving_plate_sweep':
		Cs = []
		seps = [(args.xs/args.nx)*i for i in range(1,10)]
		for sep in seps:
			args.sep = sep
			Cs.append(run_moving_plate(args))
			fn_ind += 1
		print "capacitances: ",list(Cs)
		print "seps",list(seps)
		
		#Cs = [1.0653497726674933e-08, 2.2297769915677895e-09, 1.2768072893947339e-09, 9.1002599382809441e-10, 7.1644594428365008e-10, 5.9725868741684928e-10, 5.1677583773755393e-10, 4.5896058810546686e-10, 4.1553888625802734e-10, 3.8181170895514527e-10, 3.5491456162698833e-10, 3.3300343720448786e-10, 3.1483869957879314e-10, 2.9955703041118239e-10, 2.8653925910417925e-10, 2.7533013154732755e-10, 2.6558770079274896e-10, 2.5705032146192529e-10, 2.4951451716478812e-10]
		#seps = [4e-05, 8e-05, 0.00012000000000000002, 0.00016, 0.0002, 0.00024000000000000003, 0.00028000000000000003, 0.00032, 0.00036, 0.0004, 0.00044, 0.00048000000000000007, 0.0005200000000000001, 0.0005600000000000001, 0.0006000000000000001, 0.00064, 0.00068, 0.00072, 0.00076]

		plt.figure()
		plt.plot(1e3*asarray(seps),1e9*asarray(Cs))
		plt.title('Moving plate capacitance vs. separation')
		plt.xlabel('separation distance (mm)')
		plt.ylabel('capacitance (nF)')
		plt.ylim([0,1.1*amax(1e9*asarray(Cs))])
		plt.show()