pymininec 1.2.0
Python version of the original MININEC Antenna Optimization code
License: MIT License
Requires Python: >=3.7

MININEC in Python
=================

:Author: Ralf Schlatterbeck <rsc@runtux.com>

.. |--| unicode:: U+2013   .. en dash
.. |__| unicode:: U+2013   .. en dash without spaces
    :trim:
.. |_| unicode:: U+00A0 .. Non-breaking space
    :trim:
.. |-| unicode:: U+202F .. Thin non-breaking space
    :trim:
.. |ohm| unicode:: U+02126 .. Omega
.. |numpy.linalg.solve| replace:: ``numpy.linalg.solve``
.. |scipy.integrate| replace:: ``scipy.integrate``
.. |scipy.special.ellipk| replace:: ``scipy.special.ellipk``

This is a rewrite of the original MININEC3 basic sources in Python.
Standard use-case like computation of feed impedance and far field are
implemented and are quite well tested. There is only a command line
interface.

An example of using the command-line interface uses the following
12-element Yagi/Uda antenna originally presented by Cebik [6]_ without
the resistive element loading and with slight corrections to the element
lengths to make the antenna symmetric (the errors in the display in
Cebik's article may be due to a display issue of the program used that
displays negative numbers with less accuracy than positive numbers).
You can get further help by calling ``pymininec`` with the ``--help``
option.  The command-line options for the 12-element antenna example
(also available in the file ``test/12-el.pym``) are::

    pymininec -f 148 \
    -w 22,-.51943,0.00000,0,.51943,0.00000,0,.00238 \
    -w 22,-.50165,0.22331,0,.50165,0.22331,0,.00238 \
    -w 22,-.46991,0.34215,0,.46991,0.34215,0,.00238 \
    -w 22,-.46136,0.64461,0,.46136,0.64461,0,.00238 \
    -w 22,-.46224,1.03434,0,.46224,1.03434,0,.00238 \
    -w 22,-.45989,1.55909,0,.45989,1.55909,0,.00238 \
    -w 22,-.44704,2.19682,0,.44704,2.19682,0,.00238 \
    -w 22,-.43561,2.94640,0,.43561,2.94640,0,.00238 \
    -w 22,-.42672,3.72364,0,.42672,3.72364,0,.00238 \
    -w 22,-.41783,4.53136,0,.41783,4.53136,0,.00238 \
    -w 22,-.40894,5.33400,0,.40894,5.33400,0,.00238 \
    -w 22,-.39624,6.0452,0,.39624,6.0452,0,.00238 \
    --theta=0,5,37 --phi=0,5,73 \
    --excitation-pulse=33 > 12-el.pout

When the command line options are in a file, comments (using '#' at the
start of the line) can be added.
Users on Linux can run this using (I'm sure Windows users can come up
with a similar command-line on Windows)::

    pymininec $(sed '/^#/d' test/12-el.pym) > 12-el.pout

or using grep::

    pymininec $(grep -v '#' test/12-el.pym) > 12-el.pout

This removes comments from the ``.pym`` file and passes the result as
command-line parameters to pymininec.
The resulting output file contains currents, impedance at the feedpoint
and antenna far field in dBi as tables. The output tries to reproduce
the format of the original Basic implementation of Mininec.

Usage
-----

There is a ``-h`` or equivalent ``--help`` option that gives a brief
summary of the available options. All long options can be abbreviated to
their shortest unique prefix. Values often follow an option, either as a
separate argument (with a space between option and value) or directly
connected to the option with ``=``, see the example in the introduction.

Frequency
+++++++++

The typical usage defines the frequency using the ``-f`` or
``--frequency`` option, frequency is specified in MHz. In addition it
often makes sense to define a frequency sweep: The start frequency of
the sweep is the value given by the ``-f`` option, with the options
``--frequency-increment=<increment>`` and ``--frequency-steps=<steps>``
specifying how many steps with size increment are performed.

Geometry
++++++++

Next the geometry of the antenna is defined. Currently straight wires,
Arcs (consisting of straight segments touching the arc at their ends)
and Helices (also consisting of straight segments touching the helix at
their ends) can be specified.

All lengths and coordinates are in meter, but see below for the
``--geo-scale`` option.

All geometry options take an optional first argument, the tag number.
If no tag is specified, geo objects are numbered starting with the
highest specified tag number or, if no tag is specified, one. It is
recommended to not mix geo objects with and without tag. If different
geo objects are present, arcs are numbered first, then helices, and
finally, wires.  Especially when mixing different geometry objects it is
highly recommended to not rely on auto-tagging but specify the tag
numbers explicitly.

