DJ3EI's blog

from datetime import datetime, UTC
start_time = datetime.now(UTC)
start_time.strftime("%Y-%m-%d %H:%M:%S UTC")

# The list of frequencies where SWR is to be minimized:
# I habitually use base units in my code internally
# (like m, Hz, Ω, H, F).

# Convenient human-consumable derived units are converted
# from on input and often converted back to on output.
# Here, I translate the list of frequencies in MHz to Hz:

LOW_SWR_WANTED_FS = [x*1e6 for x in (7.1, 14.2, 18.1, 24.92, 28.97)]
LOW_SWR_WANTED_FS

import math
# I'm in the habit of using 0.75 mm² stranded wire:

WIRE_AREA = 0.75e-6 # in m²
WIRE_DIAMETER = math.sqrt(WIRE_AREA/math.pi) * 2
WIRE_DIAMETER, math.pi * (WIRE_DIAMETER/2)**2

DIPOLE_HEIGHT = 10.0 # This is the height above ground as given by Brian.
Z_CABLE = 400 # This is the impedance of the feedline used for the original ZS6BKW.
Z_TX_WANTS = 50 # The ubiquitous 50 Ω we want to see.

# Brian's description mentions "city ground":
CITY_GROUND_DIEL = 3
CITY_GROUND_CONDUCTIVITY = 1e-3

from antenna_simulation_driver import run_nec2pp
import antenna_simulation_driver

antenna_simulation_count = 0

def compute_antenna_z(
    qrg: float,
    dipole_length_1: float,
    dipole_length_2: float,
    dipole_height: float = DIPOLE_HEIGHT,
    wire_diameter: float = WIRE_DIAMETER,
    ground_diel_const: float = CITY_GROUND_DIEL,
    ground_conductivity: float = CITY_GROUND_CONDUCTIVITY,
    capture_output: bool = False,
    rp_line: str = "RP 0 37 144 1003 0.0 0.0 2.5 2.5 0.0 0.0"
) -> antenna_simulation_driver.Nec2ppOutput:
    λ = 3e8 / qrg
    simulation_input = ("CM Dipole.\n"
        # end of comment:
        "CE\n"

        # tag number, number of segments, x,y,z of endpoint, x,y,z of other endpoint, wire radius
        # Somewhat arbitrary decision: One segment every 30 cm.
        f"GW 1 {math.ceil(dipole_length_1 / 0.3)} "
        f"{-dipole_length_1:.3f} .0 {dipole_height:.3f} "
        f"-0.3 .0 {dipole_height:.3f} "
        f"{wire_diameter/2:.3e}\n"

        f"GW 2 1 "
        f"-0.3 .0 {dipole_height:.3f} "
        f"0.3 .0  {dipole_height:.3f} "
        f"{wire_diameter/2:.3e}\n"

        f"GW 3 {math.ceil(dipole_length_2 / 0.3)} "
        f"0.3 .0  {dipole_height:.3f} "
        f"{dipole_length_2:.3f} .0 {dipole_height:.3f} "
        f"{wire_diameter/2:.3e}\n"

        # end of geometry
        "GE\n"

        # Ground: 2 0 0 0 finite ground with no ground-screen, dielectric constant, conductivity in mhos / m
        f"GN 2 0 0 0 {ground_diel_const} {ground_conductivity}\n"

        # copper wire
        "LD 5 1 0 0 58.1e6\n"
        "LD 5 2 0 0 58.1e6\n"
        "LD 5 3 0 0 58.1e6\n"

        # Excitation: 0 Voltage source, 1 tag number and segment number where the excitation happens,
        # the following 0 is for general sanity, the following one or two floats give the (real or complex) voltage.
        f"EX 0 2 1 0 1.0\n"

        # Frequencies: 0: linear stepping, 1: number of frequencies stepped through, 0, 0,
        # then starting frequency in MHz and stepping increment.
        # The software that reads the simulation run's output assumes we have only one frequency.
        f"FR 0 1 0 0 {qrg*1e-6:.6f} 0.00\n"
        # Actually run the simulation:
        # "RP 0 37 145 1003 0.0 0.0 2.5 2.5 0.0 0.0\n"
        f"{rp_line}\n"
        "EN\n")
    # print(simulation_input)
    global antenna_simulation_count
    antenna_simulation_count += 1 # This does not work through multiprocessing.
    return run_nec2pp(simulation_input, capture_output)

