openEMS/python/Tutorials/CRLH_Extraction.py

# -*- coding: utf-8 -*-
"""
 Tutorials / CRLH_Extraction

 Description 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()