Tags are used to refer back to the geometry objects with other
options, e.g. when specifying the feedpoint or loads on the antenna.

All geometry options take the number of segments as a parameter. This is
the number of short straight wire segments that constitute the geometry
object.

When specifying several connected geometry objects, endpoints need not
match exactly (like in the original Basic Mininec code) but must be
within 1/1000 of the shortest segment length (the same rule applies in
NEC for matching ends of geometry objects [20]_ (p. 8))

Wire
~~~~

A wire definition takes the form::

 --wire [tag,]<nseg>,<x1>,<y1>,<z1>,<x2>,<y2>,<z2>,<radius>

where (x1, y1, z1) specifies the first endpoint of the wire in
carthesian coordinates, (x2, y2, z2) defines the second endpoint of the
wire, and finally the radius is specified. The ``--wire`` option can be
abbreviated with ``-w``. The leading tag is optional.


The parameter ``nseg`` specifies the number of segments the wire is
split into for computation.

By default wires are segmented into equal-length segments. Segments
should not be longer than 1/20 the shortest wavelength λ. Mininec
typically produces results with a resonant frequency slightly too high
(so the antenna is too long), this can be improved by using a high
number of segments.

Segment Length Tapering of Wires
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

An alternative is to use segments that are not equal-length. This can be
done using the ``--taper-wire`` option. It gets 2 |__| 4 parameters. The first
is the wire tag to identify the wire that should be length-tapered. The
second parameter specifies from which end of the wire tapering should be
performed. The end can either be ``1``, ``2`` for first and second end,
respectively, and ``3`` for both ends. The third and fourth parameters
specify the minimum and maximum segments lengths, respectively. A
minimum segment length of 2.5 the radius is enforced: If segments get
shorter, the thin-wire asumption is violated and results will be wrong.
The tapering starts out with the shortest possible segment under the
constraints of minimum and maximum segment length and doubles the
segment length until the maximum is reached, it also stops making
segments longer when subsequent segments would have to be shorter than
the current segment to keep the number of segments.

Arc
~~~

A wire arc can be specified with a definition of the form::

 --arc [tag,]<nseg>,<arc-radius>,<arc-angle1><arc-angle2>,<wire-radius>

The tag is optional. The arc center is the origin. The axis is the
Y-axis. The first angle ``arc-angle1`` is the start angle,
``arc-angle2`` is the end angle. The angles are measured from the X-axis
in a left hand direction about the Y-axis in degrees.
The generated segments form a polygon *inscribed* within the arc.
If an arc in a different position or orientation is desired it can be
modified with the ``--geo-rotate`` and ``--geo-translate`` options.

Helix or Spiral
~~~~~~~~~~~~~~~

A helix can be specified with a definition of the form::

 --helix [tag,]<nseg>,<len>,<turn-len>,<wire-radius>,<xr1>,<yr1>[,<xr2>,<yr2>]

The initial tag number and the final two radii are optional.
The parameter ``len`` is the overall length of the helix. The
``turn-len`` is the length of one turn. The ``xr1`` and ``yr1``
parameters are the helix radius (from midpoint of the helix to the
midpoint of the wire) in X- and Y-direction. If the last two radii are
specified they indicate the X-direction and Y-direction radii *at the
end of the helix* while the first two radii are for the start of the
helix. The helix-radius is tapered from begin to end in this case.

