minimal_dipole3.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

# main frequency of interest
f0 = 446.1e6

# excitation frequency range
# fmin, fmax = (400e6, 500e6)
# fmin, fmax = (200e6, 600e6)
fmin, fmax = (0, 1000e6)
# fmin, fmax = (0, round(f0 * 7 * 1.15, -2))  # 3.6GHz
# fmin, fmax = (0, round(f0 * 8 * 1.15, -2))  # 4.1GHz
# fmin, fmax = (0, 7000e6)  # >15x f0 (to cover 3x, 5x, 7x, 9x, 11x, 13x harmonics)
# fmin, fmax = (0, 10000e6)

numfreq = int((fmax - fmin) / 1e5 * 2) + 1  # number of frequency points in output file
print("f0={}, fmin={}, fmax={}, numfreq={}".format(f0, fmin, fmax, numfreq))

lambda_f0 = C0 / f0 / unit
lambda_f_max = C0 / fmax / unit
print("lambda_f0={}, lambda_f_max={}".format(lambda_f0, lambda_f_max))

dipole_length = 306
# dipole_wire_radius = (dipole_length / 10 ** 2) / (2 * math.pi)
dipole_wire_radius = 1.5

# feed_resistance = 73.1
feed_resistance = 71.8
feed_resistance = 70.7

# this makes the rectangular lumped port surface fully contained within the circular dipole arm,
# but it leads to differently pitched mesh lines between dipole arm and lumped port!
# feed_radius = dipole_wire_radius / math.sqrt(2)

# this makes the circular dipole arm surface fully contained within the rectangular lumped port!
# it also leads to _same_ mesh pitch (XY plane) between dipole and lumped port!
feed_radius = dipole_wire_radius

# this results in a cube-shaped lumped port
dipole_gap = feed_radius * 2

print(
    "lambda_f0/2={}, dipole_length={}, dipole_gap={}, dipole_wire_radius={}, feed_resistance={}, feed_radius={}".format(
        lambda_f0 / 2, dipole_length, dipole_gap, dipole_wire_radius, feed_resistance, feed_radius))

sim_box = (
    np.array([1, 1, 2]) * 0.5 * lambda_f0
)
max_res = math.floor(lambda_f_max / 20)

fdtd = openEMS(NrTS=100000, EndCriteria=1e-4)  # NORMAL
# fdtd = openEMS(NrTS=100000, EndCriteria=1e-5)  # SLOWER + MORE PRECISE
# fdtd = openEMS(NrTS=100000, EndCriteria=1e-6)  # VERY SLOW + HIGHLY PRECISE

fdtd.SetGaussExcite(
    (fmin + fmax) / 2, (fmax - fmin) / 2
)

fdtd.SetBoundaryCond(["PML_8"] * 6)
# fdtd.SetBoundaryCond(["MUR"] * 6)

csx = ContinuousStructure()
fdtd.SetCSX(csx)

mesh = csx.GetGrid()
mesh.SetDeltaUnit(unit)

# setup mesh
#
k = 2  # use 2, 3, .. for finer mesh
r = 1.4
n = (2 ** k + 1) * 2
a = (n - 1) / (n / 2)
print("k={}: n={}, a={}; r={}".format(k, n, a, r))

if True:
    mesh.AddLine("x", [-sim_box[0] / 2, sim_box[0] / 2])
    mesh.AddLine("y", [-sim_box[1] / 2, sim_box[1] / 2])
    mesh.AddLine("z", [-sim_box[2] / 2, sim_box[2] / 2])

if True:
    mesh.AddLine("x", np.linspace(-feed_radius * a, feed_radius * a, n))
    mesh.AddLine("y", np.linspace(-feed_radius * a, feed_radius * a, n))
    mesh.AddLine("z", np.linspace(-feed_radius * a, feed_radius * a, n))

if True:
    mesh.AddLine(
        "z", np.linspace(-dipole_length / 2 - dipole_wire_radius * a, -dipole_length / 2 + dipole_wire_radius * a, n)
    )
    mesh.AddLine(
        "z", np.linspace(dipole_length / 2 - dipole_wire_radius * a, dipole_length / 2 + dipole_wire_radius * a, n)
    )

if True:
    mesh.SmoothMeshLines("x", max_res, ratio=r)
    mesh.SmoothMeshLines("y", max_res, ratio=r)
    mesh.SmoothMeshLines("z", max_res, ratio=r)

mesh_size = mesh.GetQtyLines("x") * mesh.GetQtyLines("y") * mesh.GetQtyLines("z")

print("Created mesh of size {} (x={}, y={}, z={})".format(mesh_size,
                                                          mesh.GetQtyLines("x"),
                                                          mesh.GetQtyLines("y"),
                                                          mesh.GetQtyLines("z")))
arm1: CSPropMetal = csx.AddMetal("arm1")
arm1.AddWire([[0, 0], [0, 0], [-dipole_gap / 2, -dipole_length / 2]], radius=dipole_wire_radius)
arm1.SetColor("#ff0000", 50)

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

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

output_dir = os.path.join(tempfile.gettempdir(), "dipole")
if not os.path.isdir(output_dir):
    os.mkdir(output_dir)
output_fn = os.path.join(output_dir, "dipole.xml")
csx.Write2XML(output_fn)
os.system('{} "{}"'.format(AppCSXCAD_BIN, output_fn))

# sys.exit(1)

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

print("Computing S11 for {} frequencies".format(numfreq))
freq = np.linspace(fmin, fmax, numfreq)
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.0) & (s11_dB == np.min(s11_dB)))[0]
if not len(idx) == 1:
    print("No resonance frequency found for far-field calculation!")
else:
    print(
        "Found resonance frequency at {} MHz with {} dB at {} Ohm for dipole length {} mm (lambda/2 {} mm)".format(
            round(freq[idx][0] / 1e6, 2),
            round(s11_dB[idx][0], 2),
            round(np.real(Zin[idx])[0], 2),
            round(dipole_length, 2),
            round(lambda_f0 / 2, 2),
        )
    )

    # create Touchstone S1P output file
    s1p_name = CSX_file = os.path.join(output_dir, "dipole.s1p")
    s1p_file = open(s1p_name, "w")
    s1p_file.write("#   Hz   S  RI   R   " + str(feed_resistance) + "\n")
    s1p_file.write("!\n")
    for index in range(0, numfreq):
        f = freq[index]
        re = np.real(s11[index])
        im = np.imag(s11[index])
        s1p_file.write(str(f) + " " + str(re) + " " + str(im) + "\n")
    s1p_file.close()
    print("Wrote Touchstone S1P output file to {}".format(s1p_name))

    # create matplotlib plot
    pyplot.figure()
    pyplot.plot(freq / 1e6, s11_dB, "k-", linewidth=2, label="$S_{11}$")
    pyplot.grid()
    pyplot.title(
        "Center-fed Lambda/2 Dipole for {} MHz\nAntenna Efficiency: S11 Reflection Coefficient vs Frequency\n".format(
            round(f0 / 1e6, 1)
        )
    )
    pyplot.ylabel("Reflection coefficient $S_{11}$ (dB)")
    pyplot.xlabel("Frequency (MHz)")
    pyplot.show()