Fast inverse transform sampling for an arbitrary probability distribution

This is a fast Python implementation of Inverse Transform Sampling for an arbitrary probability distribution. Generating 10⁶ random numbers with this crazy probability density function takes ~150 ms on a i5-3320M CPU @ 2.6 GHz, most of it from scipy.interpolate.interp1d.

	import numpy as np
	from numpy.random import random
	from scipy import interpolate
	import matplotlib.pyplot as plt
	import cProfile
	def f(x):
	# does not need to be normalized
	return np.exp(-x**2) * np.cos(3*x)**2 * (x-1)**4/np.cosh(1*x)
	def sample(g):
	x = np.linspace(-5,5,1e5)
	y = g(x) # probability density function, pdf
	cdf_y = np.cumsum(y) # cumulative distribution function, cdf
	cdf_y = cdf_y/cdf_y.max() # takes care of normalizing cdf to 1.0
	inverse_cdf = interpolate.interp1d(cdf_y,x) # this is a function
	return inverse_cdf
	def return_samples(N=1e6):
	# let's generate some samples according to the chosen pdf, f(x)
	uniform_samples = random(int(N))
	required_samples = sample(f)(uniform_samples)
	return required_samples
	cProfile.run('return_samples()')
	## plot
	x = np.linspace(-3,3,1e4)
	fig,ax = plt.subplots()
	ax.set_xlabel('x')
	ax.set_ylabel('probability density')
	ax.plot(x,f(x)/np.sum(f(x)*(x[1]-x[0])) )
	ax.hist(return_samples(1e6),bins='auto',normed=True,range=(x.min(),x.max()))
	plt.show()

view raw inverse_transform_probability_distribution.py hosted with ❤ by GitHub

Animated stripchart

This script creates a moving stripchart-style plot of a data stream. Change num_points to increase the length of the section that gets plotted. (TODO: put in a matplotlib slider widget for num_points.)

	# Monitor and plot
	from __future__ import division, print_function
	import numpy as np
	from numpy import sin, cos, random
	from scipy.constants import c,e,h,hbar,u,pi
	import matplotlib.pyplot as plt
	import matplotlib.animation as animation
	import time
	def return_values():
	time.sleep(0.2) # delay for 200 ms
	# replace the following line with a function that returns data from the labjack
	return np.random.uniform()
	# Plot and animate
	t0 = time.time()
	t = []
	t.append(time.time()-t0)
	y = []
	y.append(return_values())
	fig, ax = plt.subplots()
	t0 = time.time()
	num_points = 30 # length of trace to plot
	line, = ax.plot(t[-num_points:],y[-num_points:],lw=2)
	text_template = 'y-value = %.2g'
	text = ax.text(0.15,0.15,text_template%(y[-1]), transform=ax.transAxes, family='monospace',fontsize=20,weight='heavy')
	def animate(i):
	# add new points to time and y-value arrays
	t.append(time.time()-t0)
	y.append(return_values())
	# rescale axes to keep moving
	ax.set_xlim(t[-num_points:][0],t[-1])
	avg = np.average(y[-num_points:])
	variation = max(y[-num_points:]) - min(y[-num_points:])
	ax.set_ylim(avg - 1.2*variation/2, avg + 1.2*variation/2)
	ax.figure.canvas.draw() # this is necessary, to redraw the axes
	# refresh line and text
	line.set_xdata( t[-num_points:] )
	line.set_ydata( y[-num_points:] )
	text.set_text( text_template%(y[-1]) )
	return line, text
	def init(): #Init only required for blitting to give a clean slate.
	line.set_ydata(np.ma.array(t[-num_points:], mask=True))
	text.set_text(" ")
	return line, text
	ani = animation.FuncAnimation(fig, animate, init_func=init,
	interval=0.6, blit=False) # blit has to be False, to avoid ugly artefacts on the text
	plt.show()

view raw monitor_and_plot.py hosted with ❤ by GitHub

Linear least squares phase estimation of a sine wave

