DJ3EI's blog

import datetime
datetime.datetime.now(datetime.timezone.utc).strftime("%Y-%m-%d %H:%M:%S Z")

from math import log, pi, sqrt

c = 3e8 # speed of light in m/s

f_base = 7.1e6 # our frequency of interest in Hz
height_base = 10 # the height used throughout for the antenna, or at least the feedpoint.
wire_diameter_base = 1e-3 # meter

ω_base = 2 * pi * f_base
λ_base = c / f_base

normal_segments_count = 9

from typing import Optional

# Table for usual dielectric constant and conductivity for different grounds.
# We only use average and city.
GROUNDS = {
    "ideal": (None, None),
    "salt water": (81, 5),
    "excellent": (20, 3e-2),
    "average": (13, 5e-3),
    "industry": (5, 1e-3),
    "city": (3, 1e-3)
}

def generate_input(
    l: float,
    f: float,
    height: float,
    wire_diameter: float,
    ground_diel: Optional[float],
    ground_cond: Optional[float],
    num_segs: int,
    rp_line: str = "RP 0 37 144 1003 0.0 0.0 2.5 2.5 0.0 0.0"
) -> str:
    """Generate the input nec2++ expects.
    
    l dipole length in m.
    f frequency in Hz
    height above ground in m
    wire_diameter again in m.
    ground_diel dielectric constant.
    ground_cond ground conductivity. Omit both ground_diel and ground_cond for free space.
    num_segs Number of segments used to simulate one dipole half.
    
    """
    λ = 3e8 / f

    # Some documentation says nec2++ does not like segments shorter than 0.02 λ.
    # But using shorter segment does not seem to really hurt.
    # So let's throw caution into the wind and use 17 segments as the standard.

    input = ("CM Simple dipole antennas for various parameter studies.\n"
        # end of comment:
        "CE Here is one of many examples:\n"

        # tag number, number of segments, x,y,z of endpoint, x,y,z of other endpoint, wire radius
        f"GW 1 {2*num_segs+1}  {-l/2:.3f} .0 {height:.2f}   {l/2:.3f} .0 {height:.2f}  {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} {ground_cond}\n" \
            if ground_diel is not None or ground_cond is not None \
            else ""
          ) +

        # Copper wire
        "LD 5 1 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 1 {num_segs+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 {f*1e-6:.6f} 0.00\n"

        # Actually trigger the simulation:
        + rp_line + 
        "\nEN\n" )
    return input

# Produce a sample output:
print(generate_input(9, f_base, height_base, wire_diameter_base, *GROUNDS["average"], 17))
print()
print(generate_input(9, f_base, height_base, wire_diameter_base, None, None, normal_segments_count))

from antenna_simulation_driver import run_nec2pp, Nec2ppOutput, RadiationPatternRay

trial_sim_result = run_nec2pp(generate_input(
    10, 7.05e6, height_base, 1.6e-3, *GROUNDS["average"], normal_segments_count
))
print(
    f"Center-fed 2 x 5 m dipole {height_base} m above average ground (diel. const. 13, conductivity 5e-3 S/m)\n"
    "from copper wire of diameter 1.6 mm, at 7050 kHz.\n"
)
print("Impedance in Ω according to MMANA simulated by DK2FQ: 14.04-1001j,\n"
      "here according to NEC2: ", trial_sim_result.input_and_impedance.impedance)

# It is easy to stand tall on shoulders of giants.
# This pulls in some stuff we'll be using throughout:

import pandas as pd
import matplotlib.pyplot as plt
import scipy.optimize

# Finding free space resonant dipole.
from typing import Callable, Tuple

sim_count: int = 0

def find_resonant_dipole_parameter(
    parameter_to_ant_simulation_input: Callable[[float],str],
    parameter_min, parameter_max, xtol: float = 1e-4
) -> Tuple[float, int]:
    def null_me(parameter: float) -> float:
        global sim_count
        sim_count += 1
        return run_nec2pp(parameter_to_ant_simulation_input(parameter)).input_and_impedance.impedance.imag
    return (
        scipy.optimize.toms748(null_me, parameter_min, parameter_max, xtol = xtol, disp=True),
        sim_count
    )

freespace_length, sims = find_resonant_dipole_parameter(
    lambda l: generate_input(l, f_base, 0.0, wire_diameter_base, None, None, normal_segments_count),
    c/f_base/2*0.93, c/f_base/2
)