# Sanity check: Compare with known quarter wave dipole.
# For this dipole, MMANA calculates an impedance of 14.04-1001j ("known").
# If we see some similar impedance here, that is some indication
# that the above input script generation for NEC2 is correct:
compute_antenna_z(7.05e6, 5.0, 5.0, 10.0, 1.6e-3, 13, 5e-3).input_and_impedance.impedance

import math
SPEED_OF_LIGHT = 3e8

# Z0 (Z + j Z0 tan(2πL/λ)) / (Z0 + j Z tan(2πL/λ))
def cable_transform(z_in: complex, z_cable: complex, f: float, electric_length: float) -> complex:
    t = math.tan(2 * math.pi * electric_length * f / SPEED_OF_LIGHT)
    return z_cable * (z_in + 1j*z_cable*t) / (z_cable + 1j*z_in*t)

# Sanity check, the numbers should come out the same:
cable_transform(250+10j, 50, 10e6, SPEED_OF_LIGHT/10e6/4), 50 * 50 / (250+10j)

import scipy

# I like to use the scipy.optimize.differential_evolution algorithm,
# admittedly without checking alternatives and also without
# delving into the many parameters that could be used to control it further.

# I noticed this algorith# does not always return, as a result,
# the best candidate ever examined,
# but often only some slightly worse candiate.

# This (rather pedestrian) wrapper fixes that:
# It simply remembers the best result achived thus far
# and the candidate solution that achived it.

def optimize_wrapper(minimize_me, bounds):

    best_xs = None
    best_result = float("Infinity")

    def wrapper(xs):
        nonlocal best_xs
        nonlocal best_result
        result = minimize_me(xs)
        if result < best_result:
            best_result = result
            best_xs = list(xs)
        return result

    result = scipy.optimize.differential_evolution(wrapper, bounds)

    if result.success and best_xs is not None:
        result.x = best_xs
        result.fun = best_result

    return result

import pandas
from multiprocessing import Pool

# We only deal with electrical feedline length.
# To actually derive at physical cable length,
# you'd need to take the velocity factor into account. (We don't.)

# The feedline lengths we even consider
# are those between 1 m and 40 m.
FEEDLINE_LENGTH_MIN = 1
FEEDLINE_LENGTH_MAX = 40

# Running the simulation is much more expensive than is
# running a few experiments with different feedline lenghts.
# So after each simulation run, we optimize feedline length for best SWRs.

def find_best_feedline_length(
    fs:list[float],
    zins:list[complex],
    swrs_wanted:list[float]
) -> tuple[float, float, list[float]]:
    """Input: Arrays of frequencies and impedances (at those frequencies), so two arrays of same length.

    Output: Tupel with three slots: First: best feedline length,
    second: resulting best sum of (swr-1)² at that feedline length,
    third: array of individual SWR values at the frequencies in fs.

    You'd normally want all SWRs to be 1.0.  If so, pass an array of all 1.0 as swrs_wanted.
    But it is also possible to optimize for best fit to other SWR values.
    Pass these other SWR values as swrs_wanted. In any event, that array should also
    have the same length as do fs and zins.
    """

    # Given frequencies in fs and resulting dipole impedances at those frequencies in zins,
    # this is the function we want to minimize.
    def swr_sum(xs):
        feedline_length = xs[0]
        swrs = (
            antenna_simulation_driver.swr(cable_transform(zin, Z_CABLE, f, feedline_length), Z_TX_WANTS)
            for f, zin in zip(fs, zins)
        )
        return sum((swr_is-swr_wanted)**2 for swr_is, swr_wanted in zip(swrs, swrs_wanted))

    # Optimize feedline length for minimal overall SWRs.
    # One of these days, I should try whether
    # https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html
    # gives better results.
    o_result = optimize_wrapper(swr_sum, ((FEEDLINE_LENGTH_MIN, FEEDLINE_LENGTH_MAX), ))
    if o_result.success:
        best_feedline_length = o_result.x[0]
        best_swrs = [
            antenna_simulation_driver.swr(cable_transform(zin, Z_CABLE, f, best_feedline_length), Z_TX_WANTS)
            for f, zin in zip(fs, zins)
        ]
        best_square_sum = o_result.fun
        return best_feedline_length, best_square_sum, best_swrs
    else:
        raise RuntimeError(f"{o_result}")

