Created
November 27, 2023 11:08
-
-
Save oberstet/ce57c9579591d6d1b1f9222afa00bccc to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import os | |
import math | |
import tempfile | |
from pprint import pprint, pformat | |
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 | |
# enable NF2FF recording, computation and plotting | |
enable_nf2ff = True | |
# fire up AppCSXCAD for viewing the model before running it | |
enable_appcsxcad = False | |
# all units are in mm | |
unit = 1e-3 | |
# excitation frequency and bandwidth: 500MHz +/- 100MHz | |
# f0 = 500e6 | |
f0 = 446.1e6 | |
# f0 = 111e6 | |
# f0 = 222e6 | |
# f0 = 333e6 | |
# f0 = 444e6 | |
# f0 = 555e6 | |
# f0 = 666e6 | |
# f0 = 777e6 | |
# f0 = 888e6 | |
# f0 = 999e6 | |
# excitation bandwidth | |
fc = 0.15 * f0 # +/- ~15% => ~30% BW total < 20% BW max. for center-fed dipole | |
# length factor to apply to reach fixed point of resonance frequency | |
# being identical to excitation frequency | |
# "Found resonance frequency at 500 MHz with -44 dB at 71 Ohm" | |
# Correction factor to shorten the wavelength used from theoretical value to account | |
# for XXX | |
# | |
# Found resonance frequency at 446.1 MHz with -80.4 dB at 71.0 Ohm | |
# Dipole (lambda/2) length is 290.7 mm | |
# | |
# opt_factor = 0.8625 | |
opt_factor = 0.8651 | |
# Due to end effects a finite thickness dipole is not resonant at a length of | |
# one-half wavelength 1 2 λ {\displaystyle \ {\tfrac {1}{2}}\lambda \ } but | |
# has inductive reactance. A typical thin dipole is actually resonant (has no reactance) | |
# at a slightly shorter length around 0.475 λ , {\displaystyle \ 0.475\lambda \ ,} at | |
# which its radiation resistance is about 67 Ohms. | |
# | |
# Wallace, Richard; Andreasson, Krister (2005). Introduction to RF and | |
# Microwave Passive Components. Artech House. p. 77. ISBN 9781630810092. | |
# wave length to compute antenna length from | |
# lambda0 = round(opt_factor * C0 / 500e6 / unit) | |
lambda0 = opt_factor * C0 / f0 / unit | |
# lambda/2 dipole | |
dipole_length = lambda0 / 2 | |
# gap in between the two dipole arms (the lumped port will fill that) | |
dipole_gap = 0.0 | |
dipole_wire_radius = 0.001 | |
# dipole_wire_radius = 1.0 | |
# Radiation resistance (ohms) of center-fed half-wave dipole | |
# https://en.wikipedia.org/wiki/Radiation_resistance#Radiation_resistance_of_common_antennas | |
# feed_resistance = 73.1 | |
feed_resistance = 71.0 | |
# Radius of lumped port (dipole feed port), set to be contained / enclosed | |
# completely by the dipole wire excited | |
feed_radius = dipole_wire_radius / math.sqrt(2) | |
# feed_radius = dipole_wire_radius | |
# Overlap of lumped port (dipole feed) with the actual dipole arms excited | |
# Note: MUST be non-zero, and actually >>0, not sure .. | |
feed_overlap = 0.5 | |
# 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=100000, 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, | |
) | |
######################################################################################### | |
# setup far-field recording | |
# | |
if enable_nf2ff: | |
# wavelength of minimum/maximum frequency used (in excitation) in simulation | |
min_freq_lambda = round(C0 / (f0 - fc) / unit) | |
max_freq_lambda = round(C0 / (f0 + fc) / unit) | |
# distance of transition between near-field to far-field | |
nf_ff_transition_distance = math.ceil(min_freq_lambda / (2 * math.pi)) | |
# simulation mesh resolution for far-field | |
mesh_res_farfield = round(max_freq_lambda / 30) | |
# add the NF2FF recording box | |
start = [-nf_ff_transition_distance / 2] * 3 | |
stop = [nf_ff_transition_distance / 2] * 3 | |
nf2ff = fdtd.CreateNF2FFBox("nf2ff-box", start=start, stop=stop, opt_resolution=[mesh_res_farfield] * 3) | |
# smooth out mesh for far-field | |
# mesh.SmoothMeshLines("all", mesh_res_farfield, 1.4) | |
######################################################################################### | |
# fire up AppCSXCAD for viewing the model before running it | |
# | |
if enable_appcsxcad: | |
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) | |
# Found resonance frequency at 446.2 MHz with -42.5 dB at 71.1 Ohm | |
# Dipole (lambda/2) length is 289.8 mm | |
freq = np.linspace(f0 - fc, f0 + fc, 2001) | |
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)) | |
# Found resonance frequency at 446.1 MHz with -42.4 dB at 71.0 Ohm | |
# Dipole (lambda/2) length is 289.8 mm | |
print(s11_dB) | |
print(freq) | |
print(type(s11_dB), len(s11_dB)) | |
print(type(freq), len(freq)) | |
cutoff_db_resonance = -10.0 | |
cutoff_dbs = [cutoff_db_resonance, -15.0, -20.0, -25.0, -30.0, -40.0, -50.0] | |
cutoff_dbs_results = {} | |
idx = np.where((s11_dB < cutoff_db_resonance) & (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("=" * 80) | |
print("") | |
print( | |
"Found resonance frequency at {} MHz with {} dB at {} Ohm".format( | |
round(freq[idx][0] / 1e6, 1), | |
round(s11_dB[idx][0], 1), | |
round(np.real(Zin[idx])[0], 1), | |
) | |
) | |
print("Dipole (lambda/2) length is {} mm".format(round(dipole_length, 1))) | |
print("") | |
print("=" * 80) | |
print("") | |
cutoff_dbs_results['interest'] = {} | |
cutoff_interest_bw = round((446.2e6 - 446.0e6) / 1e6, 1) | |
for kf, f in [('lower', 446.0e6), ('center', 446.1e6), ('upper', 446.2e6)]: | |
# Calculate absolute differences | |
abs_diff = np.abs(freq - f) | |
# Find the index of the closest value | |
closest_index = np.argmin(abs_diff) | |
print( | |
"S11 at frequency {} MHz is {} dB at {} Ohm [index {}]".format( | |
round(f / 1e6, 1), | |
round(s11_dB[closest_index], 1), | |
round(np.real(Zin[closest_index]), 1), | |
closest_index, | |
) | |
) | |
cutoff_dbs_results['interest'][kf] = { | |
'idx': closest_index, | |
'freq': round(freq[closest_index] / 1e6, 1), | |
's11': round(s11_dB[closest_index], 1), | |
'r': round(np.real(Zin[closest_index]), 1), | |
'bandwidth': cutoff_interest_bw, | |
} | |
print("") | |
for cutoff_db in cutoff_dbs: | |
idx_cutoff_lower = idx[0] | |
while idx_cutoff_lower >= 0 and s11_dB[idx_cutoff_lower] < cutoff_db: | |
idx_cutoff_lower -= 1 | |
# print("cutoff_lower index: {}".format(idx_cutoff_lower)) | |
idx_cutoff_upper = idx[0] | |
while idx_cutoff_upper <= len(s11_dB) and s11_dB[idx_cutoff_upper] < cutoff_db: | |
idx_cutoff_upper += 1 | |
# print("cutoff_upper index: {}".format(idx_cutoff_upper)) | |
print("") | |
print( | |
"S11 at frequency {} MHz is {} dB at {} Ohm".format( | |
round(freq[idx_cutoff_lower] / 1e6, 1), | |
round(s11_dB[idx_cutoff_lower], 1), | |
round(np.real(Zin[idx_cutoff_lower]), 1), | |
) | |
) | |
print( | |
"S11 at frequency {} MHz is {} dB at {} Ohm".format( | |
round(freq[idx_cutoff_upper] / 1e6, 1), | |
round(s11_dB[idx_cutoff_upper], 1), | |
round(np.real(Zin[idx_cutoff_upper]), 1), | |
) | |
) | |
cutoff_dbs_results[cutoff_db] = {} | |
cutoff_dbs_results[cutoff_db]['lower'] = { | |
'idx': idx_cutoff_lower, | |
'freq': round(freq[idx_cutoff_lower] / 1e6, 1), | |
's11': round(s11_dB[idx_cutoff_lower], 1), | |
'r': round(np.real(Zin[idx_cutoff_lower]), 1), | |
} | |
cutoff_dbs_results[cutoff_db]['upper'] = { | |
'idx': idx_cutoff_upper, | |
'freq': round(freq[idx_cutoff_upper] / 1e6, 1), | |
's11': round(s11_dB[idx_cutoff_upper], 1), | |
'r': round(np.real(Zin[idx_cutoff_upper]), 1), | |
} | |
cutoff_dbs_results[cutoff_db]['bandwidth'] = round((freq[idx_cutoff_upper] - freq[idx_cutoff_lower]) / 1e6, 1) | |
print("=" * 80) | |
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( | |
"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.legend() | |
# Calculate alpha values based on s11_dB | |
min_alpha = 0.05 | |
max_alpha = 0.8 | |
alpha_values = (s11_dB - cutoff_db_resonance) / (s11_dB.min() - cutoff_db_resonance) | |
alpha_values = np.clip(alpha_values, 0, 1) # Clip values between 0 and 1 | |
alpha_values = min_alpha - alpha_values * (min_alpha - max_alpha) | |
# Add vertical lines colored based on s11_dB values | |
idx_cutoff_lower = cutoff_dbs_results[cutoff_db_resonance]['lower']['idx'] | |
idx_cutoff_upper = cutoff_dbs_results[cutoff_db_resonance]['upper']['idx'] | |
res_cutoff_interest_lower = cutoff_dbs_results['interest']['lower'] | |
res_cutoff_interest_upper = cutoff_dbs_results['interest']['upper'] | |
pprint(res_cutoff_interest_lower) | |
pprint(res_cutoff_interest_upper) | |
# index of _maximum_ of S11 at both ends of frequency band of interest | |
idx_cutoff_interest = res_cutoff_interest_lower['idx'] if res_cutoff_interest_lower['s11'] > res_cutoff_interest_upper[ | |
's11'] else res_cutoff_interest_upper['idx'] | |
print('>>' * 100, idx_cutoff_lower, idx_cutoff_upper, idx, idx_cutoff_interest) | |
for i in range(len(freq)): | |
if idx_cutoff_lower <= i <= idx_cutoff_upper: | |
if i == idx_cutoff_lower or i == idx_cutoff_upper: | |
pyplot.axvline(x=freq[i] / 1e6, color='green', alpha=1.0, zorder=2) | |
elif i == idx: | |
pyplot.axvline(x=freq[i] / 1e6, color='blue', alpha=1.0, zorder=2) | |
pyplot.axhline(y=round(s11_dB[idx][0], 1), color='blue', alpha=1.0, zorder=2) | |
elif i == idx_cutoff_interest: | |
pyplot.axhline(y=round(s11_dB[i], 1), color='blue', alpha=1.0, zorder=2, linestyle='--') | |
pyplot.text(round(freq[i] / 1e6, 1), round(s11_dB[i], 1) + 1.0, | |
"{} MHz bandwidth @ {} dB".format(res_cutoff_interest_lower['bandwidth'], round(s11_dB[i], 1)), | |
ha="center", zorder=4, fontweight='bold') | |
elif i % 2: | |
pyplot.axvline(x=freq[i] / 1e6, color='green', alpha=alpha_values[i], zorder=0, linestyle='dotted') | |
x = round(freq[idx][0] / 1e6, 1) | |
for cutoff_db in cutoff_dbs: | |
pyplot.axhline(y=cutoff_db, color='green', alpha=1.0, linestyle='dotted') | |
pyplot.text(x, cutoff_db + 1.0, | |
"{} MHz bandwidth @ {} dB".format(cutoff_dbs_results[cutoff_db]['bandwidth'], cutoff_db), | |
ha="center", zorder=4, fontweight='bold') | |
idx_cutoff_lower = cutoff_dbs_results[cutoff_db]['lower']['idx'] | |
idx_cutoff_upper = cutoff_dbs_results[cutoff_db]['upper']['idx'] | |
# Add markers with text to specific coordinates | |
markers = [ | |
{ | |
'pos': (freq[idx_cutoff_lower] / 1e6, s11_dB[idx_cutoff_lower]), | |
'offset': [-4.0, -4.0], | |
'text': " {} MHz: {} dB\n@ {} Ohm".format(round(freq[idx_cutoff_lower] / 1e6, 1), | |
round(s11_dB[idx_cutoff_lower], 1), | |
round(np.real(Zin[idx_cutoff_lower]), 1)), | |
'color': 'red', | |
'ha': 'right', | |
}, | |
{ | |
'pos': (freq[idx_cutoff_upper] / 1e6, s11_dB[idx_cutoff_upper]), | |
'offset': [4.0, -4.0], | |
'text': " {} MHz: {} dB\n@ {} Ohm".format(round(freq[idx_cutoff_upper] / 1e6, 1), | |
round(s11_dB[idx_cutoff_upper], 1), | |
round(np.real(Zin[idx_cutoff_upper]), 1)), | |
'color': 'red', | |
'ha': 'left', | |
}, | |
{ | |
'pos': (freq[idx] / 1e6, s11_dB[idx]), | |
'offset': [2.0, 2.0], | |
'text': " {} MHz: {} dB\n@ {} Ohm".format(round(freq[idx][0] / 1e6, 1), | |
round(s11_dB[idx][0], 1), | |
round(np.real(Zin[idx][0]), 1)), | |
'color': 'blue', | |
'ha': 'left', | |
}, | |
] | |
for marker in markers: | |
pyplot.scatter(marker['pos'][0], marker['pos'][1], color=marker['color'], zorder=3) | |
pyplot.text(marker['pos'][0] + marker['offset'][0], marker['pos'][1] + marker['offset'][1], marker['text'], | |
ha=marker['ha'], zorder=4, fontweight='normal') | |
######################################################################################### | |
# compute far-field from recording box and generate plots | |
# | |
if enable_nf2ff: | |
# Calculate the far field at phi=0 degrees and at phi=90 degrees | |
theta = np.arange(0., 180., 1.) | |
phi = np.arange(-180, 180, 2) | |
print("=" * 80) | |
print("\n") | |
print('Calculating the 3D far field...') | |
# https://docs.openems.de/python/openEMS/nf2ff.html#openEMS.nf2ff.nf2ff.CalcNF2FF | |
# CalcNF2FF(sim_path, freq, theta, phi, radius=1, center=[0, 0, 0], outfile=None, read_cached=False, verbose=0) | |
# theta/phi – array like – Theta/Phi angles to calculate the far-field | |
# radius – float – Radius to calculate the far-field (default is 1m) | |
nf2ff_radius = nf_ff_transition_distance | |
print("Analyzing far-field at radius {}".format(nf2ff_radius)) | |
# 1) Analyze far-field for: single center frequency of interest | |
# Works! | |
# | |
# freqs_of_interest = [f0] | |
# | |
# Result: | |
# nf2ff: Analysing far-field for 1 frequencies. | |
# Radiated power: P_rad = 1.3781441597114614e-20 W | |
# Directivity: D_max = -42.66256438449489 dBi | |
# Efficiency: nu_rad = 101.85287356184361 % | |
# Theta_HPBW = 43.0 ° | |
# 2) Analyze far-field for: lower bound, center and upper bound frequency of interest | |
# Works! | |
# | |
# freqs_of_interest = [446.0e6, 446.1e6, 446.2e6] | |
# | |
# Result: | |
# nf2ff: Analysing far-field for 3 frequencies. | |
# Radiated power: P_rad = 1.3763972766838116e-20 W | |
# Directivity: D_max = -42.662132139823896 dBi | |
# Efficiency: nu_rad = 101.76329904670399 % | |
# Theta_HPBW = 43.0 ° | |
# 3) Analyze far-field for: all frequencies with S11 at least -10 dB | |
# Works! | |
# | |
idx_cutoff_lower = cutoff_dbs_results[cutoff_db_resonance]['lower']['idx'] | |
idx_cutoff_upper = cutoff_dbs_results[cutoff_db_resonance]['upper']['idx'] | |
freqs_of_interest = freq[idx_cutoff_lower:idx_cutoff_upper + 1] | |
# | |
# Result: | |
# nf2ff: Analysing far-field for 1206 frequencies. | |
# Radiated power: P_rad = 3.313485009786238e-21 W | |
# Directivity: D_max = -42.46238526982907 dBi | |
# Efficiency: nu_rad = 24.483622446993255 % | |
# Theta_HPBW = 43.0 ° | |
# 4) Analyze far-field for: all frequencies in excitation range [f0 - fc, f0 + fc] | |
# Does NOT work! | |
# | |
# OOM killed! - "nf2ff: Analysing far-field for 2001 frequencies." | |
# | |
# freqs_of_interest = freq | |
print("Analyzing far-field for {} frequencies:\n{}".format(len(freqs_of_interest), pformat(freqs_of_interest))) | |
nf2ff_res = nf2ff.CalcNF2FF(sim_path=output_dir, freq=freqs_of_interest, theta=theta, phi=phi, radius=nf2ff_radius, | |
read_cached=True, verbose=True) | |
Dmax_dB = 10 * np.log10(nf2ff_res.Dmax[0]) | |
E_norm = 20.0 * np.log10(nf2ff_res.E_norm[0] / np.max(nf2ff_res.E_norm[0])) + 10 * np.log10(nf2ff_res.Dmax[0]) | |
theta_HPBW = theta[np.where(np.squeeze(E_norm[:, phi == 0]) < Dmax_dB - 3)[0][0]] | |
# Display power and directivity | |
print('Radiated power: P_rad = {} W'.format(nf2ff_res.Prad[0])) | |
print('Directivity: D_max = {} dBi'.format(Dmax_dB)) | |
print('Efficiency: nu_rad = {} %'.format(100 * nf2ff_res.Prad[0] / np.interp(f0, freq, feed.P_acc))) | |
print('Theta_HPBW = {} °'.format(theta_HPBW)) | |
E_norm = 20.0 * np.log10(nf2ff_res.E_norm[0] / np.max(nf2ff_res.E_norm[0])) + 10 * np.log10(nf2ff_res.Dmax[0]) | |
E_CPRH = 20.0 * np.log10(np.abs(nf2ff_res.E_cprh[0]) / np.max(nf2ff_res.E_norm[0])) + 10 * np.log10( | |
nf2ff_res.Dmax[0]) | |
E_CPLH = 20.0 * np.log10(np.abs(nf2ff_res.E_cplh[0]) / np.max(nf2ff_res.E_norm[0])) + 10 * np.log10( | |
nf2ff_res.Dmax[0]) | |
# Plot the pattern | |
pyplot.figure() | |
pyplot.plot(theta, E_norm[:, phi == 0], 'k-', linewidth=2, label='$|E|$') | |
pyplot.plot(theta, E_CPRH[:, phi == 0], 'g--', linewidth=2, label='$|E_{CPRH}|$') | |
pyplot.plot(theta, E_CPLH[:, phi == 0], 'r-.', linewidth=2, label='$|E_{CPLH}|$') | |
pyplot.grid() | |
pyplot.xlabel('Theta (deg)') | |
pyplot.ylabel('Directivity (dBi)') | |
pyplot.title('Frequency: {} GHz'.format(nf2ff_res.freq[0] / 1e9)) | |
pyplot.legend() | |
######################################################################################### | |
# show all plots | |
# | |
pyplot.show() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment