1# -*- coding: utf-8 -*-
2"""
3 Bent Patch Antenna Tutorial
4
5 Tested with
6  - python 3.4
7  - openEMS v0.0.33+
8
9 (C) 2016 Thorsten Liebig <thorsten.liebig@gmx.de>
10
11"""
12
13### Import Libraries
14import os, tempfile
15from pylab import *
16from mpl_toolkits.mplot3d import Axes3D
17
18from CSXCAD import CSXCAD
19
20from openEMS.openEMS import openEMS
21from openEMS.physical_constants import *
22
23
24### Setup the simulation
25Sim_Path = os.path.join(tempfile.gettempdir(), 'Bent_Patch')
26
27post_proc_only = False
28
29unit = 1e-3 # all length in mm
30
31f0 = 2.4e9 # center frequency, frequency of interest!
32lambda0 = round(C0/f0/unit) # wavelength in mm
33fc = 0.5e9 # 20 dB corner frequency
34
35# patch width in alpha-direction
36patch_width  = 32 # resonant length in alpha-direction
37patch_radius = 50 # radius
38patch_length = 40 # patch length in z-direction
39
40#substrate setup
41substrate_epsR   = 3.38
42substrate_kappa  = 1e-3 * 2*pi*2.45e9 * EPS0*substrate_epsR
43substrate_width  = 80
44substrate_length = 90
45substrate_thickness = 1.524
46substrate_cells = 4
47
48#setup feeding
49feed_pos   = -5.5  #feeding position in x-direction
50feed_width = 2     #feeding port width
51feed_R     = 50    #feed resistance
52
53# size of the simulation box
54SimBox_rad    = 2*100
55SimBox_height = 1.5*200
56
57### Setup FDTD parameter & excitation function
58FDTD = openEMS(CoordSystem=1, EndCriteria=1e-4) # init a cylindrical FDTD
59f0 = 2e9 # center frequency
60fc = 1e9 # 20 dB corner frequency
61FDTD.SetGaussExcite(f0, fc)
62FDTD.SetBoundaryCond(['MUR', 'MUR', 'MUR', 'MUR', 'MUR', 'MUR']) # boundary conditions
63
64### Setup the Geometry & Mesh
65# init a cylindrical mesh
66CSX = CSXCAD.ContinuousStructure(CoordSystem=1)
67FDTD.SetCSX(CSX)
68mesh = CSX.GetGrid()
69mesh.SetDeltaUnit(unit)
70
71### Setup the geometry using cylindrical coordinates
72# calculate some width as an angle in radiant
73patch_ang_width = patch_width/(patch_radius+substrate_thickness)
74substr_ang_width = substrate_width/patch_radius
75feed_angle = feed_pos/patch_radius
76
77# create patch
78patch = CSX.AddMetal('patch') # create a perfect electric conductor (PEC)
79start = [patch_radius+substrate_thickness, -patch_ang_width/2, -patch_length/2 ]
80stop  = [patch_radius+substrate_thickness,  patch_ang_width/2,  patch_length/2 ]
81patch.AddBox(priority=10, start=start, stop=stop) # add a box-primitive to the metal property 'patch'
82FDTD.AddEdges2Grid(dirs='all', properties=patch)
83
84# create substrate
85substrate = CSX.AddMaterial('substrate', epsilon=substrate_epsR, kappa=substrate_kappa  )
86start = [patch_radius                    , -substr_ang_width/2, -substrate_length/2]
87stop  = [patch_radius+substrate_thickness,  substr_ang_width/2,  substrate_length/2]
88substrate.AddBox(start=start, stop=stop)
89FDTD.AddEdges2Grid(dirs='all', properties=substrate)
90
91# save current density oon the patch
92jt_patch = CSX.AddDump('Jt_patch', dump_type=3, file_type=1)
93start = [patch_radius+substrate_thickness, -substr_ang_width/2, -substrate_length/2]
94stop  = [patch_radius+substrate_thickness, +substr_ang_width/2,  substrate_length/2]
95jt_patch.AddBox(start=start, stop=stop)
96
97# create ground
98gnd = CSX.AddMetal('gnd') # create a perfect electric conductor (PEC)
99start = [patch_radius, -substr_ang_width/2, -substrate_length/2]
100stop  = [patch_radius, +substr_ang_width/2, +substrate_length/2]
101gnd.AddBox(priority=10, start=start, stop=stop)
102FDTD.AddEdges2Grid(dirs='all', properties=gnd)
103
104# apply the excitation & resist as a current source
105start = [patch_radius                    ,  feed_angle, 0]
106stop  = [patch_radius+substrate_thickness,  feed_angle, 0]
107port = FDTD.AddLumpedPort(1 ,feed_R, start, stop, 'r', 1.0, priority=50, edges2grid='all')
108
109### Finalize the Mesh
110# add the simulation domain size
111mesh.AddLine('r', patch_radius+np.array([-20, SimBox_rad]))
112mesh.AddLine('a', [-0.75*pi, 0.75*pi])
113mesh.AddLine('z', [-SimBox_height/2, SimBox_height/2])
114
115# add some lines for the substrate
116mesh.AddLine('r', patch_radius+np.linspace(0,substrate_thickness,substrate_cells))
117
118# generate a smooth mesh with max. cell size: lambda_min / 20
119max_res = C0 / (f0+fc) / unit / 20
120max_ang = max_res/(SimBox_rad+patch_radius) # max res in radiant
121mesh.SmoothMeshLines(0, max_res, 1.4)
122mesh.SmoothMeshLines(1, max_ang, 1.4)
123mesh.SmoothMeshLines(2, max_res, 1.4)
124
125## Add the nf2ff recording box
126nf2ff = FDTD.CreateNF2FFBox()
127
128### Run the simulation
129if 0:  # debugging only
130    CSX_file = os.path.join(Sim_Path, 'bent_patch.xml')
131    if not os.path.exists(Sim_Path):
132        os.mkdir(Sim_Path)
133    CSX.Write2XML(CSX_file)
134    os.system(r'AppCSXCAD "{}"'.format(CSX_file))
135
136
137if not post_proc_only:
138    FDTD.Run(Sim_Path, verbose=3, cleanup=True)
139
140### Postprocessing & plotting
141f = np.linspace(max(1e9,f0-fc),f0+fc,401)
142port.CalcPort(Sim_Path, f)
143Zin = port.uf_tot / port.if_tot
144s11 = port.uf_ref/port.uf_inc
145s11_dB = 20.0*np.log10(np.abs(s11))
146
147figure()
148plot(f/1e9, s11_dB)
149grid()
150ylabel('s11 (dB)')
151xlabel('frequency (GHz)')
152
153P_in = 0.5*np.real(port.uf_tot * np.conj(port.if_tot)) # antenna feed power
154
155# plot feed point impedance
156figure()
157plot( f/1e6, real(Zin), 'k-', linewidth=2, label=r'$\Re(Z_{in})$' )
158grid()
159plot( f/1e6, imag(Zin), 'r--', linewidth=2, label=r'$\Im(Z_{in})$' )
160title( 'feed point impedance' )
161xlabel( 'frequency (MHz)' )
162ylabel( 'impedance ($\Omega$)' )
163legend( )
164
165
166idx = np.where((s11_dB<-10) & (s11_dB==np.min(s11_dB)))[0]
167if not len(idx)==1:
168    print('No resonance frequency found for far-field calulation')
169else:
170    f_res = f[idx[0]]
171    theta = np.arange(-180.0, 180.0, 2.0)
172    print("Calculate NF2FF")
173    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')
174
175    figure(figsize=(15, 7))
176    ax = subplot(121, polar=True)
177    E_norm = 20.0*np.log10(nf2ff_res_phi0.E_norm/np.max(nf2ff_res_phi0.E_norm)) + nf2ff_res_phi0.Dmax
178    ax.plot(np.deg2rad(theta), 10**(np.squeeze(E_norm)/20), linewidth=2, label='xz-plane')
179    ax.grid(True)
180    ax.set_xlabel('theta (deg)')
181    ax.set_theta_zero_location('N')
182    ax.set_theta_direction(-1)
183    ax.legend(loc=3)
184
185    phi = theta
186    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')
187
188    ax = subplot(122, polar=True)
189    E_norm = 20.0*np.log10(nf2ff_res_theta90.E_norm/np.max(nf2ff_res_theta90.E_norm)) + nf2ff_res_theta90.Dmax
190    ax.plot(np.deg2rad(phi), 10**(np.squeeze(E_norm)/20), linewidth=2, label='xy-plane')
191    ax.grid(True)
192    ax.set_xlabel('phi (deg)')
193    suptitle('Bent Patch Anteanna Pattern\nFrequency: {} GHz'.format(f_res/1e9), fontsize=14)
194    ax.legend(loc=3)
195
196    print( 'radiated power: Prad = {:.2e} Watt'.format(nf2ff_res_theta90.Prad[0]))
197    print( 'directivity:    Dmax = {:.1f} ({:.1f} dBi)'.format(nf2ff_res_theta90.Dmax[0], 10*np.log10(nf2ff_res_theta90.Dmax[0])))
198    print( 'efficiency:   nu_rad = {:.1f} %'.format(100*nf2ff_res_theta90.Prad[0]/real(P_in[idx[0]])))
199
200show()
201
202