dipole.py

import os
import math
import tempfile

import numpy as np
from matplotlib import pyplot

from CSXCAD import ContinuousStructure, AppCSXCAD_BIN
from CSXCAD.CSProperties import CSPropMetal

from openEMS import openEMS
from openEMS.physical_constants import C0

unit = 1e-3

# 500MHz +/- 100MHz
f0 = 500e6
fc = 100e6
lambda0 = round(C0 / f0 / unit)

# lambda/2 dipole
dipole_length = lambda0 / 2
dipole_gap = 0.0
dipole_wire_radius = 1.0

# dipole feed
feed_resistance = 73.1
feed_radius = dipole_wire_radius / math.sqrt(2)
feed_overlap = 1.0

mesh_res_dipole = 1.0
max_res = math.floor(C0 / (f0 + fc) / unit / 20)
sim_box = np.array([1, 1, 1]) * 2.0 * lambda0
# nf_ff_transition_distance = math.ceil(lambda0 / (2 * math.pi))
nf_ff_transition_distance = 2 * lambda0

output_dir = os.path.join(tempfile.gettempdir(), "dipole")
if not os.path.isdir(output_dir):
    os.mkdir(output_dir)

fdtd = openEMS(NrTS=50000, EndCriteria=1e-4)
fdtd.SetGaussExcite(f0, fc)
fdtd.SetBoundaryCond(["MUR", "MUR", "MUR", "MUR", "MUR", "MUR"])

csx = ContinuousStructure()
fdtd.SetCSX(csx)
mesh = csx.GetGrid()
mesh.SetDeltaUnit(unit)

mesh.AddLine("x", [-dipole_gap / 2 - dipole_length / 2, dipole_gap / 2 + dipole_length / 2])
mesh.SmoothMeshLines("x", mesh_res_dipole)

mesh.AddLine("x", [-sim_box[0] / 2, 0, sim_box[0] / 2])
mesh.SmoothMeshLines("x", max_res, ratio=1.4)

mesh.AddLine("y", [-sim_box[1] / 2, 0, sim_box[1] / 2])
mesh.SmoothMeshLines("y", max_res, ratio=1.4)

mesh.AddLine("z", [-sim_box[2] / 2, 0, sim_box[2] / 2])
mesh.SmoothMeshLines("z", max_res, ratio=1.4)

arm1: CSPropMetal = csx.AddMetal("arm1")
arm1.AddWire([
    [- dipole_gap / 2, -dipole_gap / 2 - dipole_length / 2],
    [0, 0],
    [0, 0]
], radius=dipole_wire_radius)
arm1.SetColor("#ff0000", 50)

arm2: CSPropMetal = csx.AddMetal("arm2")
arm2.AddWire([
    [dipole_gap / 2, dipole_gap / 2 + dipole_length / 2],
    [0, 0],
    [0, 0]
], radius=dipole_wire_radius)
arm2.SetColor("#ff0000", 50)

feed = fdtd.AddLumpedPort(
    1,
    feed_resistance,
    [- dipole_gap / 2 - feed_overlap, -feed_radius, -feed_radius],
    [dipole_gap / 2 + feed_overlap, feed_radius, feed_radius],
    "x",
    1.0,
    priority=5,
)

if False:
    output_fn = os.path.join(output_dir, "dipole.xml")
    csx.Write2XML(output_fn)
    os.system('{} "{}"'.format(AppCSXCAD_BIN, output_fn))

# mesh.SmoothMeshLines("all", max_res, 1.4)

fdtd.Run(output_dir, verbose=3, cleanup=True)

freq = np.linspace(f0 - fc, f0 + fc, 501)
feed.CalcPort(output_dir, freq)

Zin = feed.uf_tot / feed.if_tot
s11 = feed.uf_ref / feed.uf_inc
s11_dB = 20.0 * np.log10(np.abs(s11))

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 calculation!")
else:
    print("\n")
    print(
        "Found resonance frequency at {} MHz with {} dB at {} Ohm".format(
            round(freq[idx][0] / 1e6),
            round(s11_dB[idx][0]),
            round(np.real(Zin[idx])[0]),
        )
    )
    print("\n")

# plot the feed point impedance
pyplot.figure()
pyplot.plot(freq / 1e6, np.real(Zin), "k-", linewidth=2, label=r"$\Re(Z_{in})$")
pyplot.grid()
pyplot.plot(freq / 1e6, np.imag(Zin), "r--", linewidth=2, label=r"$\Im(Z_{in})$")
pyplot.title("feed point impedance")
pyplot.xlabel("frequency (MHz)")
pyplot.ylabel("impedance (Omega)")
pyplot.legend()

# plot reflection coefficient S11
pyplot.figure()
pyplot.plot(freq / 1e6, s11_dB, "k-", linewidth=2, label="$S_{11}$")
pyplot.grid()
pyplot.title("reflection coefficient $S_{11}$")
pyplot.ylabel("S-Parameter (dB)")
pyplot.xlabel("Frequency (MHz)")
pyplot.legend()

pyplot.show()