zcy
/
openEMS
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'python'

pull/29/head
Thorsten Liebig 2017-03-01 22:09:46 +01:00
parent 9741b654e0 7c0d498f56
commit 9c317dc4a8
53 changed files with 3378 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
.gitignore vendored
View File

@ -12,7 +12,6 @@ Makefile*
*.pro.user* *.pro.user*
*.user *.user
*.orig *.orig
openEMS
localPaths.pri localPaths.pri
.directory .directory
@ -22,3 +21,11 @@ CMakeCache.txt
cmake_install.cmake cmake_install.cmake
install_manifest.txt install_manifest.txt
localConfig.cmake localConfig.cmake
#python
*.pyc
*.pyo
python/**/*.cpp
!python/doc
python/doc/_build
python/doc/Tutorials/__*

14
python/README.md Normal file
View File

@ -0,0 +1,14 @@
# openEMS python interface
## Install
* Simple version:
```python
python setup.py install
```
* Extended options, e.g. for custom install path at */opt/openEMS*:
```python
python setup.py build_ext -I/opt/openEMS/include -L/opt/openEMS/lib -R/opt/openEMS/lib"
pyhton setup.py install
```
**Note:** The install command may require root on Linux, or add --user to install to ~/.local

198
python/Tutorials/Bent_Patch_Antenna.py Normal file
View File

@ -0,0 +1,198 @@
# -*- coding: utf-8 -*-
"""
Bent Patch Antenna Tutorial
Tested with
- python 3.4
- openEMS v0.0.33+
(C) 2016 Thorsten Liebig <thorsten.liebig@gmx.de>
"""
### Import Libraries
import os, tempfile
from pylab import *
from mpl_toolkits.mplot3d import Axes3D
from CSXCAD import CSXCAD
from openEMS.openEMS import openEMS
from openEMS.physical_constants import *
### Setup the simulation
Sim_Path = os.path.join(tempfile.gettempdir(), 'Bent_Patch')
post_proc_only = False
unit = 1e-3 # all length in mm
f0 = 2.4e9 # center frequency, frequency of interest!
lambda0 = round(C0/f0/unit) # wavelength in mm
fc = 0.5e9 # 20 dB corner frequency
# patch width in alpha-direction
patch_width = 32 # resonant length in alpha-direction
patch_radius = 50 # radius
patch_length = 40 # patch length in z-direction
#substrate setup
substrate_epsR = 3.38
substrate_kappa = 1e-3 * 2*pi*2.45e9 * EPS0*substrate_epsR
substrate_width = 80
substrate_length = 90
substrate_thickness = 1.524
substrate_cells = 4
#setup feeding
feed_pos = -5.5 #feeding position in x-direction
feed_width = 2 #feeding port width
feed_R = 50 #feed resistance
# size of the simulation box
SimBox_rad = 2*100
SimBox_height = 1.5*200
### Setup FDTD parameter & excitation function
FDTD = openEMS(CoordSystem=1) # init a cylindrical FDTD
f0 = 2e9 # center frequency
fc = 1e9 # 20 dB corner frequency
FDTD.SetGaussExcite(f0, fc)
FDTD.SetBoundaryCond(['MUR', 'MUR', 'MUR', 'MUR', 'MUR', 'MUR']) # boundary conditions
### Setup the Geometry & Mesh
# init a cylindrical mesh
CSX = CSXCAD.ContinuousStructure(CoordSystem=1)
FDTD.SetCSX(CSX)
mesh = CSX.GetGrid()
mesh.SetDeltaUnit(unit)
### Setup the geometry using cylindrical coordinates
# calculate some width as an angle in radiant
patch_ang_width = patch_width/(patch_radius+substrate_thickness)
substr_ang_width = substrate_width/patch_radius
feed_angle = feed_pos/patch_radius
# create patch
patch = CSX.AddMetal('patch') # create a perfect electric conductor (PEC)
start = [patch_radius+substrate_thickness, -patch_ang_width/2, -patch_length/2 ]
stop = [patch_radius+substrate_thickness, patch_ang_width/2, patch_length/2 ]
CSX.AddBox(patch, priority=10, start=start, stop=stop, edges2grid='all') # add a box-primitive to the metal property 'patch'
# create substrate
substrate = CSX.AddMaterial('substrate', epsilon=substrate_epsR, kappa=substrate_kappa )
start = [patch_radius , -substr_ang_width/2, -substrate_length/2]
stop = [patch_radius+substrate_thickness, substr_ang_width/2, substrate_length/2]
substrate.AddBox(start=start, stop=stop, edges2grid='all')
# save current density oon the patch
jt_patch = CSX.AddDump('Jt_patch', dump_type=3, file_type=1)
start = [patch_radius+substrate_thickness, -substr_ang_width/2, -substrate_length/2]
stop = [patch_radius+substrate_thickness, +substr_ang_width/2, substrate_length/2]
jt_patch.AddBox(start=start, stop=stop)
# create ground
gnd = CSX.AddMetal('gnd') # create a perfect electric conductor (PEC)
start = [patch_radius, -substr_ang_width/2, -substrate_length/2]
stop = [patch_radius, +substr_ang_width/2, +substrate_length/2]
gnd.AddBox(priority=10, start=start, stop=stop, edges2grid='all')
# apply the excitation & resist as a current source
start = [patch_radius , feed_angle, 0]
stop = [patch_radius+substrate_thickness, feed_angle, 0]
port = FDTD.AddLumpedPort(1 ,feed_R, start, stop, 'r', 1.0, priority=50, edges2grid='all')
### Finalize the Mesh
# add the simulation domain size
mesh.AddLine('r', patch_radius+np.array([-20, SimBox_rad]))
mesh.AddLine('a', [-0.75*pi, 0.75*pi])
mesh.AddLine('z', [-SimBox_height/2, SimBox_height/2])
# add some lines for the substrate
mesh.AddLine('r', patch_radius+np.linspace(0,substrate_thickness,substrate_cells))
# generate a smooth mesh with max. cell size: lambda_min / 20
max_res = C0 / (f0+fc) / unit / 20
max_ang = max_res/(SimBox_rad+patch_radius) # max res in radiant
mesh.SmoothMeshLines(0, max_res, 1.4)
mesh.SmoothMeshLines(1, max_ang, 1.4)
mesh.SmoothMeshLines(2, max_res, 1.4)
## Add the nf2ff recording box
nf2ff = FDTD.CreateNF2FFBox()
### Run the simulation
if 0: # debugging only
CSX_file = os.path.join(Sim_Path, 'bent_patch.xml')
if not os.path.exists(Sim_Path):
os.mkdir(Sim_Path)
CSX.Write2XML(CSX_file)
os.system(r'AppCSXCAD "{}"'.format(CSX_file))
if not post_proc_only:
FDTD.Run(Sim_Path, verbose=3, cleanup=True)
### Postprocessing & plotting
f = np.linspace(max(1e9,f0-fc),f0+fc,401)
port.CalcPort(Sim_Path, f)
Zin = port.uf_tot / port.if_tot
s11 = port.uf_ref/port.uf_inc
s11_dB = 20.0*np.log10(np.abs(s11))
figure()
plot(f/1e9, s11_dB)
grid()
ylabel('s11 (dB)')
xlabel('frequency (GHz)')
P_in = 0.5*np.real(port.uf_tot * np.conj(port.if_tot)) # antenna feed power
# plot feed point impedance
figure()
plot( f/1e6, real(Zin), 'k-', linewidth=2, label=r'$\Re(Z_{in})$' )
grid()
plot( f/1e6, imag(Zin), 'r--', linewidth=2, label=r'$\Im(Z_{in})$' )
title( 'feed point impedance' )
xlabel( 'frequency (MHz)' )
ylabel( 'impedance ($\Omega$)' )
legend( )
idx = np.where((s11_dB<-10) & (s11_dB==np.min(s11_dB)))[0]
if not len(idx)==1:
print('No resonance frequency found for far-field calulation')
else:
f_res = f[idx[0]]
theta = np.arange(-180.0, 180.0, 2.0)
print("Calculate NF2FF")
nf2ff_res_phi0 = nf2ff.CalcNF2FF(Sim_Path, f_res, theta, 0, center=np.array([patch_radius+substrate_thickness, 0, 0])*unit, read_cached=True, outfile='nf2ff_xz.h5')
figure(figsize=(15, 7))
ax = subplot(121, polar=True)
E_norm = 20.0*np.log10(nf2ff_res_phi0.E_norm/np.max(nf2ff_res_phi0.E_norm)) + nf2ff_res_phi0.Dmax
ax.plot(np.deg2rad(theta), 10**(np.squeeze(E_norm)/20), linewidth=2, label='xz-plane')
ax.grid(True)
ax.set_xlabel('theta (deg)')
ax.set_theta_zero_location('N')
ax.set_theta_direction(-1)
ax.legend(loc=3)
phi = theta
nf2ff_res_theta90 = nf2ff.CalcNF2FF(Sim_Path, f_res, 90, phi, center=np.array([patch_radius+substrate_thickness, 0, 0])*unit, read_cached=True, outfile='nf2ff_xy.h5')
ax = subplot(122, polar=True)
E_norm = 20.0*np.log10(nf2ff_res_theta90.E_norm/np.max(nf2ff_res_theta90.E_norm)) + nf2ff_res_theta90.Dmax
ax.plot(np.deg2rad(phi), 10**(np.squeeze(E_norm)/20), linewidth=2, label='xy-plane')
ax.grid(True)
ax.set_xlabel('phi (deg)')
suptitle('Bent Patch Anteanna Pattern\nFrequency: {} GHz'.format(f_res/1e9), fontsize=14)
ax.legend(loc=3)
print( 'radiated power: Prad = {:.2e} Watt'.format(nf2ff_res_theta90.Prad[0]))
print( 'directivity: Dmax = {:.1f} ({:.1f} dBi)'.format(nf2ff_res_theta90.Dmax[0], 10*np.log10(nf2ff_res_theta90.Dmax[0])))
print( 'efficiency: nu_rad = {:.1f} %'.format(100*nf2ff_res_theta90.Prad[0]/real(P_in[idx[0]])))
show()

239
python/Tutorials/CRLH_Extraction.py Normal file
View File

@ -0,0 +1,239 @@
# -*- coding: utf-8 -*-
"""
Tutorials / CRLH_Extraction
Describtion at:
http://openems.de/index.php/Tutorial:_CRLH_Extraction
Tested with
- python 3.4
- openEMS v0.0.34+
(C) 2016 Thorsten Liebig <thorsten.liebig@gmx.de>
"""
### Import Libraries
import os, tempfile
from pylab import *
from CSXCAD import ContinuousStructure
from openEMS import openEMS
from openEMS.physical_constants import *
### Class to represent single CRLH unit cells
class CRLH_Cells:
def __init__(self, LL, LW, Top, Bot, GLT, GLB, SL, SW, VR):
self.LL = LL # Line length
self.LW = LW # Line width
self.Top = Top # top signal height
self.Bot = Bot # bottom signal height
self.GLT = GLT # gap length top
self.GLB = GLB # gap length bottom
self.SL = SL # stub length
self.SW = SW # stub width
self.VR = VR # via radius
self.props = dict() # property dictionary
self.edge_resolution = None
def createProperties(self, CSX):
for p in ['metal_top', 'metal_bot', 'via']:
self.props[p] = CSX.AddMetal(p)
def setEdgeResolution(self, res):
self.edge_resolution = res
def createCell(self, translate = [0,0,0]):
def append_mesh(mesh1, mesh2):
for n in range(3):
if mesh1[n] is None:
mesh1[n] = mesh2[n]
elif mesh2[n] is None:
continue
else:
mesh1[n] += mesh2[n]
return mesh1
translate = array(translate)
start = [-self.LL/2 , -self.LW/2, self.Top] + translate
stop = [-self.GLT/2, self.LW/2, self.Top] + translate
box = self.props['metal_top'].AddBox(start, stop, priority=10)
mesh = box.GetGridHint('x', metal_edge_res=self.edge_resolution, down_dir=False)
append_mesh(mesh, box.GetGridHint('y', metal_edge_res=self.edge_resolution) )
start = [+self.LL/2 , -self.LW/2, self.Top] + translate
stop = [+self.GLT/2, self.LW/2, self.Top] + translate
box = self.props['metal_top'].AddBox(start, stop, priority=10)
append_mesh(mesh, box.GetGridHint('x', metal_edge_res=self.edge_resolution, up_dir=False) )
start = [-(self.LL-self.GLB)/2, -self.LW/2, self.Bot] + translate
stop = [+(self.LL-self.GLB)/2, self.LW/2, self.Bot] + translate
box = self.props['metal_bot'].AddBox(start, stop, priority=10)
append_mesh(mesh, box.GetGridHint('x', metal_edge_res=self.edge_resolution) )
start = [-self.SW/2, -self.LW/2-self.SL, self.Bot] + translate
stop = [+self.SW/2, self.LW/2+self.SL, self.Bot] + translate
box = self.props['metal_bot'].AddBox(start, stop, priority=10)
append_mesh(mesh, box.GetGridHint('xy', metal_edge_res=self.edge_resolution) )
start = [0, -self.LW/2-self.SL+self.SW/2, 0 ] + translate
stop = [0, -self.LW/2-self.SL+self.SW/2, self.Bot] + translate
self.props['via'].AddCylinder(start, stop, radius=self.VR, priority=10)
start[1] *= -1
stop [1] *= -1
self.props['via'].AddCylinder(start, stop, radius=self.VR, priority=10)
return mesh
if __name__ == '__main__':
### Setup the simulation
Sim_Path = os.path.join(tempfile.gettempdir(), 'CRLH_Extraction')
post_proc_only = False
unit = 1e-6 # specify everything in um
feed_length = 30000
substrate_thickness = [1524, 101 , 254 ]
substrate_epsr = [3.48, 3.48, 3.48]
CRLH = CRLH_Cells(LL = 14e3, LW = 4e3, GLB = 1950, GLT = 4700, SL = 7800, SW = 1000, VR = 250 , \
Top = sum(substrate_thickness), \
Bot = sum(substrate_thickness[:-1]))
# frequency range of interest
f_start = 0.8e9
f_stop = 6e9
### Setup FDTD parameters & excitation function
CSX = ContinuousStructure()
FDTD = openEMS(EndCriteria=1e-5)
FDTD.SetCSX(CSX)
mesh = CSX.GetGrid()
mesh.SetDeltaUnit(unit)
CRLH.createProperties(CSX)
FDTD.SetGaussExcite((f_start+f_stop)/2, (f_stop-f_start)/2 )
BC = {'PML_8' 'PML_8' 'MUR' 'MUR' 'PEC' 'PML_8'}
FDTD.SetBoundaryCond( ['PML_8', 'PML_8', 'MUR', 'MUR', 'PEC', 'PML_8'] )
### Setup a basic mesh and create the CRLH unit cell
resolution = C0/(f_stop*sqrt(max(substrate_epsr)))/unit /30 # resolution of lambda/30
CRLH.setEdgeResolution(resolution/4)
mesh.SetLines('x', [-feed_length-CRLH.LL/2, 0, feed_length+CRLH.LL/2])
mesh.SetLines('y', [-30000, 0, 30000])
substratelines = cumsum(substrate_thickness)
mesh.SetLines('z', [0, 20000])
mesh.AddLine('z', cumsum(substrate_thickness))
mesh.AddLine('z', linspace(substratelines[-2],substratelines[-1],4))
# create the CRLH unit cell (will define additional fixed mesh lines)
mesh_hint = CRLH.createCell()
mesh.AddLine('x', mesh_hint[0])
mesh.AddLine('y', mesh_hint[1])
# Smooth the given mesh
mesh.SmoothMeshLines('all', resolution, 1.2)
### Setup the substrate layer
substratelines = [0] + substratelines.tolist()
start, stop = mesh.GetSimArea()
for n in range(len(substrate_thickness)):
sub = CSX.AddMaterial( 'substrate_{}'.format(n), epsilon=substrate_epsr[n] )
start[2] = substratelines[n]
stop [2] = substratelines[n+1]
sub.AddBox( start, stop )
### Add the feeding MSL ports
pec = CSX.AddMetal( 'PEC' )
port = [None, None]
x_lines = mesh.GetLines('x')
portstart = [ x_lines[0], -CRLH.LW/2, substratelines[-1]]
portstop = [ -CRLH.LL/2, CRLH.LW/2, 0]
port[0] = FDTD.AddMSLPort( 1, pec, portstart, portstop, 'x', 'z', excite=-1, FeedShift=10*resolution, MeasPlaneShift=feed_length/2, priority=10)
portstart = [ x_lines[-1], -CRLH.LW/2, substratelines[-1]]
portstop = [ +CRLH.LL/2 , CRLH.LW/2, 0]
port[1] = FDTD.AddMSLPort( 2, pec, portstart, portstop, 'x', 'z', MeasPlaneShift=feed_length/2, priority=10)
### Run the simulation
if 1: # debugging only
CSX_file = os.path.join(Sim_Path, 'CRLH_Extraction.xml')
if not os.path.exists(Sim_Path):
os.mkdir(Sim_Path)
CSX.Write2XML(CSX_file)
os.system(r'AppCSXCAD "{}"'.format(CSX_file))
if not post_proc_only:
FDTD.Run(Sim_Path, verbose=3, cleanup=True)
### Post-Processing
f = linspace( f_start, f_stop, 1601 )
for p in port:
p.CalcPort( Sim_Path, f, ref_impedance = 50, ref_plane_shift = feed_length)
# calculate and plot scattering parameter
s11 = port[0].uf_ref / port[0].uf_inc
s21 = port[1].uf_ref / port[0].uf_inc
plot(f/1e9,20*log10(abs(s11)),'k-' , linewidth=2, label='$S_{11}$')
plot(f/1e9,20*log10(abs(s21)),'r--', linewidth=2, label='$S_{21}$')
grid()
legend(loc=3)
ylabel('S-Parameter (dB)')
xlabel('frequency (GHz)')
ylim([-40, 2])
### Extract CRLH parameter form ABCD matrix
A = ((1+s11)*(1-s11) + s21*s21)/(2*s21)
C = ((1-s11)*(1-s11) - s21*s21)/(2*s21) / port[1].Z_ref
Y = C
Z = 2*(A-1)/C
iZ = imag(Z)
iY = imag(Y)
fse = interp(0, iZ, f)
fsh = interp(0, iY, f)
df = f[1]-f[0]
fse_idx = np.where(f>fse)[0][0]
fsh_idx = np.where(f>fsh)[0][0]
LR = 0.5*(iZ[fse_idx]-iZ[fse_idx-1])/(2*pi*df)
CL = 1/(2*pi*fse)**2/LR
CR = 0.5*(iY[fsh_idx]-iY[fsh_idx-1])/(2*pi*df)
LL = 1/(2*pi*fsh)**2/CR
print(' Series tank: CL = {:.2f} pF, LR = {:.2f} nH -> f_se = {:.2f} GHz '.format(CL*1e12, LR*1e9, fse*1e-9))
print(' Shunt tank: CR = {:.2f} pF, LL = {:.2f} nH -> f_sh = {:.2f} GHz '.format(CR*1e12, LL*1e9, fsh*1e-9))
### Calculate analytical wave-number of an inf-array of cells
w = 2*pi*f
wse = 2*pi*fse
wsh = 2*pi*fsh
beta_calc = real(arccos(1-(w**2-wse**2)*(w**2-wsh**2)/(2*w**2/CR/LR)))
# plot
figure()
beta = -angle(s21)/CRLH.LL/unit
plot(abs(beta)*CRLH.LL*unit/pi,f*1e-9,'k-', linewidth=2, label=r'$\beta_{CRLH,\ 1\ cell}$' )
grid()
plot(beta_calc/pi,f*1e-9,'c--', linewidth=2, label=r'$\beta_{CRLH,\ \infty\ cells}$')
plot(real(port[1].beta)*CRLH.LL*unit/pi,f*1e-9,'g-', linewidth=2, label=r'$\beta_{MSL}$')
ylim([1, 6])
xlabel(r'$|\beta| p / \pi$')
ylabel('frequency (GHz)')
legend(loc=2)
show()

