Modelling a Physical Channel in the Frequency Domain

This example uses the frequency domain lyceanem.models.frequency_domain.calculate_scattering() function to predict the scattering parameters for the frequency and environment included in the model. This model allows for a very wide range of antennas and antenna arrays to be considered, but for simplicity only horn antennas will be included in this example. The simplest case would be a single source point and single receive point, rather than an aperture antenna such as a horn.

import numpy as np

Frequency and Mesh Resolution

freq = np.asarray(24.0e9)
wavelength = 3e8 / freq
mesh_resolution = 0.5 * wavelength

Setup transmitters and receivers

import lyceanem.geometry.targets as TL
import lyceanem.geometry.geometryfunctions as GF


transmit_horn_structure, transmitting_antenna_surface_coords = TL.meshedHorn(
    58e-3, 58e-3, 128e-3, 2e-3, 0.21, mesh_resolution
)

receive_horn_structure, receiving_antenna_surface_coords = TL.meshedHorn(
    58e-3, 58e-3, 128e-3, 2e-3, 0.21, mesh_resolution
)

Position Transmitter

rotate the transmitting antenna to the desired orientation, and then translate to final position. lyceanem.geometryfunctions.mesh_rotate() and lyceanem.geometryfunctions.mesh_translate() are used to achive this

rotation_vector1 = np.radians(np.asarray([90.0, 0.0, 0.0]))
rotation_vector2 = np.radians(np.asarray([0.0, 0.0, -90.0]))
rotation_vector3 = np.radians(np.asarray([0.0, 0.0, 90.0]))
transmit_horn_structure = GF.mesh_rotate(transmit_horn_structure, rotation_vector1)
transmit_horn_structure = GF.mesh_rotate(transmit_horn_structure, rotation_vector2)

transmit_horn_structure = GF.mesh_translate(
    transmit_horn_structure, np.asarray([2.529, 0, 0])
)

transmitting_antenna_surface_coords = GF.mesh_rotate(
    transmitting_antenna_surface_coords, rotation_vector1
)

transmitting_antenna_surface_coords = GF.mesh_rotate(
    transmitting_antenna_surface_coords, rotation_vector2
)

transmitting_antenna_surface_coords = GF.mesh_translate(
    transmitting_antenna_surface_coords, np.asarray([2.529, 0, 0])
)

Position Receiver

rotate the receiving horn to desired orientation and translate to final position.

receive_horn_structure = GF.mesh_rotate(receive_horn_structure, rotation_vector1)
# receive_horn_structure = GF.mesh_rotate(receive_horn_structure,rotation_vector3)
receive_horn_structure = GF.mesh_translate(
    receive_horn_structure, np.asarray([0, 1.609, 0])
)
receiving_antenna_surface_coords = GF.mesh_rotate(
    receiving_antenna_surface_coords, rotation_vector1
)
# receiving_antenna_surface_coords = GF.mesh_rotate(receiving_antenna_surface_coords,rotation_vector3)
receiving_antenna_surface_coords = GF.mesh_translate(
    receiving_antenna_surface_coords, np.asarray([0, 1.609, 0])
)

Create Scattering Plate

Create a Scattering plate a source of multipath reflections

reflectorplate, scatter_points = TL.meshedReflector(
    0.3, 0.3, 6e-3, wavelength * 0.5, sides="front"
)

position_vector = np.asarray([29e-3, 0.0, 0])
rotation_vector1 = np.radians(np.asarray([0.0, 90.0, 0.0]))
scatter_points = GF.mesh_rotate(scatter_points, rotation_vector1)
reflectorplate = GF.mesh_rotate(reflectorplate, rotation_vector1)
reflectorplate = GF.mesh_translate(reflectorplate, position_vector)
scatter_points = GF.mesh_translate(scatter_points, position_vector)

meshing reflector
args 0.3 0.3 0.006
majorsize 0.3
minorsize 0.3
thickness 0.006

Specify Reflection Angle

Rotate the scattering plate to the optimum angle for reflection from the transmitting to receiving horn

plate_orientation_angle = 45.0

rotation_vector = np.radians(np.asarray([0.0, 0.0, plate_orientation_angle]))
scatter_points = GF.mesh_rotate(scatter_points, rotation_vector)
reflectorplate = GF.mesh_rotate(reflectorplate, rotation_vector)

from lyceanem.base_classes import structures, points, antenna_structures

blockers = structures([reflectorplate, receive_horn_structure, transmit_horn_structure])
transmit_horn = antenna_structures(
    structures([transmit_horn_structure]), points([transmitting_antenna_surface_coords])
)

