Gitorious


2	2	#
3	3	*~
4	4	*.png
	5	result*
	6	*.bak


4	4	caller = sys._getframe(1).f_globals["__name__"]
5	5	exec "import %s as target" % caller
6	6
7		target.timeit = __import__("timeit")
8
9	7	modules = (
10	8	"sys",
11	9	"os",
…	…
20	20	target.matplotlib.use("Agg")
21	21	target.mpmath.mp.dps = 100
22	22
23		target.__dict__["pylab"] = __import__("pylab")
24		target.__dict__["pprint"] = __import__("pprint").pprint
	23	target.__dict__["pylab"] = __import__("pylab") # this should be eval after setting Agg backend
	24	target.__dict__["pprint"] = __import__("pprint").pprint # import function
	25
	26
25	27
26	28
27	29	# below is obsolete


		1	#!/usr/bin/env python
		2	#-- mode: python --
		3
		4	import convention
		5
		6	from fdtd import grid, mesher
		7
		8	# program start
		9
		10	data = numpy.zeros((100,100))
		11	print data
		12
		13	a_plane = grid.Plane(data)
		14	print a_plane
		15	print a_plane.__class__


		1	#!/usr/bin/env python
		2	#-- mode: python --
		3
		4	import numpy
		5
		6	class Plane(numpy.ndarray):
		7	"""
		8
		9	"""
		10	def __new__(cls, abuffer):
		11	"""
		12	create a new plane object
		13	"""
		14	abuffer = numpy.array(abuffer, dtype="float128")
		15	if len(abuffer.shape) == 2:
		16	return numpy.ndarray.__new__(cls, abuffer.shape, "float128", abuffer)
		17	else:
		18	pass
		19
		20	def __array_finalize__(self, obj):
		21	pass
		22
		23	@property
		24	def name(self, name):
		25	self.name = name
		26	return self
		27
		28
		29	p = Plane([[1,2,],[3,4]])
		30	print p
		31	print p.shape
		32	print p.dtype
		33	print p.__class__


		1	from numpy import *
		2
		3	class OneDim(ndarray):
		4	""" ONE dimension grid object """
		5
		6	def __init__(self, *args):
		7	"""
		8	"""
		9	pass
		10
		11	# {{{
		12	class Plane(ndarray):
		13	""" Two dimension grid object """
		14
		15	def __new__(cls, abuffer=None):
		16	abuffer = array(abuffer, dtype="float128")
		17	print abuffer
		18	if len(abuffer.shape) == 2:
		19	return ndarray.__new__(cls, abuffer.shape, "float128", abuffer)
		20	else:
		21	# TODO:
		22	# should raise an error here
		23	pass
		24
		25	def __init__(self, abuffer):
		26	"""
		27	Arguments:
		28	- `abuffer`:
		29	"""
		30	pass
		31	# }}}
		32
		33
		34	class Point(object):
		35	""" Grid point object used in grid array """
		36
		37	def __init__(self, properties = {}):
		38	"""
		39	pass a properties hash to initialize it
		40	Arguments:
		41	- `properties`: hash
		42	"""
		43	self._properties = properties
		44	self.eps = properties["eps"]
		45	self.mu = properties["mu"]
		46	self.sigma = properties["sigma"]
		47


5	5
6	6	from fdtd import source, grid, mesher
7	7
8		pprint(locals())
	8	from fdtd.algorithm.twodim.freespace import update_efield, update_hfield
9	9
10	10	# program start
11	11
12		p = grid.Plane(numpy.zeros((100,100)))
13		print(p)
14		p = mesher.plane(numpy.zeros((10,10)))
15		print(p)
	12	new_efield = grid.Plane(shape=(301,301), timestep=1, spacestep=2310**8)
	13	new_hfield = grid.Plane(shape=(301,301), timestep=1, spacestep=2310**8)
	14
	15
	16	for step in range(1,300):
	17	new_efield.y[151,151] += float(mpmath.sin(stepmpmath.pi/6.125))2
	18	new_efield = update_efield( new_efield, new_hfield )
	19	new_hfield.z[151,151] += float(mpmath.sin(stepmpmath.pi/6.125))2
	20	new_hfield = update_hfield( new_efield, new_hfield )
	21
	22
	23	h = numpy.array(new_hfield.z, dtype="float")
	24	for i in xrange(1,300):
	25	for j in xrange(1,300):
	26	if (i+j)%2 == 1:
	27	h[i,j] = ( h[i+1,j] + h[i-1,j] + h[i,j+1] + h[i,j-1] ) / 4
	28	pylab.imshow( h, interpolation="quadric", norm=matplotlib.colors.Normalize(-0.3,0.3,True) )
	29	pylab.colorbar()
	30	pylab.grid(True)
	31	pylab.savefig("result/fdtd_plane-%s.png" % str(step))
	32	pylab.clf()
	33
	34	d = numpy.array(new_hfield.z[0:,151], dtype="float")
	35	print d[140:159]
	36	for i in xrange(1,len(d)-1):
	37	if i%2==0:
	38	d[i] = ( d[i-1] + d[i+1] ) / 2
	39	print d[140:159]
	40	pylab.plot( d )
	41	pylab.ylim((-0.3,0.3))
	42	pylab.grid(True)
	43	pylab.savefig("result_slice/fdtd_plane_slice-%s.png" % str(step))
	44	pylab.clf()
	45
	46	os.system("open result")
	47	os.system("open result_slice")
	48	print "done"


	1		from mpmath import *


		1	import mpmath
		2	import math
		3
		4	from fdtd import grid
		5
		6	def update_efield( efield, hfield, region=None ):
		7	"""
		8	"""
		9	new_efield = grid.Plane( abuffer = efield )
		10	(xaxis,yaxis) = efield.x.shape
		11	for i in range(1,xaxis-1):
		12	for j in range(1,yaxis-1):
		13	if (i+j)%2 == 1:
		14	pass
		15	return new_efield
		16
		17
		18	def update_hfield( efield, hfield, region=None ):
		19	"""
		20	"""
		21	new_hfield = grid.Plane( abuffer = hfield )
		22	return new_hfield


		1	import mpmath
		2	import math
		3
		4	from fdtd import grid
		5
		6	def update_efield( efield, hfield ):
		7	"""
		8	Two dimension efield update equation
		9
		10	all points at ( x + y ) % 2 == 1 would be updated.
		11
		12	Arguments:
		13	- efield: previous efield
		14	- hfield: previous hfield
		15	"""
		16	new_efield = grid.Plane( abuffer=efield )
		17
		18	(xaxis,yaxis) = efield.x.shape
		19
		20	for i in range(1,xaxis-1):
		21	for j in range(1,yaxis-1):
		22	if (i+j)%2 == 1:
		23	curl_coefficient = (new_efield.timestep / (new_efield.spacestep * new_efield.eps[i,j]))
		24	new_efield.x[i,j] = efield.x[i,j] + curl_coefficient * ( hfield.z[i,j+1] - hfield.z[i,j-1] )
		25	new_efield.y[i,j] = efield.y[i,j] - curl_coefficient * ( hfield.z[i+1,j] - hfield.z[i-1,j] )
		26
		27	return new_efield
		28
		29
		30	def update_hfield( efield, hfield ):
		31	"""
		32	Two dimension hfield update equation
		33
		34	all points at ( x + y ) % 2 == 0 would be updated.
		35	"""
		36	new_hfield = grid.Plane( abuffer=hfield )
		37
		38	(xaxis,yaxis) = hfield.x.shape
		39
		40	for i in range(1,xaxis-1):
		41	for j in range(1,yaxis-1):
		42	if (i+j)%2 == 0:
		43	curl_coefficient = new_hfield.timestep / (new_hfield.spacestep * new_hfield.mu[i,j])
		44	new_hfield.z[i,j] = hfield.z[i,j] + curl_coefficient * ( efield.x[i,j+1] - efield.x[i,j-1] ) - curl_coefficient * ( efield.y[i+1,j] - efield.y[i-1,j] )
		45
		46	return new_hfield


		1	"""
		2	fdtd.grid.py
		3
		4	Implement three kinds of grid object
		5	Line as 1-D grid
		6	Plane as 2-D grid
		7	Cube as 3-D grid
		8
		9	"""
		10
		11	import numpy
		12	from scipy.constants import epsilon_0, mu_0
		13
		14	# {{{
		15	class Line(object):
		16	"""
		17	One dimension grid object
		18	"""
		19
		20	def __init__(self, step = {}):
		21	"""
		22	"""
		23
		24	def timestep(self, ):
		25	"""
		26	"""
		27	# }}}
		28
		29
		30	# {{{
		31	class Plane(object):
		32	"""
		33	Two dimension grid object
		34	"""
		35
		36	def __init__(self, abuffer=None, shape=None, eps=None, mu=None, sigmae=None, sigmah=None, timestep=None, spacestep=None):
		37	"""
		38	Give one of shape or abuffer to initialize Plane object.
		39
		40	Once abuffer is given, shape and other parameters would be obsolete.
		41
		42	Alternatively, if abuffer is None, Plane would be initialized with rest parameters.
		43
		44	If both abuffer and shape were not given, Error would be raised.
		45
		46	Arguments:
		47	- abuffer: grid buffer used to initialize plane
		48	- shape: specify plance shape
		49	- eps: primary permitivity value would be assign to whole region
		50	- mu: primary permeability value would be assign to whole region
		51	- sigmae: primary lossy coefficient of E field
		52	- sigmah: primary lossy coefficient of H field
		53	"""
		54	if abuffer and len(abuffer.x.shape) == 2:
		55	self.x = numpy.zeros_like(abuffer.x)
		56	self.y = numpy.zeros_like(abuffer.y)
		57	self.z = numpy.zeros_like(abuffer.z)
		58
		59	self.eps = abuffer.eps
		60	self.mu = abuffer.mu
		61	self.sigmae = abuffer.sigmae
		62	self.sigmah = abuffer.sigmah
		63
		64	self.timestep = abuffer.timestep
		65	self.spacestep = abuffer.spacestep
		66
		67	elif shape:
		68	self.x = numpy.zeros(shape, dtype="float128")
		69	self.y = numpy.zeros(shape, dtype="float128")
		70	self.z = numpy.zeros(shape, dtype="float128")
		71
		72	self.eps = numpy.ndarray(shape=self.x.shape, dtype="float128"); self.eps[0:,0:] = (eps or epsilon_0)
		73	self.mu = numpy.ndarray(shape=self.x.shape, dtype="float128"); self.mu[0:,0:] = (mu or mu_0)
		74	self.sigmae = numpy.ndarray(shape=self.x.shape, dtype="float128"); self.sigmae[0:,0:] = (sigmae or 0)
		75	self.sigmah = numpy.ndarray(shape=self.x.shape, dtype="float128"); self.sigmah[0:,0:] = (sigmah or 0)
		76
		77	self.timestep = (timestep or 1)
		78	self.spacestep = (spacestep or 1)
		79	pass
		80
		81	# }}}
		82
		83	# {{{
		84	class Cube(object):
		85	""" Three dimension grid object """
		86
		87	def __init__(self, abuffer):
		88	"""
		89
		90	Arguments:
		91	- `abuffer`:
		92	"""
		93	self._abuffer = abuffer
		94
		95	# }}}
		96
		97
		98
		99	# {{{
		100	class Point(object):
		101	""" Grid point object used in grid array """
		102
		103	def __init__(self, properties = {}):
		104	"""
		105	pass a properties hash to initialize it
		106	Arguments:
		107	- `properties`: hash
		108	"""
		109	self._properties = properties
		110	self.eps = properties["eps"]
		111	self.mu = properties["mu"]
		112	self.sigma = properties["sigma"]
		113
		114	# }}}


		1	from numpy import *
		2	from fdtd import grid
		3
		4	def line(abuffer):
		5	"""
		6
		7	Arguments:
		8	- `abuffer`:
		9	"""
		10	return grid.Line(abuffer)
		11
		12
		13	def plane(abuffer):
		14	"""
		15
		16	Arguments:
		17	- `abuffer`:
		18	"""
		19	return grid.Plane(abuffer)
		20
		21	def cube(abuffer):
		22	"""
		23
		24	Arguments:
		25	- `abuffer`:
		26	"""
		27	return grid.Cube(abuffer)


1	1	#!/usr/bin/env python
2	2	#-- mode: python --
3
4
5	3	import convention
6	4
7	5	from fdtd.source import *
8	6
9	7	# program start
	8	t = numpy.arange(-10, 10, 0.01) # x axis range
10	9
11	10	# plot polynomial
12		t = numpy.arange(-10, 10, 0.001)
13
14	11	pylab.plot(t, map( polynomial_pulse, t ) )
15	12	pylab.plot(t, map( polynomial_pulse_p, t ) )
16	13	pylab.plot(t, map( polynomial_pulse_pp, t ) )
17	14
18	15	pylab.xlim(-2, 2)
19		# pylab.ylim(-0.5, 1.5)
20	16	pylab.grid(True)
21	17	pylab.savefig("poly.png")
22		pylab.gcf().clear()
	18	pylab.clf()
23	19	os.system("open poly.png")
24	20
25	21	# plot gaussian
26	22	one = numpy.ones_like(t)
27	23
28		pylab.plot(t, map( gaussian, t, one2, one3 ) )
29		pylab.plot(t, map( gaussian_p, t, one2, one3 ) )
	24	pylab.plot(t, [ gaussian(x,2,3) for x in t ])
	25	pylab.plot(t, [ gaussian_p(x,2,3,) for x in t ])
30	26
31	27	pylab.grid(True)
32	28	pylab.savefig("gaussian.png")
33		pylab.gcf().clear()
	29	pylab.clf()
34	30	os.system("open gaussian.png")
	31
	32	# simple fourier transform a step function
	33	s = numpy.zeros(1000)
	34	s[:100] = 1
	35	pylab.plot(s)
	36	pylab.plot(scipy.fft(s))
	37	pylab.grid(True)
	38	pylab.savefig("simple_fft.png")
	39	pylab.clf()
	40	os.system("open simple_fft.png")
	41
	42	# plot fourier transform of gaussian
	43
	44	g = numpy.array([ gaussian(x,2,3) for x in t ], dtype="float")
	45
	46	pylab.plot(t, scipy.fft(g))
	47
	48	pylab.savefig("gaussian_p_fft.png")
	49	os.system("open gaussian_p_fft.png")

Commit b607d4677360269a7ca8f7d82786454dd3712945

Gitorious