191
python/Tutorials/Helical_Antenna.py Normal file
View File

@ -0,0 +1,191 @@
# -*- coding: utf-8 -*-
"""
Helical Antenna Tutorial
Tested with
- python 3.4
- openEMS v0.0.33+
(C) 2015-2016 Thorsten Liebig <thorsten.liebig@gmx.de>
"""
### Import Libraries
import os, tempfile
from pylab import *
from CSXCAD import CSXCAD
from openEMS import openEMS
from openEMS.physical_constants import *
### Setup the simulation
Sim_Path = os.path.join(tempfile.gettempdir(), 'Helical_Ant')
post_proc_only = False
unit = 1e-3 # all length in mm
f0 = 2.4e9 # center frequency, frequency of interest!
lambda0 = round(C0/f0/unit) # wavelength in mm
fc = 0.5e9 # 20 dB corner frequency
Helix_radius = 20 # --> diameter is ~ lambda/pi
Helix_turns = 10 # --> expected gain is G ~ 4 * 10 = 40 (16dBi)
Helix_pitch = 30 # --> pitch is ~ lambda/4
Helix_mesh_res = 3
gnd_radius = lambda0/2
# feeding
feed_heigth = 3
feed_R = 120 #feed impedance
# size of the simulation box
SimBox = array([1, 1, 1.5])*2.0*lambda0
### Setup FDTD parameter & excitation function
FDTD = openEMS(EndCriteria=1e-4)
FDTD.SetGaussExcite( f0, fc )
FDTD.SetBoundaryCond( ['MUR', 'MUR', 'MUR', 'MUR', 'MUR', 'PML_8'] )
### Setup Geometry & Mesh
CSX = CSXCAD.ContinuousStructure()
FDTD.SetCSX(CSX)
mesh = CSX.GetGrid()
mesh.SetDeltaUnit(unit)
max_res = floor(C0 / (f0+fc) / unit / 20) # cell size: lambda/20
# create helix mesh
mesh.AddLine('x', [-Helix_radius, 0, Helix_radius])
mesh.SmoothMeshLines('x', Helix_mesh_res)
# add the air-box
mesh.AddLine('x', [-SimBox[0]/2-gnd_radius, SimBox[0]/2+gnd_radius])
# create a smooth mesh between specified fixed mesh lines
mesh.SmoothMeshLines('x', max_res, ratio=1.4)
# copy x-mesh to y-direction
mesh.SetLines('y', mesh.GetLines('x'))
# create helix mesh in z-direction
mesh.AddLine('z', [0, feed_heigth, Helix_turns*Helix_pitch+feed_heigth])
mesh.SmoothMeshLines('z', Helix_mesh_res)
# add the air-box
mesh.AddLine('z', [-SimBox[2]/2, max(mesh.GetLines('z'))+SimBox[2]/2 ])
# create a smooth mesh between specified fixed mesh lines
mesh.SmoothMeshLines('z', max_res, ratio=1.4)
### Create the Geometry
## * Create the metal helix using the wire primitive.
## * Create a metal gorund plane as cylinder.
# create a perfect electric conductor (PEC)
helix_metal = CSX.AddMetal('helix' )
ang = linspace(0,2*pi,21)
coil_x = Helix_radius*cos(ang)
coil_y = Helix_radius*sin(ang)
coil_z = ang/2/pi*Helix_pitch
Helix_x=np.array([])
Helix_y=np.array([])
Helix_z=np.array([])
zpos = feed_heigth
for n in range(Helix_turns-1):
Helix_x = r_[Helix_x, coil_x]
Helix_y = r_[Helix_y, coil_y]
Helix_z = r_[Helix_z ,coil_z+zpos]
zpos = zpos + Helix_pitch
p = np.array([Helix_x, Helix_y, Helix_z])
helix_metal.AddCurve(p)
# create ground circular ground
gnd = CSX.AddMetal( 'gnd' ) # create a perfect electric conductor (PEC)
# add a box using cylindrical coordinates
start = [0, 0, -0.1]
stop = [0, 0, 0.1]
gnd.AddCylinder(start, stop, radius=gnd_radius)
# apply the excitation & resist as a current source
start = [Helix_radius, 0, 0]
stop = [Helix_radius, 0, feed_heigth]
port = FDTD.AddLumpedPort(1 ,feed_R, start, stop, 'z', 1.0, priority=5)
# nf2ff calc
nf2ff = FDTD.CreateNF2FFBox(opt_resolution=[lambda0/15]*3)
### Run the simulation
if 0: # debugging only
CSX_file = os.path.join(Sim_Path, 'helix.xml')
if not os.path.exists(Sim_Path):
os.mkdir(Sim_Path)
CSX.Write2XML(CSX_file)
os.system(r'AppCSXCAD "{}"'.format(CSX_file))
if not post_proc_only:
FDTD.Run(Sim_Path, verbose=3, cleanup=True)
### Postprocessing & plotting
freq = linspace( f0-fc, f0+fc, 501 )
port.CalcPort(Sim_Path, freq)
Zin = port.uf_tot / port.if_tot
s11 = port.uf_ref / port.uf_inc
## Plot the feed point impedance
figure()
plot( freq/1e6, real(Zin), 'k-', linewidth=2, label=r'$\Re(Z_{in})$' )
grid()
plot( freq/1e6, imag(Zin), 'r--', linewidth=2, label=r'$\Im(Z_{in})$' )
title( 'feed point impedance' )
xlabel( 'frequency (MHz)' )
ylabel( 'impedance ($\Omega$)' )
legend( )
## Plot reflection coefficient S11
figure()
plot( freq/1e6, 20*log10(abs(s11)), 'k-', linewidth=2 )
grid()
title( 'reflection coefficient $S_{11}$' )
xlabel( 'frequency (MHz)' )
ylabel( 'reflection coefficient $|S_{11}|$' )
### Create the NFFF contour
## * calculate the far field at phi=0 degrees and at phi=90 degrees
theta = arange(0.,180.,1.)
phi = arange(-180,180,2)
disp( 'calculating the 3D far field...' )
nf2ff_res = nf2ff.CalcNF2FF(Sim_Path, f0, theta, phi, read_cached=True, verbose=True )
Dmax_dB = 10*log10(nf2ff_res.Dmax[0])
E_norm = 20.0*log10(nf2ff_res.E_norm[0]/np.max(nf2ff_res.E_norm[0])) + 10*log10(nf2ff_res.Dmax[0])
theta_HPBW = theta[ np.where(squeeze(E_norm[:,phi==0])<Dmax_dB-3)[0][0] ]
## * Display power and directivity
print('radiated power: Prad = {} W'.format(nf2ff_res.Prad[0]))
print('directivity: Dmax = {} dBi'.format(Dmax_dB))
print('efficiency: nu_rad = {} %'.format(100*nf2ff_res.Prad[0]/interp(f0, freq, port.P_acc)))
print('theta_HPBW = {} °'.format(theta_HPBW))
E_norm = 20.0*log10(nf2ff_res.E_norm[0]/np.max(nf2ff_res.E_norm[0])) + 10*log10(nf2ff_res.Dmax[0])
E_CPRH = 20.0*log10(np.abs(nf2ff_res.E_cprh[0])/np.max(nf2ff_res.E_norm[0])) + 10*log10(nf2ff_res.Dmax[0])
E_CPLH = 20.0*log10(np.abs(nf2ff_res.E_cplh[0])/np.max(nf2ff_res.E_norm[0])) + 10*log10(nf2ff_res.Dmax[0])
## * Plot the pattern
figure()
plot(theta, E_norm[:,phi==0],'k-' , linewidth=2, label='$|E|$')
plot(theta, E_CPRH[:,phi==0],'g--', linewidth=2, label='$|E_{CPRH}|$')
plot(theta, E_CPLH[:,phi==0],'r-.', linewidth=2, label='$|E_{CPLH}|$')
grid()
xlabel('theta (deg)')
ylabel('directivity (dBi)')
title('Frequency: {} GHz'.format(nf2ff_res.freq[0]/1e9))
legend()
show()

123
python/Tutorials/MSL_NotchFilter.py Normal file
View File

