Моделирование антенн на Python при помощи PyNEC

29 марта 2021

NEC2++ — это движок для моделирования антенн, совместимый с NEC2. Движок написан на C++ и распространяется под лицензией GPL. Также предусмотрены байндинги к языку Python. О них далее и пойдет речь.

Для экспериментов нам понадобятся следующие пакеты:

pip3 install PyNEC==1.7.3.4 matplotlib

Во избежание конфликтов между пакетами и чтобы не засорять систему, рекомендую использовать virtualenv. PyNEC указан конкретной версии по той причине, что в самой последней версии на момент написания статьи был открытый баг. Он делал использование пакета невозможным. К моменту, когда вы будете читать эти строки, баг может быть устранен.

Моделировать будем антенну inverted-V на диапазон 20 метров, установленную на 10-и метровой мачте. У нас будет несколько скриптов, анализирующих разные свойства антенны. Общий код положим в файл common.py:

import PyNEC
from math import pi, sin

def free_space(nec):
    #  0 - no ground plane
    #  1 - ground plane with current expansion
    # -1 - ground plane without current expansion
    nec.geometry_complete(0)

def good_ground(nec):
    conductivity = 0.0303 # S/m
    rel_dielectric_constant = 20.0
    nec.geometry_complete(1)
    nec.gn_card(0, 0, rel_dielectric_constant, conductivity,
                0, 0, 0, 0)

def average_ground(nec):
    conductivity = 0.005 # S/m
    rel_dielectric_constant = 13.0
    nec.geometry_complete(1)
    nec.gn_card(0, 0, rel_dielectric_constant, conductivity,
                0, 0, 0, 0)

def poor_ground(nec):
    conductivity = 0.001 # S/m
    rel_dielectric_constant = 3.0
    nec.geometry_complete(1)
    nec.gn_card(0, 0, rel_dielectric_constant, conductivity,
                0, 0, 0, 0)

def geometry(length=5.22, angle=100, height=10, nr_segments=21,
             ground=average_ground):
    nec = PyNEC.nec_context()
    nec.set_extended_thin_wire_kernel(False)

    wire_radius = 0.0005 # 0.5 mm
    length_ratio = 1.0
    radius_ratio = 1.0

    eps = 0.1
    b = length*sin(pi*angle/2/180)
    a = length*sin(pi*(90 - angle/2)/180)

    geo = nec.get_geometry()
    # center
    geo.wire(1, 1, 0, -eps, height, 0, eps, height,
             wire_radius, length_ratio, radius_ratio)
    # left
    geo.wire(2, nr_segments, 0, -b, height-a, 0, -eps, height,
             wire_radius, length_ratio, radius_ratio)
    # right
    geo.wire(3, nr_segments, 0, b, height-a, 0, eps, height,
             wire_radius, length_ratio, radius_ratio)

    ground(nec)

    # voltage excitation in the center of the antenna
    segment_nr = 1
    voltage = 1.0+0j
    ex_voltage_excitation = 0
    nec.ex_card(ex_voltage_excitation, segment_nr, 1, 1,
                voltage.real, voltage.imag, 0, 0, 0, 0)
    return nec

Здесь основная функция — это geometry. Она возвращает объект nec_context. Объект имеет N-ое количество методов-оберток к движку, а также хранит стейт с информацией о свойствах земли, проводах, из которых состоит антенна, и что именно в ней нужно промоделировать.

Не могут не бросаться в глаза странные имена методов, вроде gn_card и ex_card. Эта терминология позаимствована из движка NEC2. Движок принимает на вход файл, содержащий строки в определенном формате. Файл в терминологии NEC2 называется колодой (deck), а строки в нем — картами (card). Нужно понимать, что история NEC2 уходит корнями во времена, когда программисты писали на Fortran и работали с перфокартами. Отсюда и такая терминология.

Так вот, допустим, мы хотим описать землю. Для этого в файл нужно добавить строчку GN (какие-то аргументы). Говорят, что это GN-карта. Описание всех возможных карт, ровно как и их аргументов, можно найти в документации на NEC2. Пример колоды есть на Википедии. Еще в cocoaNEC можно открыть вкладку «Cards» и посмотреть, в какую колоду был странслирован код на NC.

Обладая этой информацией, код перестает казаться таким уж запутанным. Здесь мы объявили четыре варианта земли. Аргументы для GN-карт были банально подсмотрены в cocoaNEC. Далее в функции geometry мы рисуем inverted-V из трех проводов (geo.wire), подстилаем под ними землю, и говорим, что антенна питается в первый сегмент провода с меткой «1» (nec.ex_card).

Имея модель антенны, какая идея первой приходит на ум? Конечно же, построить график КСВ:

#!/usr/bin/env python3

import matplotlib.pyplot as plt
import PyNEC
from common import *

def vswr(z, z0):
    Gamma =  abs((z - z0) / (z + z0))
    return float((1 + Gamma) / (1 - Gamma))

nec = geometry()

system_impedance = 50
count = 100
start_freq = 13.8
stop_freq = 14.6
step_size = (stop_freq - start_freq) / count
ifrq_linear_step = 0
nec.fr_card(ifrq_linear_step, count, start_freq, step_size)
nec.xq_card(0) # execute simulation

freqs = []
vswrs = []
for idx in range(0, count):
    ipt = nec.get_input_parameters(idx)
    z = ipt.get_impedance()
    freqs.append(ipt.get_frequency() / 1000000)
    vswrs.append(vswr(z, system_impedance))

dpi = 80
plt.figure(dpi = dpi, figsize = (512 / dpi, 384 / dpi) )
plt.plot(freqs, vswrs)
plt.title("SWR vs Frequency")
plt.xlabel("Frequency, MHz")
plt.ylabel("SWR")
plt.axis([start_freq, stop_freq, 1.0, 3.0])
plt.grid(True)
plt.savefig("swr.png")