def zin_and_eff_from_args(args: tuple[float, float, float, ...]) -> tuple[complex, float]:
    """Helper function to run the simulation:
    Given arguments as our compute_antenna eats it (at least three are required),
    return that antenna's complex impedance and its radiation efficiency.
    """
    result = compute_antenna_z(*args)
    return result.input_and_impedance.impedance, result.average_gain.average_power_gain / 2

def optimize(wire_diameter: float,
             swrs_wanted:list[float] = [1.0 for f in LOW_SWR_WANTED_FS],
             min_half_dipole_length: float = 3,
             max_half_dipole_length: float = 50
            ):
    """The main optimization function."""
    global antenna_simulation_count
    best_found_thus_far = float("Infinity")

    def minimize_me(xs):
        global antenna_simulation_count
        nonlocal best_found_thus_far
        # This "outer optimization" only optimizes the dipole length.
        dipole_length_1 = xs[0]

        # Arguments for one "volley" of antenna simulations.
        # The "volley" consists of one antenna simulation for each frequency in LOW_SWR_WANTED_FS.
        argss = [
            (f, dipole_length_1, dipole_length_1, DIPOLE_HEIGHT, wire_diameter) for f in LOW_SWR_WANTED_FS
        ]
        # Run the simulations of one "volley" in parallel simultaneously, in separate processes.
        # This way, the work is distributed over the cores of a multicore computer.
        with Pool() as pool:
            zin_eff_s = list(pool.map(zin_and_eff_from_args, argss))
            # The global count is not accessible from separate processes, so we increment it here:
            antenna_simulation_count += len(argss)

        # Now that we have the impedances at various frequencies,
        # run a separate optimization to find the best common feedline length:
        feedline_length, swr_square_sum, swrs = \
            find_best_feedline_length(LOW_SWR_WANTED_FS, [zin_eff[0] for zin_eff in zin_eff_s], swrs_wanted)
        # This is the value we'll return:
        result = swr_square_sum
        # To keep the operator entertained while waiting for the optimization result,
        # output the sequence of best results found thus far:
        if result < best_found_thus_far:
            best_found_thus_far = result
            print(
                f"dipole_length_1: {dipole_length_1:.3f} m\n"
                f"dipole_length_2: {dipole_length_1:.3f} m\n"
                f"feedl_length:    {feedline_length:.3f} m\n"
                f"result:          {result:.5f}\n"
                f"SWRs:            {[round(swr,3) for swr in swrs]}\n"
                f"effs:            {[zin_eff[1] for zin_eff in zin_eff_s]}\n"
                f"simcount:        {antenna_simulation_count}\n"
            )
        return result

    # We assume that the dipole length we are looking for
    # will be somewhere in the range of 6 to 100 meters.
    o_result = optimize_wrapper(minimize_me, ((min_half_dipole_length, max_half_dipole_length),))

    if o_result.success:
        best_dipole_length_1 = o_result.x[0]
        best_argss = [
                (f, best_dipole_length_1, best_dipole_length_1) for f in LOW_SWR_WANTED_FS
        ]
        # Do the calculations one more time
        # (we did that same calculation before as part of the optimization,
        # but did not bother to keep the results):
        with Pool() as pool:
            best_zin_eff_s = list(pool.map(zin_and_eff_from_args, best_argss))
            antenna_simulation_count += len(best_argss)
        best_feedline_length, best_swr_square_sum, best_swrs = \
            find_best_feedline_length(
                LOW_SWR_WANTED_FS,
                [best_zin_eff[0] for best_zin_eff in best_zin_eff_s],
                swrs_wanted
            )
        print(
            "Best antenne found:\n"
            f"dipole_length_1: {best_dipole_length_1:.3f} m\n"
            f"dipole_length_2: {best_dipole_length_1:.3f} m\n"
            f"feedl_length:    {best_feedline_length:.3f} m\n"
            f"result:          {best_swr_square_sum:.5f}\n"
            f"SWRs, 3 digits:  {[round(swr,3) for swr in best_swrs]}\n"
            f"SWRs, 2 digits:  {[round(swr,2) for swr in best_swrs]}\n"
            f"efficiencies:    {[best_zin_eff[1] for best_zin_eff in best_zin_eff_s]}\n"
            f"simcount:        {antenna_simulation_count}\n"
        )
        return (best_dipole_length_1, best_dipole_length_1, best_feedline_length)
    else:
        raise RuntimeError(f"{o_result}")