@ -0,0 +1,123 @@
# -*- coding: utf-8 -*-
"""
Microstrip Notch Filter Tutorial
Describtion at:
http://openems.de/doc/openEMS/Tutorials.html#microstrip-notch-filter
Tested with
- python 3.4
- openEMS v0.0.34+
(C) 2016 Thorsten Liebig <thorsten.liebig@gmx.de>
"""
### Import Libraries
import os, tempfile
from pylab import *
from CSXCAD import ContinuousStructure
from openEMS import openEMS
from openEMS.physical_constants import *
### Setup the simulation
Sim_Path = os.path.join(tempfile.gettempdir(), 'NotchFilter')
post_proc_only = False
unit = 1e-6 # specify everything in um
MSL_length = 50000
MSL_width = 600
substrate_thickness = 254
substrate_epr = 3.66
stub_length = 12e3
f_max = 7e9
### Setup FDTD parameters & excitation function
FDTD = openEMS()
FDTD.SetGaussExcite( f_max/2, f_max/2 )
FDTD.SetBoundaryCond( ['PML_8', 'PML_8', 'MUR', 'MUR', 'PEC', 'MUR'] )
### Setup Geometry & Mesh
CSX = ContinuousStructure()
FDTD.SetCSX(CSX)
mesh = CSX.GetGrid()
mesh.SetDeltaUnit(unit)
resolution = C0/(f_max*sqrt(substrate_epr))/unit/50 # resolution of lambda/50
third_mesh = array([2*resolution/3, -resolution/3])/4
## Do manual meshing
mesh.AddLine('x', 0)
mesh.AddLine('x', MSL_width/2+third_mesh)
mesh.AddLine('x', -MSL_width/2-third_mesh)
mesh.SmoothMeshLines('x', resolution/4)
mesh.AddLine('x', [-MSL_length, MSL_length])
mesh.SmoothMeshLines('x', resolution)
mesh.AddLine('y', 0)
mesh.AddLine('y', MSL_width/2+third_mesh)
mesh.AddLine('y', -MSL_width/2-third_mesh)
mesh.SmoothMeshLines('y', resolution/4)
mesh.AddLine('y', [-15*MSL_width, 15*MSL_width+stub_length])
mesh.AddLine('y', (MSL_width/2+stub_length)+third_mesh)
mesh.SmoothMeshLines('y', resolution)
mesh.AddLine('z', linspace(0,substrate_thickness,5))
mesh.AddLine('z', 3000)
mesh.SmoothMeshLines('z', resolution)
## Add the substrate
substrate = CSX.AddMaterial( 'RO4350B', epsilon=substrate_epr)
start = [-MSL_length, -15*MSL_width, 0]
stop = [+MSL_length, +15*MSL_width+stub_length, substrate_thickness]
substrate.AddBox(start, stop )
## MSL port setup
port = [None, None]
pec = CSX.AddMetal( 'PEC' )
portstart = [ -MSL_length, -MSL_width/2, substrate_thickness]
portstop = [ 0, MSL_width/2, 0]
port[0] = FDTD.AddMSLPort( 1, pec, portstart, portstop, 'x', 'z', excite=-1, FeedShift=10*resolution, MeasPlaneShift=MSL_length/3, priority=10)
portstart = [MSL_length, -MSL_width/2, substrate_thickness]
portstop = [0 , MSL_width/2, 0]
port[1] = FDTD.AddMSLPort( 2, pec, portstart, portstop, 'x', 'z', MeasPlaneShift=MSL_length/3, priority=10 )
## Filter-Stub Definition
start = [-MSL_width/2, MSL_width/2, substrate_thickness]
stop = [ MSL_width/2, MSL_width/2+stub_length, substrate_thickness]
pec.AddBox(start, stop, priority=10 )
### Run the simulation
if 0: # debugging only
CSX_file = os.path.join(Sim_Path, 'notch.xml')
if not os.path.exists(Sim_Path):
os.mkdir(Sim_Path)
CSX.Write2XML(CSX_file)
os.system(r'AppCSXCAD "{}"'.format(CSX_file))
if not post_proc_only:
FDTD.Run(Sim_Path, verbose=3, cleanup=True)
### Post-processing and plotting
f = linspace( 1e6, f_max, 1601 )
for p in port:
p.CalcPort( Sim_Path, f, ref_impedance = 50)
s11 = port[0].uf_ref / port[0].uf_inc
s21 = port[1].uf_ref / port[0].uf_inc
plot(f/1e9,20*log10(abs(s11)),'k-',linewidth=2 , label='$S_{11}$')
grid()
plot(f/1e9,20*log10(abs(s21)),'r--',linewidth=2 , label='$S_{21}$')
legend()
ylabel('S-Parameter (dB)')
xlabel('frequency (GHz)')
ylim([-40, 2])
show()

126
python/Tutorials/RCS_Sphere.py Normal file
View File

@ -0,0 +1,126 @@
# -*- coding: utf-8 -*-
"""
Tutorials / radar cross section of a metal sphere
Tested with
- python 3.4
- openEMS v0.0.34+
(C) 2016 Thorsten Liebig <thorsten.liebig@gmx.de>
"""
### Import Libraries
import os, tempfile
from pylab import *
from CSXCAD import ContinuousStructure
from openEMS import openEMS
from openEMS.physical_constants import *
from openEMS.ports import UI_data
### Setup the simulation
Sim_Path = os.path.join(tempfile.gettempdir(), 'RCS_Sphere')
post_proc_only = False
unit = 1e-3 # all length in mm
sphere_rad = 200
inc_angle = 0 #incident angle (to x-axis) in deg
# size of the simulation box
SimBox = 1200
PW_Box = 750
### Setup FDTD parameters & excitation function
FDTD = openEMS(EndCriteria=1e-5)
f_start = 50e6 # start frequency
f_stop = 1000e6 # stop frequency
f0 = 500e6
FDTD.SetGaussExcite( 0.5*(f_start+f_stop), 0.5*(f_stop-f_start) )
FDTD.SetBoundaryCond( ['PML_8', 'PML_8', 'PML_8', 'PML_8', 'PML_8', 'PML_8'] )
### Setup Geometry & Mesh
CSX = ContinuousStructure()
FDTD.SetCSX(CSX)
mesh = CSX.GetGrid()
mesh.SetDeltaUnit(unit)
#create mesh
mesh.SetLines('x', [-SimBox/2, 0, SimBox/2])
mesh.SmoothMeshLines('x', C0 / f_stop / unit / 20) # cell size: lambda/20
mesh.SetLines('y', mesh.GetLines('x'))
mesh.SetLines('z', mesh.GetLines('x'))
### Create a metal sphere and plane wave source
sphere_metal = CSX.AddMetal( 'sphere' ) # create a perfect electric conductor (PEC)
sphere_metal.AddSphere(priority=10, center=[0, 0, 0], radius=sphere_rad)
# plane wave excitation
k_dir = [cos(inc_angle), sin(inc_angle), 0] # plane wave direction
E_dir = [0, 0, 1] # plane wave polarization --> E_z
pw_exc = CSX.AddExcitation('plane_wave', exc_type=10, exc_val=E_dir)
pw_exc.SetPropagationDir(k_dir)
pw_exc.SetFrequency(f0)
start = np.array([-PW_Box/2, -PW_Box/2, -PW_Box/2])
stop = -start
pw_exc.AddBox(start, stop)
# nf2ff calc
nf2ff = FDTD.CreateNF2FFBox()
### Run the simulation
if 0: # debugging only
CSX_file = os.path.join(Sim_Path, 'RCS_Sphere.xml')
if not os.path.exists(Sim_Path):
os.mkdir(Sim_Path)
CSX.Write2XML(CSX_file)
os.system(r'AppCSXCAD "{}"'.format(CSX_file))
if not post_proc_only:
FDTD.Run(Sim_Path, verbose=3, cleanup=True)
### Postprocessing & plotting
# get Gaussian pulse stength at frequency f0
ef = UI_data('et', Sim_Path, freq=f0)
Pin = 0.5*norm(E_dir)**2/Z0 * abs(ef.ui_f_val[0])**2
#
nf2ff_res = nf2ff.CalcNF2FF(Sim_Path, f0, 90, arange(-180, 180.1, 2))
RCS = 4*pi/Pin[0]*nf2ff_res.P_rad[0]
fig = figure()
ax = fig.add_subplot(111, polar=True)
ax.plot( nf2ff_res.phi, RCS[0], 'k-', linewidth=2 )
ax.grid(True)
# calculate RCS over frequency
freq = linspace(f_start,f_stop,100)
ef = UI_data( 'et', Sim_Path, freq ) # time domain/freq domain voltage
Pin = 0.5*norm(E_dir)**2/Z0 * abs(np.array(ef.ui_f_val[0]))**2
nf2ff_res = nf2ff.CalcNF2FF(Sim_Path, freq, 90, 180+inc_angle, outfile='back_nf2ff.h5')
back_scat = np.array([4*pi/Pin[fn]*nf2ff_res.P_rad[fn][0][0] for fn in range(len(freq))])
figure()
plot(freq/1e6,back_scat, linewidth=2)
grid()
xlabel('frequency (MHz)')
ylabel('RCS ($m^2$)')
title('radar cross section')
figure()
semilogy(sphere_rad*unit/C0*freq,back_scat/(pi*sphere_rad*unit*sphere_rad*unit), linewidth=2)
ylim([10^-2, 10^1])
grid()
xlabel('sphere radius / wavelength')
ylabel('RCS / ($\pi a^2$)')
title('normalized radar cross section')
show()

125
python/Tutorials/Rect_Waveguide.py Normal file
View File

@ -0,0 +1,125 @@
# -*- coding: utf-8 -*-
"""
Rectangular Waveguide Tutorial
Describtion at:
http://openems.de/doc/openEMS/Tutorials.html#rectangular-waveguide
Tested with
- python 3.4
- openEMS v0.0.34+
(C) 2015-2016 Thorsten Liebig <thorsten.liebig@gmx.de>
"""
### Import Libraries
import os, tempfile
from pylab import *
from CSXCAD import ContinuousStructure
from openEMS import openEMS
from openEMS.physical_constants import *
### Setup the simulation
Sim_Path = os.path.join(tempfile.gettempdir(), 'Rect_WG')
post_proc_only = False
unit = 1e-6; #drawing unit in um
# waveguide dimensions
# WR42
a = 10700; #waveguide width
b = 4300; #waveguide heigth
length = 50000;
# frequency range of interest
f_start = 20e9;
f_0 = 24e9;
f_stop = 26e9;
lambda0 = C0/f_0/unit;
#waveguide TE-mode definition
TE_mode = 'TE10';
#targeted mesh resolution
mesh_res = lambda0/30
### Setup FDTD parameter & excitation function
FDTD = openEMS(NrTS=1e4);
FDTD.SetGaussExcite(0.5*(f_start+f_stop),0.5*(f_stop-f_start));
# boundary conditions
FDTD.SetBoundaryCond([0, 0, 0, 0, 3, 3]);
### Setup geometry & mesh
CSX = ContinuousStructure()
FDTD.SetCSX(CSX)
mesh = CSX.GetGrid()
mesh.SetDeltaUnit(unit)
mesh.AddLine('x', [0, a])
mesh.AddLine('y', [0, b])
mesh.AddLine('z', [0, length])
## Apply the waveguide port
ports = []
start=[0, 0, 10*mesh_res];
stop =[a, b, 15*mesh_res];
mesh.AddLine('z', [start[2], stop[2]])
ports.append(FDTD.AddRectWaveGuidePort( 0, start, stop, 'z', a*unit, b*unit, TE_mode, 1))
start=[0, 0, length-10*mesh_res];
stop =[a, b, length-15*mesh_res];
mesh.AddLine('z', [start[2], stop[2]])
ports.append(FDTD.AddRectWaveGuidePort( 1, start, stop, 'z', a*unit, b*unit, TE_mode))
mesh.SmoothMeshLines('all', mesh_res, ratio=1.4)
### Define dump box...
Et = CSX.AddDump('Et', file_type=0, sub_sampling=[2,2,2])
start = [0, 0, 0];
stop = [a, b, length];
Et.AddBox(start, stop);
### Run the simulation
if 0: # debugging only
CSX_file = os.path.join(Sim_Path, 'rect_wg.xml')
if not os.path.exists(Sim_Path):
os.mkdir(Sim_Path)
CSX.Write2XML(CSX_file)
os.system(r'AppCSXCAD "{}"'.format(CSX_file))
if not post_proc_only:
FDTD.Run(Sim_Path, verbose=3, cleanup=True)
### Postprocessing & plotting
freq = linspace(f_start,f_stop,201)
for port in ports:
port.CalcPort(Sim_Path, freq)
s11 = ports[0].uf_ref / ports[0].uf_inc
s21 = ports[1].uf_ref / ports[0].uf_inc
ZL = ports[0].uf_tot / ports[0].if_tot
ZL_a = ports[0].ZL # analytic waveguide impedance
## Plot s-parameter
figure()
plot(freq*1e-6,20*log10(abs(s11)),'k-',linewidth=2, label='$S_{11}$')
grid()
plot(freq*1e-6,20*log10(abs(s21)),'r--',linewidth=2, label='$S_{21}$')
legend();
ylabel('S-Parameter (dB)')
xlabel(r'frequency (MHz) $\rightarrow$')
## Compare analytic and numerical wave-impedance
figure()
plot(freq*1e-6,real(ZL), linewidth=2, label='$\Re\{Z_L\}$')
grid()
plot(freq*1e-6,imag(ZL),'r--', linewidth=2, label='$\Im\{Z_L\}$')
plot(freq*1e-6,ZL_a,'g-.',linewidth=2, label='$Z_{L, analytic}$')
ylabel('ZL $(\Omega)$')
xlabel(r'frequency (MHz) $\rightarrow$')
legend()
show()

151
python/Tutorials/Simple_Patch_Antenna.py Normal file
View File

@ -0,0 +1,151 @@
# -*- coding: utf-8 -*-
"""
Created on Fri Dec 18 20:56:53 2015
@author: thorsten
"""
### Import Libraries
import os, tempfile
from pylab import *
from CSXCAD import ContinuousStructure
from openEMS import openEMS
from openEMS.physical_constants import *
### General parameter setup
Sim_Path = os.path.join(tempfile.gettempdir(), 'Simp_Patch')
post_proc_only = False
# patch width (resonant length) in x-direction
patch_width = 32 #
# patch length in y-direction
patch_length = 40
#substrate setup
substrate_epsR = 3.38
substrate_kappa = 1e-3 * 2*pi*2.45e9 * EPS0*substrate_epsR
substrate_width = 60
substrate_length = 60
substrate_thickness = 1.524
substrate_cells = 4
#setup feeding
feed_pos = -6 #feeding position in x-direction
feed_R = 50 #feed resistance
# size of the simulation box
SimBox = np.array([200, 200, 150])
# setup FDTD parameter & excitation function
f0 = 2e9 # center frequency
fc = 1e9 # 20 dB corner frequency
### FDTD setup
## * Limit the simulation to 30k timesteps
## * Define a reduced end criteria of -40dB
FDTD = openEMS(NrTS=30000, EndCriteria=1e-4)
FDTD.SetGaussExcite( f0, fc )
FDTD.SetBoundaryCond( ['MUR', 'MUR', 'MUR', 'MUR', 'MUR', 'MUR'] )
CSX = ContinuousStructure()
FDTD.SetCSX(CSX)
mesh = CSX.GetGrid()
mesh.SetDeltaUnit(1e-3)
mesh_res = C0/(f0+fc)/1e-3/20
### Generate properties, primitives and mesh-grid
#initialize the mesh with the "air-box" dimensions
mesh.AddLine('x', [-SimBox[0]/2, SimBox[0]/2])
mesh.AddLine('y', [-SimBox[1]/2, SimBox[1]/2] )
mesh.AddLine('z', [-SimBox[2]/3, SimBox[2]*2/3] )
# create patch
patch = CSX.AddMetal( 'patch' ) # create a perfect electric conductor (PEC)
start = [-patch_width/2, -patch_length/2, substrate_thickness]
stop = [ patch_width/2 , patch_length/2, substrate_thickness]
patch.AddBox(priority=10, start=start, stop=stop) # add a box-primitive to the metal property 'patch'
FDTD.AddEdges2Grid(dirs='xy', properties=patch, metal_edge_res=mesh_res/2)
# create substrate
substrate = CSX.AddMaterial( 'substrate', epsilon=substrate_epsR, kappa=substrate_kappa)
start = [-substrate_width/2, -substrate_length/2, 0]
stop = [ substrate_width/2, substrate_length/2, substrate_thickness]
substrate.AddBox( priority=0, start=start, stop=stop )
# add extra cells to discretize the substrate thickness
mesh.AddLine('z', linspace(0,substrate_thickness,substrate_cells+1))
# create ground (same size as substrate)
gnd = CSX.AddMetal( 'gnd' ) # create a perfect electric conductor (PEC)
start[2]=0
stop[2] =0
gnd.AddBox(start, stop, priority=10)
FDTD.AddEdges2Grid(dirs='xy', properties=gnd)
# apply the excitation & resist as a current source
start = [feed_pos, 0, 0]
stop = [feed_pos, 0, substrate_thickness]
port = FDTD.AddLumpedPort(1, feed_R, start, stop, 'z', 1.0, priority=5, edges2grid='xy')
mesh.SmoothMeshLines('all', mesh_res, 1.4)
# Add the nf2ff recording box
nf2ff = FDTD.CreateNF2FFBox()
### Run the simulation
if 0: # debugging only
CSX_file = os.path.join(Sim_Path, 'simp_patch.xml')
if not os.path.exists(Sim_Path):
os.mkdir(Sim_Path)
CSX.Write2XML(CSX_file)
os.system(r'AppCSXCAD "{}"'.format(CSX_file))
if not post_proc_only:
FDTD.Run(Sim_Path, verbose=3, cleanup=True)
### Post-processing and plotting
f = np.linspace(max(1e9,f0-fc),f0+fc,401)
port.CalcPort(Sim_Path, f)
s11 = port.uf_ref/port.uf_inc
s11_dB = 20.0*np.log10(np.abs(s11))
figure()
plot(f/1e9, s11_dB, 'k-', linewidth=2, label='$S_{11}$')
grid()
legend()
ylabel('S-Parameter (dB)')
xlabel('Frequency (GHz)')
idx = np.where((s11_dB<-10) & (s11_dB==np.min(s11_dB)))[0]
if not len(idx)==1:
print('No resonance frequency found for far-field calulation')
else:
f_res = f[idx[0]]
theta = np.arange(-180.0, 180.0, 2.0)
phi = [0., 90.]
nf2ff_res = nf2ff.CalcNF2FF(Sim_Path, f_res, theta, phi, center=[0,0,1e-3])
figure()
E_norm = 20.0*np.log10(nf2ff_res.E_norm[0]/np.max(nf2ff_res.E_norm[0])) + nf2ff_res.Dmax[0]
plot(theta, np.squeeze(E_norm[:,0]), 'k-', linewidth=2, label='xz-plane')
plot(theta, np.squeeze(E_norm[:,1]), 'r--', linewidth=2, label='yz-plane')
grid()
ylabel('Directivity (dBi)')
xlabel('Theta (deg)')
title('Frequency: {} GHz'.format(f_res/1e9))
legend()
Zin = port.uf_tot/port.if_tot
figure()
plot(f/1e9, np.real(Zin), 'k-', linewidth=2, label='$\Re\{Z_{in}\}$')
plot(f/1e9, np.imag(Zin), 'r--', linewidth=2, label='$\Im\{Z_{in}\}$')
grid()
legend()
ylabel('Zin (Ohm)')
xlabel('Frequency (GHz)')
show()