This snippet extracts the phase out of a noisy sinusoidal signal. The algorithm is a clever linear least squares (LLS) scheme, courtesy of the IEEE waveform measurement standard for characterizing ADCs and DACs. It works in spite of the fact that the sine wave is non-linear -- which just underscores the point that the linear bit of LLS is linearity in the parameters, not the fit function. It works quite a bit faster than the Levenberg-Marquardt nonlinear least squares routine.

 

	# Least squares fitting
	# IEEE algorithm version for known frequency
	from __future__ import division
	import numpy as np
	from numpy import pi
	from numpy import fft
	from numpy import linalg as lg
	import pylab as plt
	import matplotlib.gridspec as gridspec
	## 3-parameter sine fit, with linear least squares
	dt = 1e-2
	T = np.random.uniform(3,15) # record length
	A = 1.0 # sine amplitude
	noise_amplitude = np.random.uniform(0,20)
	f = np.random.uniform(50,80) # sine frequency
	N = int(T/dt) + 1
	t = np.arange(0,T,dt)
	phase = np.random.uniform(0,2*pi)
	noise = np.random.normal(0,noise_amplitude,N) # noise vector
	y = A * np.sin(2*pi*f*t + phase) + noise
	y = np.matrix(y).T
	D = np.matrix( [np.cos(2*pi*f*t), np.sin(2*pi*f*t), np.ones(N)] ).T
	x_opt = lg.inv(D.T * D) * D.T * y # parameter vector
	a,b,c = x_opt # = A * np.sin(phase), A * np.cos(phase), np.average(y)
	residuals = y - D * x_opt
	S = lg.norm(residuals)**2
	x_cov = np.asscalar( S/(N - len(x_opt)) ) * lg.inv(D.T * D) # covariance matrix
	phi = np.arctan2(a,b)%pi # phase estimate
	delta_a, delta_b = np.sqrt(x_cov[0,0]), np.sqrt(x_cov[1,1])
	delta_phi = np.cos(phi)**2 * np.abs(a/b) * np.sqrt( (delta_a/a)**2 + (delta_b/b)**2 )
	def corr(x,y): return np.dot(x,y)/np.sqrt( np.dot(x,x) * np.dot(y,y) )
	def correlation_length(x):
	l = 0
	while abs(corr( x,np.roll(x,l) )) > 0.5: l = l+1
	return l
	N_eff = N/correlation_length(noise)
	delta_phi_estimated = np.sqrt(2) * np.std(residuals)/np.sqrt( (a**2 + b**2) * N_eff )
	phase_error = phi%pi - phase%pi
	print "SNR = ", round( lg.norm(x_opt[:-1])/np.std(residuals) ,3)
	print
	print "phase error = ",round(phase_error,3), "("+`round(delta_phi,3)`+")", "(fit)"
	print "predicted uncertainty = ", "("+`round(delta_phi_estimated,3)`+")"
	residuals, y, y_reconstructed = np.array(residuals), np.array(y), np.array( D * x_opt )
	## Plotting
	fig = plt.figure()
	gs = gridspec.GridSpec(3,6)
	ax1 = plt.subplot(gs[0:2,:-1])
	ax2 = plt.subplot(gs[2:,:-1])
	ax3 = plt.subplot(gs[2:,-1])
	ax3.xaxis.set_visible(False)
	ax3.yaxis.set_visible(False)
	ax1.plot(t,y,'b',lw=2)
	ax1.plot(t,y_reconstructed,'g',lw=2)
	ax1.set_ylabel("y")
	ax1.grid(axis='y')
	ax2.plot(t,residuals, 'g')
	ax2.set_ylabel("residuals")
	ax2.set_xlabel("Time, t [s]")
	n, bins, patches = ax3.hist(residuals, 50, normed=True, \
	orientation='horizontal', facecolor='green', \
	histtype='stepfilled')
	plt.show()

view raw sine_fitting - known frequency.py hosted with ❤ by GitHub

Exit fluff. Enter IEP !

After dinking around with getting iPython notebooks to work without trouble on my linux laptop, I've finally hit the limit of my patience for the sluggishness and foppish foofery of iPython notebooks. I tried to revert to my old gedit + Dreampie combination, but having been spoiled by executing code inline, I looked around for other alternatives. And found an amazing option: IEP 

The features I like include

• syntax highlighting and the simple, uncluttered interface
• organized file structure for projects
• easy integration with matplotlib (and not having to switch back and forth with %pylab tags to make code run on different machines)
• being able to run selected code snippets with just Alt + Enter !

One word for it. Liberating! 

Arduino gets its Due