best_dipole_length_1, best_dipole_length_2, best_feedline_length = \
    optimize(WIRE_DIAMETER)

intermediate_time = datetime.now(UTC)
intermediate_time.strftime("%Y-%m-%d %H:%M:%S UTC")
intermediate_duration = (intermediate_time - start_time).total_seconds()
print(f"{antenna_simulation_count} simulations thus far in {intermediate_duration:.1f} s, "
      f"so {intermediate_duration*1e3 / antenna_simulation_count:.1f} ms/simulation")

# Calculate deviation in percent

def deviation_in_percent(actual, should):
    return 100 * (actual-should) / should

deviation_in_percent(14.39, 14.25), deviation_in_percent(13.33, 13.3)

# Let us re-calculate the SWR with the values given by the author,
# using our software.

recalculated_swrs = []

for f in LOW_SWR_WANTED_FS:
    z_in = compute_antenna_z(f, 14.25, 14.25).input_and_impedance.impedance
    z_antenna = cable_transform(z_in, Z_CABLE, f, 13.3)
    recalculated_swr = antenna_simulation_driver.swr(z_antenna)
    recalculated_swrs.append(recalculated_swr)
    print(f"At {f*1e-6:6.3f} MHz, ZS6BKW has SWR {recalculated_swr:.2f}")

# Another optimisation run, this time to find an antenna
# that best reproduces the SWR values as published by Brian:

def find_wire_diameter_brian_used() -> float:
    brians_swrs = [1.1, 1.1, 1.3, 1.4, 1.4]
    global antenna_simulation_count
    best_found_thus_far = float("Infinity")

    def minimize_me(xs):
        global antenna_simulation_count
        nonlocal best_found_thus_far
        # This "outer optimization" only optimizes the dipole length.
        dipole_length_1 = xs[0]
        wire_diameter = xs[1]

        # Run one "volley" of antenna simulations,
        # which consists of one antenna simulation for each frequency in LOW_SWR_WANTED_FS.
        argss = [
            (f, dipole_length_1, dipole_length_1, DIPOLE_HEIGHT, wire_diameter) for f in LOW_SWR_WANTED_FS
        ]
        # Run the simulations of one "volley" in parallel, in separate processes.
        # This way, the work is distributed over the cores of a multicore computer.
        with Pool() as pool:
            zin_eff_s = list(pool.map(zin_and_eff_from_args, argss))
            # The global count is not accessible from separate processes, so we increment it here:
            antenna_simulation_count += len(argss)

        # Now that we have the impedances at various frequencies,
        # run a separate optimization to find the best common feedline length:
        feedline_length, swr_square_sum, swrs = \
            find_best_feedline_length(
                LOW_SWR_WANTED_FS,
                [zin_eff[0] for zin_eff in zin_eff_s],
                brians_swrs
            )
        # This is the value we'll return:
        result = swr_square_sum
        # To keep the operator entertained, output the best result found thus far:
        if result < best_found_thus_far:
            print(
                f"wire_diameter:   {wire_diameter*1e3:.2f} mm\n"
                f"dipole_length_1: {dipole_length_1:.3f} m\n"
                f"dipole_length_2: {dipole_length_1:.3f} m\n"
                f"feedl_length:    {feedline_length:.3f} m\n"
                f"result:          {result:.5f}\n"
                f"SWRs:            {[round(swr,3) for swr in swrs]}\n"
                f"effs:            {[zin_eff[1] for zin_eff in zin_eff_s]}\n"
                f"simcount:        {antenna_simulation_count}\n"
            )
            best_found_thus_far = result
        return result

    # We assume that the dipole length we are looking for will be somewhere
    # in the range of 6 to 100 meters.  What we give to the optimizer
    # is one half of the dipole length.
    o_result = optimize_wrapper(minimize_me, ((13.5, 15.0),(0.5e-3, 5e-3),))
    if o_result.success:
        best_dipole_length_1 = o_result.x[0]
        best_wire_diameter = o_result.x[1]
        best_argss = [
            (f, best_dipole_length_1, best_dipole_length_1, DIPOLE_HEIGHT, best_wire_diameter)
            for f in LOW_SWR_WANTED_FS
        ]
        with Pool() as pool:
            best_zin_eff_s = list(pool.map(zin_and_eff_from_args, best_argss))
            antenna_simulation_count += len(best_argss)

        best_feedline_length, best_swr_square_sum, best_swrs = \
            find_best_feedline_length(
                LOW_SWR_WANTED_FS,
                [best_zin_eff[0] for best_zin_eff in best_zin_eff_s],
                brians_swrs
            )

        print(
            "Best antenne found:\n"
            f"wire_diameger:   {best_wire_diameter*1e3:.2f} mm\n"
            f"dipole_length_1: {best_dipole_length_1:.3f} m\n"
            f"dipole_length_2: {best_dipole_length_1:.3f} m\n"
            f"feedl_length:    {best_feedline_length:.3f} m\n"
            f"result:          {best_swr_square_sum:.5f}\n"
            f"SWRs, 3 digits:  {[round(swr,3) for swr in best_swrs]}\n"
            f"SWRs, 2 digits:  {[round(swr,2) for swr in best_swrs]}\n"
            f"effs:            {[best_zin_eff[1] for best_zin_eff in best_zin_eff_s]}\n"
            f"simcount:        {antenna_simulation_count}\n"
        )
        return (best_wire_diameter, best_dipole_length_1, best_dipole_length_1, best_feedline_length)
    else:
        raise RuntimeError(f"{o_result}")

