python: add tutorials (ported from Matlab)

Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>
This commit is contained in:
Thorsten Liebig 2016-08-28 21:42:00 +02:00
parent 8a2eb2df33
commit 956125e360
4 changed files with 573 additions and 0 deletions

View File

@ -0,0 +1,191 @@
# -*- coding: utf-8 -*-
"""
Tutorials / helical antenna
Tested with
- python 3.4
- openEMS v0.0.33+
(C) 2015-2016 Thorsten Liebig <thorsten.liebig@gmx.de>
"""
import os, tempfile
from pylab import *
from CSXCAD import CSXCAD
from openEMS.openEMS import openEMS
from openEMS.physical_constants import *
Sim_Path = os.path.join(tempfile.gettempdir(), 'Helical_Ant')
post_proc_only = False
## setup the simulation
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 CSXCAD 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 helix using the wire primitive
helix_metal = CSX.AddMetal('helix' ) # create a perfect electric conductor (PEC)
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])
CSX.AddCurve(helix_metal, 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]
CSX.AddCylinder(gnd, 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)
if 0: # debugging only
CSX_file = os.path.join(Sim_Path, 'helix.xml')
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 & do the plots
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 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}|$' )
## NFFF contour plots ####################################################
## 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.CalcNF2FF(Sim_Path, f0, theta, phi, read_cached=True, verbose=True )
#
Dmax_dB = 10*log10(nf2ff.Dmax[0])
E_norm = 20.0*log10(nf2ff.E_norm[0]/np.max(nf2ff.E_norm[0])) + 10*log10(nf2ff.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.Prad[0]))
print('directivity: Dmax = {} dBi'.format(Dmax_dB))
print('efficiency: nu_rad = {} %'.format(100*nf2ff.Prad[0]/interp(f0, freq, port.P_acc)))
print('theta_HPBW = {} °'.format(theta_HPBW))
##
E_norm = 20.0*log10(nf2ff.E_norm[0]/np.max(nf2ff.E_norm[0])) + 10*log10(nf2ff.Dmax[0])
E_CPRH = 20.0*log10(np.abs(nf2ff.E_cprh[0])/np.max(nf2ff.E_norm[0])) + 10*log10(nf2ff.Dmax[0])
E_CPLH = 20.0*log10(np.abs(nf2ff.E_cplh[0])/np.max(nf2ff.E_norm[0])) + 10*log10(nf2ff.Dmax[0])
##
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.freq[0]/1e9))
legend()
### dump to vtk
# TODO
show()

View File

@ -0,0 +1,114 @@
# -*- coding: utf-8 -*-
"""
Tutorials / MSL_NotchFilter
Tested with
- python 3.4
- openEMS v0.0.33+
(C) 2016 Thorsten Liebig <thorsten.liebig@gmx.de>
"""
import os, tempfile
from pylab import *
from CSXCAD import CSXCAD
from openEMS import openEMS
from openEMS.physical_constants import *
post_proc_only = False
Sim_Path = os.path.join(tempfile.gettempdir(), 'NotchFilter')
## setup the simulation ###################################################
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.openEMS()
FDTD.SetGaussExcite( f_max/2, f_max/2 )
FDTD.SetBoundaryCond( ['PML_8', 'PML_8', 'MUR', 'MUR', 'PEC', 'MUR'] )
## setup CSXCAD geometry & mesh ###########################################
CSX = CSXCAD.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
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)
## 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]
CSX.AddBox( substrate, start, stop )
## MSL port
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
start = [-MSL_width/2, MSL_width/2, substrate_thickness]
stop = [ MSL_width/2, MSL_width/2+stub_length, substrate_thickness]
CSX.AddBox( pec, start, stop, priority=10 )
if 0: # debugging only
CSX_file = os.path.join(Sim_Path, 'notch.xml')
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( 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()

View File

@ -0,0 +1,124 @@
# -*- coding: utf-8 -*-
"""
Created on Tue Dec 29 15:40:23 2015
Tutorials / Rect_Waveguide
Describtion at:
http://openems.de/index.php/Tutorial:_Rectangular_Waveguide
Tested with
- python 3.4
- openEMS v0.0.33+
(C) 2015-2016 Thorsten Liebig <thorsten.liebig@gmx.de>
"""
import os, tempfile
from pylab import *
from CSXCAD import CSXCAD
from openEMS.openEMS import openEMS
from openEMS.physical_constants import *
Sim_Path = os.path.join(tempfile.gettempdir(), 'Rect_WG')
## setup the simulation ###################################################
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 CSXCAD geometry & mesh
CSX = CSXCAD.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];
CSX.AddBox(Et, start, stop);
## Write openEMS compatoble xml-file ######################################
if 0: # debugging only
CSX_file = os.path.join(Sim_Path, 'rect_wg.xml')
CSX.Write2XML(CSX_file)
os.system(r'AppCSXCAD "{}"'.format(CSX_file))
if not post_proc_only:
FDTD.Run(Sim_Path, verbose=3, cleanup=True)
### postproc ###############################################################
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()