After a few years away from Arduinos, I started thinking about them a few days ago as a possible solution for a laser feedback controller project. In the past, I have been unimpressed with the Arduino's ADC, and with its outputs: either PWMs that need to be post-filtered, or klugey DAC chips. And so I was delighted to see the specs of the new(ish) Arduino Due: 12 bit analog ADCs sampling (as fast as ~50 kHz), and real 12 bit DACs right on the board. Finally, something with enough juice to do useful things in the lab!

To celebrate, here's a simple arbitrary waveform generator, taking advantage of the Due Timer library.

	// Arduino Due simple synthesizer
	// arbitrary waveforms on DAC0
	// using Timer library to run off internal clock
	#include <DueTimer.h>
	int led = 13;
	int sensorPin = A0;
	int triggerPin = 4;
	int outputValue = 0;
	int s = 0;
	long samplePeriod = 5; // us. Needs to be >= 4 us
	int *waveform;
	int sine20[] = {2047, 2680, 3250, 3703, 3994, 4095, 3994, 3703, 3250, 2680, 2047, 1414, 844, 391, 100, 0, 100, 391, 844, 1414, -1};
	int cosine100[] = {4095, 4090, 4078, 4058, 4030, 3994, 3951, 3900, 3841, 3776, 3703, 3625, 3540, 3449, 3352, 3250, 3144, 3033, 2919, 2801, 2680, 2556, 2431, 2304, 2176, 2047, 1918, 1790, 1663, 1538, 1414, 1293, 1175, 1061, 950, 844, 742, 645, 554, 469, 391, 318, 253, 194, 143, 100, 64, 36, 16, 4, 0, 4, 16, 36, 64, 100, 143, 194, 253, 318, 391, 469, 554, 645, 742, 844, 950, 1061, 1175, 1293, 1414, 1538, 1663, 1790, 1918, 2047, 2176, 2304, 2431, 2556, 2680, 2801, 2919, 3033, 3144, 3250, 3352, 3449, 3540, 3625, 3703, 3776, 3841, 3900, 3951, 3994, 4030, 4058, 4078, 4090,-1};
	int sine100[] = {2047, 2176, 2304, 2431, 2556, 2680, 2801, 2919, 3033, 3144, 3250, 3352, 3449, 3540, 3625, 3703, 3776, 3841, 3900, 3951, 3994, 4030, 4058, 4078, 4090, 4095, 4090, 4078, 4058, 4030, 3994, 3951, 3900, 3841, 3776, 3703, 3625, 3540, 3449, 3352, 3250, 3144, 3033, 2919, 2801, 2680, 2556, 2431, 2304, 2176, 2047, 1918, 1790, 1663, 1538, 1414, 1293, 1175, 1061, 950, 844, 742, 645, 554, 469, 391, 318, 253, 194, 143, 100, 64, 36, 16, 4, 0, 4, 16, 36, 64, 100, 143, 194, 253, 318, 391, 469, 554, 645, 742, 844, 950, 1061, 1175, 1293, 1414, 1538, 1663, 1790, 1918,-1};
	int square[] = {0,4095,-1};
	int sawtooth512[] = {0, 8, 16, 24, 32, 40, 48, 56, 64, 72, 80, 88, 96, 104, 112, 120, 128, 136, 144, 152, 160, 168, 176, 184, 192, 200, 208, 216, 224, 232, 240, 248, 256, 264, 272, 280, 288, 296, 304, 312, 320, 328, 336, 344, 352, 360, 368, 376, 384, 392, 400, 408, 416, 424, 432, 440, 448, 456, 464, 472, 480, 488, 496, 504, 512, 520, 528, 536, 544, 552, 560, 568, 576, 584, 592, 600, 608, 616, 624, 632, 640, 648, 656, 664, 672, 680, 688, 696, 704, 712, 720, 728, 736, 744, 752, 760, 768, 776, 784, 792, 800, 808, 816, 824, 832, 840, 848, 856, 864, 872, 880, 888, 896, 904, 912, 920, 928, 936, 944, 952, 960, 968, 976, 984, 992, 1000, 1008, 1016, 1024, 1032, 1040, 1048, 1056, 1064, 1072, 1080, 1088, 1096, 1104, 1112, 1120, 1128, 1136, 1144, 1152, 1160, 1168, 1176, 1184, 1192, 1200, 1208, 1216, 1224, 1232, 1240, 1248, 1256, 1264, 1272, 1280, 1288, 1296, 1304, 1312, 1320, 1328, 1336, 1344, 1352, 1360, 1368, 1376, 1384, 1392, 1400, 1408, 1416, 1424, 1432, 1440, 1448, 1456, 1464, 1472, 1480, 1488, 1496, 1504, 1512, 1520, 1528, 1536, 1544, 1552, 1560, 1568, 1576, 1584, 1592, 1600, 1608, 1616, 1624, 1632, 1640, 1648, 1656, 1664, 1672, 1680, 1688, 1696, 1704, 1712, 1720, 1728, 1736, 1744, 1752, 1760, 1768, 1776, 1784, 1792, 1800, 1808, 1816, 1824, 1832, 1840, 1848, 1856, 1864, 1872, 1880, 1888, 1896, 1904, 1912, 1920, 1928, 1936, 1944, 1952, 1960, 1968, 1976, 1984, 1992, 2000, 2008, 2016, 2024, 2032, 2040, 2048, 2056, 2064, 2072, 2080, 2088, 2096, 2104, 2112, 2120, 2128, 2136, 2144, 2152, 2160, 2168, 2176, 2184, 2192, 2200, 2208, 2216, 2224, 2232, 2240, 2248, 2256, 2264, 2272, 2280, 2288, 2296, 2304, 2312, 2320, 2328, 2336, 2344, 2352, 2360, 2368, 2376, 2384, 2392, 2400, 2408, 2416, 2424, 2432, 2440, 2448, 2456, 2464, 2472, 2480, 2488, 2496, 2504, 2512, 2520, 2528, 2536, 2544, 2552, 2560, 2568, 2576, 2584, 2592, 2600, 2608, 2616, 2624, 2632, 2640, 2648, 2656, 2664, 2672, 2680, 2688, 2696, 2704, 2712, 2720, 2728, 2736, 2744, 2752, 2760, 2768, 2776, 2784, 2792, 2800, 2808, 2816, 2824, 2832, 2840, 2848, 2856, 2864, 2872, 2880, 2888, 2896, 2904, 2912, 2920, 2928, 2936, 2944, 2952, 2960, 2968, 2976, 2984, 2992, 3000, 3008, 3016, 3024, 3032, 3040, 3048, 3056, 3064, 3072, 3080, 3088, 3096, 3104, 3112, 3120, 3128, 3136, 3144, 3152, 3160, 3168, 3176, 3184, 3192, 3200, 3208, 3216, 3224, 3232, 3240, 3248, 3256, 3264, 3272, 3280, 3288, 3296, 3304, 3312, 3320, 3328, 3336, 3344, 3352, 3360, 3368, 3376, 3384, 3392, 3400, 3408, 3416, 3424, 3432, 3440, 3448, 3456, 3464, 3472, 3480, 3488, 3496, 3504, 3512, 3520, 3528, 3536, 3544, 3552, 3560, 3568, 3576, 3584, 3592, 3600, 3608, 3616, 3624, 3632, 3640, 3648, 3656, 3664, 3672, 3680, 3688, 3696, 3704, 3712, 3720, 3728, 3736, 3744, 3752, 3760, 3768, 3776, 3784, 3792, 3800, 3808, 3816, 3824, 3832, 3840, 3848, 3856, 3864, 3872, 3880, 3888, 3896, 3904, 3912, 3920, 3928, 3936, 3944, 3952, 3960, 3968, 3976, 3984, 3992, 4000, 4008, 4016, 4024, 4032, 4040, 4048, 4056, 4064, 4072, 4080, 40880, 8, 16, 24, 32, 40, 48, 56, 64, 72, 80, 88, 96, 104, 112, 120, 128, 136, 144, 152, 160, 168, 176, 184, 192, 200, 208, 216, 224, 232, 240, 248, 256, 264, 272, 280, 288, 296, 304, 312, 320, 328, 336, 344, 352, 360, 368, 376, 384, 392, 400, 408, 416, 424, 432, 440, 448, 456, 464, 472, 480, 488, 496, 504, 512, 520, 528, 536, 544, 552, 560, 568, 576, 584, 592, 600, 608, 616, 624, 632, 640, 648, 656, 664, 672, 680, 688, 696, 704, 712, 720, 728, 736, 744, 752, 760, 768, 776, 784, 792, 800, 808, 816, 824, 832, 840, 848, 856, 864, 872, 880, 888, 896, 904, 912, 920, 928, 936, 944, 952, 960, 968, 976, 984, 992, 1000, 1008, 1016, 1024, 1032, 1040, 1048, 1056, 1064, 1072, 1080, 1088, 1096, 1104, 1112, 1120, 1128, 1136, 1144, 1152, 1160, 1168, 1176, 1184, 1192, 1200, 1208, 1216, 1224, 1232, 1240, 1248, 1256, 1264, 1272, 1280, 1288, 1296, 1304, 1312, 1320, 1328, 1336, 1344, 1352, 1360, 1368, 1376, 1384, 1392, 1400, 1408, 1416, 1424, 1432, 1440, 1448, 1456, 1464, 1472, 1480, 1488, 1496, 1504, 1512, 1520, 1528, 1536, 1544, 1552, 1560, 1568, 1576, 1584, 1592, 1600, 1608, 1616, 1624, 1632, 1640, 1648, 1656, 1664, 1672, 1680, 1688, 1696, 1704, 1712, 1720, 1728, 1736, 1744, 1752, 1760, 1768, 1776, 1784, 1792, 1800, 1808, 1816, 1824, 1832, 1840, 1848, 1856, 1864, 1872, 1880, 1888, 1896, 1904, 1912, 1920, 1928, 1936, 1944, 1952, 1960, 1968, 1976, 1984, 1992, 2000, 2008, 2016, 2024, 2032, 2040, 2048, 2056, 2064, 2072, 2080, 2088, 2096, 2104, 2112, 2120, 2128, 2136, 2144, 2152, 2160, 2168, 2176, 2184, 2192, 2200, 2208, 2216, 2224, 2232, 2240, 2248, 2256, 2264, 2272, 2280, 2288, 2296, 2304, 2312, 2320, 2328, 2336, 2344, 2352, 2360, 2368, 2376, 2384, 2392, 2400, 2408, 2416, 2424, 2432, 2440, 2448, 2456, 2464, 2472, 2480, 2488, 2496, 2504, 2512, 2520, 2528, 2536, 2544, 2552, 2560, 2568, 2576, 2584, 2592, 2600, 2608, 2616, 2624, 2632, 2640, 2648, 2656, 2664, 2672, 2680, 2688, 2696, 2704, 2712, 2720, 2728, 2736, 2744, 2752, 2760, 2768, 2776, 2784, 2792, 2800, 2808, 2816, 2824, 2832, 2840, 2848, 2856, 2864, 2872, 2880, 2888, 2896, 2904, 2912, 2920, 2928, 2936, 2944, 2952, 2960, 2968, 2976, 2984, 2992, 3000, 3008, 3016, 3024, 3032, 3040, 3048, 3056, 3064, 3072, 3080, 3088, 3096, 3104, 3112, 3120, 3128, 3136, 3144, 3152, 3160, 3168, 3176, 3184, 3192, 3200, 3208, 3216, 3224, 3232, 3240, 3248, 3256, 3264, 3272, 3280, 3288, 3296, 3304, 3312, 3320, 3328, 3336, 3344, 3352, 3360, 3368, 3376, 3384, 3392, 3400, 3408, 3416, 3424, 3432, 3440, 3448, 3456, 3464, 3472, 3480, 3488, 3496, 3504, 3512, 3520, 3528, 3536, 3544, 3552, 3560, 3568, 3576, 3584, 3592, 3600, 3608, 3616, 3624, 3632, 3640, 3648, 3656, 3664, 3672, 3680, 3688, 3696, 3704, 3712, 3720, 3728, 3736, 3744, 3752, 3760, 3768, 3776, 3784, 3792, 3800, 3808, 3816, 3824, 3832, 3840, 3848, 3856, 3864, 3872, 3880, 3888, 3896, 3904, 3912, 3920, 3928, 3936, 3944, 3952, 3960, 3968, 3976, 3984, 3992, 4000, 4008, 4016, 4024, 4032, 4040, 4048, 4056, 4064, 4072, 4080, 4088,-1};
	int random64[] = {1964, 1366, 646, 1394, 1480, 1632, 1456, 1356, 1257, 1953, 862, 384, 2021, 1974, 45, 37, 1844, 1222, 1435, 77, 1061, 1128, 376, 160, 1198, 365, 1315, 658, 1837, 1800, 656, 239, 691, 823, 446, 355, 1870, 840, 354, 1979, 1830, 1727, 1029, 1642, 1453, 344, 2, 1343, 471, 1227, 803, 1273, 311, 1191, 1919, 1513, 1365, 536, 338, 467, 1175, 237, 1301, 1375, -1};
	void setup() {
	// Pick a waveform
	waveform = sine20;
	analogWriteResolution(12);
	Timer0.attachInterrupt(wave).start(samplePeriod);
	// setting time intervals is more reliable than frequency, because of integer math
	// Timer1.attachInterrupt(trigger).start(20*samplePeriod);
	}
	void wave(){
	outputValue = *(waveform + s);
	if (outputValue==-1){ // wrap around and repeat the waveform
	s = 0;
	outputValue = *(waveform + s);
	}
	analogWrite(DAC0, outputValue); // write the selected waveform on DAC1
	s++;
	}
	void trigger(){
	digitalWrite(triggerPin, HIGH);
	delayMicroseconds(50);
	digitalWrite(triggerPin, LOW);
	}
	void loop(){}