repro_swr_wire_diameter, repro_swr_dipole_length_1, repro_swr_dipole_length_2, repro_swr_feedline_length = \
    find_wire_diameter_brian_used()

intermediate_time = datetime.now(UTC)
intermediate_time.strftime("%Y-%m-%d %H:%M:%S UTC")
intermediate_duration = (intermediate_time - start_time).total_seconds()
print(f"{antenna_simulation_count} simulations thus far in {intermediate_duration} s, "
      f"so {intermediate_duration*1e3 / antenna_simulation_count} ms/simulation")

# Try to optimize for thick wire:
optimize(2e-3)
pass

# Try to optimize for thick wire, this time even giving a head start:
optimize(1.97e-3, min_half_dipole_length = 14.2, max_half_dipole_length = 14.5)

intermediate_time = datetime.now(UTC)
intermediate_time.strftime("%Y-%m-%d %H:%M:%S UTC")
intermediate_duration = (intermediate_time - start_time).total_seconds()
print(f"{antenna_simulation_count} simulations thus far in {intermediate_duration} s, "
      f"so {intermediate_duration*1e3 / antenna_simulation_count} ms/simulation")

def draw_diagrams(from_f, to_f, dipole_length_1, dipole_length_2, feedline_length):
    global antenna_simulation_count
    num_of_points = 500
    fs = [from_f + i * (to_f - from_f) / num_of_points for i in range (0,num_of_points+1)]
    argss = [(f, dipole_length_1, dipole_length_2) for f in fs]
    with Pool() as pool:
        z_and_eff_s = list(pool.map(zin_and_eff_from_args, argss))
        antenna_simulation_count += len(argss)

    swrs = [
        antenna_simulation_driver.swr(
            cable_transform(z_and_eff[0], Z_CABLE, f, feedline_length),
            Z_TX_WANTS
        )
        for z_and_eff, f in zip(z_and_eff_s, fs)
    ]
        

    df = pandas.DataFrame(
        {
            # "efficiency": [z_and_eff[1] for z_and_eff in z_and_eff_s],
            "swr": swrs,            
        },
        index = [f * 1e-6 for f in fs]
    )
    print(f"Frequency of best SWR and that SWR:\n{df[["swr"]].loc[df[["swr"]].idxmin()]}")
    df.plot(
        figsize = (12,12),
        subplots = True,
        title = f"SWR of ZS6BKW antenna.",
        xlabel = "Frequency / MHz",
        ylabel = "SWR",
        ylim = (0.9, 5.5),
        grid = True
    )
    return fs, swrs