View File

@ -0,0 +1,144 @@
# -*- coding: utf-8 -*-
"""
Created on Fri Dec 18 20:56:53 2015
@author: thorsten
"""
import os, tempfile
from pylab import *
from CSXCAD import CSXCAD
from openEMS.openEMS import openEMS
from openEMS.physical_constants import *
Sim_Path = os.path.join(tempfile.gettempdir(), 'Simp_Patch')
post_proc_only = True
# patch width in x-direction
patch_width = 32 # resonant length
# 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 = openEMS(NrTS=30000, EndCriteria=1e-4)
FDTD.SetGaussExcite( f0, fc )
FDTD.SetBoundaryCond( ['MUR', 'MUR', 'MUR', 'MUR', 'MUR', 'MUR'] )
mesh_res = C0/(f0+fc)/1e-3/20
CSX = CSXCAD.ContinuousStructure()
FDTD.SetCSX(CSX)
#initialize the mesh with the "air-box" dimensions
mesh = {}
mesh['x'] = [-SimBox[0]/2, feed_pos, SimBox[0]/2]
mesh['y'] = [-SimBox[1]/2, SimBox[1]/2]
mesh['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]
pb=CSX.AddBox(patch, priority=10, start=start, stop=stop) # add a box-primitive to the metal property 'patch'
edge_mesh = np.array([-1/3.0, 2/3.0]) * mesh_res/2
mesh['x'] = r_[mesh['x'], start[0]-edge_mesh, stop[0]+edge_mesh]
mesh['y'] = r_[mesh['y'], start[1]-edge_mesh, stop[1]+edge_mesh]
## create substrate
substrate = CSX.AddMaterial( 'substräte', Epsilon=substrate_epsR, Kappa=substrate_kappa)
start = [-substrate_width/2, -substrate_length/2, 0]
stop = [ substrate_width/2, substrate_length/2, substrate_thickness]
sb=CSX.AddBox( substrate, priority=0, start=start, stop=stop )
# add extra cells to discretize the substrate thickness
mesh['z'] = r_[mesh['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
gb=CSX.AddBox(gnd, start, stop, priority=10)
mesh['x'] = r_[mesh['x'], start[0], stop[0]]
mesh['y'] = r_[mesh['y'], start[1], stop[1]]
## 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)
mesh['x'] = r_[mesh['x'], start[0]]
mesh['y'] = r_[mesh['y'], start[1]]
CSX.DefineGrid(mesh, unit=1e-3, smooth_mesh_res=mesh_res)
nf2ff = FDTD.CreateNF2FFBox()
if 0: # debugging only
CSX_file = os.path.join(Sim_Path, 'simp_patch.xml')
CSX.Write2XML(CSX_file)
os.system(r'AppCSXCAD "{}"'.format(CSX_file))
if not post_proc_only:
FDTD.Run(Sim_Path, verbose=3, cleanup=True)
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)
grid()
ylabel('s11 (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.CalcNF2FF(Sim_Path, f_res, theta, phi, center=[0,0,1e-3], read_cached=False )
figure()
E_norm = 20.0*np.log10(nf2ff.E_norm[0]/np.max(nf2ff.E_norm[0])) + nf2ff.Dmax[0]
plot(theta, np.squeeze(E_norm[:,0]), label='xz-plane')
plot(theta, np.squeeze(E_norm[:,1]), 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))
plot(f/1e9, np.imag(Zin))
grid()
ylabel('Zin (Ohm)')
xlabel('frequency (GHz)')
show()