9
python/doc/Tutorials/Antenna_Tutorials.rst Normal file
View File

@ -0,0 +1,9 @@
Antennas
--------
.. toctree::
:maxdepth: 1
Simple_Patch_Antenna
Helical_Antenna
Bent_Patch_Antenna

35
python/doc/Tutorials/Bent_Patch_Antenna.rst Normal file
View File

@ -0,0 +1,35 @@
Bent Patch Antenna
==================
* Setup & Simulate a bent patch antenna using a cylindrical mesh
Introduction
-------------
**This tutorial covers:**
* Setup of a Bent Patch Antenna (see for comparison: :ref:`simple_patch_antenna`)
* setup of a *cylindrical FDTD mesh*.
* Calculate the S-Parameter and input impedance
* Calculate far-field pattern 2D/3D
Python Script
-------------
Get the latest version `from git <http://www.openems.de/gitweb/?p=openEMS.git;a=blob_plain;f=matlab/Tutorials/Bent_Patch_Antenna.m;hb=refs/heads/master>`_.
.. include:: ./__Bent_Patch_Antenna.txt
Images
-------------
.. figure:: images/Bent_Patch.png
:width: 49%
:alt: alternate text
3D view of the Bent Patch Antenna (AppCSXCAD)
.. figure:: images/Bent_Patch_Pattern.png
:width: 80%
:alt: Farfield pattern
Farfield pattern on an xy- and xz-plane

40
python/doc/Tutorials/CRLH_Extraction.rst Normal file
View File

@ -0,0 +1,40 @@
CRLH Parameter Extraction
=========================
* Setup a composite-right/left-handed (CRLH) unit cell and extract the equivalent circuit parameter.
Introduction
-------------
**This tutorial covers:**
* Setup a feeding mircostrip line & port
* Apply an inhomogeneous mesh used for improved accuracy and simulation speed
* Use an internal clss to setup a CRLH unit cell
* Use the port voltages and currents to extract the unit cell equivalent circuit parameter
.. figure:: images/CRLH_cell.png
:width: 80%
:alt: CRLH unit cell with feeding MSL.
CRLH unit cell with feeding MSL.
Python Script
-------------
Get the latest version `from git <http://www.openems.de/gitweb/?p=openEMS.git;a=blob_plain;f=matlab/Tutorials/Bent_Patch_Antenna.m;hb=refs/heads/master>`_.
.. include:: ./__CRLH_Extraction.txt
Images
-------------
.. figure:: images/CRLH_Spara.png
:width: 80%
:alt: CRLH cell S-parameter
CRLH cell S-parameter
.. figure:: images/CRLH_dispersion.png
:width: 80%
:alt: CRLH unit cell dispersion diagram
CRLH unit cell dispersion diagram

32
python/doc/Tutorials/Helical_Antenna.rst Normal file
View File

@ -0,0 +1,32 @@
Helical Antenna
===============
Introduction
-------------
**This tutorial covers:**
* setup of a helix using the wire primitive
* setup a lumped feeding port (R_in = 120 Ohms)
* adding a near-field to far-field (nf2ff) box using an efficient subsampling
* calculate the S-Parameter of the antenna
* calculate and plot the far-field pattern
Python Script
-------------
Get the latest version `from git <http://www.openems.de/gitweb/?p=openEMS.git;a=blob_plain;f=matlab/Tutorials/Helical_Antenna.m;hb=refs/heads/master>`_.
.. include:: ./__Helical_Antenna.txt
Images
-------------
.. figure:: images/Helix_Ant.png
:width: 49%
:alt: alternate text
3D view of the Helical Antenna (AppCSXCAD)
.. figure:: images/Helix_Ant_Pattern.png
:width: 49%
:alt: alternate text
Far-Field pattern showing a right-handed circular polarization.

10
python/doc/Tutorials/Intro_Tutorials.rst Normal file
View File

@ -0,0 +1,10 @@
.. _intro_tutorials:
Introductional Tutorials
------------------------
.. toctree::
Rect_Waveguide
RCS_Sphere

27
python/doc/Tutorials/MSL_NotchFilter.rst Normal file
View File

@ -0,0 +1,27 @@
Microstrip Notch Filter
=======================
* A straight MSL line with a open-ended stub to create a simple microwave filter.
Introduction
-------------
**This tutorial covers:**
* Setup a mircostrip line (MSL) and MSL port
* Apply an inhomogeneous mesh used for improved accuracy and simulation speed
* Calculate the S-Parameter of the filter
Python Script
-------------
Get the latest version `from git <http://openems.de/gitweb/?p=openEMS.git;a=blob_plain;f=matlab/Tutorials/MSL_NotchFilter.m;hb=HEAD>`_.
.. include:: ./__MSL_NotchFilter.txt
Images
-------------
.. figure:: images/Notch_Filter_SPara.png
:width: 49%
:alt: S-Parameter over Frequency
S-Parameter over Frequency

10
python/doc/Tutorials/MicroWave_Tutorials.rst Normal file
View File

@ -0,0 +1,10 @@
.. _microwave_tutorials:
Micro Wave Tutorials
--------------------
.. toctree::
MSL_NotchFilter
CRLH_Extraction

32
python/doc/Tutorials/RCS_Sphere.rst Normal file
View File

@ -0,0 +1,32 @@
Metal Sphere Radar Cross Section
================================
* A 3D simulation demonstrating a the total-field/scattered-field approach on a metallic sphere with a RCS (radar cross section) calculation.
Introduction
-------------
**This tutorial covers:**
* The total-field/scattered-field approach
* Calculation of a radar cross section (RCS)
Python Script
-------------
Get the latest version `from git <https://raw.githubusercontent.com/thliebig/openEMS/master/python/Tutorials/RCS_Sphere.py>`_.
.. include:: ./__RCS_Sphere.txt
Images
-------------
.. figure:: images/RCS_pattern.png
:width: 49%
:alt: Radar cross section pattern
Radar cross section pattern
.. figure:: images/RCS_norm.png
:width: 49%
:alt: normalized radar cross section
Normalized radar cross Section over normalized wavelength

27
python/doc/Tutorials/Rect_Waveguide.rst Normal file
View File

@ -0,0 +1,27 @@
Rectangular Waveguide
=====================
* A simple rectangular waveguide, showing the openEMS mode profile capabilities.
Introduction
-------------
**This tutorial covers:**
* Setup a mode profile excitation
* Create voltage and current probes using the mode profile
* Calculate the waveguide impedance and S-Parameter
Python Script
-------------
Get the latest version `from git <http://openems.de/gitweb/?p=openEMS.git;a=blob_plain;f=matlab/Tutorials/Rect_Waveguide.m;hb=HEAD>`_.
.. include:: ./__Rect_Waveguide.txt
Images
-------------
.. figure:: images/Rect_WG_SPara.png
:width: 49%
:alt: S-Parameter over Frequency
S-Parameter over Frequency

42
python/doc/Tutorials/Simple_Patch_Antenna.rst Normal file
View File

@ -0,0 +1,42 @@
.. _simple_patch_antenna:
Simple Patch Antenna
====================
Introduction
------------
A simple patch antenna for 2.4 GHz.
**This tutorial covers:**
* Setup a patch, substrate and ground.
* Setup of a lumped feeding port.
* Adding a near-field to far-field (nf2ff) recording box.
* Calculate the S-Parameter of the antenna.
* Calculate and plot the far-field pattern
Python Script
-------------
Get the latest version `from git <http://www.openems.de/gitweb/?p=openEMS.git;a=blob_plain;f=matlab/Tutorials/Simple_Patch_Antenna.m;hb=refs/heads/master>`_.
.. include:: ./__Simple_Patch_Antenna.txt
Images
------
.. figure:: images/Simp_Patch_S11.png
:width: 49%
:alt: S11 over Frequency
S-Parameter over Frequency
.. figure:: images/Simp_Patch_Zin.png
:width: 49%
:alt: Input Impedance
Antenna Input Impedance
.. figure:: images/Simp_Patch_Pattern.png
:width: 49%
:alt: Farfield pattern
Farfield pattern for the xy- and yz-plane.

BIN
python/doc/Tutorials/images/Bent_Patch.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 346 KiB

BIN
python/doc/Tutorials/images/Bent_Patch_Pattern.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 125 KiB

BIN
python/doc/Tutorials/images/Bent_Patch_SPara.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 44 KiB

BIN
python/doc/Tutorials/images/CRLH_Spara.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 42 KiB

BIN
python/doc/Tutorials/images/CRLH_cell.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 157 KiB

BIN
python/doc/Tutorials/images/CRLH_dispersion.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 57 KiB

BIN
python/doc/Tutorials/images/Helix_Ant.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 304 KiB

BIN
python/doc/Tutorials/images/Helix_Ant_Pattern.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 62 KiB

BIN
python/doc/Tutorials/images/Notch_Filter_SPara.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 48 KiB

BIN
python/doc/Tutorials/images/RCS_norm.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 42 KiB

BIN
python/doc/Tutorials/images/RCS_pattern.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 60 KiB

BIN
python/doc/Tutorials/images/Rect_WG_SPara.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 28 KiB

BIN
python/doc/Tutorials/images/Simp_Patch_Pattern.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 53 KiB

BIN
python/doc/Tutorials/images/Simp_Patch_S11.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 28 KiB

BIN
python/doc/Tutorials/images/Simp_Patch_Zin.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 39 KiB

12
python/doc/Tutorials/index.rst Normal file
View File

@ -0,0 +1,12 @@
.. _tutorials:
#########
Tutorials
#########
.. toctree::
:maxdepth: 2
Intro_Tutorials
MicroWave_Tutorials
Antenna_Tutorials

297
python/doc/conf.py Normal file
View File