The helix option can also be used to specify a spiral -- no ideal flat
spiral can be specified because the number of turns is determined from
the ratio of ``len` to ``turn-len`` but a very small ``len`` can be
used to create an almost-flat spiral.

The ``len`` and ``turn-len`` parameters may be negative. The sign does
not indicate that the helix grows in negative direction -- in fact, the
helix *always* starts at Z=0 and grows in a positive Z direction. The
signs determine where the helix starts and the winding direction. A
right-handed helix is produced when both signs are the same. The helix
starts at Y=0, X=``xr1`` when the ``len`` parameter is positive and at
X=0, Y=``yr1`` if it is negative.

This reproduces the behavior of the unofficial NEC ``GH`` card (``GH``
was an unofficial addition to NEC and is present in most
implementations, e.g., in ``nec2c``).

A Note on Modeling Helical Antennas
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

When modeling helical antennas, the first fraction of a turn and the
wire thickness have great influence on the feedpoint impedance. This
effect has been noted by other modelers [17]_. This means it is not
usually a good idea to start the helix on the ground plane directly and
locate the feedpoint on the ground plane. It is better to have a short
stub at the feedpoint or for a center-fed helix have a short wire from
the feedpoint on the ground plane to the start of the helix (typically
this wire has the same angle as the pitch angle of the helix [18]_).
When directly feeding the helix on the ground plane, depending on the
wire thickness, very small (non-physical) impedances can result because
this violates an asumption of the simulator that the wire does not lie
partially in the ground plane.

The effect that the feedpoint impedance gets lower when the wire is
thicker near the ground plane, is a physical effect, however. Kraus (the
inventor of the helical antenna) recommends to use this effect for
lowering the impedance for matching to a 50 |ohm| coax feed [19]_. He
uses a metal strip bonded to the helix conductor and alternatively
suggests to tune the spacing between conductor and ground plane when the
helix is not fed at the center.

Geometry Transformation
+++++++++++++++++++++++

Sometimes it is necessary to modify parts of the geometry. Three
geometry transformation options are available. To rotate part or all of
an antenna the ``--geo-rotate`` option is used. It gets 4 |__| 5 comma
separated parameters. The first is a numeric key for sorting geo
transformations: The order of transformations matters, so it is
necessary to specify the order. The next three options are the rotations
around X- Y- and Z-axis. An optional fifth parameter specifies the geo
object tag (e.g. wire tag) to rotate. If no tag is given the whole
antenna is rotated. If more than one rotation is non-zero, the X-axis
rotation is performed first, then the Y-axis rotation and finally the
Z-axis rotation.

The ``--geo-translate`` option again gets 4 |__| 5 comma separated
parameters. The first is again a sort key. The next three parameters
specify displacement in X- Y- and Z-direction. Finally again a tag can
be specified to define the geometry object to translate. If left out the
whole antenna is moved. This is often useful for modifying the height
above ground of an antenna: The whole antenna can be shifted up without
having to edit all the Z-components of all geometry elements.

Finally the ``--geo-scale`` option scales all geometry parameters
(including the wire radius) by a given factor. The factor is the first
parameter, an optional second parameter again gives a geometry tag. If
the tag is omitted the whole antenna is scaled. The scaling is always
applied last so that the ``--geo-translate`` option applies to the
original lengths.

An example is in ``test/vdipole-rot-trans.pym``: This has the geo
transformation options::

    --geo-rotate=1,0,0,90
    --geo-rotate=2,90,0,0
    --geo-translate=3,0,0,7.33

This first rotates the antenna around the Z-axis by 90° (sort-key 1),
then around the X-axis by 90° (sort-key 2), and finally the whole
antenna is shifted up by 7.33m (sort-key 3). Note that in this case we
cannot combine the rotation around Z- and X-axes into a single
``--geo-rotate`` option because this would rotate first around the
X-axis which would get a different result than first rotating around the
Z-axis.

Pulses
++++++

Mininec uses the concept of a *pulse* for defining where feed voltages
and loads apply. Think of a pulse as the point between two segments.
This means that at the *end* of wires (unless a second or third wire
connects there) there is no pulse. So a single wire consisting of 3
segments contains only 2 pulses, or generally a wire with ``N`` segments
contains ``N-1`` pulses. Pulses are automatically numbered starting with 1.

When a new wire is defined joining the endpoint of an already-existing
wire which has no connections yet, the pulse at the wire junction is
"owned" by the new wire: It becomes the first pulse on the new wire.

If more than two wires join at a coordinate, it is not a good idea to
allocate a feedpoint or load to that pulse: The feedpoint or load
would be only between two of the three or more wires. In such a case it
is better to insert a small length of wire where the feedpoint or load
is placed as in the following picture.

.. figure:: https://raw.githubusercontent.com/schlatterbeck/pymininec/master/feed.png
    :align: center


Feedpoint(s)
++++++++++++

For an antenna at least one feedpoint needs to be defined. This is done
using the ``--excitation-pulse`` option. The pulse number is either
absolute over *all* pulses of the antenna or a comma-separated sequence
of two values can be specified where the first is the pulse number
*relative to* the wire specified by the second number, the wire tag.
By default the excitation voltage is 1V but this can be changed by
specifying a ``--excitation-voltage`` option which gets a complex number
in volts. If multiple feedpoints are defined this is done by multiple
``--excitation-pulse`` and ``--excitation-voltage`` options.

Lumped Loads
++++++++++++

Adding loads to an antenna structure is a two-step process. In the first
step the loads are defined. In the second step they are attached to
pulses.

The easiest load type is specified with the ``-l`` or ``--load`` option.
It gets a complex number as argument. Note that this simple load type
does not change with frequency. Simple loads are sorted first when
attaching loads.

Laplace loads are the most general type of load. With it combinations of
`lumped element loads`_ can be modeled. In a combination of serial and
parallel lumped components, an inductance is modeled with L*s, a
capacitance is modeled with 1/(C*s) and a resistance with R. After
analyzing a complex circuit, a polynomial of s results in the numerator
and denominator of a fraction. The denominator is specified with the
``--laplace-load-a`` option and the numerator with the
``--laplace-load-b`` option. Both take a comma-separated list of real
numbers, representing the coefficients of the polynomial in increasing
order of exponentials. Laplace loads are sorted last when attaching
loads.

Another load type that is internally based on laplace loads is specified
with the ``--rlc-load`` option. It gets three parameters, the resistance
in Ohm, the inductance in Henry and the capacitance in Farad. A zero
in the Farad position indicates a short instead of a capacitance.  All
three lumped components are connected in series.  RLC loads are sorted
second when attaching loads.

Finally trap loads |--| which are also based on laplace loads internally
|--| allow the specification of traps in an antenna. They
consist of a serial connection of a resistor with an inductance
(modeling the non-zero resistance of a real inductance) in parallel with
a capacitance. The ``--trap-load`` option gets three comma-separated
real numbers, the resistance, the inductance, and the capacitance in
Ohm, Henry, and Farad, respectively. Trap loads are sorted third when
attaching loads.

Loads are attached to pulses with the ``--attach-load`` option. The
option takes 2 |__| 3 comma separated parameters. The first is the load
index. The load indeces are computed by iterating over all simple loads,
then all RLC loads, then all trap loads and finally all laplace loads
assigning them a load index starting with one.

Distributed Load on Wires
+++++++++++++++++++++++++

Non-ideal wires have distributed conductivity. With the option
``--skin-effect-conductivity`` distributed conductivity of a wire can be
specified. Alternatively the ``--skin-effect-resistivity`` option can be
used if the resistivity of the wire is known. Both option get one or two
parameters. The first parameter is the conductivity or resistivity,
respectively. The second optional parameter specifies the geometry
(e.g. wire) tag. If no tag is given, the skin effect load is attached to
*all* geometry objects.

Wires can have insulation. The effect of insulation on the distributed
impedance of a wire is modeled with the ``--insulation-load`` option. It
gets 2 |__| 3 parameters. The first parameter specifies the radius of the
wire *including* insulation. The second specifies the relative
permittivity of the insulation. The third optional parameter specifies
the geometry (e.g. wire) tag. If no tag is given the insulation load is
attached to all wires.

At most one insulation load and at most one skin effect load can be
specified per wire.

Ground and Radials
++++++++++++++++++

Ground can be specified with the ``--medium`` option. If not given, free
space is asumed. Multiple ``--medium`` options can be specified in which
case the subsequent media are either concentric around the first ground
or linearly allocated in X-direction. The ``--boundary`` option
specifies if the media are concentric (``--boundary=circular``) or in
X-direction (``--boundary=linear``) the default is a linear boundary.
The ``--medium`` option gets 3 |__| 4 comma-separated parameters, the
permittivity (dielectric constant), the conductivity, and the height.
If the first three are zero, ideal ground is asumed. With ideal ground
only a single ``--medium`` option is allowed.

The fourth parameter gives the width of the ground (the distance to the
next medium), this is a length in X-direction for linear boundary and a
radius for circular boundary. The fourth parameter is not used for the
last ``--medium`` option.

Note that you typically want *negative* heights for media further out,
this allows modelling of summits. Mininec *allows* the specification of
*higher* grounds but the results will be questionable as no reflection
at the higher ground is modelled. The first medium must always be at
height zero.

For the first medium, radials can be specified. Radials are allowed only
for non-ideal ground. The option ``--radial-count`` gives the number of
radials. The option ``--radial-radius`` gives the radius of the
radial-wires. Specifying radials will automatically select circular
boundary. The length of the radials is defined by the distance to the
next medium. So with radials at least two ``--medium`` options are
required.

Specifying What is Computed
+++++++++++++++++++++++++++

With the ``--option`` option it can be specified what outputs are
computed and printed. This option can be specified multiple times.
It can take the arguments ``far-field``, ``far-field-absolute``,
``near-field``, and ``none``. When ``none`` is specified as the only
option, only currents and feed point impedance are printed.

The ``far-field`` options selects printing of the far field in dBi.
The ``far-field-absolute`` option selects printing of the far field in
V/m. This option can be modified by specifying a different power level
using the `--ff-power`` option and the ``--ff-distance`` option to
specify the distance in radial direction of the far field measurement
point. Far field measurements are taken at elevation and azimuth angles
specified with the ``--theta`` and ``--phi`` options, respectively.
The elevation angle theta is measured from the zenith while the azimuth
angle phi is measured from the X-axis. Both, the ``--theta`` and the
``--phi`` option take tree comma-separated arguments: The start angle,
the angle increment, and the number of angles. By default, theta is
"0,10,10", so it runs from the zenith to ground in 10 degree steps. The
default for phi is "0,10,37", so it runs around the azimuth circle in 10
degree steps, computing the 0° and 360° on the X-axis value twice. This
is needed for some 3d-plotting tools for plotting a closed surface for
the 3d gain pattern.