Visualise the Scene Geometry

import pyvista as pv

pl = pv.Plotter()
pl.add_mesh(pv.from_meshio(scatter_points), scalars="Area")
pl.add_mesh(pv.from_meshio(reflectorplate), color="grey")
pl.add_mesh(pv.from_meshio(transmitting_antenna_surface_coords), scalars="Area")
pl.add_mesh(pv.from_meshio(receiving_antenna_surface_coords), scalars="Area")
pl.add_mesh(pv.from_meshio(receive_horn_structure), color="green")
pl.add_mesh(pv.from_meshio(transmit_horn_structure), color="green")
pl.add_axes()
pl.show()

# Specify desired Transmit Polarisation
# --------------------------------------
# The transmit polarisation has a significant effect on the channel characteristics. In this example the transmit
# horn will be vertically polarised, (e-vector aligned with the y direction)

desired_E_axis = np.zeros((1, 3), dtype=np.float32)
desired_E_axis[0, 1] = 1.0

Static Scene

Interactive Scene

Frequency Domain Scattering

Once the arrangement of interest has been setup, lyceanem.models.frequency_domain.calculate_scattering() can be called, using raycasting to calculate the scattering parameters based upon the inputs. The scattering parameter determines how many reflections will be considered. A value of 0 would mean only line of sight contributions will be calculated, with 1 including single reflections, and 2 including double reflections as well.

import lyceanem.models.frequency_domain as FD

Ex, Ey, Ez = FD.calculate_scattering(
    aperture_coords=transmitting_antenna_surface_coords,
    sink_coords=receiving_antenna_surface_coords,
    antenna_solid=blockers,
    desired_E_axis=transmit_horn.excitation_function(
        desired_e_vector=desired_E_axis, transmit_power=0.25
    ),
    scatter_points=scatter_points,
    wavelength=wavelength,
    scattering=0,
    project_vectors=False,
    beta=(2 * np.pi) / wavelength,
)
Excuda, Eycuda, Ezcuda = FD.calculate_scattering(
    aperture_coords=transmitting_antenna_surface_coords,
    sink_coords=receiving_antenna_surface_coords,
    antenna_solid=blockers,
    desired_E_axis=transmit_horn.excitation_function(
        desired_e_vector=desired_E_axis, transmit_power=0.25
    ),
    scatter_points=scatter_points,
    wavelength=wavelength,
    scattering=0,
    project_vectors=False,
    beta=(2 * np.pi) / wavelength,
    cuda=True,
)
print("sumdiff", (np.sum((Ex - Excuda))))
print("sumdiff", (np.sum((Ey - Eycuda))))
print("sumdiff", (np.sum((Ez - Ezcuda))))

C:\Users\lycea\miniconda3\envs\CudaDevelopment\Lib\site-packages\lyceanem\electromagnetics\empropagation.py:3719: ComplexWarning: Casting complex values to real discards the imaginary part
  uvn_axes[2, :] = point_vector