view raw sine_test.ino hosted with ❤ by GitHub

Impedance of transmission lines

This code calculates the characteristic impedance of transmission lines of specified cross section, by calculating the electric potential through relaxation on a grid. You'll need to be able to run blitz from scipy.weave.

Here is an example of what it can do.

	# Transmission line characteristic impedance
	# Amar, 12 July 2012
	# version 3
	# All lengths in cm, E-fields in V/cm, time in ns
	import pylab as plt
	import numpy as np
	from numpy import pi
	from scipy.weave import blitz
	import cPickle, time
	eta0 = 376.730313461 # Ohms
	"""
	# Blitz test:
	a,b = 1,0
	print blitz("b = a+a")
	"""
	################### Computation grid ##########################
	size = 50.0 # mm
	numberOfGridBits = 10
	numberOfGridPoints = 2**numberOfGridBits
	dx = dy = 2*size/numberOfGridPoints
	x = np.arange(-size/2.,size/2.,dx)
	y = np.arange(-size/2.,size/2.,dy)
	X,Y = np.meshgrid(x,y)
	N = len(x)
	################## Electrodes #################################
	class CircularElectrode:
	def __init__(self,x0,y0,D,V0,innie=True):
	self.x0, self.y0, self.V0 = x0,y0,V0
	self.D = D
	if innie: self.mask = ( ((X-self.x0)**2 + (Y-self.y0)**2) <= self.D**2/4. )
	else: self.mask = ( ((X-self.x0)**2 + (Y-self.y0)**2) >= self.D**2/4. )
	#electrodes = [ CircularElectrode(0,0,.24425,V0,innie=True), CircularElectrode(0,0,.5625,0,innie=False) ] # GR900
	class RectangularElectrode:
	def __init__(self,x0,y0,w,h,V0,innie=True):
	self.x0, self.y0, self.V0 = x0,y0,V0
	self.w,self.h = w,h
	if innie: self.mask = np.logical_and( np.abs(X-self.x0) <= self.w/2., np.abs(Y-self.y0) <= self.h/2. )
	else: self.mask = np.logical_or( np.abs(X-self.x0) >= self.w/2., np.abs(Y-self.y0) >= self.h/2. )
	# Note this subtleness between logical_and vs. logical_or for innies vs. outies
	################# Calculation routine ########################
	expr = "V[step:-step, step:-step] = ((V[0:-2*step, step:-step] + V[2*step:, step:-step]) + (V[step:-step,0:-2*step] + V[step:-step, 2*step:]))/4.0"
	expr2 = "V = V_mask_inv*V + V_mask*V_electrode"
	def calculateVoltage(electrodes, V0=1.0):
	"""Calculation of voltage for an electrode configuration"""
	V_electrode = np.zeros(X.shape)
	for electrode in electrodes: V_electrode[electrode.mask] = electrode.V0
	V_mask = sum([electrode.mask for electrode in electrodes])
	V_mask_inv = np.logical_not(V_mask)
	V = V_mask_inv*np.average(V_electrode) # initialized
	print "Preconditioning ...",
	for i in xrange(4,numberOfGridBits):
	step = 2**numberOfGridBits/2**i
	for j in xrange(int(5e1)):
	blitz(expr, check_size=0)
	blitz(expr2,check_size=0)
	print "done"
	step = 1.0
	for i in xrange(int(4e3)):
	blitz(expr, check_size=0)
	blitz(expr2,check_size=0)
	V_unref = V.copy()
	for i in xrange(int(1e3)):
	blitz(expr, check_size=0)
	blitz(expr2,check_size=0)
	V_ref = V.copy()
	fractionalChange = np.sqrt( np.sum((V_ref - V_unref)**2)/np.sum(V_unref**2) )
	# Electric field
	ESquared = np.zeros(X.shape)
	ESquared[1:-1, 1:-1] = ((V[2:, 1:-1] - V[0:-2, 1:-1])/(2*dx))**2 + ((V[1:-1, 2:] - V[1:-1,0:-2])/(2*dy))**2
	# Characteristic impedance
	Z0 = eta0/( np.sum(ESquared)*dx*dy ) * V0**2
	return Z0, V, (fractionalChange, V_ref-V_unref) # characteristic impedance, Ohms
	##################### Test problem #############################
	V0 = 1.0
	outerConductor = CircularElectrode(0,0,46.06,0,innie=False)
	innerConductor = CircularElectrode(0,0,20.00,V0,innie=True)
	electrodes = [ outerConductor, innerConductor ]
	Z0,V,convergence = calculateVoltage(electrodes,V0)
	print "Z0 = ", round(Z0,2),
	print "converged to ", convergence[0]
	V0 = 2.0
	outerConductor = RectangularElectrode(0,0,49.6,49.6,0,innie=False)
	innerConductor = RectangularElectrode(0,0,19.8,19.8,V0,innie=True)
	electrodes = [ outerConductor, innerConductor ]
	Z0,V,convergence = calculateVoltage(electrodes,V0)
	print "Z0 = ", round(Z0,2),
	print "converged to ", convergence[0]
	######### Plots ##############################
	# Electric field meshgrid
	Ex,Ey = np.zeros(X.shape),np.zeros(X.shape)
	Ex[1:-1, 1:-1] = (V[2:, 1:-1] - V[0:-2, 1:-1])/(2*dx)
	Ey[1:-1, 1:-1] = (V[1:-1, 2:] - V[1:-1,0:-2])/(2*dy)
	E = np.sqrt( Ex**2 + Ey**2 )
	plt.figure()
	plt.title(r"Voltage (V)")
	plt.xlabel("x (mm)")
	plt.ylabel("y (mm)")
	C1 = plt.contourf(X,Y,V,30)
	plt.colorbar(C1)
	p.figure()
	p.title(r"Electric field, $\mathcal{E}_x$ (V/cm)")
	p.xlabel("x (mm)")
	p.ylabel("y (mm)")
	C2 = p.contourf(X,Y,E,30)
	p.colorbar(C2)
	##################### Batch problem setup ######################
	outerConductor = RectangularElectrode(0,0,30,50,0,innie=False)
	V0 = 1.0
	w = np.arange(5,25,2)
	h = np.arange(1,15,2)
	W,H = np.meshgrid(w,h)
	Z = np.zeros(W.shape)
	for i in xrange(len(w)):
	for j in xrange(len(h)):
	innerConductor = RectangularElectrode(0,-15,w[i],h[j],V0,innie=True)
	electrodes = [ outerConductor, innerConductor ]
	Z0,V,junk = calculateVoltage(electrodes)
	Z[j,i] = Z0
	print w[i],h[j],round(Z0)
	p.figure()
	p.title(r"Characteristic impedance ($\Omega$)")
	p.xlabel("Center conductor width, w (mm)")
	p.ylabel("Center conductor height, h (mm)")
	C = p.contourf(W,H,Z,30)
	p.colorbar(C)
	p.figure()
	p.title("Characteristic impedance ($\Omega$)\n"+"Inner conductor @ "+`(innerConductor.x0,innerConductor.y0)`)
	p.xlabel("Center conductor width, w (mm)")
	p.ylabel("Center conductor height, h (mm)")
	C = p.contour(W,H,Z,np.arange(0,100,5), linewidth=10, colors='b')
	p.clabel(C, inline=1, fontsize=12, linewidth=5)

view raw characteristicImpedanceTEM.py hosted with ❤ by GitHub