The ``--option=near-field`` specifies printing of the near field.
This also needs specification of the ``--near-field`` option which gets
9 comma-separated parameters: The first three define the start (x, y, z)
coordinate of near-field measurements, the next three define the
increment of far field measurements and the last three define the number
of increments in each direction. With the ``--nf-power`` option it
is possible to modify the power level for the near field computation.

Without any ``--option``, far field is printed if no near field
options are present.

Miscellaneous Options
+++++++++++++++++++++

With the option ``--output-cmdline`` the given command-line options can
be printed. This is useful for tests and when using the API: All options
can be written out to reproduce the current settings. The option takes a
file name as an argument.

With the option ``--output-basic-input``, input for the original Mininec
code in Basic can be printed. The Basic code uses prompts to ask the
user for input. With this option the complete user input can be
generated. Running the Basic code with Yabasi_, the user input can be
fed into the Basic program with the ``-i`` option which is useful for
comparing the values computed by pymininec to the values computed by the
original Basic code. The option takes a file name as an argument.
There are different versions of the Basic code which have slightly
different user-input prompts. The version of the Basic code for which to
generate input can be specified with the ``--mininec-version`` option.
This option currently takes the values 9, 12, or 13. The default is the
widely-available version 9.

With the ``-T`` or ``--timing`` option, printing of runtimes of various
parts of the computation is requested. The option takes no arguments.