@ -0,0 +1,297 @@
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#
# openEMS documentation build configuration file, created by
# sphinx-quickstart on Thu Sep 8 20:42:18 2016.
#
# This file is execfile()d with the current directory set to its
# containing dir.
#
# Note that not all possible configuration values are present in this
# autogenerated file.
#
# All configuration values have a default; values that are commented out
# serve to show the default.
import sys
import os
import sphinx_rtd_theme
# If extensions (or modules to document with autodoc) are in another directory,
# add these directories to sys.path here. If the directory is relative to the
# documentation root, use os.path.abspath to make it absolute, like shown here.
#sys.path.insert(0, os.path.abspath('.'))
# -- General configuration ------------------------------------------------
# If your documentation needs a minimal Sphinx version, state it here.
#needs_sphinx = '1.0'
# Add any Sphinx extension module names here, as strings. They can be
# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
# ones.
extensions = [
'sphinx.ext.autodoc',
'sphinx.ext.intersphinx',
'sphinx.ext.todo',
'sphinx.ext.mathjax',
'numpydoc',
'sphinx.ext.autosummary',
]
# Add any paths that contain templates here, relative to this directory.
templates_path = ['_templates']
# The suffix(es) of source filenames.
# You can specify multiple suffix as a list of string:
# source_suffix = ['.rst', '.md']
source_suffix = '.rst'
# The encoding of source files.
#source_encoding = 'utf-8-sig'
# The master toctree document.
master_doc = 'index'
# General information about the project.
project = 'openEMS'
copyright = '2016, Thorsten Liebig'
author = 'Thorsten Liebig'
# The version info for the project you're documenting, acts as replacement for
# |version| and |release|, also used in various other places throughout the
# built documents.
#
# The short X.Y version.
version = '0.0.34'
# The full version, including alpha/beta/rc tags.
release = '0.0.34'
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.
#
# This is also used if you do content translation via gettext catalogs.
# Usually you set "language" from the command line for these cases.
language = None
# There are two options for replacing |today|: either, you set today to some
# non-false value, then it is used:
#today = ''
# Else, today_fmt is used as the format for a strftime call.
#today_fmt = '%B %d, %Y'
# List of patterns, relative to source directory, that match files and
# directories to ignore when looking for source files.
exclude_patterns = ['_build']
# The reST default role (used for this markup: `text`) to use for all
# documents.
#default_role = None
# If true, '()' will be appended to :func: etc. cross-reference text.
#add_function_parentheses = True
# If true, the current module name will be prepended to all description
# unit titles (such as .. function::).
#add_module_names = True
# If true, sectionauthor and moduleauthor directives will be shown in the
# output. They are ignored by default.
#show_authors = False
# The name of the Pygments (syntax highlighting) style to use.
pygments_style = 'sphinx'
# A list of ignored prefixes for module index sorting.
#modindex_common_prefix = []
# If true, keep warnings as "system message" paragraphs in the built documents.
#keep_warnings = False
# If true, `todo` and `todoList` produce output, else they produce nothing.
todo_include_todos = True
# -- Options for HTML output ----------------------------------------------
# The theme to use for HTML and HTML Help pages. See the documentation for
# a list of builtin themes.
html_theme = 'sphinx_rtd_theme'
# Theme options are theme-specific and customize the look and feel of a theme
# further. For a list of options available for each theme, see the
# documentation.
#html_theme_options = {}
# Add any paths that contain custom themes here, relative to this directory.
#html_theme_path = []
# The name for this set of Sphinx documents. If None, it defaults to
# "<project> v<release> documentation".
#html_title = None
# A shorter title for the navigation bar. Default is the same as html_title.
#html_short_title = None
# The name of an image file (relative to this directory) to place at the top
# of the sidebar.
#html_logo = None
# The name of an image file (relative to this directory) to use as a favicon of
# the docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32
# pixels large.
#html_favicon = None
# Add any paths that contain custom static files (such as style sheets) here,
# relative to this directory. They are copied after the builtin static files,
# so a file named "default.css" will overwrite the builtin "default.css".
html_static_path = ['_static']
# Add any extra paths that contain custom files (such as robots.txt or
# .htaccess) here, relative to this directory. These files are copied
# directly to the root of the documentation.
#html_extra_path = []
# If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
# using the given strftime format.
#html_last_updated_fmt = '%b %d, %Y'
# If true, SmartyPants will be used to convert quotes and dashes to
# typographically correct entities.
#html_use_smartypants = True
# Custom sidebar templates, maps document names to template names.
#html_sidebars = {}
# Additional templates that should be rendered to pages, maps page names to
# template names.
#html_additional_pages = {}
# If false, no module index is generated.
#html_domain_indices = True
# If false, no index is generated.
#html_use_index = True
# If true, the index is split into individual pages for each letter.
#html_split_index = False
# If true, links to the reST sources are added to the pages.
#html_show_sourcelink = True
# If true, "Created using Sphinx" is shown in the HTML footer. Default is True.
#html_show_sphinx = True
# If true, "(C) Copyright ..." is shown in the HTML footer. Default is True.
#html_show_copyright = True
# If true, an OpenSearch description file will be output, and all pages will
# contain a <link> tag referring to it. The value of this option must be the
# base URL from which the finished HTML is served.
#html_use_opensearch = ''
# This is the file name suffix for HTML files (e.g. ".xhtml").
#html_file_suffix = None
# Language to be used for generating the HTML full-text search index.
# Sphinx supports the following languages:
# 'da', 'de', 'en', 'es', 'fi', 'fr', 'h', 'it', 'ja'
# 'nl', 'no', 'pt', 'ro', 'r', 'sv', 'tr'
#html_search_language = 'en'
# A dictionary with options for the search language support, empty by default.
# Now only 'ja' uses this config value
#html_search_options = {'type': 'default'}
# The name of a javascript file (relative to the configuration directory) that
# implements a search results scorer. If empty, the default will be used.
#html_search_scorer = 'scorer.js'
# Output file base name for HTML help builder.
htmlhelp_basename = 'openEMSdoc'
# -- Options for LaTeX output ---------------------------------------------
latex_elements = {
# The paper size ('letterpaper' or 'a4paper').
#'papersize': 'letterpaper',
# The font size ('10pt', '11pt' or '12pt').
#'pointsize': '10pt',
# Additional stuff for the LaTeX preamble.
#'preamble': '',
# Latex figure (float) alignment
#'figure_align': 'htbp',
}
# Grouping the document tree into LaTeX files. List of tuples
# (source start file, target name, title,
# author, documentclass [howto, manual, or own class]).
latex_documents = [
(master_doc, 'openEMS.tex', 'openEMS Documentation',
'Thorsten Liebig', 'manual'),
]
# The name of an image file (relative to this directory) to place at the top of
# the title page.
#latex_logo = None
# For "manual" documents, if this is true, then toplevel headings are parts,
# not chapters.
#latex_use_parts = False
# If true, show page references after internal links.
#latex_show_pagerefs = False
# If true, show URL addresses after external links.
#latex_show_urls = False
# Documents to append as an appendix to all manuals.
#latex_appendices = []
# If false, no module index is generated.
#latex_domain_indices = True
# -- Options for manual page output ---------------------------------------
# One entry per manual page. List of tuples
# (source start file, name, description, authors, manual section).
man_pages = [
(master_doc, 'openems', 'openEMS Documentation',
[author], 1)
]
# If true, show URL addresses after external links.
#man_show_urls = False
# -- Options for Texinfo output -------------------------------------------
# Grouping the document tree into Texinfo files. List of tuples
# (source start file, target name, title, author,
# dir menu entry, description, category)
texinfo_documents = [
(master_doc, 'openEMS', 'openEMS Documentation',
author, 'openEMS', 'One line description of project.',
'Miscellaneous'),
]
# Documents to append as an appendix to all manuals.
#texinfo_appendices = []
# If false, no module index is generated.
#texinfo_domain_indices = True
# How to display URL addresses: 'footnote', 'no', or 'inline'.
#texinfo_show_urls = 'footnote'
# If true, do not generate a @detailmenu in the "Top" node's menu.
#texinfo_no_detailmenu = False
numpydoc_show_class_members = False
# Example configuration for intersphinx: refer to the Python standard library.
intersphinx_mapping = {'CSXCAD': ('http://openems.de/doc/CSXCAD/', None)}

63
python/doc/convert_tutorials.py Executable file
View File

@ -0,0 +1,63 @@
#!/usr/bin/python3
# -*- coding: utf-8 -*-
"""
Created on Sat Sep 10 17:12:53 2016
@author: thorsten
"""
import os
import glob
DOC_DIR = os.path.dirname(__file__)
ROOT_DIR = os.path.join(DOC_DIR, '..')
def main():
in_path = os.path.join(ROOT_DIR, 'Tutorials')
fns = glob.glob(os.path.join(in_path, '*.py'))
for fn in fns:
bn = os.path.basename(fn)
out_fn = os.path.join(DOC_DIR, 'Tutorials', '__' + bn.replace('.py', '.txt'))
in_code_block = False
in_ignore_block = False
out_fh = open(out_fn, 'w')
for line in open(fn, 'r'):
if in_ignore_block==False and line.startswith('"""'):
in_ignore_block = True
in_code_block = False
continue
elif in_ignore_block==True and line.startswith('"""'):
in_ignore_block = False
in_code_block = False
continue
elif in_ignore_block==True:
in_code_block = False
continue
elif line.startswith('# -*-'):
continue
elif not line.startswith('##'):
if not in_code_block:
if len(line.strip())==0:
continue
out_fh.write('\n.. code-block:: python\n\n')
in_code_block = True
out_fh.write(' ' + line)
elif line.startswith('###'):
if in_code_block:
out_fh.write('\n')
in_code_block = False
line = line.replace('#','').strip()
out_fh.write('**' + line + '**\n\n')
# out_fh.write('"'*len(line) + '\n')
elif line.startswith('##'):
if in_code_block:
out_fh.write('\n')
in_code_block = False
out_fh.write(line.replace('#','').strip() + '\n')
out_fh.close()
if __name__ == '__main__':
main()

24
python/doc/index.rst Normal file
View File

@ -0,0 +1,24 @@
.. openEMS documentation master file, created by
sphinx-quickstart on Thu Sep 8 20:42:18 2016.
You can adapt this file completely to your liking, but it should at least
contain the root `toctree` directive.
Welcome to openEMS's documentation!
===================================
Contents:
.. toctree::
:maxdepth: 3
Tutorials/index
openEMS_API
Indices and tables
==================
* :ref:`genindex`
* :ref:`modindex`
* :ref:`search`

16
python/doc/nf2ff.rst Normal file
View File

@ -0,0 +1,16 @@
.. _nf2ff:
NF2FF
-----
.. automodule:: openEMS.nf2ff
NF2FF
-----
.. autoclass:: nf2ff
:members:
NF2FF Results
-----------------
.. autoclass:: nf2ff_results
:members:

8
python/doc/openEMS.rst Normal file
View File

@ -0,0 +1,8 @@
.. _openems:
openEMS
-------
.. automodule:: openEMS
:members: openEMS
:undoc-members:

11
python/doc/openEMS_API.rst Normal file
View File

@ -0,0 +1,11 @@
.. _openems_api:
openEMS Python Interface
========================
.. toctree::
openEMS
ports
nf2ff

37
python/doc/ports.rst Normal file
View File

@ -0,0 +1,37 @@
.. _ports:
Ports
-----
.. automodule:: openEMS.ports
Port (Base Class)
-----------------
.. autoclass:: Port
:members:
:show-inheritance:
Lumped Port
-----------
.. autoclass:: LumpedPort
:members:
:show-inheritance:
MSL Port
--------
.. autoclass:: MSLPort
:members:
:show-inheritance:
Waveguide Port
--------------
.. autoclass:: WaveguidePort
:members:
:show-inheritance:
Rect Waveguide Port
-------------------
.. autoclass:: RectWGPort
:members:
:show-inheritance:

4
python/openEMS/__init__.py Normal file
View File

@ -0,0 +1,4 @@
# -*- coding: utf-8 -*-
#
# Shortcut openEMS import
from openEMS.openEMS import openEMS

49
python/openEMS/_nf2ff.pxd Normal file
View File

@ -0,0 +1,49 @@
# -*- coding: utf-8 -*-
#
# Copyright (C) 2015,20016 Thorsten Liebig (Thorsten.Liebig@gmx.de)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
from libcpp.string cimport string
from libcpp.vector cimport vector
from libcpp.complex cimport complex
from libcpp cimport bool
cimport cython.numeric
cdef extern from "openEMS/nf2ff.h":
cdef cppclass cpp_nf2ff "nf2ff":
cpp_nf2ff(vector[float] freq, vector[float] theta, vector[float] phi, vector[float] center, unsigned int numThreads) except +
bool AnalyseFile(string E_Field_file, string H_Field_file)
void SetRadius(float radius)
void SetPermittivity(vector[float] permittivity);
void SetPermeability(vector[float] permeability);
void SetMirror(int _type, int _dir, float pos);
double GetTotalRadPower(size_t f_idx)
double GetMaxDirectivity(size_t f_idx)
complex[double]** GetETheta(size_t f_idx)
complex[double]** GetEPhi(size_t f_idx)
double** GetRadPower(size_t f_idx)
bool Write2HDF5(string filename)
void SetVerboseLevel(int level)
cdef class _nf2ff:
cdef cpp_nf2ff *thisptr

59
python/openEMS/_nf2ff.pyx Normal file
View File

@ -0,0 +1,59 @@
# -*- coding: utf-8 -*-
#
# Copyright (C) 2015,20016 Thorsten Liebig (Thorsten.Liebig@gmx.de)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
cimport _nf2ff
import numpy as np
import os
from CSXCAD.Utilities import CheckNyDir
cdef class _nf2ff:
def __cinit__(self, freq, theta, phi, center, numThreads=0, **kw):
if type(freq) in [float, int]:
freq = list(float(freq))
if type(theta) in [float, int]:
theta = list(float(theta))
if type(phi) in [float, int]:
phi = list(float(phi))
self.thisptr = new cpp_nf2ff(freq, theta, phi, center, numThreads)
if 'verbose' in kw:
self.SetVerboseLevel(kw['verbose'])
del kw['verbose']
assert len(kw)==0, 'Unknown keyword(s): {}'.format(kw)
def AnalyseFile(self, e_file, h_file):
assert os.path.exists(e_file)
assert os.path.exists(h_file)
return self.thisptr.AnalyseFile(e_file.encode('UTF-8'), h_file.encode('UTF-8'))
def SetMirror(self, mirr_type, ny, pos):
if mirr_type<=0:
return
assert mirr_type<3
ny = CheckNyDir(ny)
self.thisptr.SetMirror(mirr_type, ny, pos)
def SetRadius(self, radius):
self.thisptr.SetRadius(radius)
def Write2HDF5(self, filename):
return self.thisptr.Write2HDF5(filename.encode('UTF-8'))
def SetVerboseLevel(self, level):
self.thisptr.SetVerboseLevel(level)

77
python/openEMS/automesh.py Normal file
View File

@ -0,0 +1,77 @@
# -*- coding: utf-8 -*-
"""
Created on Sun Feb 19 20:29:25 2017
@author: thorsten
"""
import sys
import numpy as np
from CSXCAD import CSPrimitives
from CSXCAD.Utilities import CheckNyDir, GetMultiDirs
def mesh_hint_from_primitive(primitive, dirs, **kw):
if primitive.GetType() is CSPrimitives.POINT:
return mesh_hint_from_point(primitive, dirs, **kw)
if primitive.GetType() is CSPrimitives.BOX:
return mesh_hint_from_box(primitive, dirs, **kw)
else:
return None
def mesh_hint_from_point(point, dirs, **kw):
""" mesh_hint_from_point(point, dirs)
Get a grid hint for the coordinates of the point.
:param dirs: str -- 'x','y','z' or 'xy', 'yz' or 'xyz' or 'all'
:returns: (3,) list of mesh hints
"""
hint = [None, None, None]
coord = point.GetCoord()
for ny in GetMultiDirs(dirs):
hint[ny] = [coord[ny],]
return hint
def mesh_hint_from_box(box, dirs, **kw):
""" mesh_hint_from_box(box, dirs, metal_edge_res=None, **kw)
Get a grid hint for the edges of the given box with an an optional 2D metal
edge resolution.
:param dirs: str -- 'x','y','z' or 'xy', 'yz' or 'xyz' or 'all'
:param metal_edge_res: float -- 2D flat edge resolution
:returns: (3,) list of mesh hints
"""
metal_edge_res = kw.get('metal_edge_res', None)
up_dir = kw.get('up_dir' , True)
down_dir = kw.get('down_dir', True)
if metal_edge_res is None:
mer = 0
else:
mer = np.array([-1.0, 2.0])/3 * metal_edge_res
if box.HasTransform():
sys.stderr.write('FDTD::automesh: Warning, cannot add edges to grid with transformations enabled\n')
return
hint = [None, None, None]
start = np.fmin(box.GetStart(), box.GetStop())
stop = np.fmax(box.GetStart(), box.GetStop())
for ny in GetMultiDirs(dirs):
hint[ny] = []
if metal_edge_res is not None and stop[ny]-start[ny]>metal_edge_res:
if down_dir:
hint[ny].append(start[ny]-mer[0])
hint[ny].append(start[ny]-mer[1])
if up_dir:
hint[ny].append(stop[ny]+mer[0])
hint[ny].append(stop[ny]+mer[1])
elif stop[ny]-start[ny]:
if down_dir:
hint[ny].append(start[ny])
if up_dir:
hint[ny].append(stop[ny])
else:
hint[ny].append(start[ny])
return hint

210
python/openEMS/nf2ff.py Normal file
View File