C:\Users\lycea\miniconda3\envs\CudaDevelopment\Lib\site-packages\lyceanem\electromagnetics\empropagation.py:3736: ComplexWarning: Casting complex values to real discards the imaginary part
  uvn_axes[0, :] = np.cross(local_axes[2, :], point_vector) / np.linalg.norm(
C:\Users\lycea\miniconda3\envs\CudaDevelopment\Lib\site-packages\lyceanem\electromagnetics\empropagation.py:3758: ComplexWarning: Casting complex values to real discards the imaginary part
  uvn_axes[1, :] = np.cross(point_vector, uvn_axes[0, :]) / np.linalg.norm(
C:\Users\lycea\miniconda3\envs\CudaDevelopment\Lib\site-packages\numba_cuda\numba\cuda\dispatcher.py:693: NumbaPerformanceWarning: Grid size 40 will likely result in GPU under-utilization due to low occupancy.
  warn(NumbaPerformanceWarning(msg))
scatter_source_sink
sumdiff (-7.478353599761232e-08+4.723203539624922e-09j)
sumdiff (4.6496070964463726e-08-1.828645673300738e-09j)
sumdiff (-0.00036213591731276695-0.00097871361311673j)

Examine Scattering

The resultant scattering is decomposed into the Ex,Ey,Ez components at the receiving antenna, by itself this is not that interesting, so for this example we will rotate the reflector back, and then create a loop to step the reflector through different angles from 0 to 90 degrees in 1 degree steps.

angle_values = np.linspace(0, 90, 181)
angle_increment = np.diff(angle_values)[0]
responsex = np.zeros((len(angle_values)), dtype="complex")
responsey = np.zeros((len(angle_values)), dtype="complex")
responsez = np.zeros((len(angle_values)), dtype="complex")

plate_orientation_angle = -45.0

rotation_vector = np.radians(np.asarray([0.0, 0.0, plate_orientation_angle + 0.0]))
scatter_points = GF.mesh_rotate(scatter_points, rotation_vector)
reflectorplate = GF.mesh_rotate(reflectorplate, rotation_vector)


for angle_inc in range(len(angle_values)):
    rotation_vector = np.radians(np.asarray([0.0, 0.0, angle_values[angle_inc]]))
    scatter_points_temp = GF.mesh_rotate(scatter_points, rotation_vector)
    reflectorplate_temp = GF.mesh_rotate(reflectorplate, rotation_vector)
    blockers = structures(
        [reflectorplate_temp, receive_horn_structure, transmit_horn_structure]
    )

    # Scattered Path

    Ex, Ey, Ez = FD.calculate_scattering(
        aperture_coords=transmitting_antenna_surface_coords,
        sink_coords=scatter_points_temp,
        antenna_solid=blockers,
        desired_E_axis=transmit_horn.excitation_function(
            desired_e_vector=desired_E_axis, transmit_power=0.25
        ),
        scatter_points=scatter_points_temp,
        wavelength=wavelength,
        scattering=0,
        project_vectors=False,
        beta=(2 * np.pi) / wavelength,
    )
    scattered_field = np.array([Ex, Ey, Ez]).transpose()

    Ex2, Ey2, Ez2 = FD.calculate_scattering(
        aperture_coords=scatter_points_temp,
        sink_coords=receiving_antenna_surface_coords,
        antenna_solid=blockers,
        desired_E_axis=scattered_field,
        scatter_points=scatter_points_temp,
        wavelength=wavelength,
        scattering=0,
        project_vectors=False,
        beta=(2 * np.pi) / wavelength,
    )

    # Line of Sight Path

    Ex3, Ey3, Ez3 = FD.calculate_scattering(
        aperture_coords=transmitting_antenna_surface_coords,
        sink_coords=receiving_antenna_surface_coords,
        antenna_solid=blockers,
        desired_E_axis=transmit_horn.excitation_function(
            desired_e_vector=desired_E_axis, transmit_power=0.25
        ),
        scatter_points=scatter_points_temp,
        wavelength=wavelength,
        scattering=0,
        project_vectors=False,
        beta=(2 * np.pi) / wavelength,
    )
    responsex[angle_inc] = np.sum(
        (Ex2 + Ex3) * receiving_antenna_surface_coords.point_data["Area"]
    )
    responsey[angle_inc] = np.sum(
        (Ey2 + Ey3) * receiving_antenna_surface_coords.point_data["Area"]
    )
    responsez[angle_inc] = np.sum(
        (Ez2 + Ez3) * receiving_antenna_surface_coords.point_data["Area"]
    )

C:\Users\lycea\miniconda3\envs\CudaDevelopment\Lib\site-packages\numba_cuda\numba\cuda\dispatcher.py:693: NumbaPerformanceWarning: Grid size 107 will likely result in GPU under-utilization due to low occupancy.
  warn(NumbaPerformanceWarning(msg))

Plot Normalised Response

Using matplotlib, plot the results

import matplotlib.pyplot as plt

normalised_max = np.nanmax(
    np.array(
        [
            np.nanmax(20 * np.log10(np.abs(responsex))),
            np.nanmax(20 * np.log10(np.abs(responsey))),
            np.nanmax(20 * np.log10(np.abs(responsez))),
        ]
    )
)
ExdB = 20 * np.log10(np.abs(responsex)) - normalised_max
EydB = 20 * np.log10(np.abs(responsey)) - normalised_max
EzdB = 20 * np.log10(np.abs(responsez)) - normalised_max

fig, ax = plt.subplots()
ax.plot(angle_values - 45, ExdB, label="Ex")
ax.plot(angle_values - 45, EydB, label="Ey")
ax.plot(angle_values - 45, EzdB, label="Ez")
plt.xlabel("$\\theta_{N}$ (degrees)")
plt.ylabel("Normalised Level (dB)")
ax.set_ylim(-40.0, 0)
ax.set_xlim(np.min(angle_values) - 45, np.max(angle_values) - 45)
ax.set_xticks(np.linspace(np.min(angle_values) - 45, np.max(angle_values) - 45, 19))
ax.set_yticks(np.linspace(-40, 0.0, 21))
legend = ax.legend(loc="upper right", shadow=True)
plt.grid()
plt.show()