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

f0 = 446.1e6
fmin = 200e6
fmax = 600e6

# NOTES:
# 1) we must increase the frequency binning since we want to hit a fixed point at 100kHz resolution
# 2) it seems 2x, that is 50kHz binning removes most artifacts, which makes sense considering likely rounding
#    happening within simulation/frequency-binning)
# 3) I'm not sure about the +1. my guess would be the binning for some reason wants a symmetric array with
#    a center point and 2 array halves left/right to center point of equal length. is that correct?
#
numfreq = int((fmax - fmin) / 1e5 * 2) + 1  # number of frequency points in output file
assert numfreq == 8001

lambda_f0 = C0 / f0 / unit / 2
lambda_f_min = C0 / fmin / unit
lambda_f_max = C0 / fmax / unit

# NOTES:
# 1) hitting a fixed point is incredibly sensitive to changes in many aspects,
#    most importantly of course the dipole length corrected for reactance of finite conductor (dipole arms)
#
# Found resonance frequency at 446.1 MHz with -73.3 dB at 70.1 Ohm for dipole length 299.4624 mm (lambda/2 336.0148598968841 mm)
#
# QUESTIONS:
# 1) why does the simulation produce _different_ results when re-rerunning without changing _anything_? this is
#    highly irritating to me, it doesn't increase my trust in the results, does it mean results are
#    a question of good luck? in my eyes, being able to produce repeatable results is absolutely essential
#    for any simulation, or in general, any experiment? I thought that would be "scientific" rather than "gambling"?
#
# Found resonance frequency at 446.0 MHz with -69.8 dB at 70.1 Ohm for dipole length 299.4624 mm (lambda/2 336.0148598968841 mm)
# Found resonance frequency at 446.1 MHz with -67.8 dB at 70.2 Ohm for dipole length 299.4624 mm (lambda/2 336.0148598968841 mm)
# Found resonance frequency at 446.3 MHz with -64.9 dB at 70.2 Ohm for dipole length 299.4624 mm (lambda/2 336.0148598968841 mm)
# Found resonance frequency at 446.1 MHz with -68.6 dB at 70.2 Ohm for dipole length 299.4624 mm (lambda/2 336.0148598968841 mm)
# Found resonance frequency at 446.1 MHz with -73.3 dB at 70.1 Ohm for dipole length 299.4624 mm (lambda/2 336.0148598968841 mm)
# Found resonance frequency at 446.1 MHz with -73.3 dB at 70.1 Ohm for dipole length 299.4624 mm (lambda/2 336.0148598968841 mm)
# Found resonance frequency at 446.0 MHz with -70.5 dB at 70.1 Ohm for dipole length 299.4624 mm (lambda/2 336.0148598968841 mm)
#
#    Note how in above results I absolutely changed nothing and merely re-run the simulation. both the resonance
#    frequency and the impedance, as well as S11 changes in the results.
#
#    wtf? I want to optimize for a _fixed_ point, and I want results that are repeatable (identical) regardless of
#    when/how often I run. how do I ensure that? I openEMS even designed to allow that? are other EM simulation tools
#    also exhibiting such non-determinism? I'm deeply puzzled:(
#
# 2) should the whole EM simulation be run multiple times without changing anything (no geo, no mesh, nothing), and
#    then average/min/max the _results_ to derive final actual, robust results? IOW: "sample whole simulation" would
#    be like "run experiment", and the do some basic statistics over those runs? could be done of course, I just need
#    to know if that is what I am _supposed_ to do - or if I goofed up sth else and the former would be pointless;)
#
dipole_length = 299.4624
dipole_gap = 1.0
dipole_wire_radius = 2.0

feed_resistance = 70.1
feed_radius = dipole_wire_radius / math.sqrt(2)

# QUESTIONS:
# 1) why is that using lambda_f_max, not lambda_f_min?
# 2) why is that setting 2x size for Z-axis?
#
sim_box = (
    np.array([1, 1, 2]) * 0.5 * lambda_f_max
)  # half wavelength distance to boundary is fine for S-parameter simulation
max_res = math.floor(lambda_f_max / 20)
assert max_res == 24

fdtd = openEMS(NrTS=100000, EndCriteria=1e-4)

# QUESTIONS:
# 1) will this work when using this with an excitation frequency range starting at DC (fmin == 0Hz)?
#
fdtd.SetGaussExcite(
    (fmin + fmax) / 2, (fmax - fmin) / 2
)  # we want the wide spectrum spread of the GAUSSIAN pulse for S-parameter calculation!

# QUESTIONS:
# 1) why is PML_8 chosen vs MUR, and what is that anyway?
#
fdtd.SetBoundaryCond(["PML_8"] * 6)
# fdtd.SetBoundaryCond(["MUR"] * 6)

csx = ContinuousStructure()
fdtd.SetCSX(csx)

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

mesh.AddLine("z", np.linspace(-dipole_gap / 2, dipole_gap / 2, 5))
mesh.AddLine(
    "z", np.linspace(-dipole_length / 2 - 5 * dipole_wire_radius, -dipole_length / 2 + 5 * dipole_wire_radius, 11)
)
mesh.AddLine(
    "z", np.linspace(dipole_length / 2 - 5 * dipole_wire_radius, dipole_length / 2 + 5 * dipole_wire_radius, 11)
)

# QUESTIONS:
# 1) apparently, fine meshing in X-axis (and Y-axis couple of lines below) is REQUIRED (that was my main mistake),
#    however:
#
#           24 == max_res > dipole_wire_radius == 2 (!)
#
#    so how does this mesh up the dipole wire anyways? the mesh resolution applied is 10x larger than the
#    geometry of the object being meshed. I am confused;)
#
assert max_res == 24 and dipole_wire_radius == 2 and max_res > dipole_wire_radius
mesh.AddLine("x", np.linspace(-dipole_wire_radius / 2, -dipole_wire_radius / 2, 5))
mesh.AddLine("x", [-sim_box[0] / 2, 0, sim_box[0] / 2])
mesh.SmoothMeshLines("x", max_res, ratio=1.4)

mesh.AddLine("y", np.linspace(-dipole_wire_radius / 2, -dipole_wire_radius / 2, 5))
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([[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))

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, 1),
            round(s11_dB[idx][0], 1),
            round(np.real(Zin[idx])[0], 1),
            dipole_length,
            lambda_f0,
        )
    )

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