Import GIT HEAD of openEMS sub-project (with Python interface)

1 files changed, 239 insertions, 0 deletions

diff --git a/openEMS/python/Tutorials/CRLH_Extraction.py b/openEMS/python/Tutorials/CRLH_Extraction.py
new file mode 100644
index 0000000..4e1850e
--- /dev/null
+++ b/openEMS/python/Tutorials/CRLH_Extraction.py
@@ -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()
\ No newline at end of file