Measuring Timings
-----------------

Starting with version 1 there is a command-line option ``-T`` which
outputs computation timings on the standard error output. This was used
for measuring the results of recent vectorization of computations.
Speedups are roughly:

- About a factor of 50 for computation of the impedance matrix.
  So we're down from around 23 seconds to 0.44 seconds for a 12 element
  Yagi/Uda with 22 segments per element.
- About a factor of 200 for computation of the far field.
  So we're down from around 19 seconds to 0.09 seconds for the 12
  element Yagi/Uda with 5° resolution of azimuth and elevation angles.
  Even the computation of a 1° resolution takes below 2s for this
  antenna.
- About a factor of 5 for near-field computations. This could be further
  improved by batching the near field coordinates in chunks. I'm
  currently not using near-field computations much, so further
  improvements will wait until I have more need...

Plotting
--------

The output tables produced by ``pymininec``
are not very useful to get an idea of the far field behaviour of
an antenna. The companion program `plot-antenna`_ used to be bundled
with ``pymininec`` but was moved to its own project. You can currently
plot elevation and azimuth diagram of an antenna, a 3D-plot, the
geometry and VSWR. All either as a standalone program (using matplotlib)
or exported as HTML to the browser (using plotly).

Test Coverage and Code Quality
------------------------------

This section contains some notes on code quality and recent
improvements.

Test Coverage: making Sure it is Consistent with Original Mininec
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

There are several tests against the `original Basic source code`_, for
the test cases see the subdirectory ``test``. One of the test cases is
a simple 7MHz wire dipole with half the wavelength and 10 segments.
In one case the wire is 0.01m (1cm) thick, we use such a thick wire to
make the mininec code work harder because it cannot use the thin wire
assumptions. Another test is for the thin wire case. Also added are the
inverted-L and the T antenna from the original Mininec reports. All
these may also serve as examples.  Tests statement coverage is currently
at 100%.

