1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 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()
|