swr_fs, swrs = draw_diagrams(6e6, 30e6, best_dipole_length_1, best_dipole_length_2, best_feedline_length)

def find_swr_change(min_f, max_f, swr_boundary = 2) -> float:
    """Find and return a frequency f between min_f and max_f where SWR is precisely swr_boundary.
    
    This analyses the antenna given by the global variables
    best_dipole_length_1, best_dipole_length_2 and best_feedline_length.
    """
    def find_my_0(f: float) -> float:
        z_in = compute_antenna_z(f, best_dipole_length_1, best_dipole_length_2).input_and_impedance.impedance
        z = cable_transform(z_in, Z_CABLE, f, best_feedline_length)
        swr = antenna_simulation_driver.swr(z, Z_TX_WANTS)
        return swr - 2
    frequency_of_change = scipy.optimize.toms748(find_my_0, min_f, max_f, disp=True, rtol=1e-10, xtol=100)
    return frequency_of_change


def find_bandwidths(fs: list[float], swrs: list[float]) -> list[tuple[float, float]]:
    """Go through a scan of swrs and find the intervals of SWR < 2, as tuples of f_min, f_max."""
    swr_boundary = 2
    if swrs[0] <= swr_boundary:
        raise RuntimeError("This starts unexpected")
    if swrs[-1] <= swr_boundary:
        raise RuntimeError("This ends unexpected")
    good_bands: list[tuple[float, float]] = []
    good_band_start = None
    for i in range(0, len(fs) - 1):
        if swrs[i+1] < 2 < swrs[i]:
            if good_band_start is None:
                good_band_start = find_swr_change(fs[i], fs[i+1], swr_boundary)
            else:
                raise RuntimeError("Unexpected")
        elif swrs[i] < 2 < swrs[i+1]:
            if good_band_start is not None:
                good_band_end = find_swr_change(fs[i], fs[i+1], swr_boundary)
                good_bands.append((good_band_start, good_band_end))
                good_band_start = None
            else:
                raise RuntimeError("Unexpected")
    return good_bands

good_sw_bands = find_bandwidths(swr_fs, swrs)

for from_f, to_f in good_sw_bands:
    print(
        f"{from_f*1e-6:5.2f} - {to_f*1e-6:5.2f} MHz, center {0.5e-6*(from_f+to_f):5.2f} MHz, "
        f"bandwidth {(to_f - from_f)*1e-3:3.0f} kHz"
    ) 

sixmeter_fs, sixmeter_swrs = \
    draw_diagrams(
        49.7e6, 52.3e6, best_dipole_length_1, best_dipole_length_2, best_feedline_length
    )

good_sixmeter_bands = find_bandwidths(sixmeter_fs, sixmeter_swrs)

for from_f, to_f in good_sixmeter_bands:
    print(
        f"{from_f*1e-6:5.2f} - {to_f*1e-6:5.2f} MHz, center {0.5e-6*(from_f+to_f):5.2f} MHz, "
        f"bandwidth {(to_f - from_f)*1e-3:3.0f} kHz"
    ) 

sixmeter_result = compute_antenna_z(
    51e6, best_dipole_length_1, best_dipole_length_2, rp_line="RP 0 91 360 1001 0.0 0.0 1 1 0.0 0.0"
)

def find_best_direction(sim_result: antenna_simulation_driver.Nec2ppOutput) -> list[tuple[float, float, float]]:
    best_theta, best_phi, best_power_gain = None, None, -float("Infinity")

    results = []
    
    for ray in sim_result.radiation_pattern.rays:
        if best_power_gain < ray.power_gain_total:
            best_theta = ray.theta
            best_phi = ray.phi
            best_power_gain = ray.power_gain_total
            results = [(best_theta, best_phi, best_power_gain)]
        elif best_power_gain == ray.power_gain_total:
            results.append((ray.theta, ray.phi, ray.power_gain_total))
    
    return results
for some_theta, some_phi, some_power_gain in find_best_direction(sixmeter_result):
    print(f"{90-some_theta:5.1f}° above horizon, {some_phi:5.1f}° from wire, {some_power_gain:5.2f} dB")

def calculate_ideal_dipole():

    def freespace_dipole_simulation(qrg: float) -> antenna_simulation_driver.Nec2ppOutput:

        dipole_length_1 = 10
        dipole_length_2 = 10
        dipole_height = 0
        wire_diameter = 0.001

        # rp_line="RP 0 181 361 1003 0.0 0.0 1.0 1.0 0.0 0.0"
        rp_line="RP 0 361 721 1001 0.0 0.0 0.5 0.5 0.0 0.0"

        simulation_input = ("CM Dipole.\n"
        # end of comment:
        "CE\n"

        # tag number, number of segments, x,y,z of endpoint, x,y,z of other endpoint, wire radius
        f"GW 1 {math.ceil(dipole_length_1 / 0.3)} "
        f"{-dipole_length_1:.3f} .0 {dipole_height:.3f} "
        f"-0.3 .0 {dipole_height:.3f} "
        f"{wire_diameter/2:.3e}\n"

        f"GW 2 1 "
        f"-0.3 .0 {dipole_height:.3f} "
        f"0.3 .0  {dipole_height:.3f} "
        f"{wire_diameter/2:.3e}\n"

        f"GW 3 {math.ceil(dipole_length_2 / 0.3)} "
        f"0.3 .0  {dipole_height:.3f} "
        f"{dipole_length_2:.3f} .0 {dipole_height:.3f} "
        f"{wire_diameter/2:.3e}\n"

        # end of geometry
        "GE\n"

        # free space
        "GN -1\n"

        # perfect ground
        # "GN 1\n"

        # bad city ground
        # f"GN 2 0 0 0 {CITY_GROUND_DIEL} {CITY_GROUND_CONDUCTIVITY}\n"

        # ideal wire
        "LD -1\n"

        # copper wire
        # "LD 5 1 0 0 58.1e6\n"
        # "LD 5 2 0 0 58.1e6\n"
        # "LD 5 3 0 0 58.1e6\n"

        # Excitation: 0 Voltage source, 1 tag number and segment number where the excitation happens,
        # the following 0 is for general sanity, the following one or two floats give the (real or complex) voltage.
        f"EX 0 2 1 0 1.0\n"

        # Frequencies: 0: linear stepping, 1: number of frequencies stepped through, 0, 0,
        # then starting frequency in MHz and stepping increment.
        # The software that reads the simulation run's output assumes we have only one frequency.
        f"FR 0 1 0 0 {qrg*1e-6:.6f} 0.00\n"
        # Actually run the simulation:
        # "RP 0 37 145 1003 0.0 0.0 2.5 2.5 0.0 0.0\n"
        f"{rp_line}\n"
        "EN\n")
        # print(simulation_input)
        global antenna_simulation_count
        antenna_simulation_count += 1 # This does not work through multiprocessing.
        return run_nec2pp(simulation_input, True)

    # Find frequency of resonance:
    frequency_of_resonance  = scipy.optimize.toms748(
        lambda f: freespace_dipole_simulation(f).input_and_impedance.impedance.imag,
        7.3e6,
        7.6e6,
        disp=True
    )
    print(f"Found resonance at {frequency_of_resonance*1e-3:.3f} kHz for 20 m dipole of 1 mm diameter wire.")
    return freespace_dipole_simulation(frequency_of_resonance)

ideal_dipole_data = calculate_ideal_dipole()

print(ideal_dipole_data.power_budget)
print(ideal_dipole_data.average_gain)
print(ideal_dipole_data.input_and_impedance)

best_gain = max((ray.power_gain_total for ray in ideal_dipole_data.radiation_pattern.rays))

print(f"Dipole gain found to be {best_gain:.2f} dBi by simulation, should have been 2.15 dBi.")

# with open("dipole.out", "w") as dipol_f:
#     dipol_f.write(result.raw_output)

adjusted_best_gain = best_gain + \
    10 * math.log(1 / ideal_dipole_data.average_gain.average_power_gain) / math.log(10)
adjusted_best_gain

end_time = datetime.now(UTC)
end_time.strftime("%Y-%m-%d %H:%M:%S UTC")
duration = (end_time - start_time).total_seconds()
print(
    f"{antenna_simulation_count} simulations in {duration} s, "
    f"so {duration*1e3 / antenna_simulation_count} ms/simulation"
)