There was a line that is flagged as not covered by the ``pytest``
framework if the Python version is below 3.10. This is a ``continue``
statement in ``compute_impedance_matrix`` near the end (as of this
writing line 1388). This is a bug in Python in versions below 3.10:
When setting a breakpoint in the python debugger on the continue
statement, the breakpoint is never reached although the continue
statement is correctly executed. A workaround would be to put a dummy
assignment before the continue statement and verify the test coverage
now reports the continue statement as covered.
I've `reported this as a bug in the pytest project`_ and `as a bug in
python`_, the bugs are closed now because Python3.9 does no longer get
maintenance.

For all the test examples it was carefully verified that the results are
close to the original results in Basic (see `Running examples in Basic`_
to see how you can run the original Basic code in the 21th century).

Some of the differences were due to rounding errors in the single
precision implementation in Basic compared to a double precision
implementation in Python. This happened because pcbasic_ uses single
precision float as the default. When running the Basic code with my
Basic interpreter Yabasi_, the default is double precision float.

I'm using numeric code from `numpy`_ where possible to speed up
computation, e.g. solving the impedance matrix is done using
|numpy.linalg.solve|_ instead of a line-by-line translation from Basic.
You can verify the differences yourself. In the ``test`` directory there
are input files with extension ``.mini`` which are intended (after
conversion to carriage-return convention when using pcbasic_, Yabasi_
can use the ``.mini`` files directly) to be used as input to the
original Basic code. The output of the Basic code is in files with the
extension ``.bout`` while the output of the Python code is in files
with the extension ``.pout``. The ``.pout`` files are compared in the
regression tests. The ``.pym`` files in the ``test`` directory are the
command-line arguments to recreate the ``.pout`` files with
``mininec.py``. An uppercase ``.Bout`` extension is used for output
generated with Yabasi_ where the distinction matters.

Near Field Accuracy
~~~~~~~~~~~~~~~~~~~

In his thesis [5]_, Zeineddin investigates numerical instabilities when
comparing near and far field. He solves this by doing certain
computations for the near field in double precision arithmetics.
I've tried to replicate these experiments and the numerical
instabilities are reproduceable in the Basic version. In the Python
version the instabilities are not present (because everything is in
double precision). But the absolute field values computed in Python are
lower than the ones reported by Zeineddin (and the Basic code *does*
reproduce Zeineddins values). I have not verified if these differences
result from Zeineddin using an earlier version of the Mininec Basic
code.

It doesn't look like there is a problem in the computations of the
currents in the Python code, the computed currents are lower than in
Basic which leads to lower field values. But the computed impedance
matrix when comparing both versions has very low error, see the test
``test_matrix_fill_ohio_example`` in ``test/test_mininec.py`` and the
routine ``plot_z_errors`` to plot the errors (in percent) in
``test/ohio.py``. Compared to the values computed by NEC [5]_, the Basic
code produces slightly higher values for near and far field while the
Python code produces slightly lower values than NEC. I've not tried to
simulate this myself in NEC yet.

You can find the files in
``test/ohio*`` (the thesis was at Ohio University). This time there is a
python script ``ohio.py`` to compute the near and far field values
without recomputing the impedance matrix. This script can show the near
and far field values in a plot and the difference in a second plot.
There are two distances for which these are computed, so the code
produces four plots. There is a second script to plot the Basic near and
far field differences ``plot_bas_ohio.py``.

Code Quality Before Vectorization
+++++++++++++++++++++++++++++++++

Before the vectorization this was the state of the code:

The current Python code is still hard to understand |--| it's the
result of a line-by-line translation from Basic, especially where I
didn't (yet) understand the intention of the code. The same holds for
Variable names which might not (yet) reflect the intention of the code.
I *did* move things like computation of the angle of a complex number,
or the computation of the absolute value, or multiplication/division of
complex numbers to the corresponding complex arithmetic in python where
I detected the pattern.

So the *de-spaghettification* was not successful in some parts of the
code yet :-) My notes from the reverse-engineering can be found in the
file ``basic-notes.txt`` which has explanations of some of the variables
used in mininec and some sub routines with descriptions (mostly taken
from ``REM`` statements) of the Basic code.

The code is also still quite slow: An example of a 12 element Yagi/Uda
antenna used in modeling examples by Cebik [6]_ takes about 50 seconds
on my PC (this has 264 segments, more than the original Mininec ever
supported) when I'm using 5 degree increments for theta and phi angles
and about 11 minutes (!) for 1 degree angles. The reason is that
everything currently is implemented (like in Basic) as nested loops.
This could (and should) be changed to use vector and matrix operations
in `numpy`_. In the inner loop of the matrix fill operation there are
several integrals computed using `gaussian quadrature`_ or a numeric
solution to an `elliptic integral`_. These are now implemented using
methods (or at least constants in the case of `gaussian quadrature`_)
from |scipy.integrate|_ and |scipy.special.ellipk|_.

Code Quality After Vectorization
++++++++++++++++++++++++++++++++

Before beginning the vectorization I've changed the implicit pulse
computations (this used a very complicated indexing schema to access
pulse information) to an explicit data structure in
``mininec/pulse.py``. This improved understandability of the code
considerably (so that I was able to refactor it further to vectorize
computations).

The current version still has obscure variable names from the Basic
implementations and in many cases it is not clear what intermediate
values in computations mean. Since Basic does not have complex numbers,
the semantics of computations can only be guessed. I hope to improve
on this when I get a version of [2]_ |--| the version available as
ADA181682_ contains many completely unreadable pages. If you have a
source of that report with better readability, let me know!

Multiple Inverted-V Example
+++++++++++++++++++++++++++

An old `web-page from 1998 by Dr. Carol F. Milazzo, KP4MD`_ has examples
of antennas simulated with Mininec. The first of these examples is three
crossed inverted-V (one of which has loading inductors to boost the
effective length). The simulation results of pymininec are in the
ballpark of the Mininec-based *NEC4WIN* which was used by KP4MD. But it
looks like *NEC4WIN* might use what it prints as "Diam." as the radius
of the wire (see Fig. 1 in the website) as the radius (see Antenna Model
Files in the Appendix). At least if this format is inherited from NEC
the last column of the wire definition would hold the radius and this
interpretation of the format also is more consistent with the simulation
results of Pymininec. The following table shows the original data
compared to using half of the diameter in the original model in
Pymininec ("Pymininec r") and the diameter as the radius (Pymininec 2r).
When using the (supposed) diameter for the radius, the output data
matches better to the website data.

+---------------+----------------+--------------+--------------+--------------+
| Frequency     |                | Original     | Pymininec r  | Pymininec 2r |
+---------------+----------------+--------------+--------------+--------------+
| 7MHz          | Gain Azimuth   | -2.42 dBi    | -2.52 dBi    | -2.49 dBi    |
|               +----------------+--------------+--------------+--------------+
|               | Gain Elevation |  7.21 dBi    |  7.21 dBi    |  7.21 dBi    |
|               +----------------+--------------+--------------+--------------+
|               | Impedance      | 38.74 +6.77j | 38.82 -3.66j | 39.28 +1.49j |
+---------------+----------------+--------------+--------------+--------------+
| 14MHz         | Gain Azimuth   |  4.33 dBi    |  4.60 dBi    |  4.37 dBi    |
|               +----------------+--------------+--------------+--------------+
|               | Gain Elevation |  7.23 dBi    |  7.73 dBi    |  7.38 dBi    |
|               +----------------+--------------+--------------+--------------+
|               | Impedance      | 46.16 -326j  | 31.86 -307j  | 43.00 -313j  |
+---------------+----------------+--------------+--------------+--------------+

All of KP4MD's examples have been converted to Pymininec and are available as
``inve802B.pym``, ``hloop40-14.pym``, ``hloop40-7.pym``,
``vloop20.pym``, and ``lzh20.pym`` in the ``test`` directory. Only the
``inve802B.pym`` (with the inverted-Vs) uses the diameter in the
original example as the radius in Pymininec, all others use half of the
value in the original example (which is supposed to be the diameter) as
the radius. But most examples match better to the values computed by
KP4MD when doubling the radius.

The Other Edge of The Sword
+++++++++++++++++++++++++++

There are some new tests that check the feedpoint impedance against
known computations from the literature. In particular an old article by
Roy Lewallen [8]_ with the same title as this section.

The column "Python" is from pymininec, the column "Basic
Yabasi" is the original Basic implementation run with my Basic
interpreter Yabasi_. The column "Basic pcbasic" uses the pcbasic_
interpreter.

Note that the "Bent Dipole" is bent horizontally (not an inverted V),
all wire ends are the same height. I have not been able so far to
reproduce the results of the special segmentation scheme that uses only
14 segements with the same results as indicated in the article (see then
entry ``14*`` for the bent dipole). When trying to reproduce it exactly
the imaginary part is much lower (more capacity). The segmentation
scheme is also not very good: In mininec adjacent segment should only
have a factor of 2 in length, not more. The segmentation special scheme
has a jump of factor 5, maybe this makes it numerically instable so that
we get much different results with double precision float.

For the bent dipole I've made three more experiments: One with tapering
from both ends (entry ``14t2``) and two with tapering from one end (entry
``14t1`` and ``14t1l``). Example ``14t1`` has no limit on segment length
while entry ``14t1l`` enforces a minimum segment length of 1/200 lambda.
In all the cases where tapering is from one end, the end with the
feedpoint has the smallest segment length. None of these experiments
comes close to the 14 segment experiment in the paper.

Straight Dipole
~~~~~~~~~~~~~~~

+------+----------------+----------------+----------------+----------------+
| Segs | Lewallen       | Python         | Basic Yabasi   | Basic pcbasic  |
+------+----------------+----------------+----------------+----------------+
| 10   | 74.073+20.292j | 74.074+20.298j | 74.074+20.298j | 74.074+20.300j |
+------+----------------+----------------+----------------+----------------+
| 20   | 75.870+21.877j | 75.872+21.897j | 75.872+21.897j | 75.872+21.897j |
+------+----------------+----------------+----------------+----------------+
| 30   | 76.573+23.218j | 76.567+23.169j | 76.567+23.169j | 76.572+23.203j |
+------+----------------+----------------+----------------+----------------+
| 40   | 76.972+24.053j | 76.972+24.052j | 76.972+24.052j | 76.973+24.068j |
+------+----------------+----------------+----------------+----------------+
| 50   | 77.222+24.517j | 77.240+24.647j | 77.240+24.647j |                |
+------+----------------+----------------+----------------+----------------+

Bent Dipole
~~~~~~~~~~~

+-------+----------------+----------------+----------------+----------------+
| Segs  | Lewallen       | Python         | Basic Yabasi   | Basic pcbasic  |
+-------+----------------+----------------+----------------+----------------+
| 10    | 11.509-76.933j | 11.498-77.045j | 11.498-77.045j | 11.498-77.044j |
+-------+----------------+----------------+----------------+----------------+
| 20    | 11.751-53.812j | 11.740-53.929j | 11.740-53.929j | 11.740-53.932j |
+-------+----------------+----------------+----------------+----------------+
| 30    | 11.819-46.934j | 11.808-47.068j | 11.808-47.068j | 11.808-47.055j |
+-------+----------------+----------------+----------------+----------------+
| 40    | 11.848-43.783j | 11.837-43.893j | 11.837-43.893j | 11.838-43.858j |
+-------+----------------+----------------+----------------+----------------+
| 50    | 11.861-41.988j | 11.851-42.107j | 11.851-42.107j |                |
+-------+----------------+----------------+----------------+----------------+
| 14*   | 11.312-43.119j | 11.104-47.879j |                                 |
+-------+----------------+----------------+---------------------------------+
| 14t1  |                | 10.859-42.486j |                                 | 
+-------+----------------+----------------+---------------------------------+
| 14t1l |                | 11.118-46.593j |                                 | 
+-------+----------------+----------------+---------------------------------+
| 14t2  |                | 11.314-45.659j |                                 | 
+-------+----------------+----------------+---------------------------------+


Running the Tests
+++++++++++++++++

You can run the tests with::

  python3 -m pytest test

If coverage should be reported this becomes::

  python3 -m pytest --cov mininec test

For a more detailed coverage report use::

  python3 -m pytest --cov-report term-missing --cov mininec test

This will show a detailed report of the lines that are not covered by
tests.

Skin Effect Loads
-----------------

[This section uses math in ReStructuredText which may not render
correctly on all platforms. In particular, `Github has an open issue`_
on this for more than a decade now. The formulas *are*
`supported on pypi`_]

To support skin effect loads on geometry objects (e.g. wires) we need to
compute the internal impedance of a segment. The `Wikipedia article
on skin effect`_ has the following formula for the internal impedance
per unit length:

.. math::
    \newcommand{\Int}{{\mathrm\scriptscriptstyle int}}
    \newcommand{\ber}{\mathop{\mathrm{ber}}\nolimits}
    \newcommand{\bei}{\mathop{\mathrm{bei}}\nolimits}

.. math::
    Z_\Int = \frac{k\rho}{2\pi r}\frac{J_0 (kr)}{J_1 (kr)}

where

.. math::
    k = \sqrt{\frac{-j\omega\mu}{\rho}}

and :math:`r` is the radius, :math:`J_v` are the Bessel functions of
the first kind of order :math:`v`. :math:`Z_\Int` is the impedance *per
unit le