@ -0,0 +1,210 @@
# -*- coding: utf-8 -*-
#
# Copyright (C) 2015,20016 Thorsten Liebig (Thorsten.Liebig@gmx.de)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
import os
import numpy as np
import h5py
from openEMS import _nf2ff
from openEMS import utilities
class nf2ff:
"""
Create an nf2ff recording box. The nf2ff can either record in time-domain
or frequency-domain. Further more certain directions and boundary condition
mirroring can be enabled/disabled.
:param name: str -- Name for this recording box.
:param start/stop: (3,) array -- Box start/stop coordinates.
:param directions: (6,) bool array -- Enable/Disables directions.
:param mirror: (6,) int array -- 0 (Off), 1 (PEC) or 2 (PMC) boundary mirroring
:param frequency: array like -- List of frequencies (FD-domain recording)
"""
def __init__(self, CSX, name, start, stop, **kw):
self.CSX = CSX
self.name = name
self.start = start
self.stop = stop
self.freq = None
self.theta = None
self.phi = None
self.center = None
self.directions = [True]*6 # all directions by default
if 'directions' in kw:
self.directions = kw['directions']
del kw['directions']
assert len(self.directions)==6
self.mirror = [0]*6
if 'mirror' in kw:
self.mirror = kw['mirror']
del kw['mirror']
assert len(self.mirror)==6
self.dump_type = 0 # default Et/Ht
self.dump_mode = 1 # default cell interpolated
self.freq = None # broadband recording by defualt
if 'frequency' in kw:
self.freq = kw['frequency']
del kw['frequency']
self.dump_type = 10 # Ef/Hf
if np.isscalar(self.freq):
self.freq = [self.freq]
self.e_file = '{}_E'.format(self.name)
self.h_file = '{}_H'.format(self.name)
self.e_dump = CSX.AddDump(self.e_file, dump_type=self.dump_type , dump_mode=self.dump_mode, file_type=1, **kw)
self.h_dump = CSX.AddDump(self.h_file, dump_type=self.dump_type+1, dump_mode=self.dump_mode, file_type=1, **kw)
if self.freq is not None:
self.e_dump.SetFrequency(self.freq)
self.h_dump.SetFrequency(self.freq)
# print(self.directions)
for ny in range(3):
pos = 2*ny
if self.directions[pos]:
l_start = np.array(start)
l_stop = np.array(stop)
l_stop[ny] = l_start[ny]
self.e_dump.AddBox(l_start, l_stop)
self.h_dump.AddBox(l_start, l_stop)
if self.directions[pos+1]:
l_start = np.array(start)
l_stop = np.array(stop)
l_start[ny] = l_stop[ny]
self.e_dump.AddBox(l_start, l_stop)
self.h_dump.AddBox(l_start, l_stop)
def CalcNF2FF(self, sim_path, freq, theta, phi, radius=1, center=[0,0,0], outfile=None, read_cached=False, verbose=0):
""" CalcNF2FF(sim_path, freq, theta, phi, center=[0,0,0], outfile=None, read_cached=True, verbose=0):
Calculate the far-field after the simulation is done.
:param sim_path: str -- Simulation path
:param freq: array like -- list of frequency for transformation
:param theta/phi: array like -- Theta/Phi angles to calculate the far-field
:param radius: float -- Radius to calculate the far-field (default is 1m)
:param center: (3,) array -- phase center, must be inside the recording box
:param outfile: str -- File to save results in. (defaults to recording name)
:param read_cached: bool -- enable/disable read already existing results (default off)
:param verbose: int -- set verbose level (default 0)
:returns: nf2ff_results class instance
"""
if np.isscalar(freq):
freq = [freq]
self.freq = freq
if np.isscalar(theta):
theta = [theta]
self.theta = theta
if np.isscalar(phi):
phi = [phi]
self.phi = phi
self.center = center
if outfile is None:
fn = os.path.join(sim_path, self.name + '.h5')
else:
fn = os.path.join(sim_path, outfile)
if not read_cached or not os.path.exists(fn):
nfc = _nf2ff._nf2ff(self.freq, np.deg2rad(theta), np.deg2rad(phi), center, verbose=verbose)
for ny in range(3):
nfc.SetMirror(self.mirror[2*ny] , ny, self.start[ny])
nfc.SetMirror(self.mirror[2*ny+1], ny, self.stop[ny])
nfc.SetRadius(radius)
for n in range(6):
fn_e = os.path.join(sim_path, self.e_file + '_{}.h5'.format(n))
fn_h = os.path.join(sim_path, self.h_file + '_{}.h5'.format(n))
if os.path.exists(fn_e) and os.path.exists(fn_h):
assert nfc.AnalyseFile(fn_e, fn_h)
nfc.Write2HDF5(fn)
result = nf2ff_results(fn)
if result.phi is not None:
assert np.abs((result.r-radius)/radius)<1e-6, 'Radius does not match. Did you read an invalid chached result? Try "read_cached=False"'
assert utilities.Check_Array_Equal(np.rad2deg(result.theta), self.theta, 1e-4), 'Theta array does not match. Did you read an invalid chached result? Try "read_cached=False"'
assert utilities.Check_Array_Equal(np.rad2deg(result.phi), self.phi, 1e-4), 'Phi array does not match. Did you read an invalid chached result? Try "read_cached=False"'
assert utilities.Check_Array_Equal(result.freq, self.freq, 1e-6, relative=True), 'Frequency array does not match. Did you read an invalid chached result? Try "read_cached=False"'
return result
class nf2ff_results:
"""
nf2ff result class containing all results obtained by the nf2ff calculation.
Usueally returned from nf2ff.CalcNF2FF
Available attributes:
* `fn` : file name
* `theta`: theta angles
* `phi` : phi angles
* `r` : radius
* `freq` : frequencies
* `Dmax` : directivity over frequency
* `Prad` : total radiated power over frequency
* `E_theta` : theta component of electric field over frequency/theta/phi
* `E_phi` : phi component of electric field over frequency/theta/phi
* `E_norm` : abs component of electric field over frequency/theta/phi
* `E_cprh` : theta component of electric field over frequency/theta/phi
* `E_cplh` : theta component of electric field over frequency/theta/phi
* `P_rad` : radiated power (S) over frequency/theta/phi
"""
def __init__(self, fn):
self.fn = fn
h5_file = h5py.File(fn, 'r')
mesh_grp = h5_file['Mesh']
self.phi = np.array(mesh_grp['phi'])
self.theta = np.array(mesh_grp['theta'])
self.r = np.array(mesh_grp['r'])
data = h5_file['nf2ff']
self.freq = np.array(data.attrs['Frequency'])
self.Dmax = np.array(data.attrs['Dmax'])
self.Prad = np.array(data.attrs['Prad'])
THETA, PHI = np.meshgrid(self.theta, self.phi, indexing='ij')
cos_phi = np.cos(PHI)
sin_phi = np.sin(PHI)
self.E_theta = []
self.E_phi = []
self.P_rad = []
self.E_norm = []
self.E_cprh = []
self.E_cplh = []
for n in range(len(self.freq)):
E_theta = np.array(h5_file['/nf2ff/E_theta/FD/f{}_real'.format(n)]) + 1j*np.array(h5_file['/nf2ff/E_theta/FD/f{}_imag'.format(n)])
E_theta = np.swapaxes(E_theta, 0, 1)
E_phi = np.array(h5_file['/nf2ff/E_phi/FD/f{}_real'.format(n)]) + 1j*np.array(h5_file['/nf2ff/E_phi/FD/f{}_imag'.format(n)])
E_phi = np.swapaxes(E_phi, 0, 1)
self.P_rad .append(np.swapaxes(np.array(h5_file['/nf2ff/P_rad/FD/f{}'.format(n)]), 0, 1))
self.E_theta.append(E_theta)
self.E_phi .append(E_phi)
self.E_norm .append(np.sqrt(np.abs(E_theta)**2 + np.abs(E_phi)**2))
self.E_cprh .append((cos_phi+1j*sin_phi) * (E_theta+1j*E_phi)/np.sqrt(2.0))
self.E_cplh .append((cos_phi-1j*sin_phi) * (E_theta-1j*E_phi)/np.sqrt(2.0))

60
python/openEMS/openEMS.pxd Normal file
View File

@ -0,0 +1,60 @@
# -*- coding: utf-8 -*-
#
# Copyright (C) 2015,20016 Thorsten Liebig (Thorsten.Liebig@gmx.de)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
from libcpp.string cimport string
from libcpp cimport bool
from CSXCAD.CSXCAD cimport _ContinuousStructure, ContinuousStructure
cdef extern from "openEMS/openems.h":
cdef cppclass _openEMS "openEMS":
_openEMS() except +
void SetNumberOfTimeSteps(unsigned int val)
void SetCSX(_ContinuousStructure* csx)
void SetEndCriteria(double val)
void SetOverSampling(int val)
void SetCellConstantMaterial(bool val)
void SetCylinderCoords(bool val)
void SetupCylinderMultiGrid(string val)
void SetTimeStepMethod(int val)
void SetTimeStep(double val)
void SetTimeStepFactor(double val)
void SetMaxTime(double val)
void Set_BC_Type(int idx, int _type)
int Get_BC_Type(int idx)
void Set_BC_PML(int idx, unsigned int size)
int Get_PML_Size(int idx)
void Set_Mur_PhaseVel(int idx, double val)
void SetGaussExcite(double f0, double fc)
void SetVerboseLevel(int level)
int SetupFDTD()
void RunFDTD()
@staticmethod
void WelcomeScreen()
cdef class openEMS:
cdef _openEMS *thisptr
cdef readonly ContinuousStructure __CSX # hold a C++ instance which we're wrapping

440
python/openEMS/openEMS.pyx Normal file
View File

