diff --git a/docs/examples/01_aperture_projection.py b/docs/examples/01_aperture_projection.py deleted file mode 100644 index 9799548..0000000 --- a/docs/examples/01_aperture_projection.py +++ /dev/null @@ -1,166 +0,0 @@ -#!/usr/bin/env python -# coding: utf-8 - -""" -Calculating Antenna Array Performance Envelope using Aperture Projection -========================================================================== -This is a demonstration using the aperture projection function in the context of a conformal antenna array mounted upon -an unmanned aerial vehicle. - -Aperture Projection as a technique is based upon Hannan's formulation of the gain of an aperture based upon its surface -area and the freuqency of interest. This is defined in terms of the maximum gain :math:`G_{max}`, the effective area of -the aperture :math:`A_{e}`, and the wavelength of interest :math:`\lambda`. - -.. math:: - G_{max}=\dfrac{4 \pi A_{e}}{\lambda^{2}} - -While this has been in common use since the 70s, as a formula it is limited to planar surfaces, and only providing the -maximum gain in the boresight direction for that surface. - -Aperture projection as a function is based upon the rectilinear projection of the aperture into the farfield. This can -then be used with Hannan's formula to predict the maximum achievable directivity for all farfield directions of -interest. - -As this method is built into a raytracing environment, the maximum performance for an aperture on any platform can also -be predicted using the :func:`lyceanem.models.frequency_domain.aperture_projection` function. - -""" -import copy - -import numpy as np - -# %% -# Setting Farfield Resolution and Wavelength -# ------------------------------------------- -# LyceanEM uses Elevation and Azimuth to record spherical coordinates, ranging from -180 to 180 degrees in azimuth, -# and from -90 to 90 degrees in elevation. In order to launch the aperture projection function, the resolution in -# both azimuth and elevation is requried. -# In order to ensure a fast example, 37 points have been used here for both, giving a total of 1369 farfield points. -# -# The wavelength of interest is also an important variable for antenna array analysis, so we set it now for 10GHz, -# an X band aperture. - -az_res = 37 -elev_res = 37 -wavelength = 3e8 / 10e9 - -# %% -# Geometries -# ------------------------ -# In order to make things easy to start, an example geometry has been included within LyceanEM for a UAV, and the -# open3d trianglemesh structures can be accessed by importing the data subpackage -import lyceanem.tests.reflectordata as data - -body, array, _ = data.exampleUAV(10e9) - - - -# %% -## .. image:: ../_static/open3d_structure.png - -# crop the inner surface of the array trianglemesh (not strictly required, as the UAV main body provides blocking to -# the hidden surfaces, but correctly an aperture will only have an outer face. -surface_array = copy.deepcopy(array) -surface_array.cells[0].data = np.asarray(array.cells[0].data)[: (array.cells[0].data).shape[0] // 2, :] - -surface_array.cell_data["normals"] = np.array(array.cell_data["normals"])[: (array.cells[0].data).shape[0] // 2] - -# %% -# Structures -# -------------- -# LyceanEM uses a class named 'structures' to store and maniuplate joined 3D solids. Currently all that is implemented -# is the class itself, and methods to allow translation and rotation of the trianglemesh solids. A structure can be -# passed to the models to provide the environment to be considered as blockers. -# structures are created by calling the class, and passing it a list of the open3d trianglemesh structures to be added. -from lyceanem.base_classes import structures - -blockers = structures([body]) - -# %% -# Aperture Projection -# ----------------------- -# Aperture Projection is imported from the frequency domain models, requiring the aperture of interest, wavelength to -# be considered, and the azimuth and elevation ranges. The function then returns the directivity envelope as a numpy -# array of floats, and an open3d point cloud with points and colors corresponding to the directivity envelope of the -# provided aperture, scaling from yellow at maximum to dark purple at minimum. -from lyceanem.models.frequency_domain import aperture_projection - -directivity_envelope, pcd = aperture_projection( - surface_array, - environment=blockers, - wavelength=wavelength, - az_range=np.linspace(-180.0, 180.0, az_res), - elev_range=np.linspace(-90.0, 90.0, elev_res), -) -# %% -# Open3D Visualisation -# ------------------------ -# The resultant maximum directivity envelope is provided as both a numpy array of directivities for each angle, but -# also as an open3d point cloud. This allows easy visualisation using :func:`open3d.visualization.draw_geometries`. -# %% - - -# %% -# .. image:: ../_static/open3d_results_rendering.png - - -# Maximum Directivity -print( - "Maximum Directivity of {:3.1f} dBi".format( - np.max(10 * np.log10(directivity_envelope)) - ) -) - -# %% -# Plotting the Output -# ------------------------ -# While the open3d visualisation is very intuitive for examining the results of the aperture projection, it is -# difficult to consider the full 3D space, and cannot be included in documentation in this form. However, matplotlib -# can be used to generate contour plots with 3dB contours to give a more systematic understanding of the resultant -# maximum directivity envelope. - -import matplotlib.pyplot as plt - -# set directivity limits on the closest multiple of 5 -plot_max = ((np.ceil(np.nanmax(10 * np.log10(directivity_envelope))) // 5.0) + 1) * 5 -azmesh, elevmesh = np.meshgrid( - np.linspace(-180.0, 180.0, az_res), np.linspace(-90, 90, elev_res) -) -fig, ax = plt.subplots(constrained_layout=True) -origin = "lower" - -levels = np.linspace(plot_max - 40, plot_max, 81) -CS = ax.contourf( - azmesh, - elevmesh, - 10 * np.log10(directivity_envelope), - levels, - origin=origin, - extend="both", -) -cbar = fig.colorbar(CS) -cbar.ax.set_ylabel("Directivity (dBi)") -cbar.set_ticks(np.linspace(plot_max - 40, plot_max, 9)) -cbar.ax.set_yticklabels(np.linspace(plot_max - 40, plot_max, 9).astype("str")) -levels2 = np.linspace( - np.nanmax(10 * np.log10(directivity_envelope)) - 60, - np.nanmax(10 * np.log10(directivity_envelope)), - 21, -) -CS4 = ax.contour( - azmesh, - elevmesh, - 10 * np.log10(directivity_envelope), - levels2, - colors=("k",), - linewidths=(2,), - origin=origin, -) -ax.set_ylim(-90, 90) -ax.set_xlim(-180.0, 180) -ax.set_xticks(np.linspace(-180, 180, 13)) -ax.set_yticks(np.linspace(-90, 90, 13)) -ax.set_xlabel("Azimuth (degrees)") -ax.set_ylabel("Elevation (degrees)") -ax.set_title("Maximum Directivity Envelope") -fig.show() diff --git a/docs/examples/02_coherently_polarised_array.py b/docs/examples/02_coherently_polarised_array.py deleted file mode 100644 index 4a85de2..0000000 --- a/docs/examples/02_coherently_polarised_array.py +++ /dev/null @@ -1,120 +0,0 @@ -#!/usr/bin/env python -# coding: utf-8 -""" -Modelling a Coherently Polarised Aperture -====================================================== - -This example uses the frequency domain :func:`lyceanem.models.frequency_domain.calculate_farfield` function to predict -the farfield pattern for a linearly polarised aperture. This could represent an antenna array without any beamforming -weights. - - -""" -import copy - -import numpy as np -import meshio -# %% -# Setting Farfield Resolution and Wavelength -# ------------------------------------------- -# LyceanEM uses Elevation and Azimuth to record spherical coordinates, ranging from -180 to 180 degrees in azimuth, -# and from -90 to 90 degrees in elevation. In order to launch the aperture projection function, the resolution in -# both azimuth and elevation is requried. -# In order to ensure a fast example, 37 points have been used here for both, giving a total of 1369 farfield points. -# -# The wavelength of interest is also an important variable for antenna array analysis, so we set it now for 10GHz, -# an X band aperture. - -az_res = 181 -elev_res = 181 -wavelength = 3e8 / 10e9 - -# %% -# Geometries -# ------------------------ -# In order to make things easy to start, an example geometry has been included within LyceanEM for a UAV, and the -# :class:`open3d.geometry.TriangleMesh` structures can be accessed by importing the data subpackage -import lyceanem.tests.reflectordata as data - -body, array, source_coords = data.exampleUAV(10e9) - - - -# %% -## .. image:: ../_static/open3d_structure.png - -# crop the inner surface of the array trianglemesh (not strictly required, as the UAV main body provides blocking to -# the hidden surfaces, but correctly an aperture will only have an outer face. -surface_array = copy.deepcopy(array) -surface_array.cells[0].data = np.asarray(array.cells[0].data)[: (array.cells[0].data).shape[0] // 2, :] - -surface_array.cell_data["normals"] = np.array(array.cell_data["normals"])[: (array.cells[0].data).shape[0] // 2] - -from lyceanem.base_classes import structures - -blockers = structures([body, array]) - -from lyceanem.models.frequency_domain import calculate_farfield - - - - -import pyvista as pv - - -source_points = surface_array.points - - -# %% -# .. image:: ../_static/sourcecloudfromshapeuav.png - -# %% -# Drawbacks of :func:`lyceanem.geometry.geometryfunctions.sourcecloudfromshape` -# ------------------------------------------------------------------------------ -# As can be seen by comparing the two source point sets, :func:`lyceanem.geometry.geometryfunctions.sourcecloudfromshape` -# has a significant drawback when used for complex sharply curved antenna arrays, as the poisson disk sampling method -# does not produce consistently spaced results. - -desired_E_axis = np.zeros((1, 3), dtype=np.float32) -desired_E_axis[0, 2] = 1.0 - -Etheta, Ephi = calculate_farfield( - source_coords, - blockers, - desired_E_axis, - az_range=np.linspace(-180, 180, az_res), - el_range=np.linspace(-90, 90, elev_res), - wavelength=wavelength, - farfield_distance=20, - project_vectors=True, -) - -# %% -# Storing and Manipulating Antenna Patterns -# --------------------------------------------- -# The resultant antenna pattern can be stored in :class:`lyceanem.base.antenna_pattern` as it has been modelled as one -# distributed aperture, the advantage of this class is the integrated display, conversion and export functions. It is -# very simple to define, and save the pattern, and then display with a call -# to :func:`lyceanem.base.antenna_pattern.display_pattern`. This produces 3D polar plots which can be manipulated to -# give a better view of the whole pattern, but if contour plots are required, then this can also be produced by passing -# plottype='Contour' to the function. - -from lyceanem.base_classes import antenna_pattern - -UAV_Static_Pattern = antenna_pattern( - azimuth_resolution=az_res, elevation_resolution=elev_res -) -UAV_Static_Pattern.pattern[:, :, 0] = Etheta -UAV_Static_Pattern.pattern[:, :, 0] = Ephi - -UAV_Static_Pattern.display_pattern() - -# %% -# .. image:: ../_static/sphx_glr_02_coherently_polarised_array_001.png -# .. image:: ../_static/sphx_glr_02_coherently_polarised_array_002.png - -UAV_Static_Pattern.display_pattern(plottype="Contour") - -# %% -# .. image:: ../_static/sphx_glr_02_coherently_polarised_array_003.png -# .. image:: ../_static/sphx_glr_02_coherently_polarised_array_004.png diff --git a/docs/examples/04_time_domain_channel_modelling.py b/docs/examples/04_time_domain_channel_modelling.py deleted file mode 100644 index ac98bb3..0000000 --- a/docs/examples/04_time_domain_channel_modelling.py +++ /dev/null @@ -1,303 +0,0 @@ -#!/usr/bin/env python -# coding: utf-8 -""" -Modelling a Physical Channel in the Time Domain -====================================================== - -This example uses the frequency domain :func:`lyceanem.models.time_domain.calculate_scattering` function to -predict the time domain response for a given excitation signal 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 -import meshio -# %% -# Frequency and Mesh Resolution -# ------------------------------ -# -sampling_freq = 60e9 -model_time = 1e-7 -num_samples = int(model_time * (sampling_freq)) - -# simulate receiver noise -bandwidth = 8e9 -kb = 1.38065e-23 -receiver_impedence = 50 -thermal_noise_power = 4 * kb * 293.15 * receiver_impedence * bandwidth -noise_power = -80 # dbw -mean_noise = 0 - -model_freq = 16e9 -wavelength = 3e8 / model_freq - -# %% -# 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, wavelength * 0.5 -) -receive_horn_structure, receiving_antenna_surface_coords = TL.meshedHorn( - 58e-3, 58e-3, 128e-3, 2e-3, 0.21, wavelength * 0.5 -) - -# %% -# Position Transmitter -# ---------------------- -# rotate the transmitting antenna to the desired orientation, and then translate to final position. -# :func:`lyceanem.geometry.geometryfunctions.open3drotate` allows both the center of rotation to be defined, and -# ensures the right syntax is used for Open3d, as it was changed from 0.9.0 to 0.10.0 and onwards. -# -rotation_vector1 = np.radians(np.asarray([90.0, 0.0, 0.0])) -rotation_vector2 = 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.translate_mesh(transmit_horn_structure, np.asarray([2.695, 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.translate_mesh(transmitting_antenna_surface_coords, np.asarray([2.695, 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.translate_mesh(receive_horn_structure, np.asarray([0, 1.427, 0])) -receiving_antenna_surface_coords = GF.mesh_rotate(receiving_antenna_surface_coords, rotation_vector1) -receiving_antenna_surface_coords = GF.translate_mesh(receiving_antenna_surface_coords, np.asarray([0, 1.427, 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.translate_mesh(reflectorplate, position_vector) -scatter_points = GF.translate_mesh(scatter_points, position_vector) - - -# %% -# 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 - -blockers = structures([reflectorplate, receive_horn_structure, transmit_horn_structure]) - -# %% -# Visualise the Scene Geometry -# ------------------------------ -# Use open3d function :func:`open3d.visualization.draw_geometries` to visualise the scene and ensure that all the -# relavent sources and scatter points are correct. Point normal vectors can be displayed by pressing 'n' while the -# window is open. - - -# %% -# .. image:: ../_static/03_frequency_domain_channel_model_picture_01.png - -# %% -# 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 z direction) - -desired_E_axis = np.zeros((1, 3), dtype=np.float32) -desired_E_axis[0, 1] = 1.0 - -# %% -# Time Domain Scattering -# ---------------------------- -# -import scipy.signal as sig -import lyceanem.models.time_domain as TD -from lyceanem.base_classes import structures - - -angle_values = np.linspace(0, 90, 91) -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 + angle_increment]) -) -scatter_points = GF.mesh_rotate(scatter_points, rotation_vector) -reflectorplate = GF.mesh_rotate(reflectorplate, rotation_vector) - -from tqdm import tqdm - -wake_times = np.zeros((len(angle_values))) -Ex = np.zeros((len(angle_values), num_samples)) -Ey = np.zeros((len(angle_values), num_samples)) -Ez = np.zeros((len(angle_values), num_samples)) - -for angle_inc in tqdm(range(len(angle_values))): - rotation_vector = np.radians(np.asarray([0.0, 0.0, angle_increment])) - scatter_points = GF.mesh_rotate(scatter_points, rotation_vector) - reflectorplate = GF.mesh_rotate(reflectorplate, rotation_vector) - blockers = structures( - [reflectorplate, transmit_horn_structure, receive_horn_structure] - ) - pulse_time = 5e-9 - output_power = 0.01 # dBwatts - powerdbm = 10 * np.log10(output_power) + 30 - v_transmit = ((10 ** (powerdbm / 20)) * receiver_impedence) ** 0.5 - output_amplitude_rms = v_transmit / (1 / np.sqrt(2)) - output_amplitude_peak = v_transmit - - desired_E_axis = np.zeros((3), dtype=np.float32) - desired_E_axis[2] = 1.0 - noise_volts_peak = (10 ** (noise_power / 10) * receiver_impedence) * 0.5 - - excitation_signal = output_amplitude_rms * sig.chirp( - np.linspace(0, pulse_time, int(pulse_time * sampling_freq)), - model_freq - bandwidth, - pulse_time, - model_freq, - method="linear", - phi=0, - vertex_zero=True, - ) + np.random.normal(mean_noise, noise_volts_peak, int(pulse_time * sampling_freq)) - ( - Ex[angle_inc, :], - Ey[angle_inc, :], - Ez[angle_inc, :], - wake_times[angle_inc], - ) = TD.calculate_scattering( - transmitting_antenna_surface_coords, - receiving_antenna_surface_coords, - excitation_signal, - blockers, - desired_E_axis, - scatter_points=scatter_points, - wavelength=wavelength, - scattering=1, - elements=False, - sampling_freq=sampling_freq, - num_samples=num_samples, - ) - - noise_volts = np.random.normal(mean_noise, noise_volts_peak, num_samples) - Ex[angle_inc, :] = Ex[angle_inc, :] + noise_volts - Ey[angle_inc, :] = Ey[angle_inc, :] + noise_volts - Ez[angle_inc, :] = Ez[angle_inc, :] + noise_volts - - -# %% -# Plot Normalised Response -# ---------------------------- -# Using matplotlib, plot the results - - -import matplotlib.pyplot as plt - -time_index = np.linspace(0, model_time * 1e9, num_samples) -time, anglegrid = np.meshgrid(time_index[:1801], angle_values - 45) -norm_max = np.nanmax( - np.array( - [ - np.nanmax(10 * np.log10((Ex ** 2) / receiver_impedence)), - np.nanmax(10 * np.log10((Ey ** 2) / receiver_impedence)), - np.nanmax(10 * np.log10((Ez ** 2) / receiver_impedence)), - ] - ) -) - -fig2, ax2 = plt.subplots(constrained_layout=True) -origin = "lower" -# Now make a contour plot with the levels specified, -# and with the colormap generated automatically from a list -# of colors. - -levels = np.linspace(-80, 0, 41) - -CS = ax2.contourf( - anglegrid, - time, - 10 * np.log10((Ez[:, :1801] ** 2) / receiver_impedence) - norm_max, - levels, - origin=origin, - extend="both", -) -cbar = fig2.colorbar(CS) -cbar.ax.set_ylabel("Received Power (dBm)") - -ax2.set_ylim(0, 30) -ax2.set_xlim(-45, 45) - -ax2.set_xticks(np.linspace(-45, 45, 7)) -ax2.set_yticks(np.linspace(0, 30, 16)) - -ax2.set_xlabel("Rotation Angle (degrees)") -ax2.set_ylabel("Time of Flight (ns)") -ax2.set_title("Received Power vs Time for rotating Plate (24GHz)") - -from scipy.fft import fft, fftfreq -import scipy - -xf = fftfreq(len(time_index), 1 / sampling_freq)[: len(time_index) // 2] -input_signal = excitation_signal * (output_amplitude_peak) -inputfft = fft(input_signal) -input_freq = fftfreq(120, 1 / sampling_freq)[:60] -freqfuncabs = scipy.interpolate.interp1d(input_freq, np.abs(inputfft[:60])) -freqfuncangle = scipy.interpolate.interp1d(input_freq, np.angle(inputfft[:60])) -newinput = freqfuncabs(xf[1600]) * np.exp(freqfuncangle(xf[1600])) -Exf = fft(Ex) -Eyf = fft(Ey) -Ezf = fft(Ez) - -# %% -# .. image:: ../_static/sphx_glr_04_time_domain_channel_modelling_001.png - -# %% -# Frequency Specific Results -# ------------------------------- -# The time of flight plot is useful to displaying the output of the model, giving a understanding about what is -# physically happening in the channel, but to get an idea of the behaviour in the frequency domain we need to use a -# fourier transform to move from time and voltages to frequency. - -s21x = 20 * np.log10(np.abs(Exf[:, 1600] / newinput)) -s21y = 20 * np.log10(np.abs(Eyf[:, 1600] / newinput)) -s21z = 20 * np.log10(np.abs(Ezf[:, 1600] / newinput)) -tdangles = np.linspace(-45, 45, 91) -fig, ax = plt.subplots() -ax.plot(tdangles, s21x - np.max(s21z), label="Ex") -ax.plot(tdangles, s21y - np.max(s21z), label="Ey") -ax.plot(tdangles, s21z - np.max(s21z), label="Ez") -plt.xlabel("$\\theta_{N}$ (degrees)") -plt.ylabel("Normalised Level (dB)") -ax.set_ylim(-60.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(-60, 0.0, 21)) -legend = ax.legend(loc="upper right", shadow=True) -plt.grid() -plt.title("$S_{21}$ at 16GHz") -plt.show() - -# %% -# .. image:: ../_static/sphx_glr_04_time_domain_channel_modelling_002.png diff --git a/docs/examples/05_array_beamforming.py b/docs/examples/05_array_beamforming.py deleted file mode 100644 index afc705f..0000000 --- a/docs/examples/05_array_beamforming.py +++ /dev/null @@ -1,113 +0,0 @@ -#!/usr/bin/env python -# coding: utf-8 -""" -Array Beamforming -====================================================== - -This example uses the frequency domain :func:`lyceanem.models.frequency_domain.calculate_farfield` function to predict -the farfield patterns for a linearly polarised aperture with multiple elements. This is then beamformed to all farfield points using multiple open loop beamforming algorithms to attemp to 'map' out the acheivable beamforming for the antenna array using :func:`lyceanem.electromagnetics.beamforming.MaximumDirectivityMap`. - -The Steering Efficiency can then be evaluated using :func:`lyceanem.electromagnetics.beamforming.Steering_Efficiency` for the resultant achieved beamforming. - - -""" -import numpy as np - -# %% -# Setting Farfield Resolution and Wavelength -# ------------------------------------------- -# LyceanEM uses Elevation and Azimuth to record spherical coordinates, ranging from -180 to 180 degrees in azimuth, -# and from -90 to 90 degrees in elevation. In order to launch the aperture projection function, the resolution in -# both azimuth and elevation is requried. -# In order to ensure a fast example, 37 points have been used here for both, giving a total of 1369 farfield points. -# -# The wavelength of interest is also an important variable for antenna array analysis, so we set it now for 10GHz, -# an X band aperture. - -az_res = 181 -elev_res = 37 -wavelength = 3e8 / 10e9 - -# %% -# Geometries -# ------------------------ -# In order to make things easy to start, an example geometry has been included within LyceanEM for a UAV, and the -# :class:`open3d.geometry.TriangleMesh` structures can be accessed by importing the data subpackage -import lyceanem.tests.reflectordata as data - -body, array, source_coords = data.exampleUAV(10e9) - - - -# %% -# .. image:: ../_static/UAVArraywithPoints.png - - -from lyceanem.base_classes import structures - -blockers = structures([body, array]) - -# %% -# Model Farfield Array Patterns -# ------------------------------- -# The same function is used to predict the farfield pattern of each element in the array, but the variable 'elements' -# is set as True, instructing the function to return the antenna patterns as 3D arrays arranged with axes element, -# elevation points, and azimuth points. These can then be beamformed using the desired beamforming algorithm. LyceanEM -# currently includes two open loop algorithms for phase weights :func:`lyceanem.electromagnetics.beamforming.EGCWeights`, -# and :func:`lyceanem.electromagnetics.beamforming.WavefrontWeights` -from lyceanem.models.frequency_domain import calculate_farfield - -desired_E_axis = np.zeros((1, 3), dtype=np.float32) -desired_E_axis[0, 2] = 1.0 - -Etheta, Ephi = calculate_farfield( - source_coords, - blockers, - desired_E_axis, - az_range=np.linspace(-180, 180, az_res), - el_range=np.linspace(-90, 90, elev_res), - wavelength=wavelength, - farfield_distance=20, - elements=True, - project_vectors=True, -) - - -from lyceanem.electromagnetics.beamforming import MaximumDirectivityMap - -az_range = np.linspace(-180, 180, az_res) -el_range = np.linspace(-90, 90, elev_res) -directivity_map = MaximumDirectivityMap( - Etheta, Ephi, source_coords, wavelength, az_range, el_range -) - -from lyceanem.electromagnetics.beamforming import PatternPlot - -az_mesh, elev_mesh = np.meshgrid(az_range, el_range) - -PatternPlot( - directivity_map[:, :, 2], az_mesh, elev_mesh, logtype="power", plottype="Contour" -) - -# %% -# .. image:: ../_static/sphx_glr_05_array_beamforming_001.png - -from lyceanem.electromagnetics.beamforming import Steering_Efficiency - -setheta, sephi, setot = Steering_Efficiency( - directivity_map[:, :, 0], - directivity_map[:, :, 1], - directivity_map[:, :, 2], - np.radians(np.diff(el_range)[0]), - np.radians(np.diff(az_range)[0]), - 4 * np.pi, -) - -print("Steering Effciency of {:3.1f}%".format(setot)) - - -print( - "Maximum Directivity of {:3.1f} dBi".format( - np.max(10 * np.log10(directivity_map[:, :, 2])) - ) -) diff --git a/docs/examples/06_farfield_polarisation.py b/docs/examples/06_farfield_polarisation.py index e799ac9..a8c5cdd 100644 --- a/docs/examples/06_farfield_polarisation.py +++ b/docs/examples/06_farfield_polarisation.py @@ -11,8 +11,7 @@ """ import numpy as np -import open3d as o3d - +import meshio # %% # Setting Farfield Resolution and Wavelength # ------------------------------------------- @@ -34,14 +33,13 @@ from lyceanem.base_classes import points,structures,antenna_structures -aperture_coords=o3d.geometry.PointCloud() point1=np.asarray([0.0,0,0]).reshape(1,3) normal1=np.asarray([0,0,1.0]).reshape(1,3) -aperture_coords.points=o3d.utility.Vector3dVector(point1) -aperture_coords.normals=o3d.utility.Vector3dVector(normal1) +aperture_coords = meshio.Mesh(points=point1, cells=[]) +aperture_coords.point_data['normals']=normal1 aperture=points([aperture_coords]) blockers=structures([None]) -point_antenna=antenna_structures(blockers, aperture) +point_antenna=antenna_structures(None, aperture) from lyceanem.models.frequency_domain import calculate_farfield @@ -52,8 +50,8 @@ desired_E_axis = np.zeros((1, 3), dtype=np.complex64) desired_E_axis[0, 0] = 1.0 Etheta, Ephi = calculate_farfield( - point_antenna.export_all_points(), - point_antenna.export_all_structures(), + aperture_coords, + None, point_antenna.excitation_function(desired_e_vector=desired_E_axis), az_range=np.linspace(-180, 180, az_res), el_range=np.linspace(-90, 90, elev_res), @@ -83,8 +81,8 @@ desired_E_axis = np.zeros((1, 3), dtype=np.complex64) desired_E_axis[0, 1] = 1.0 Etheta, Ephi = calculate_farfield( - point_antenna.export_all_points(), - point_antenna.export_all_structures(), + aperture.export_all_points(), + blockers.export_all_structures(), point_antenna.excitation_function(desired_e_vector=desired_E_axis), az_range=np.linspace(-180, 180, az_res), el_range=np.linspace(-90, 90, elev_res), @@ -108,8 +106,8 @@ desired_E_axis = np.zeros((1, 3), dtype=np.complex64) desired_E_axis[0, 2] = 1.0 Etheta, Ephi = calculate_farfield( - point_antenna.export_all_points(), - point_antenna.export_all_structures(), + aperture.export_all_points(), + blockers.export_all_structures(), point_antenna.excitation_function(desired_e_vector=desired_E_axis), az_range=np.linspace(-180, 180, az_res), el_range=np.linspace(-90, 90, elev_res), @@ -134,7 +132,7 @@ desired_E_axis = np.zeros((1, 3), dtype=np.complex64) desired_E_axis[0, 0] = 1.0 Etheta, Ephi = calculate_farfield( - point_antenna.export_all_points(), + aperture.export_all_points(), point_antenna.export_all_structures(), point_antenna.excitation_function(desired_e_vector=desired_E_axis), az_range=np.linspace(-180, 180, az_res), @@ -152,7 +150,7 @@ desired_E_axis = np.zeros((1, 3), dtype=np.complex64) desired_E_axis[0, 1] = 1.0 Etheta, Ephi = calculate_farfield( - point_antenna.export_all_points(), + aperture.export_all_points(), point_antenna.export_all_structures(), point_antenna.excitation_function(desired_e_vector=desired_E_axis), az_range=np.linspace(-180, 180, az_res), @@ -170,7 +168,7 @@ desired_E_axis = np.zeros((1, 3), dtype=np.complex64) desired_E_axis[0, 2] = 1.0 Etheta, Ephi = calculate_farfield( - point_antenna.export_all_points(), + aperture.export_all_points(), point_antenna.export_all_structures(), point_antenna.excitation_function(desired_e_vector=desired_E_axis), az_range=np.linspace(-180, 180, az_res), diff --git a/lyceanem/base_classes.py b/lyceanem/base_classes.py index f41342d..8c7dceb 100644 --- a/lyceanem/base_classes.py +++ b/lyceanem/base_classes.py @@ -178,15 +178,17 @@ def export_points(self, point_index=None): if item == 0: points = np.array(self.points[item].points) points = np.append(points, self.points[item].points, axis=0) - point_data = np.array(len(self.points[0].point_data), points.shape[0]) + point_data = {} for key in self.points[0].point_data.keys(): - point_data[key] = np.array(points.shape[0]) - pointssofar = 0 for item in range(len(self.points)): - point_data_element = np.array(self.point_data[item][key]) - point_data[key][pointssofar:(pointssofar+point_data_element.shape[0])] = point_data_element - pointssofar+= point_data_element.shape[0] - combinded_points = meshio.Mesh(points, point_data=point_data) + + point_data_element = np.array(self.points[item].point_data[key]) + + if item == 0: + point_data[key] = point_data_element + else: + point_data[key] = np.append(point_data[key], point_data_element, axis=0) + combinded_points = meshio.Mesh(points, cells = [], point_data=point_data) combinded_points = GF.mesh_transform(combinded_points, self.pose, False) return combinded_points @@ -200,15 +202,15 @@ def export_points(self, point_index=None): points = np.append(points, self.points[item].points, axis=0) point_data = {} for key in self.points[0].point_data.keys(): - point_data[key] = np.array(points.shape[0]) - pointssofar = 0 for item in point_index: - point_data_element = np.array(self.point_data[item][key]) - point_data[key][pointssofar:(pointssofar+point_data_element.shape[0])] = point_data_element - pointssofar+= point_data_element.shape[0] + point_data_element = np.array(self.points[item].point_data[key]) + if item == 0: + point_data[key] = point_data_element + else: + point_data[key] = np.append(point_data[key], point_data_element, axis=0) - combinded_points = meshio.Mesh(points, point_data=point_data) + combinded_points = meshio.Mesh(points,cells = [], point_data=point_data) combinded_points = GF.mesh_transform(combinded_points, self.pose, False)