Здесь FR-карта задает интервал частот, на котором будет анализироваться антенна, а XQ-карта запускает моделирование. В итоге мы получаем данные для 100 частот с шагом 8 кГц. Для каждой частоты определяется входное сопротивление антенны. Из него вычисляется КСВ, график которого нам и нужен. График рисуется при помощи библиотеки Matplotlib.

Результат:

График КСВ, построенный на Python при помощи PyNEC

С тем же успехом входное сопротивление антенны можно отобразить на диаграмме Вольперта-Смита. Заинтересованным читателям предлагается модифицировать скрипт соответствующим образом в качестве упражнения.

Казалось бы, на КСВ можно остановиться, но уж нет! Раз есть модель антенны, то нарисуем полярную диаграмму ее усиления под углом 33° к горизонту на середине телеграфного участка:

#!/usr/bin/env python3

import matplotlib.pyplot as plt
import PyNEC
from common import *

nec = geometry()

freq_mhz = 14.035
ifrq_linear_step = 0
nec.fr_card(ifrq_linear_step, 1, freq_mhz, 0)

elevation = 33
theta_angle = 90 - elevation

nec.rp_card(calc_mode=0, n_theta=1, n_phi=360,
    output_format=1, normalization=0, D=0, A=0, # XNDA = 1000
    theta0=theta_angle, phi0=0, delta_theta=0, delta_phi=1,
    radial_distance=5000, gain_norm=0)
nec.xq_card(0) # execute simulation

rpt = nec.get_radiation_pattern(0)
gain = rpt.get_gain_tot()

gm = max(gain)
xs = [2*pi*float(phi)/360 for phi in range(gain.size)]
ys = [-36 if g-gm < -36 else g-gm for g in gain]

dpi = 80
fig = plt.figure(dpi = dpi, figsize = (512 / dpi, 384 / dpi) )
ax = plt.subplot(111, projection='polar')
ax.set_title('Elevation {} deg'.format(elevation))
ax.set_xlabel('Max gain: {:.2f} dBi'.format(gm))
ax.plot(xs, ys, linestyle = 'solid', color='red')
ax.set_rmax(0)
ax.set_rticks([-6*i for i in range(0,7)])
ax.set_yticklabels([''] + [str(-6*i) for i in range(1,7)])
ax.set_theta_direction(-1)
ax.set_rlabel_position(-90)
ax.set_thetagrids(range(0, 360, 15))
ax.grid(True)

fig.savefig('elevation.png')

FR-карта в этот раз задает одну, фиксированную частоту. За диаграмму направленности отвечает RP-карта. В этой карте нужно указать начальные значения углов phi и theta в сферической системе координат, с каким шагом и сколько раз они будут обновляться, плюс некоторые другие параметры. Заметьте, что зениту соответствует theta = 0. В итоге для всех пар (phi, theta) получаем усиление антенны в dBi. Его остается только отобразить на полярной диаграмме. Сталкиваться с последней нам уже доводилось:

Полярная диаграмма усиления антенны, построенная с помощью Python и PyNEC

И для полноты картины построим аналогичную диаграмму для фиксированного азимута и переменной элевации:

#!/usr/bin/env python3

import matplotlib.pyplot as plt
import PyNEC
from common import *

nec = geometry()

freq_mhz = 14.035
ifrq_linear_step = 0
nec.fr_card(ifrq_linear_step, 1, freq_mhz, 0)

azimuth = 0
nec.rp_card(calc_mode=0, n_theta=181, n_phi=1,
    output_format=1, normalization=0, D=0, A=0, # XNDA = 1000
    theta0=-90, phi0=azimuth, delta_theta=1, delta_phi=0,
    radial_distance=5000, gain_norm=0)
nec.xq_card(0) # execute simulation

rpt = nec.get_radiation_pattern(0)
gain = rpt.get_gain_tot()

gm = max(gain)
xs = [2*pi*float(90-theta)/360 for theta in range(gain.size)]
ys = [-36 if g-gm < -36 else g-gm for g in gain]

dpi = 80
fig = plt.figure(dpi = dpi, figsize = (512 / dpi, 384 / dpi) )
ax = plt.subplot(111, projection='polar')
ax.set_title('Azimuth {} deg'.format(azimuth))
ax.set_xlabel('Max gain: {:.2f} dBi'.format(gm))
ax.set_thetamin(-90) # set the limits
ax.set_thetamax(90)
ax.set_theta_offset(2*pi*90/360)
ax.plot(xs, ys, linestyle = 'solid', color='red')
ax.set_rmax(0)
ax.set_rticks([-6*i for i in range(0,7)])
ax.set_yticklabels([''] + [str(-6*i) for i in range(1,7)])
ax.set_thetagrids(range(-90, 91, 15))
ax.set_theta_direction(-1)
ax.grid(True)

fig.savefig('azimuth.png')

Результат:

Диаграмма направленности антенны, построенная с помощью Python и PyNEC

Диаграмма выглядит немного непривычно. Дело в том, что масштаб на ней, как и на предыдущей диаграмме, линейный. Но в книгах и в антенных моделировщиках мы привыкли видеть логарифмический масштаб. Заинтересованным читателям предлагается устранить этот дефект в качестве упражнения. Также в качестве упражнения предлагается нарисовать диаграмму направленности антенны в 3D. Полную версию исходников вы найдете на GitHub.

При помощи PyNEC могут быть решены разные задачи. Например, его можно использовать для оптимизации антенн по заданному критерию. Кроме того, библиотека может послужить основой для кроссплатформенного антенного моделировщика.

Метки: , , , .