@ -0,0 +1,440 @@
# -*- coding: utf-8 -*-
#
# Copyright (C) 2015,20016 Thorsten Liebig (Thorsten.Liebig@gmx.de)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
import os, sys, shutil
import numpy as np
cimport openEMS
from . import ports, nf2ff, automesh
from CSXCAD.Utilities import GetMultiDirs
cdef class openEMS:
""" openEMS
This class is the main control class for the FDTD options and setup and
to run the final simulation.
Examples
--------
>>> CSX = CSXCAD.ContinuousStructure()
>>>
>>> grid = CSX.GetGrid()
>>> grid.SetLines('x', np.arange(-50,50,1))
>>> grid.SetLines('y', np.arange(-50,50,1))
>>> grid.SetLines('z', np.arange(-2,2.1,1))
>>> grid.SetDeltaUnit(1e-3)
>>>
>>> FDTD = openEMS(NrTS=1e4, EndCriteria=1e-4)
>>>
>>> FDTD.SetCSX(CSX)
>>> FDTD.SetBoundaryCond(['PML_8', 'PML_8', 'PML_8', 'PML_8', 'PEC', 'PEC'])
>>> FDTD.SetGaussExcite(0, 10e9)
>>>
>>> FDTD.AddLumpedPort(port_nr=1, R=50, start=[10, 0, -2], stop=[10, 0, 2], p_dir='z', excite=1)
>>>
>>> FDTD.Run(sim_path='/tmp/test')
:param NrTS: max. number of timesteps to simulate (e.g. default=1e9)
:param EndCriteria: end criteria, e.g. 1e-5, simulations stops if energy has decayed by this value (<1e-4 is recommended, default=1e-5)
:param MaxTime: max. real time in seconds to simulate
:param OverSampling: nyquist oversampling of time domain dumps
:param CoordSystem: choose coordinate system (0 Cartesian, 1 Cylindrical)
:param MultiGrid: define a cylindrical sub-grid radius
:param TimeStep: force to use a given timestep (dangerous!)
:param TimeStepFactor: reduce the timestep by a given factor (>0 to <=1)
:param TimeStepMethod: 1 or 3 chose timestep method (1=CFL, 3=Rennigs (default))
:param CellConstantMaterial: set to 1 to assume a material is constant inside a cell (material probing in cell center)
"""
@staticmethod
def WelcomeScreen():
"""
Show the openEMS welcome screen.
"""
_openEMS.WelcomeScreen()
def __cinit__(self, *args, **kw):
self.thisptr = new _openEMS()
self.__CSX = None
if 'NrTS' in kw:
self.SetNumberOfTimeSteps(kw['NrTS'])
del kw['NrTS']
else:
self.SetNumberOfTimeSteps(1e9)
if 'EndCriteria' in kw:
self.SetEndCriteria(kw['EndCriteria'])
del kw['EndCriteria']
if 'MaxTime' in kw:
self.SetMaxTime(kw['MaxTime'])
del kw['MaxTime']
if 'OverSampling' in kw:
self.SetOverSampling(kw['OverSampling'])
del kw['OverSampling']
if 'CoordSystem' in kw:
self.SetCoordSystem(kw['CoordSystem'])
del kw['CoordSystem']
if 'TimeStep' in kw:
self.SetTimeStep(kw['TimeStep'])
del kw['TimeStep']
if 'TimeStepFactor' in kw:
self.SetTimeStepFactor(kw['TimeStepFactor'])
del kw['TimeStepFactor']
if 'TimeStepMethod' in kw:
self.SetTimeStepMethod(kw['TimeStepMethod'])
del kw['TimeStepMethod']
if 'CellConstantMaterial' in kw:
self.SetCellConstantMaterial(kw['CellConstantMaterial'])
del kw['CellConstantMaterial']
if 'MultiGrid' in kw:
self.SetMultiGrid(kw['MultiGrid'])
del kw['MultiGrid']
assert len(kw)==0, 'Unknown keyword arguments: "{}"'.format(kw)
def __dealloc__(self):
del self.thisptr
if self.__CSX is not None:
self.__CSX.thisptr = NULL
def SetNumberOfTimeSteps(self, val):
""" SetNumberOfTimeSteps(val)
Set the number of timesteps. E.g. 5e4 (default is 1e9)
"""
self.thisptr.SetNumberOfTimeSteps(val)
def SetEndCriteria(self, val):
""" SetEndCriteria(val)
Set the end critera value. E.g. 1e-6 for -60dB
"""
self.thisptr.SetEndCriteria(val)
def SetOverSampling(self, val):
""" SetOverSampling(val)
Set the time domain signal oversampling as multiple of the Nyquist-rate.
"""
self.thisptr.SetOverSampling(val)
def SetCellConstantMaterial(self, val):
""" SetCellConstantMaterial(val)
Set cell material averaging to assume constant material inside each primary cell. (Advanced option)
:param val: bool -- Enable or Disable (default disabled)
"""
self.thisptr.SetCellConstantMaterial(val)
def SetCoordSystem(self, val):
""" SetCoordSystem(val)
Set the coordinate system. 0 --> Cartesian (default), 1 --> cylindrical
"""
assert (val==0 or val==1), 'SetCoordSystem: Invalid coordinate system'
if val==0:
self.thisptr.SetCylinderCoords(False)
elif val==1:
self.thisptr.SetCylinderCoords(True)
def SetMultiGrid(self, radii):
""" SetMultiGrid(radii)
Define radii at which a cylindrical multi grid should be defined.
:param radii: array like, multigrid radii
See Also
--------
openEMS.SetCylinderCoords
"""
assert len(radii)>0, 'SetMultiGrid: invalid multi grid definition'
grid_str = ','.join(['{}'.format(x) for x in radii])
self.thisptr.SetupCylinderMultiGrid(grid_str.encode('UTF-8'))
def SetCylinderCoords(self):
""" SetCylinderCoords()
Enable use of cylindircal coordinates.
See Also
--------
openEMS.SetMultiGrid
"""
self.thisptr.SetCylinderCoords(True)
def SetTimeStepMethod(self, val):
""" SetTimeStepMethod(val)
Set the time step calculation method. (Advanced option)
Options:
* 1: CFL criteria
* 3: Advanced Rennings criteria (default)
:param val: int -- 1 or 3 (See above)
"""
self.thisptr.SetTimeStepMethod(val)
def SetTimeStep(self, val):
""" SetTimeStep(val)
Set/force the timestep. (Advanced option)
It is highly recommended to not use this method! You may use the
SetTimeStepFactor instead to reduce the time step if necessary!
"""
self.thisptr.SetTimeStep(val)
def SetTimeStepFactor(self, val):
""" SetTimeStepFactor(val)
Set a time step factor (>0..1) to increase FDTD stability.
:param val: float -- >0..1
"""
self.thisptr.SetTimeStepFactor(val)
def SetMaxTime(self, val):
""" SetMaxTime(val)
Set max simulation time for a max. number of timesteps.
"""
self.thisptr.SetMaxTime(val)
def SetGaussExcite(self, f0, fc):
""" SetGaussExcite(f0, fc)
Set a Gaussian pulse as excitation signal.
:param f0: float -- Center frequency in Hz.
:param fc: float -- -20dB bandwidth in Hz.
"""
self.thisptr.SetGaussExcite(f0, fc)
def SetBoundaryCond(self, BC):
""" SetBoundaryCond(BC)
Set the boundary conditions for all six FDTD directions.
Options:
* 0 or 'PEC' : perfect electric conductor (default)
* 1 or 'PMC' : perfect magnetic conductor, useful for symmetries
* 2 or 'MUR' : simple MUR absorbing boundary conditions
* 3 or 'PML-8' : PML absorbing boundary conditions
:param BC: (8,) array or list -- see options above
"""
assert len(BC)==6
for n in range(len(BC)):
if type(BC[n])==int:
self.thisptr.Set_BC_Type(n, BC[n])
continue
if BC[n] in ['PEC', 'PMC', 'MUR']:
self.thisptr.Set_BC_Type(n, ['PEC', 'PMC', 'MUR'].index(BC[n]))
continue
if BC[n].startswith('PML_'):
size = int(BC[n].strip('PML_'))
self.thisptr.Set_BC_PML(n, size)
continue
raise Exception('Unknown boundary condition')
def AddLumpedPort(self, port_nr, R, start, stop, p_dir, excite=0, **kw):
""" AddLumpedPort(port_nr, R, start, stop, p_dir, excite=0, **kw)
Add a lumped port wit the given values and location.
See Also
--------
openEMS.ports.LumpedPort
"""
assert self.__CSX is not None, 'AddLumpedPort: CSX is not set!'
port = ports.LumpedPort(self.__CSX, port_nr, R, start, stop, p_dir, excite, **kw)
edges2grid = kw.get('edges2grid', None)
if edges2grid is not None:
grid = self.__CSX.GetGrid()
for n in GetMultiDirs(edges2grid):
grid.AddLine(n, start[n])
if start[n] != stop[n]:
grid.AddLine(n, stop[n])
return port
def AddWaveGuidePort(self, port_nr, start, stop, p_dir, E_func, H_func, kc, excite=0, **kw):
""" AddWaveGuidePort(self, port_nr, start, stop, p_dir, E_func, H_func, kc, excite=0, **kw)
Add a arbitrary waveguide port.
See Also
--------
openEMS.ports.WaveguidePort
"""
assert self.__CSX is not None, 'AddWaveGuidePort: CSX is not set!'
return ports.WaveguidePort(self.__CSX, port_nr, start, stop, p_dir, E_func, H_func, kc, excite, **kw)
def AddRectWaveGuidePort(self, port_nr, start, stop, p_dir, a, b, mode_name, excite=0, **kw):
""" AddRectWaveGuidePort(port_nr, start, stop, p_dir, a, b, mode_name, excite=0, **kw)
Add a rectilinear waveguide port.
See Also
--------
openEMS.ports.RectWGPort
"""
assert self.__CSX is not None, 'AddRectWaveGuidePort: CSX is not set!'
return ports.RectWGPort(self.__CSX, port_nr, start, stop, p_dir, a, b, mode_name, excite, **kw)
def AddMSLPort(self, port_nr, metal_prop, start, stop, prop_dir, exc_dir, excite=0, **kw):
""" AddMSLPort(port_nr, metal_prop, start, stop, prop_dir, exc_dir, excite=0, **kw)
Add a microstrip transmission line port.
See Also
--------
openEMS.ports.MSLPort
"""
assert self.__CSX is not None, 'AddMSLPort: CSX is not set!'
return ports.MSLPort(self.__CSX, port_nr, metal_prop, start, stop, prop_dir, exc_dir, excite, **kw)
def CreateNF2FFBox(self, name='nf2ff', start=None, stop=None, **kw):
""" CreateNF2FFBox(name='nf2ff', start=None, stop=None, **kw)
Create a near-field to far-field box.
This method will automatically adept the recording box to the current
FDTD grid and boundary conditions.
Notes
-----
* Make sure the mesh grid and all boundary conditions are finially defined.
See Also
--------
openEMS.nf2ff.nf2ff
"""
assert self.__CSX is not None, 'CreateNF2FFBox: CSX is not set!'
directions = [True]*6
mirror = [0]*6
BC_size = [0]*6
BC_type = [0]*6
for n in range(6):
BC_type[n] = self.thisptr.Get_BC_Type(n)
if BC_type[n]==0:
directions[n]= False
mirror[n] = 1 # PEC mirror
elif BC_type[n]==1:
directions[n]= False
mirror[n] = 2 # PMC mirror
elif BC_type[n]==2:
BC_size[n] = 2
elif BC_type[n]==3:
BC_size[n] = self.thisptr.Get_PML_Size(n)+1
if start is None or stop is None:
grid = self.__CSX.GetGrid()
assert grid.IsValid(), 'Error::CreateNF2FFBox: Grid is invalid'
start = np.zeros(3)
stop = np.zeros(3)
for n in range(3):
l = grid.GetLines(n)
BC_type = self.thisptr.Get_BC_Type(2*n)
assert len(l)>(BC_size[2*n]+BC_size[2*n+1]), 'Error::CreateNF2FFBox: not enough lines in some direction'
start[n] = l[BC_size[2*n]]
stop[n] = l[-1*BC_size[2*n+1]-1]
return nf2ff.nf2ff(self.__CSX, name, start, stop, directions=directions, mirror=mirror, **kw)
def SetCSX(self, ContinuousStructure CSX):
""" SetCSX(CSX)
Set the CSXCAD Continuous Structure for CAD data handling.
See Also
--------
CSXCAD.ContinuousStructure
"""
self.__CSX = CSX
self.thisptr.SetCSX(CSX.thisptr)
def GetCSX(self):
return self.__CSX
def AddEdges2Grid(self, dirs, primitives=None, properties=None, **kw):
""" AddEdges2Grid(primitives, dirs, **kw)
Add the edges of the given primitives to the FDTD grid.
:param dirs: primitives -- one or more primitives
:param dirs: str -- 'x','y','z' or 'xy', 'yz' or 'xyz' or 'all'
"""
csx = self.GetCSX()
if csx is None:
raise Exception('AddEdges2Grid: Unable to access CSX!')
prim_list = []
if primitives is not None and type(primitives) is not list:
prim_list.append(primitives)
elif primitives is not None:
prim_list += primitives
if properties is not None and type(properties) is not list:
prim_list += properties.GetAllPrimitives()
elif primitives is not None:
for prop in properties:
prim_list += prop.GetAllPrimitives()
grid = csx.GetGrid()
for prim in prim_list:
hint = automesh.mesh_hint_from_primitive(prim, dirs, **kw)
if hint is None:
continue
for n in range(3):
if hint[n] is None:
continue
grid.AddLine(n, hint[n])
def Run(self, sim_path, cleanup=False, setup_only=False, verbose=None):
""" Run(sim_path, cleanup=False, setup_only=False, verbose=None)
Run the openEMS FDTD simulation.
:param sim_path: str -- path to run in and create result data
:param cleanup: bool -- remove exisiting sim_path to cleanup old results
:param setup_only: bool -- only perform FDTD setup, do not run simulation
:param verbose: int -- set the openEMS verbosity level 0..3
"""
if cleanup and os.path.exists(sim_path):
shutil.rmtree(sim_path)
os.mkdir(sim_path)
if not os.path.exists(sim_path):
os.mkdir(sim_path)
cwd = os.getcwd()
os.chdir(sim_path)
if verbose is not None:
self.thisptr.SetVerboseLevel(verbose)
assert os.getcwd() == sim_path
_openEMS.WelcomeScreen()
cdef int EC
EC = self.thisptr.SetupFDTD()
if EC!=0:
print('Run: Setup failed, error code: {}'.format(EC))
if setup_only or EC!=0:
return EC
self.thisptr.RunFDTD()

26
python/openEMS/physical_constants.py Normal file
View File

@ -0,0 +1,26 @@
# -*- coding: utf-8 -*-
#
# Copyright (C) 2015,20016 Thorsten Liebig (Thorsten.Liebig@gmx.de)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
import numpy as np
C0 = 299792458 # m/s
MUE0 = 4e-7*np.pi # N/A^2
EPS0 = 1/(MUE0*C0**2) # F/m
# free space wave impedance
Z0 = np.sqrt(MUE0/EPS0) # Ohm

433
python/openEMS/ports.py Normal file
View File

@ -0,0 +1,433 @@
# -*- coding: utf-8 -*-
#
# Copyright (C) 2015,20016 Thorsten Liebig (Thorsten.Liebig@gmx.de)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
import os
import numpy as np
from CSXCAD.Utilities import CheckNyDir
from openEMS import utilities
from openEMS.physical_constants import *
class UI_data:
def __init__(self, fns, path, freq, signal_type='pulse', **kw):
self.path = path
if type(fns)==str:
fns = [fns]
self.fns = fns
if np.isscalar(freq):
freq = [freq]
self.freq = freq
self.ui_time = []
self.ui_val = []
self.ui_f_val = []
for fn in fns:
tmp = np.loadtxt(os.path.join(path, fn),comments='%')
self.ui_time.append(tmp[:,0])
self.ui_val.append(tmp[:,1])
self.ui_f_val.append(utilities.DFT_time2freq(tmp[:,0], tmp[:,1], freq, signal_type=signal_type))
# Port Base-Class
class Port:
"""
The port base class.
:param CSX: Continuous Structure
:param port_nr: int -- port number
:param R: float -- port reference impedance, e.g. 50 (Ohms)
:param start, stop: (3,) array -- Start/Stop box coordinates
:param p_dir: int -- port direction
:param excite: float -- port excitation amplitude
:param priority: int -- priority of all contained primtives
:param PortNamePrefix: str -- a prefix for all ports-names
:param delay: float -- a positiv delay value to e.g. emulate a phase shift
"""
def __init__(self, CSX, port_nr, start, stop, excite, **kw):
self.CSX = CSX
self.number = port_nr
self.excite = excite
self.start = np.array(start, np.float)
self.stop = np.array(stop, np.float)
self.Z_ref = None
self.U_filenames = []
self.I_filenames = []
self.priority = 0
if 'priority' in kw:
self.priority = kw['priority']
self.prefix = ''
if 'PortNamePrefix' in kw:
self.prefix = kw['PortNamePrefix']
self.delay = 0
if 'delay' in kw:
self.delay = kw['delay']
self.lbl_temp = self.prefix + 'port_{}' + '_{}'.format(self.number)
def ReadUIData(self, sim_path, freq, signal_type ='pulse'):
self.u_data = UI_data(self.U_filenames, sim_path, freq, signal_type )
self.uf_tot = 0
self.ut_tot = 0
for n in range(len(self.U_filenames)):
self.uf_tot += self.u_data.ui_f_val[n]
self.ut_tot += self.u_data.ui_val[n]
self.i_data = UI_data(self.I_filenames, sim_path, freq, signal_type )
self.if_tot = 0
self.it_tot = 0
for n in range(len(self.U_filenames)):
self.if_tot += self.i_data.ui_f_val[n]
self.it_tot += self.i_data.ui_val[n]
def CalcPort(self, sim_path, freq, ref_impedance=None, ref_plane_shift=None, signal_type='pulse'):
self.ReadUIData(sim_path, freq, signal_type)
if ref_impedance is not None:
self.Z_ref = ref_impedance
assert self.Z_ref is not None
if ref_plane_shift is not None:
assert hasattr(self, 'beta')
shift = ref_plane_shift
if self.measplane_shift:
shift -= self.measplane_shift
shift *= self.CSX.GetGrid().GetDeltaUnit()
phase = np.real(self.beta)*shift
uf_tot = self.uf_tot * np.cos(-phase) + 1j * self.if_tot * self.Z_ref * np.sin(-phase)
if_tot = self.if_tot * np.cos(-phase) + 1j * self.uf_tot / self.Z_ref * np.sin(-phase)
self.uf_tot = uf_tot
self.if_tot = if_tot
self.uf_inc = 0.5 * ( self.uf_tot + self.if_tot * self.Z_ref )
self.if_inc = 0.5 * ( self.if_tot + self.uf_tot / self.Z_ref )
self.uf_ref = self.uf_tot - self.uf_inc
self.if_ref = self.if_inc - self.if_tot
if type(self.Z_ref) == float:
self.ut_inc = 0.5 * ( self.ut_tot + self.it_tot * self.Z_ref )
self.it_inc = 0.5 * ( self.it_tot + self.ut_tot / self.Z_ref )
self.ut_ref = self.ut_tot - self.ut_inc
self.it_ref = self.it_inc - self.it_tot
# calc some more port parameter
# incoming power
self.P_inc = 0.5*np.real(self.uf_inc*np.conj(self.if_inc))
# reflected power
self.P_ref = 0.5*np.real(self.uf_ref*np.conj(self.if_ref))
# accepted power (incoming - reflected)
self.P_acc = 0.5*np.real(self.uf_tot*np.conj(self.if_tot))
class LumpedPort(Port):
"""
The lumped port.
See Also
--------
Port
"""
def __init__(self, CSX, port_nr, R, start, stop, exc_dir, excite=0, **kw):
super(LumpedPort, self).__init__(CSX, port_nr=port_nr, start=start, stop=stop, excite=excite, **kw)
self.R = R
self.exc_ny = CheckNyDir(exc_dir)
self.direction = np.sign(self.stop[self.exc_ny]-self.start[self.exc_ny])
assert self.start[self.exc_ny]!=self.stop[self.exc_ny], 'LumpedPort: start and stop may not be identical in excitation direction'
if self.R > 0:
lumped_R = CSX.AddLumpedElement(self.lbl_temp.format('resist'), ny=self.exc_ny, caps=True, R=self.R)
elif self.R==0:
lumped_R = CSX.AddMetal(self.lbl_temp.format('resist'))
lumped_R.AddBox(self.start, self.stop, priority=self.priority)
if excite!=0:
exc_vec = np.zeros(3)
exc_vec[self.exc_ny] = -1*self.direction*excite
exc = CSX.AddExcitation(self.lbl_temp.format('excite'), exc_type=0, exc_val=exc_vec, delay=self.delay)
exc.AddBox(self.start, self.stop, priority=self.priority)
self.U_filenames = [self.lbl_temp.format('ut'), ]
u_start = 0.5*(self.start+self.stop)
u_start[self.exc_ny] = self.start[self.exc_ny]
u_stop = 0.5*(self.start+self.stop)
u_stop[self.exc_ny] = self.stop[self.exc_ny]
u_probe = CSX.AddProbe(self.U_filenames[0], p_type=0, weight=-1*self.direction)
u_probe.AddBox(u_start, u_stop)
self.I_filenames = [self.lbl_temp.format('it'), ]
i_start = np.array(self.start)
i_start[self.exc_ny] = 0.5*(self.start[self.exc_ny]+self.stop[self.exc_ny])
i_stop = np.array(self.stop)
i_stop[self.exc_ny] = 0.5*(self.start[self.exc_ny]+self.stop[self.exc_ny])
i_probe = CSX.AddProbe(self.I_filenames[0], p_type=1, weight=self.direction, norm_dir=self.exc_ny)
i_probe.AddBox(i_start, i_stop)
def CalcPort(self, sim_path, freq, ref_impedance=None, ref_plane_shift=None, signal_type='pulse'):
if ref_impedance is None:
self.Z_ref = self.R
if ref_plane_shift is not None:
Warning('A lumped port does not support a reference plane shift! Ignoring...')
super(LumpedPort, self).CalcPort(sim_path, freq, ref_impedance, ref_plane_shift, signal_type)
class MSLPort(Port):
"""
The microstrip transmission line port.
:param prop_dir: int/str -- direction of propagation
See Also
--------
Port
"""
def __init__(self, CSX, port_nr, metal_prop, start, stop, prop_dir, exc_dir, excite=0, **kw):
super(MSLPort, self).__init__(CSX, port_nr=port_nr, start=start, stop=stop, excite=excite, **kw)
self.exc_ny = CheckNyDir(exc_dir)
self.prop_ny = CheckNyDir(prop_dir)
self.direction = np.sign(stop[self.prop_ny]-start[self.prop_ny])
self.upside_down = np.sign(stop[self.exc_ny] -start[self.exc_ny])
assert (self.start!=self.stop).all()
# assert stop[self.prop_ny]!=start[self.prop_ny], 'port length in propergation direction may not be zero!'
# assert stop[self.exc_ny] !=start[self.exc_ny], 'port length in propergation direction may not be zero!'
assert self.exc_ny!=self.prop_ny
self.feed_shift = 0
if 'FeedShift' in kw:
self.feed_shift = kw['FeedShift']
self.measplane_shift = 0.5*np.abs(self.start[self.prop_ny]-self.stop[self.prop_ny])
if 'MeasPlaneShift' in kw:
self.measplane_shift = kw['MeasPlaneShift']
self.measplane_pos = self.start[self.prop_ny] + self.measplane_shift*self.direction
self.feed_R = np.inf
if 'Feed_R' in kw:
self.feed_R = kw['Feed_R']
# add metal msl-plane
MSL_start = np.array(self.start)
MSL_stop = np.array(self.stop)
MSL_stop[self.exc_ny] = MSL_start[self.exc_ny]
metal_prop.AddBox(MSL_start, MSL_stop, priority=self.priority )
mesh = CSX.GetGrid()
prop_lines = mesh.GetLines(self.prop_ny)
assert len(prop_lines)>5, 'At least 5 lines in propagation direction required!'
meas_pos_idx = np.argmin(np.abs(prop_lines-self.measplane_pos))
if meas_pos_idx==0:
meas_pos_idx=1
if meas_pos_idx>=len(prop_lines)-1:
meas_pos_idx=len(prop_lines)-2
self.measplane_shift = np.abs(self.start[self.prop_ny]-prop_lines[meas_pos_idx])
prope_idx = np.array([meas_pos_idx-1, meas_pos_idx, meas_pos_idx+1], np.int)
if self.direction<0:
prope_idx = np.flipud(prope_idx)
u_prope_pos = prop_lines[prope_idx]
self.U_filenames = []
self.U_delta = np.diff(u_prope_pos)
suffix = ['A', 'B', 'C']
for n in range(len(prope_idx)):
u_start = 0.5*(self.start+self.stop)
u_stop = 0.5*(self.start+self.stop)
u_start[self.prop_ny] = u_prope_pos[n]
u_stop[self.prop_ny] = u_prope_pos[n]
u_start[self.exc_ny] = self.start[self.exc_ny]
u_stop[self.exc_ny] = self.stop [self.exc_ny]
u_name = self.lbl_temp.format('ut') + suffix[n]
self.U_filenames.append(u_name)
u_probe = CSX.AddProbe(u_name, p_type=0, weight=self.upside_down)
u_probe.AddBox(u_start, u_stop)
i_prope_pos = u_prope_pos[0:2] + np.diff(u_prope_pos)/2.0
self.I_filenames = []
self.I_delta = np.diff(i_prope_pos)
i_start = np.array(self.start)
i_stop = np.array(self.stop)
i_stop[self.exc_ny] = self.start[self.exc_ny]
for n in range(len(i_prope_pos)):
i_start[self.prop_ny] = i_prope_pos[n]
i_stop[self.prop_ny] = i_prope_pos[n]
i_name = self.lbl_temp.format('it') + suffix[n]
self.I_filenames.append(i_name)
i_probe = CSX.AddProbe(i_name, p_type=1, weight=self.direction, norm_dir=self.prop_ny)
i_probe.AddBox(i_start, i_stop)
if excite!=0:
excide_pos_idx = np.argmin(np.abs(prop_lines-(self.start[self.prop_ny] + self.feed_shift*self.direction)))
exc_start = np.array(self.start)
exc_stop = np.array(self.stop)
exc_start[self.prop_ny] = prop_lines[excide_pos_idx]
exc_stop [self.prop_ny] = prop_lines[excide_pos_idx]
exc_vec = np.zeros(3)
exc_vec[self.exc_ny] = -1*self.upside_down*excite
exc = CSX.AddExcitation(self.lbl_temp.format('excite'), exc_type=0, exc_val=exc_vec, delay=self.delay)
exc.AddBox(exc_start, exc_stop, priority=self.priority)
if self.feed_R>=0 and not np.isinf(self.feed_R):
R_start = np.array(self.start)
R_stop = np.array(self.stop)
R_stop [self.prop_ny] = R_start[self.prop_ny]
if self.feed_R==0:
metal_prop.AddBox(R_start, R_stop)
else:
lumped_R = CSX.AddLumpedElement(self.lbl_temp.format('resist'), ny=self.exc_ny, caps=True, R=self.feed_R)
lumped_R.AddBox(R_start, R_stop)
def ReadUIData(self, sim_path, freq, signal_type ='pulse'):
self.u_data = UI_data(self.U_filenames, sim_path, freq, signal_type )
self.uf_tot = self.u_data.ui_f_val[1]
self.i_data = UI_data(self.I_filenames, sim_path, freq, signal_type )
self.if_tot = 0.5*(self.i_data.ui_f_val[0]+self.i_data.ui_f_val[1])
unit = self.CSX.GetGrid().GetDeltaUnit()
Et = self.u_data.ui_f_val[1]
dEt = (self.u_data.ui_f_val[2] - self.u_data.ui_f_val[0]) / (np.sum(np.abs(self.U_delta)) * unit)
Ht = self.if_tot # space averaging: Ht is now defined at the same pos as Et
dHt = (self.i_data.ui_f_val[1] - self.i_data.ui_f_val[0]) / (np.abs(self.I_delta[0]) * unit)
beta = np.sqrt( - dEt * dHt / (Ht * Et) )
beta[np.real(beta) < 0] *= -1 # determine correct sign (unlike the paper)
self.beta = beta
# determine ZL
self.Z_ref = np.sqrt(Et * dEt / (Ht * dHt))
class WaveguidePort(Port):
"""
Base class for any waveguide port.
See Also
--------
Port, RectWGPort
"""
def __init__(self, CSX, port_nr, start, stop, exc_dir, E_WG_func, H_WG_func, kc, excite=0, **kw):
super(WaveguidePort, self).__init__(CSX, port_nr=port_nr, start=start, stop=stop, excite=excite, **kw)
self.exc_ny = CheckNyDir(exc_dir)
self.ny_P = (self.exc_ny+1)%3
self.ny_PP = (self.exc_ny+2)%3
self.direction = np.sign(stop[self.exc_ny]-start[self.exc_ny])
self.ref_index = 1
assert not (self.excite!=0 and stop[self.exc_ny]==start[self.exc_ny]), 'port length in excitation direction may not be zero if port is excited!'
self.kc = kc
self.E_func = E_WG_func
self.H_func = H_WG_func
if excite!=0:
e_start = np.array(start)
e_stop = np.array(stop)
e_stop[self.exc_ny] = e_start[self.exc_ny]
e_vec = np.ones(3)
e_vec[self.exc_ny]=0
exc = CSX.AddExcitation(self.lbl_temp.format('excite'), exc_type=0, exc_val=e_vec, delay=self.delay)
exc.SetWeightFunction([str(x) for x in self.E_func])
exc.AddBox(e_start, e_stop, priority=self.priority)
# voltage/current planes
m_start = np.array(start)
m_stop = np.array(stop)
m_start[self.exc_ny] = m_stop[self.exc_ny]
self.measplane_shift = np.abs(stop[self.exc_ny] - start[self.exc_ny])
self.U_filenames = [self.lbl_temp.format('ut'), ]
u_probe = CSX.AddProbe(self.U_filenames[0], p_type=10, mode_function=self.E_func)
u_probe.AddBox(m_start, m_stop)
self.I_filenames = [self.lbl_temp.format('it'), ]
i_probe = CSX.AddProbe(self.I_filenames[0], p_type=11, weight=self.direction, mode_function=self.H_func)
i_probe.AddBox(m_start, m_stop)
def CalcPort(self, sim_path, freq, ref_impedance=None, ref_plane_shift=None, signal_type='pulse'):
k = 2.0*np.pi*freq/C0*self.ref_index
self.beta = np.sqrt(k**2 - self.kc**2)
self.ZL = k * Z0 / self.beta #analytic waveguide impedance
if ref_impedance is None:
self.Z_ref = self.ZL
super(WaveguidePort, self).CalcPort(sim_path, freq, ref_impedance, ref_plane_shift, signal_type)
class RectWGPort(WaveguidePort):
"""
Rectangular waveguide port.
:param a,b: float -- Width/Height of rectangular waveguide port
See Also
--------
Port, WaveguidePort
"""
def __init__(self, CSX, port_nr, start, stop, exc_dir, a, b, mode_name, excite=0, **kw):
Port.__init__(self, CSX, port_nr, start, stop, excite=0, **kw)
self.exc_ny = CheckNyDir(exc_dir)
self.ny_P = (self.exc_ny+1)%3
self.ny_PP = (self.exc_ny+2)%3
self.WG_size = [a, b]
self.WG_mode = mode_name
assert len(self.WG_mode)==4, 'Invalid mode definition'
self.unit = self.CSX.GetGrid().GetDeltaUnit()
if self.WG_mode.startswith('TE'):
self.TE = True
self.TM = False
else:
self.TE = False
self.TM = True
self.M = float(self.WG_mode[2])
self.N = float(self.WG_mode[3])
assert self.TE, 'Currently only TE-modes are supported! Mode found: {}'.format(self.WG_mode)
# values by David M. Pozar, Microwave Engineering, third edition
a = self.WG_size[0]
b = self.WG_size[1]
xyz = 'xyz'
if self.start[self.ny_P]!=0:
name_P = '({}-{})'.format(xyz[self.ny_P], self.start[self.ny_P])
else:
name_P = xyz[self.ny_P]
if self.start[self.ny_PP]!=0:
name_PP = '({}-{})'.format(xyz[self.ny_P], self.start[self.ny_P])
else:
name_PP = xyz[self.ny_P]
kc = np.sqrt((self.M*np.pi/a)**2 + (self.N*np.pi/b)**2)
a /= self.unit
b /= self.unit
E_func = [0,0,0]
H_func = [0,0,0]
if self.N>0:
E_func[self.ny_P] = '{}*cos({}*{})*sin({}*{})'.format(self.N/b , self.M*np.pi/a, name_P, self.N*np.pi/b, name_PP)
if self.M>0:
E_func[self.ny_PP] = '{}*sin({}*{})*cos({}*{})'.format(-1*self.M/a, self.M*np.pi/a, name_P, self.N*np.pi/b, name_PP)
if self.M>0:
H_func[self.ny_P] = '{}*sin({}*{})*cos({}*{})'.format(self.M/a, self.M*np.pi/a, name_P, self.N*np.pi/b, name_PP)
if self.N>0:
H_func[self.ny_PP] = '{}*cos({}*{})*sin({}*{})'.format(self.N/b, self.M*np.pi/a, name_P, self.N*np.pi/b, name_PP)
super(RectWGPort, self).__init__(CSX, port_nr=port_nr, start=start, stop=stop, exc_dir=exc_dir, E_WG_func=E_func, H_WG_func=H_func, kc=kc, excite=excite, **kw)

66
python/openEMS/utilities.py Normal file
View File

@ -0,0 +1,66 @@
# -*- coding: utf-8 -*-
#
# Copyright (C) 2015,20016 Thorsten Liebig (Thorsten.Liebig@gmx.de)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
import numpy as np
def DFT_time2freq( t, val, freq, signal_type='pulse'):
assert len(t)==len(val)
assert len(freq)>0
f_val = np.zeros(len(freq))*1j
for n_f in range(len(freq)):
f_val[n_f] = np.sum( val*np.exp( -1j*2*np.pi*freq[n_f] * t ) )
if signal_type == 'pulse':
f_val *= t[1]-t[0]
elif signal_type == 'periodic':
f_val /= len(t)
else:
raise Exception('Unknown signal type: "{}"'.format(signal_type))
return 2*f_val # single-sided spectrum
def Check_Array_Equal(a,b, tol, relative=False):
a = np.array(a)
b = np.array(b)
if a.shape!=b.shape:
return False
if tol==0:
return (a==b).all()
if relative:
d = np.abs((a-b)/a)
else:
d = np.abs((a-b))
return np.max(d)<tol
if __name__=="__main__":
import pylab as plt
t = np.linspace(0,2,201)
s = np.sin(2*np.pi*2*t)
plt.plot(t,s)
f = np.linspace(0,3,101)
sf = DFT_time2freq(t, s, f, 'periodic')
plt.figure()
plt.plot(f, np.abs(sf))
plt.show()

47
python/setup.py Normal file
View File

@ -0,0 +1,47 @@
# -*- coding: utf-8 -*-
"""
Created on Sun Dec 13 23:48:22 2015
@author: thorsten
"""
from distutils.core import setup
from distutils.extension import Extension
from Cython.Build import cythonize
import os, sys
ROOT_DIR = os.path.dirname(__file__)
sys.path.append(os.path.join(ROOT_DIR,'..','..','CSXCAD','python'))
extensions = [
Extension("*", [os.path.join(os.path.dirname(__file__), "openEMS","*.pyx")],
language="c++", # generate C++ code
libraries = ['CSXCAD','openEMS', 'nf2ff']),
]
setup(
name="openEMS",
version = '0.0.33',
description = "Python interface for the openEMS FDTD library",
classifiers = [
'Development Status :: 3 - Alpha',
'Intended Audience :: Developers',
'Intended Audience :: Information Technology',
'Intended Audience :: Science/Research',
'License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)',
'Programming Language :: Python',
'Topic :: Scientific/Engineering',
'Topic :: Software Development :: Libraries :: Python Modules',
'Operating System :: POSIX :: Linux',
'Operating System :: Microsoft :: Windows',
],
author = 'Thorsten Liebig',
author_email = 'Thorsten.Liebig@gmx.de',
maintainer = 'Thorsten Liebig',
maintainer_email = 'Thorsten.Liebig@gmx.de',
url = 'http://openEMS.de',
packages=["openEMS", ],
package_data={'openEMS': ['*.pxd']},
ext_modules = cythonize(extensions)
)