freespace_length, sim_count

import math

# I gauge DX-performance of an antenna by the EIRP as follows:
# Assuming the reflecting ionosphere to be at a certain height.
ionosphere_height = 400e3
# I want to reach a faraway station at a certain distance.
distance_faraway_station = 2000e3
# I calculate the required altitude angle (above horizon, initially in radians):
def calculate_altitude_angle(ionosphere_height: float, distance_faraway_station: float) -> float:
    from math import pi, sqrt, cos, acos
    earth_radius = 6371e3
    # Consider this triangle:
    # line a: connects the transmitting antenna with the earth's center
    # line b: line the radiation ray (for simplicity, assumed to be straight)
    #         travels until it hits the ionosphere (for simplicity, assumed to be a point)
    # line c: line connecting that point in the ionosphere with the earth's center
    a = earth_radius
    c = earth_radius + ionosphere_height

    # Angle at earth center, in radians:
    beta = distance_faraway_station / 2 / earth_radius

    # Law of cosines:
    b = sqrt(a**2 + c**2 - 2*a*c*cos(beta))

    # Angle between a and b:
    gamma = acos((a**2 + b**2 - c**2)/(2*a*b))

    # Altitude (angle between ray and line to horizon)
    return gamma - pi/2

altitude_angle = calculate_altitude_angle(ionosphere_height, distance_faraway_station)
altitude_angle_degree = altitude_angle * 180 / math.pi
altitude_angle_degree    

def loss_db(loss: float) -> float:
    return -10 * math.log(loss) / math.log(10)

loss_db(0.5)

from dataclasses import dataclass

@dataclass
class ResultsAtHeight:
    height: float
    resonant_length: float
    impedance: complex
    loss_dB: float
    gain_up: float
    gain_dx: float
    sims_needed: int

def calc_for_height(
    height: float,
    ground_diel: Optional[float],
    ground_cond: Optional[float],
    altitude_angle_degree_wanted: float,
) -> ResultsAtHeight:
    resonant_length, sims_for_resonant_length = find_resonant_dipole_parameter(
        lambda l: generate_input(
            l,
            f_base,
            height,
            wire_diameter_base,
            ground_diel,
            ground_cond,
            normal_segments_count
        ), freespace_length * 0.9, freespace_length * 1.1
    )
    
    antenna_at_height_result = run_nec2pp(
        generate_input(
            resonant_length,
            f_base,
            height,
            wire_diameter_base,
            ground_diel,
            ground_cond,
            normal_segments_count,
            rp_line="RP 0 91 361 1001 0.0 0.0 1 1 0.0 0.0"
        )
    )
    # with open("dipole.out","w") as dipf:
    #     dipf.write(antenna_at_height_result.raw_output)

    result = ResultsAtHeight(
        height = height,
        resonant_length = resonant_length,
        impedance = antenna_at_height_result.input_and_impedance.impedance,
        loss_dB = loss_db(antenna_at_height_result.average_gain.average_power_gain / 2),
        gain_up = antenna_at_height_result.radiation_pattern.closest_ray(90, 0).power_gain_total,
        gain_dx = antenna_at_height_result.radiation_pattern.closest_ray(90, 90 - altitude_angle_degree).power_gain_total,
        sims_needed = sims_for_resonant_length + 1
    )
    return result

# calc_for_height(10, *GROUNDS["average"], 10)

# Prepare data:
import itertools
import multiprocessing

heights =  list(itertools.chain((0.1 * h for h in range(20,100)), (10 + 0.25 * h for h in range(0, 361))))

def calculate_height_avg_ground(h: float) -> ResultsAtHeight:
    return calc_for_height(h, *GROUNDS["average"], altitude_angle_degree)

def calculate_height_city_ground(h: float) -> ResultsAtHeight:
    return calc_for_height(h, *GROUNDS["city"], altitude_angle_degree)

start_ts = datetime.datetime.now(datetime.timezone.utc)

with multiprocessing.Pool() as pool:
    avg_data = pool.map(calculate_height_avg_ground, heights)
    city_data = pool.map(calculate_height_city_ground, heights)
    
end_ts = datetime.datetime.now(datetime.timezone.utc)
len(heights)

import os
number_of_simulations = \
    sum((d.sims_needed for d in avg_data)) + sum((d.sims_needed for d in city_data))
seconds_avg = (end_ts-start_ts).total_seconds()
print(
    f"{number_of_simulations} simulations done in {seconds_avg:.1f} seconds,"
    f"so {number_of_simulations/seconds_avg:.1f} sims/s or "
    f"{seconds_avg*1e3/number_of_simulations:.1f} ms per sim "
    f"(on a machine claiming {os.process_cpu_count()} CPUs)."
)

df = pd.DataFrame(
    {
        "resonant length avg": [100 * r.resonant_length / freespace_length for r in avg_data],
        "resonant length city": [100 * r.resonant_length / freespace_length for r in city_data],
        "loss / dB avg": [r.loss_dB for r in avg_data],
        "loss / dB city": [r.loss_dB for r in city_data],
        "impedance avg": [r.impedance.real for r in avg_data],
        "impedance city": [r.impedance.real for r in city_data],
        "gain NVIS dB avg": [r.gain_up for r in avg_data],
        "gain NVIS dB city": [r.gain_up for r in city_data],
        "gain DX dB avg": [r.gain_dx for r in avg_data],
        "gain DX dB city": [r.gain_dx for r in city_data],
    },
    index=heights
)
selected_heights = [(i <= 100 and i%10 == 0) or i%20 == 0 for i in range(0, len(df))]
df.iloc[selected_heights].round(2)

# Resonace length as a function of height,
# in percent of resonant length of the same wire in free space:

ax = df[["resonant length avg", "resonant length city"]].plot(
    figsize=(12,8),
    grid=True,
    subplots=False,
    lw=7,
    fontsize=21
)
ax.set_title("Resonant dipole length, in percent of resonant length in free space", fontsize=21)
ax.set_xlabel("Height above ground / m", fontsize=21)
ax.set_ylabel("%", fontsize=21)
pass

ax = df[["loss / dB avg", "loss / dB city"]].plot(
    figsize=(12,8),
    grid=True,
    subplots=False,
    lw=7,
    fontsize=21
)
ax.set_title("Overall dipole loss (wire + ground).", fontsize=21)
ax.set_xlabel("Height above ground / m", fontsize=21)
ax.set_ylabel("Loss / dB", fontsize=21)
df[["loss / dB avg", "loss / dB city"]].iloc[selected_heights].round(1)

# Impedance plot.
# As for each height, we first adjust the dipole length for resonance,
# the imaginary path is always 0, so Z is real.

ax = df[["impedance avg", "impedance city"]].plot(
    figsize=(12,8),
    grid=True,
    subplots=False,
    lw=7,
    fontsize=21
)
ax.set_title("Impedance of resonant half-wave dipole", fontsize=21)
ax.set_xlabel("Height above ground / m", fontsize=21)
ax.set_ylabel("Z / Ω", fontsize=21)
df[["impedance avg", "impedance city"]].iloc[selected_heights].round(1)

# How does the gain change with resonant dipole height?
# We are interested in a 2000 km QSO, which needs a fairly low radiation angle above horizon:

ax = df[["gain DX dB avg", "gain DX dB city"]].plot(
    figsize=(12,8),
    grid=True,
    subplots=False,
    lw=7,
    fontsize=21
)
ax.set_title(f"Gain for a {distance_faraway_station*1e-3:.0f} km QSO "
    f"(radiation {altitude_angle_degree:.0f}° above horizon).", fontsize=21)
ax.set_xlabel("Height above ground / m", fontsize=21)
ax.set_ylabel("Gain dBi", fontsize=21)
df[["gain DX dB avg", "gain DX dB city"]].iloc[selected_heights].round(1)

# How does the gain change with resonant dipole height?
# Let us try NVIS this time, so gain straight up:

ax = df[["gain NVIS dB avg", "gain NVIS dB city"]].plot(
    figsize=(12,8),
    grid=True,
    subplots=False,
    lw=7,
    fontsize=21
)
ax.set_title("Gain for NVIS QSOs (radiation straight up).", fontsize=21)
ax.set_xlabel("Height above ground / m", fontsize=21)
ax.set_ylabel("Gain dB EIRP", fontsize=21)
df[["gain NVIS dB avg", "gain NVIS dB city"]].iloc[selected_heights].round(1)

datetime.datetime.now(datetime.timezone.utc).strftime("%Y-%m-%d %H:%M:%S Z")