Source code for lyceanem.electromagnetics.beamforming

import cmath

import matplotlib.animation as animation
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d.art3d as art3d
import numpy as np
import pylab as pl
import scipy.stats
from matplotlib import cm
from matplotlib.patches import Wedge
from numba import cuda, float32, njit, prange
from scipy.spatial import distance

from ..raycasting import rayfunctions as RF
from ..utility import math_functions


[docs]def Steering_Efficiency( Dtheta, Dphi, Dtot, first_dimension_angle, second_dimension_angle, angular_coverage ): """ Calculate Steering Efficiency for the provided pattern, in radians Parameters ---------- Dtheta : numpy 2D array of floats or complex DESCRIPTION. Dphi : numpy 2D array of floats or complex DESCRIPTION. Dtot : numpy 2D array of floats or complex DESCRIPTION. angular coverage : float the total angular coverage to be considered, should be $4\pi$ steradians Returns ------- setheta : float steering efficiency in Dtheta sephi : float steering efficiency in Dphi setot : float steering efficiency in Dtotal """ with np.errstate(divide="ignore"): a_index = 10 * np.log10(np.abs(Dtheta)) >= ( 10 * np.log10(np.nanmax(np.abs(Dtheta))) - 3 ) b_index = 10 * np.log10(np.abs(Dphi)) >= ( 10 * np.log10(np.nanmax(np.abs(Dphi))) - 3 ) tot_index = 10 * np.log10(np.abs(Dtot)) >= ( 10 * np.log10(np.nanmax(np.abs(Dtot))) - 3 ) # setheta = (np.sum(a_index) / (Dtheta.shape[0] * Dtheta.shape[1])) * 100 # sephi = (np.sum(b_index) / (Dphi.shape[0] * Dphi.shape[1])) * 100 # setot = (np.sum(tot_index) / (Dtot.shape[0] * Dtot.shape[1])) * 100 setheta = ( (np.sum(a_index) * (first_dimension_angle * second_dimension_angle)) / angular_coverage ) * 100 sephi = ( (np.sum(b_index) * (first_dimension_angle * second_dimension_angle)) / angular_coverage ) * 100 setot = ( (np.sum(tot_index) * (first_dimension_angle * second_dimension_angle)) / angular_coverage ) * 100 return setheta, sephi, setot
[docs]@njit(cache=True, nogil=True) def WavefrontWeights(source_coords, steering_vector, wavelength): """ calculate the weights for a given set of element coordinates, wavelength, and steering vector (cartesian) """ weights = np.zeros((source_coords.shape[0]), dtype=np.complex64) # calculate distances of coords from steering_vector by using it to calculate arbitarily distant point # dist=distance.cdist(source_coords,(steering_vector*1e9)) # _,_,_,dist=calc_dv(source_coords,(steering_vector*1e9)) dist = np.zeros((source_coords.shape[0]), dtype=np.float32) dist = np.sqrt( np.abs( (source_coords[:, 0] - steering_vector.ravel()[0] * 1e9) ** 2 + (source_coords[:, 1] - steering_vector.ravel()[1] * 1e9) ** 2 + (source_coords[:, 2] - steering_vector.ravel()[2] * 1e9) ** 2 ) ) # RF.fast_calc_dv(source_coords,target,dv,dist) dist = dist - np.min(dist) # calculate required time delays, and then convert to phase delays delays = dist / scipy.constants.speed_of_light weights[:] = np.exp( 1j * 2 * np.pi * (scipy.constants.speed_of_light / wavelength) * delays ) return weights
def ArbitaryCoherenceWeights(source_coords, target_coord, wavelength): """ Generate Wavefront coherence weights based upon the desired wavelength and the coordinates of the target point """ weights = np.zeros((len(source_coords), 1), dtype=np.complex64) # calculate distances of coords from steering_vector by using it to calculate arbitarily distant point dist = distance.cdist(source_coords, target_coord) dist = dist - np.min(dist) # calculate required time delays, and then convert to phase delays delays = dist / scipy.constants.speed_of_light weights[:] = np.exp( 1j * 2 * np.pi * (scipy.constants.speed_of_light / wavelength) * delays ) return weights def TimeDelayWeights(source_coords, steering_vector, magnitude=1.0, maximum_delay=10): """ Generate time delay weights to focus on a target coordinate, with delays in nanoseconds """ weights = np.zeros((len(source_coords)), dtype=np.complex64) # calculate distances of coords from steering_vector by using it to calculate arbitarily distant point dist = np.sqrt( np.abs( (source_coords[:, 0] - steering_vector.ravel()[0] * 1e9) ** 2 + (source_coords[:, 1] - steering_vector.ravel()[1] * 1e9) ** 2 + (source_coords[:, 2] - steering_vector.ravel()[2] * 1e9) ** 2 ) ) dist = dist - np.min(dist) # dist is now relative distance from target vector, want to delay closest maximum, then have farthest be delayed zero # calculate required time delays, and then convert to nanoseconds, stored as a complex number with the magnitude weights delays = (dist / scipy.constants.speed_of_light) * 1e9 delays = np.max(delays) - delays weights[:] = magnitude + delays * 1j return weights def TimeDelayBeamform(excitation_function, weights, sampling_rate): """ The time delay beamform function takes an n by 2 or n by 3 array, and applies the supplied time delay weights to each by rolling each slice of the array by the required number of sampling intervals. Only positive time delays should be applied, but if positive and negative values are required for weighting, a constant value should be applied so that no delay is less than 0ns. Parameters ---------- excitation_function : float excitation function for time domain antenna array weights : complex magnitude is the real part, while the complex part is the time_delay in ns sampling_rate : float sampling rate in Hz, used to calculate the shifts to align the excitation function. Returns ------- beamformed_function """ # extract delays and convert to seconds, then integrer shifts delay = (sampling_rate * (np.imag(weights) * 1e-9)).astype(int).ravel() if len(excitation_function.shape) == 2: magnitudes = np.real(weights).reshape(-1, 1) elif len(excitation_function.shape) == 3: magnitudes = np.real(weights).reshape(-1, 1, 1) elif len(excitation_function.shape) == 4: magnitudes = np.real(weights).reshape(-1, 1, 1, 1) for row in range(excitation_function.shape[0]): excitation_function = shift_slice(excitation_function, row, delay[row]) excitation_function = excitation_function * magnitudes return excitation_function def shift_slice(array, row, shift): """ Parameters ---------- array : TYPE DESCRIPTION. row : TYPE DESCRIPTION. shift : TYPE DESCRIPTION. Returns ------- array : TYPE DESCRIPTION. """ if len(array.shape) == 2: array[row, :] = np.roll(array[row, :], shift) if shift > 0: array[row, :shift] = 0 elif shift < 0: array[row, -shift:] = 0 if len(array.shape) == 3: array[row, :, :] = np.roll(array[row, :, :], shift) if shift > 0: array[row, :, :shift] = 0 elif shift < 0: array[row, :, -shift:] = 0 if len(array.shape) == 4: array[row, :, :, :] = np.roll(array[row, :, :, :], shift) if shift > 0: array[row, :, :, :shift] = 0 elif shift < 0: array[row, :, :, -shift:] = 0 return array
[docs]@njit(cache=True, nogil=True) def EGCWeights( Etheta, Ephi, command_angles, polarization_switch="Etheta", az_range=np.linspace(-180.0, 180.0, 19), elev_range=np.linspace(-90.0, 90.0, 19), ): """ calculate the equal gain combining weights for a given set of element coordinates, wavelength, and command angles (az,elev) """ weights = np.zeros((Etheta.shape[0]), dtype=np.complex64) az_index = np.argmin(np.abs(az_range - command_angles[0])) elev_index = np.argmin(np.abs(elev_range - command_angles[1])) if polarization_switch == "Etheta": angle_vector = -np.angle(Etheta[:, elev_index, az_index].astype(np.complex64)) else: angle_vector = -np.angle(Ephi[:, elev_index, az_index].astype(np.complex64)) weights[:] = np.exp(1j * angle_vector) return weights
def OAMWeights(x, y, mode): """ generate OAM mode weights, based upon the radial angle of each element. """ # assumed array is x directed angles = np.arctan2(x, y) weights = np.zeros((len(x)), dtype=np.complex64) weights = np.exp(1j * mode * angles) return weights def OAMFourier( Ex, Ey, Ez, coordinates, prime_vector, mode_limit, first_dimension, second_dimension, coord_format="AzEl", ): """ producing mode index, mode coefficiencts, and mode probabilities with the co and crosspolar (Etheta,Ephi), but can probably be done the same for Ex,Ey,Ez """ # establish coordinate set, in this case theta and phi, but would work just as well with elevation and azimuth, assume that array is propagating in the z+ direction. if coord_format == "xyz": mode_index, mode_coefficients, mode_probabilites = OAMFourierCartesian( Ex, Ey, Ez, coordinates, mode_limit, first_dimension, second_dimension ) elif coord_format == "AzEl": mode_index, mode_coefficients, mode_probabilites = OAMFourierSpherical( Ex, Ey, Ez, coordinates, mode_limit, first_dimension, second_dimension ) return mode_index, mode_coefficients, mode_probabilites def OAMFourierCartesian(Ex, Ey, Ez, coordinates, mode_limit): """ assume propagation is in the Ez dimension """ mode_index = np.linspace(-mode_limit, mode_limit, mode_limit * 2 + 1) mode_coefficients = np.zeros((mode_index.shape[0], 3), dtype=np.complex64) az, el, r = math_functions.cart2sph( coordinates[:, 0], coordinates[:, 1], coordinates[:, 2] ) # a coefficient of mode m, at angle theta is defined in terms of the # integral of the phi dimension in the range 0 to 2pi. for oam_m in range(len(mode_index)): mode_coefficients[oam_m, 0] = (1 / (2 * np.pi)) * np.sum( Ex.ravel() * np.exp(1j * mode_index[oam_m] * az), axis=0 ) mode_coefficients[oam_m, 1] = (1 / (2 * np.pi)) * np.sum( Ey.ravel() * np.exp(1j * mode_index[oam_m] * az), axis=0 ) mode_coefficients[oam_m, 2] = (1 / (2 * np.pi)) * np.sum( Ez.ravel() * np.exp(1j * mode_index[oam_m] * az), axis=0 ) powers = np.sum(np.abs(mode_coefficients**2), axis=0) mode_probabilities = np.zeros((mode_index.shape[0], 3), dtype=np.float32) mode_probabilities[:, 0] = np.abs( (1 / np.sum(powers)) * (mode_coefficients[:, 0] ** 2) ) mode_probabilities[:, 1] = np.abs( (1 / np.sum(powers)) * (mode_coefficients[:, 1] ** 2) ) mode_probabilities[:, 2] = np.abs( (1 / np.sum(powers)) * (mode_coefficients[:, 2] ** 2) ) return mode_index, mode_coefficients, mode_probabilities def OAMFourierSpherical(Ex, Ey, Ez, coordinates, mode_limit, az_range, elev_range): """ assume probagation is in the Ez dimension """ mode_index = np.linspace(-mode_limit, mode_limit, mode_limit * 2 + 1) mode_coefficients = np.zeros((mode_index.shape[0], len(elev_range), 3)) # a coefficient of mode m, at angle theta is defined in terms of the # integral of the azimuth dimension in the range -pi to pi. for oam_m in range(mode_index): # mode_coefficients(oam_m,:,1)=(1/(2*pi))*sum(CoPolar(:,1:theta_limit).*exp(-sqrt(-1)*mode_index(oam_m)*p'*ones(1,theta_limit))); # mode_coefficients(oam_m,:,2)=(1/(2*pi))*sum(CrossPolar(:,1:theta_limit).*exp(-sqrt(-1)*mode_index(oam_m)*p'*ones(1,theta_limit))); mode_coefficients[oam_m, :, 0] = (1 / (2 * np.pi)) * np.sum( Ex[:, :] * np.exp(-1j * mode_index[oam_m]), axis=0 ) # copolar_power=sum(sum(abs(mode_coefficients(:,:,1)).^2)); # crosspolar_power=sum(sum(abs(mode_coefficients(:,:,2)).^2)); mode_probabilities = np.zeros((mode_index.shape[0], 2), dtype=np.float32) # mode_prob(:,1)=(1/sum([copolar_power,crosspolar_power]))*sum(abs(mode_coefficients(:,:,1)').^2); # mode_prob(:,2)=(1/sum([copolar_power,crosspolar_power]))*sum(abs(mode_coefficients(:,:,2)').^2); return mode_index, mode_coefficients, mode_probabilities # @njit(parallel=True,cache=True, nogil=True)
[docs]def MaximumDirectivityMap( Etheta, Ephi, source_coords, wavelength, az_range=np.linspace(-180.0, 180.0, 19), elev_range=np.linspace(-90.0, 90.0, 19), az_index=None, elev_index=None, forming="Total", total_solid_angle=(4 * np.pi), phase_resolution=24, ): """ Uses wavefront beamsteering, and equal gain combining algorithms to steer the antenna array to each possible command angle in the farfield, mapping out the maximum achieved directivity at the command angle for each command angle set. Parameters ---------- Etheta : 3D numpy array The Etheta polarisation farfield patterns, arranged in terms of the number of elements, azimuth resolution, and elevation resolution Ephi : 3D numpy array The Ephi polarisation farfield patterns, arranged in terms of the number of elements, azimuth resolution, and elevation resolution source_coords : :class:`open3d.geometry.PointCloud` The source coordinates of each element, corresponding to the order of element patterns in Etheta and Ephi. Units should be m wavelength : float The wavelength of interest az_range : 1D numpy array of float The azimuth values for the farfield mesh, arranged from smallest to largest elev_range : 1D numpy array of float The elevation values for the farfield mesh, arranged from smallest to largest az_index : 1D array of int optional parameter, can be specified as an index of the azimuth values of interest via indexing, defaults to [None], which ensures all values are covered elev_index : 1D array of int optional parameter, can be specified as an index of the elevation values of interest via indexing, defaults to [None], which ensures all values are covered forming : str Which polarisation should be beamformed, the default is [Total], beamforming the total directivity pattern, avoiding issues with elements which have a strongly $E\\theta$ or $E\\phi$ pattern. This can also be set to [Etheta] or [Ephi] total_solid_angle : float the total solid angle covered by the farfield patterns, this defaults to $4\\pi$ for a full spherical pattern phase_resolution : int the desired phase resolution of the beamforming architecture in bits. Default is [24], which means no practical truncation will occur. If beam mapping at a single resolution is required, then this can be set between 2 and 24. If multiple values are required, it may be more efficient to use :func:`lyceanem.electromagnetics.beamforming.MaximumDirectivityMapDiscrete`, which allows a 1D array of resolutions to be supplied, and produces a maximum directivity map for each. Returns ------- directivity_map : 3D numpy array of float The achieved maximum directivity map. At each point the directivity corresponds to the achieved directivity at that command angle. Arranged as elev axis, azimuth axis, Dtheta,Dphi,Dtot """ if elev_index == None: # if no elev index is provided then generate for all possible values (assumes every elevation point is of interest) elev_index = np.linspace(0, len(elev_range) - 1, len(elev_range)).astype(int) if az_index == None: # if no az index is provided then generate for all possible values (assumes every azimuth point is of interest) az_index = np.linspace(0, len(az_range) - 1, len(az_range)).astype(int) az_res = len(az_index) elev_res = len(elev_index) source_points = np.asarray(source_coords.points) directivity_map = np.zeros((elev_res, az_res, 3)) command_angles = np.zeros((2), dtype=np.float32) for az_inc in range(az_res): for elev_inc in range(elev_res): command_angles[0] = az_range[az_index[az_inc]] command_angles[1] = elev_range[elev_index[elev_inc]] steering_vector = np.zeros((1, 3)) ( steering_vector[0, 0], steering_vector[0, 1], steering_vector[0, 2], ) = math_functions.sph2cart( np.radians(command_angles[0]), np.radians(command_angles[1]), 1 ) WS_weights = WavefrontWeights(source_points, steering_vector, wavelength) EGC_weights = EGCWeights( Etheta, Ephi, command_angles, az_range=az_range, elev_range=elev_range ) EGC_weights2 = EGCWeights( Etheta, Ephi, command_angles, az_range=az_range, elev_range=elev_range, polarization_switch="Ephi", ) if phase_resolution <= 12: WS_weights = WeightTruncation(WS_weights, phase_resolution) EGC_weights = WeightTruncation(EGC_weights, phase_resolution) EGC_weights2 = WeightTruncation(EGC_weights2, phase_resolution) Ethetabeamformed = np.sum( EGC_weights.reshape(Etheta.shape[0], 1, 1) * Etheta, axis=0 ) Ephibeamformed = np.sum( EGC_weights.reshape(Etheta.shape[0], 1, 1) * Ephi, axis=0 ) Ethetabeamformed2 = np.sum( EGC_weights2.reshape(Etheta.shape[0], 1, 1) * Etheta, axis=0 ) Ephibeamformed2 = np.sum( EGC_weights2.reshape(Etheta.shape[0], 1, 1) * Ephi, axis=0 ) Ethetabeamformed3 = np.sum( WS_weights.reshape(Etheta.shape[0], 1, 1) * Etheta, axis=0 ) Ephibeamformed3 = np.sum( WS_weights.reshape(Etheta.shape[0], 1, 1) * Ephi, axis=0 ) Dtheta, Dphi, Dtot, _ = directivity_transform( Ethetabeamformed, Ephibeamformed, az_range=az_range, elev_range=elev_range, ) Dtheta2, Dphi2, Dtot2, _ = directivity_transform( Ethetabeamformed2, Ephibeamformed2, az_range=az_range, elev_range=elev_range, total_solid_angle=total_solid_angle, ) Dtheta3, Dphi3, Dtot3, _ = directivity_transform( Ethetabeamformed3, Ephibeamformed3, az_range=az_range, elev_range=elev_range, total_solid_angle=total_solid_angle, ) if forming == "Total": comparitor = np.asarray( [ Dtot[elev_inc, az_inc], Dtot2[elev_inc, az_inc], Dtot3[elev_inc, az_inc], ] ) elif forming == "Etheta": comparitor = np.asarray( [ Dtheta[elev_inc, az_inc], Dtheta2[elev_inc, az_inc], Dtheta3[elev_inc, az_inc], ] ) elif forming == "Ephi": comparitor = np.asarray( [ Dphi[elev_inc, az_inc], Dphi2[elev_inc, az_inc], Dphi3[elev_inc, az_inc], ] ) if np.any(np.isnan(comparitor)): print("error") best_index = np.where(comparitor == np.max(comparitor))[0][0] if best_index == 0: directivity_map[elev_inc, az_inc, 0] = Dtheta[elev_inc, az_inc] directivity_map[elev_inc, az_inc, 1] = Dphi[elev_inc, az_inc] directivity_map[elev_inc, az_inc, 2] = Dtot[elev_inc, az_inc] elif best_index == 1: directivity_map[elev_inc, az_inc, 0] = Dtheta2[elev_inc, az_inc] directivity_map[elev_inc, az_inc, 1] = Dphi2[elev_inc, az_inc] directivity_map[elev_inc, az_inc, 2] = Dtot2[elev_inc, az_inc] elif best_index == 2: directivity_map[elev_inc, az_inc, 0] = Dtheta3[elev_inc, az_inc] directivity_map[elev_inc, az_inc, 1] = Dphi3[elev_inc, az_inc] directivity_map[elev_inc, az_inc, 2] = Dtot3[elev_inc, az_inc] return directivity_map
# @njit(parallel=False, cache=True, nogil=True)
[docs]def MaximumDirectivityMapDiscrete( Etheta, Ephi, source_coords, wavelength, az_range=np.linspace(-180.0, 180.0, 19), elev_range=np.linspace(-180.0, 180.0, 19), az_index=None, elev_index=None, forming="Total", total_solid_angle=(4 * np.pi), phase_resolution=np.asarray([24]), ): """ Uses wavefront beamsteering, and equal gain combining algorithms to steer the antenna array to each possible command angle in the farfield, mapping out the maximum achieved directivity at the command angle for each command angle set. Parameters ---------- Etheta : 3D numpy array The $E\theta$ polarisation farfield patterns, arranged in terms of the number of elements, azimuth resolution, and elevation resolution Ephi : 3D numpy array The $E\phi$ polarisation farfield patterns, arranged in terms of the number of elements, azimuth resolution, and elevation resolution source_coords : :class:`open3d.geometry.PointCloud` The source coordinates of each element, corresponding to the order of element patterns in $E\theta$ and $E\phi$. Units should be m wavelength : float The wavelength of interest az_res : int Azimuth resolution elev_res : int Elevation resolution az_range : 1D numpy array of float The azimuth values for the farfield mesh, arranged from smallest to largest elev_range : 1D numpy array of float The elevation values for the farfield mesh, arranged from smallest to largest forming : str Which polarisation should be beamformed, the default is [Total], beamforming the total directivity pattern, avoiding issues with elements which have a strongly $E\theta$ or $E\phi$ pattern. This can also be set to [Etheta] or [Ephi] total_solid_angle : float the total solid angle covered by the farfield patterns, this defaults to $4\pi$ for a full spherical pattern phase_resolution : 1D numpy array of int the desired phase resolution of the beamforming architecture in bits. Default is [24], which means no practical truncation will occur. If beam mapping is desired at a single resolution is required, then this can be set between 2 and 24, if more than one resolution value is required, then a 1D array of values can be specified. resolutions to be supplied, and produces a maximum directivity map for each. Returns ------- directivity_map : 4D numpy array of float The achieved maximum directivity map (a 3D numpy array for each phase resolution) . For each phase resolution for each point the directivity corresponds to the achieved directivity at that command angle. """ if elev_index == None: # if no elev index is provided then generate for all possible values (assumes every elevation point is of interest) elev_index = np.linspace(0, len(elev_range) - 1, len(elev_range)).astype(int) if az_index == None: # if no az index is provided then generate for all possible values (assumes every azimuth point is of interest) az_index = np.linspace(0, len(az_range) - 1, len(az_range)).astype(int) az_res = len(az_index) elev_res = len(elev_index) source_points = np.asarray(source_coords.points) directivity_map = np.zeros( (elev_res, az_res, 3, phase_resolution.shape[0]), dtype=np.float32 ) command_angles = np.zeros((2), dtype=np.float32) for az_inc in range(az_res): for elev_inc in range(elev_res): inc_res = 0 for res_inc in range(phase_resolution.shape[0]): resolution = phase_resolution[res_inc] command_angles[0] = az_range[az_inc] command_angles[1] = elev_range[elev_inc] steering_vector = np.zeros((1, 3)) ( steering_vector[0, 0], steering_vector[0, 1], steering_vector[0, 2], ) = math_functions.sph2cart( np.radians(command_angles[0]), np.radians(command_angles[1]), 1 ) WS_weights = WavefrontWeights( source_points, steering_vector, wavelength ) EGC_weights = EGCWeights( Etheta, Ephi, command_angles, az_range=az_range, elev_range=elev_range, ) EGC_weights2 = EGCWeights( Etheta, Ephi, command_angles, az_range=az_range, elev_range=elev_range, polarization_switch="Ephi", ) WS_weights = WeightTruncation(WS_weights, resolution) EGC_weights = WeightTruncation(EGC_weights, resolution) EGC_weights2 = WeightTruncation(EGC_weights2, resolution) Ethetabeamformed = np.sum( EGC_weights.reshape(Etheta.shape[0], 1, 1) * Etheta, axis=0 ) Ephibeamformed = np.sum( EGC_weights.reshape(Etheta.shape[0], 1, 1) * Ephi, axis=0 ) Ethetabeamformed2 = np.sum( EGC_weights2.reshape(Etheta.shape[0], 1, 1) * Etheta, axis=0 ) Ephibeamformed2 = np.sum( EGC_weights2.reshape(Etheta.shape[0], 1, 1) * Ephi, axis=0 ) Ethetabeamformed3 = np.sum( WS_weights.reshape(Etheta.shape[0], 1, 1) * Etheta, axis=0 ) Ephibeamformed3 = np.sum( WS_weights.reshape(Etheta.shape[0], 1, 1) * Ephi, axis=0 ) Dtheta, Dphi, Dtot, _ = directivity_transform( Ethetabeamformed, Ephibeamformed, az_range=az_range, elev_range=elev_range, ) Dtheta2, Dphi2, Dtot2, _ = directivity_transform( Ethetabeamformed2, Ephibeamformed2, az_range=az_range, elev_range=elev_range, total_solid_angle=total_solid_angle, ) Dtheta3, Dphi3, Dtot3, _ = directivity_transform( Ethetabeamformed3, Ephibeamformed3, az_range=az_range, elev_range=elev_range, total_solid_angle=total_solid_angle, ) if forming == "Total": comparitor = np.asarray( [ Dtot[elev_inc, az_inc], Dtot2[elev_inc, az_inc], Dtot3[elev_inc, az_inc], ] ) elif forming == "Etheta": comparitor = np.asarray( [ Dtheta[elev_inc, az_inc], Dtheta2[elev_inc, az_inc], Dtheta3[elev_inc, az_inc], ] ) elif forming == "Ephi": comparitor = np.asarray( [ Dphi[elev_inc, az_inc], Dphi2[elev_inc, az_inc], Dphi3[elev_inc, az_inc], ] ) if np.any(np.isnan(comparitor)): print("error") best_index = np.where(comparitor == np.max(comparitor))[0][0] if best_index == 0: directivity_map[elev_inc, az_inc, 0, inc_res] = Dtheta[ elev_inc, az_inc ] directivity_map[elev_inc, az_inc, 1, inc_res] = Dphi[ elev_inc, az_inc ] directivity_map[elev_inc, az_inc, 2, inc_res] = Dtot[ elev_inc, az_inc ] elif best_index == 1: directivity_map[elev_inc, az_inc, 0, inc_res] = Dtheta2[ elev_inc, az_inc ] directivity_map[elev_inc, az_inc, 1, inc_res] = Dphi2[ elev_inc, az_inc ] directivity_map[elev_inc, az_inc, 2, inc_res] = Dtot2[ elev_inc, az_inc ] elif best_index == 2: directivity_map[elev_inc, az_inc, 0, inc_res] = Dtheta3[ elev_inc, az_inc ] directivity_map[elev_inc, az_inc, 1, inc_res] = Dphi3[ elev_inc, az_inc ] directivity_map[elev_inc, az_inc, 2, inc_res] = Dtot3[ elev_inc, az_inc ] inc_res += 1 return directivity_map
@njit(cache=True, nogil=True) def MaximumfieldMapDiscrete( Etheta, Ephi, source_coords, wavelength, az_res, elev_res, az_range=np.linspace(-180.0, 180.0, 19), elev_range=np.linspace(-180.0, 180.0, 19), forming="Total", total_solid_angle=(4 * np.pi), phase_resolution=[24], ): efield_map = np.zeros( (elev_res, az_res, 3, len(phase_resolution)), dtype=np.complex64 ) command_angles = np.zeros((2), dtype=np.float32) for az_inc in prange(az_res): for elev_inc in range(elev_res): inc_res = 0 for resolution in phase_resolution: command_angles[0] = az_range[az_inc] command_angles[1] = elev_range[elev_inc] steering_vector = np.zeros((1, 3)) ( steering_vector[0, 0], steering_vector[0, 1], steering_vector[0, 2], ) = math_functions.sph2cart( np.radians(command_angles[0]), np.radians(command_angles[1]), 1 ) WS_weights = WavefrontWeights( source_coords, steering_vector, wavelength ) EGC_weights = EGCWeights( Etheta, Ephi, command_angles, az_range=az_range, elev_range=elev_range, ) EGC_weights2 = EGCWeights( Etheta, Ephi, command_angles, az_range=az_range, elev_range=elev_range, polarization_switch="Ephi", ) WS_weights = WeightTruncation(WS_weights, resolution) EGC_weights = WeightTruncation(EGC_weights, resolution) EGC_weights2 = WeightTruncation(EGC_weights2, resolution) Ethetabeamformed = np.sum( EGC_weights.reshape(Etheta.shape[0], 1, 1) * Etheta, axis=0 ) Ephibeamformed = np.sum( EGC_weights.reshape(Etheta.shape[0], 1, 1) * Ephi, axis=0 ) Ethetabeamformed2 = np.sum( EGC_weights2.reshape(Etheta.shape[0], 1, 1) * Etheta, axis=0 ) Ephibeamformed2 = np.sum( EGC_weights2.reshape(Etheta.shape[0], 1, 1) * Ephi, axis=0 ) Ethetabeamformed3 = np.sum( WS_weights.reshape(Etheta.shape[0], 1, 1) * Etheta, axis=0 ) Ephibeamformed3 = np.sum( WS_weights.reshape(Etheta.shape[0], 1, 1) * Ephi, axis=0 ) Etotal = Ethetabeamformed**2 + Ephibeamformed**2 Etotal2 = Ethetabeamformed2**2 + Ephibeamformed2**2 Etotal3 = Ethetabeamformed3**2 + Ephibeamformed3**2 # Dtheta,Dphi,Dtot,_=directivity_transform(Ethetabeamformed,Ephibeamformed,az_range=az_range,elev_range=elev_range) # Dtheta2,Dphi2,Dtot2,_=directivity_transform(Ethetabeamformed2,Ephibeamformed2,az_range=az_range,elev_range=elev_range,total_solid_angle=total_solid_angle) # Dtheta3,Dphi3,Dtot3,_=directivity_transform(Ethetabeamformed3,Ephibeamformed3,az_range=az_range,elev_range=elev_range,total_solid_angle=total_solid_angle) if forming == "Total": comparitor = np.asarray( [ Etotal[elev_inc, az_inc], Etotal2[elev_inc, az_inc], Etotal3[elev_inc, az_inc], ] ) if forming == "Etheta": comparitor = np.asarray( [ Ethetabeamformed[elev_inc, az_inc], Ethetabeamformed2[elev_inc, az_inc], Ethetabeamformed3[elev_inc, az_inc], ] ) elif forming == "Ephi": comparitor = np.asarray( [ Ephibeamformed[elev_inc, az_inc], Ephibeamformed2[elev_inc, az_inc], Ephibeamformed3[elev_inc, az_inc], ] ) if np.any(np.isnan(comparitor)): print("error") best_index = np.where(comparitor == np.max(comparitor))[0][0] if best_index == 0: efield_map[elev_inc, az_inc, 0, inc_res] = Ethetabeamformed[ elev_inc, az_inc ] efield_map[elev_inc, az_inc, 1, inc_res] = Ephibeamformed[ elev_inc, az_inc ] efield_map[elev_inc, az_inc, 2, inc_res] = Etotal[elev_inc, az_inc] elif best_index == 1: efield_map[elev_inc, az_inc, 0, inc_res] = Ethetabeamformed2[ elev_inc, az_inc ] efield_map[elev_inc, az_inc, 1, inc_res] = Ephibeamformed2[ elev_inc, az_inc ] efield_map[elev_inc, az_inc, 2, inc_res] = Etotal2[elev_inc, az_inc] elif best_index == 2: efield_map[elev_inc, az_inc, 0, inc_res] = Ethetabeamformed3[ elev_inc, az_inc ] efield_map[elev_inc, az_inc, 1, inc_res] = Ephibeamformed3[ elev_inc, az_inc ] efield_map[elev_inc, az_inc, 2, inc_res] = Etotal3[elev_inc, az_inc] inc_res += 1 return efield_map @njit(cache=True, nogil=True) def directivity_transform( Etheta, Ephi, az_range=np.linspace(-180.0, 180.0, 19), elev_range=np.linspace(-90.0, 90.0, 19), total_solid_angle=(4 * np.pi), ): # transform Etheta and Ephi data into antenna directivity # directivity is defined in terms of the power radiated in a specific direction, over the average radiated power # power per unit solid angle Dmax = np.zeros((3), dtype=np.float32) Umax = np.zeros((3), dtype=np.float32) Utheta = np.abs(Etheta**2) Uphi = np.abs(Ephi**2) Utotal = np.abs(Etheta**2) + np.abs(Ephi**2) power_sum = 0.0 temp_vector = np.zeros((4), dtype=np.float32) temp_vector2 = np.zeros((4), dtype=np.float32) for elinc in range(len(elev_range) - 1): for azinc in range(len(az_range) - 1): temp_vector[0] = Utheta[elinc, azinc] temp_vector[1] = Utheta[elinc + 1, azinc] temp_vector[2] = Utheta[elinc, azinc + 1] temp_vector[3] = Utheta[elinc + 1, azinc + 1] temp_vector2[0] = Uphi[elinc, azinc] temp_vector2[1] = Uphi[elinc + 1, azinc] temp_vector2[2] = Uphi[elinc, azinc + 1] temp_vector2[3] = Uphi[elinc + 1, azinc + 1] r1 = np.mean(temp_vector) r2 = np.mean(temp_vector2) power_sum = ( power_sum + r1 * np.abs( np.sin(np.radians(az_range[azinc] + az_range[azinc + 1]) / 2.0) ) + r2 * np.abs( np.sin(np.radians(az_range[azinc] + az_range[azinc + 1]) / 2.0) ) ) Umax[0] = np.nanmax(Utheta) Umax[1] = np.nanmax(Uphi) Umax[2] = np.nanmax(Utotal) Uav = ( power_sum * ( (np.radians(az_range[1] - az_range[0])) * (np.radians(elev_range[1] - elev_range[0])) ) / total_solid_angle ) Dmax = Umax / Uav Dtheta = Utheta / Uav Dphi = Uphi / Uav Dtot = Utotal / Uav return Dtheta, Dphi, Dtot, Dmax @njit(cache=True, nogil=True) def WeightTruncation(weights, resolution): quant = 2**resolution levels = np.zeros((weights.shape), dtype=np.complex64) # shift numpy complex angles from -pi to pi, into 0 to 2pi, then quantise for `quant' levels levels = np.round( ((quant - 1) * ((np.angle(weights) + np.pi) / (2 * np.pi))), 0, levels ) / (quant - 1) new_weights = np.abs(weights) * np.exp(1j * ((levels * 2 * np.pi) - np.pi)) return new_weights def AnimatedPlot( dimension1, dimension2, dimension3, data, pattern_min=-40.0, pattern_max=0.0, dimension1label=None, dimension2label=None, colorbarlabel=None, title_text=None, ticknum=9, fps=5, save_location=None, ): def animate(i): # title_text='Sampled Phase {:.2f} Wavelengths from the Transmitting Antenna'.format(dimension3[i]) ax.clear() ax.contourf( dimension1, dimension2, data[:, :, i], levels, cmap="viridis", origin=origin ) ax.contour( dimension1, dimension2, data[:, :, i], levels, colors=("k",), origin=origin ) ax.set_xlim([np.min(dimension1), np.max(dimension1)]) ax.set_ylim([np.min(dimension2), np.max(dimension2)]) ax.grid() # ax.set_xticks(np.linspace(-180, 180, 9)) # ax.set_yticks(np.linspace(-90, 90.0, 13)) if dimension1label != None: ax.set_xlabel(dimension1label) if dimension2label != None: ax.set_ylabel(dimension2label) if title_text != None: ax.set_title(title_text.format(dimension3[i])) plt.rcParams["figure.figsize"] = [7.00, 3.50] plt.rcParams["figure.autolayout"] = True ticknum = 9 fig, ax = plt.subplots(constrained_layout=True) origin = "lower" # pattern_min = -np.pi / 2 # pattern_max = np.pi / 2 levels = np.linspace(pattern_min, pattern_max, 73) CS = ax.contourf( dimension1, dimension2, data[:, :, 0], levels, cmap="viridis", origin=origin ) cbar = fig.colorbar(CS, ticks=np.linspace(pattern_min, pattern_max, ticknum)) # bar_label="Sampled Phase Angle (Radians)" c_label_values = np.linspace(pattern_min, pattern_max, ticknum) c_labels = np.char.mod("%.2f", c_label_values) cbar.set_ticklabels(c_labels.tolist()) ax.set_xlim([np.min(dimension1), np.max(dimension1)]) ax.set_ylim([np.min(dimension2), np.max(dimension2)]) # ax.set_xticks(np.linspace(-180, 180, 9)) # ax.set_yticks(np.linspace(-90, 90.0, 13)) ax.set_xlabel("x ($\lambda$)") ax.set_ylabel("y ($\lambda$)") # setup for 3dB contours # contournum = np.ceil((pattern_max - pattern_min) / 3).astype(int) # levels2 = np.linspace(-contournum * 3, plot_max, contournum + 1) # title_text=None ax.grid() if dimension1label != None: ax.set_xlabel(dimension1label) if dimension2label != None: ax.set_ylabel(dimension2label) if colorbarlabel != None: cbar.ax.set_ylabel(colorbarlabel) if title_text != None: ax.set_title(title_text.format(dimension3[0])) ani = animation.FuncAnimation(fig, animate, 100, interval=50, blit=False) plt.show() if save_location != None: # f = r"C:/Users/lycea/Documents/10-19 Research Projects/farfieldanimation.gif" writergif = animation.PillowWriter(fps=fps) ani.save(save_location, writer=writergif)
[docs]def PatternPlot( data, az, elev, pattern_min=-40, plot_max=0.0, plottype="Polar", logtype="amplitude", ticknum=6, title_text=None, ): """ Plot the relavent 3D data in relative power (dB) or normalised directivity (dBi) Parameters ----------- data : 2D array of floats or complex the data to plot az : 2D array of floats the azimuth angles for each datapoint in [data] in degrees elev : 2D array of floats the elevation angles for each datapoint in [data] in degrees pattern_min : float the desired scale minimum in dB, default is [-40] plot_max : float the desired scale maximum in dB, default is [0] plottype : str the plot type, either [Polar], [Cartesian-Surf], or [Contour]. The default is [Polar] logtype : str the type of data being considered, either [amplitude] or [power], to ensure the correct logarithm is used, default is [amplitude] ticknum : int the number of ticks on the colorbar, default is [6] title_text : str the graph title, defaults to [None] Returns -------- None """ # condition data data = np.abs(data) # calculate log profile if logtype == "power": logdata = 10 * np.log10(data) bar_label = "Relative Power (dB)" else: logdata = 20 * np.log10(data) if plot_max == 0.0: logdata -= np.nanmax(logdata) bar_label = "Normalised Directivity (dBi)" else: bar_label = "Directivity (dBi)" logdata[logdata <= pattern_min] = pattern_min if plottype == "Polar": norm_log = (logdata - pattern_min) / np.abs(pattern_min) sinks = np.zeros((len(np.ravel(az)), 3), dtype=np.float32) sinks[:, 0], sinks[:, 1], sinks[:, 2] = RF.azeltocart( np.ravel(az), np.ravel(elev), np.ravel(norm_log) ) dist = np.sqrt( sinks[:, 0].reshape(az.shape) ** 2 + sinks[:, 1].reshape(az.shape) ** 2 + sinks[:, 2].reshape(az.shape) ** 2 ) dist_max = np.max(dist) my_col = cm.viridis(dist / dist_max) fig = plt.figure() ax = fig.add_subplot(111, projection="3d") V = np.array([[1.1, 0, 0], [0, 1.1, 0], [0, 0, 1.1]], dtype=np.float32) origin = np.zeros((3, 3), dtype=np.float32) # origin point offset = np.array([0.8, 0.8, 0.8], dtype=np.float32).reshape(1, 3) ax.quiver( origin[0, 0] - offset, origin[0, 0] - offset, origin[0, 0] - offset, V[0, 0], V[1, 0], V[2, 0], color=["red"], ) ax.quiver( origin[0, 1] - offset, origin[0, 1] - offset, origin[0, 1] - offset, V[0, 1], V[1, 1], V[2, 1], color=["green"], ) ax.quiver( origin[0, 2] - offset, origin[0, 2] - offset, origin[0, 2] - offset, V[0, 2], V[1, 2], V[2, 2], color=["blue"], ) plot_handle = ax.plot_surface( sinks[:, 0].reshape(az.shape), sinks[:, 1].reshape(az.shape), sinks[:, 2].reshape(az.shape), facecolors=my_col, linewidth=0, antialiased=False, clim=[0, 1], ) ax.set_xlim([-1, 1]) ax.set_ylim([-1, 1]) ax.set_zlim([-1, 1]) plt.axis("off") # plot_handle.set_clim([0,1]) cbar = fig.colorbar( plot_handle, ticks=np.linspace(0, 1.0, ticknum), extend="both" ) cbar.ax.set_ylabel(bar_label) c_labels = np.linspace(pattern_min, plot_max, ticknum).astype("str") cbar.set_ticklabels(c_labels.tolist()) p = Wedge((0, 0), 1.01, 0, 360, width=0.0001, color="gray") ax.add_patch(p) art3d.pathpatch_2d_to_3d(p, z=0, zdir="x") p = Wedge((0, 0), 1.01, 0, 360, width=0.0001, color="gray") ax.add_patch(p) art3d.pathpatch_2d_to_3d(p, z=0, zdir="y") p = Wedge((0, 0), 1.01, 0, 360, width=0.0001, color="gray") ax.add_patch(p) art3d.pathpatch_2d_to_3d(p, z=0, zdir="z") ax.view_init(elev=45.0, azim=-45) if title_text != None: ax.set_title(title_text) elif plottype == "Cartesian-Surf": fig = plt.figure() ax = fig.add_subplot(projection="3d") ax.plot_surface(az, elev, logdata, cmap="viridis", edgecolor="none") ax.set_xlim([np.min(az), np.max(az)]) ax.set_ylim([np.min(elev), np.max(elev)]) ax.set_zlim([pattern_min, plot_max]) ax.set_zticks(np.linspace(pattern_min, plot_max, ticknum)) ax.set_xticks(np.linspace(-180, 180, 9)) ax.set_yticks(np.linspace(-90, 90.0, 5)) ax.set_xlabel("Azimuth (degrees)") ax.set_ylabel("Elevation (degrees)") ax.set_zlabel(bar_label) if title_text != None: ax.set_title(title_text) elif plottype == "Contour": fig, ax = plt.subplots(constrained_layout=True) origin = "lower" levels = np.linspace(pattern_min, plot_max, ticknum * 10) CS = ax.contourf(az, elev, logdata, levels, cmap="viridis", origin=origin) cbar = fig.colorbar(CS, ticks=np.linspace(pattern_min, plot_max, ticknum)) cbar.ax.set_ylabel(bar_label) c_labels = np.linspace(pattern_min, plot_max, ticknum).astype("str") cbar.set_ticklabels(c_labels.tolist()) ax.set_xlim([np.min(az), np.max(az)]) ax.set_ylim([np.min(elev), np.max(elev)]) ax.set_xticks(np.linspace(-180, 180, 9)) ax.set_yticks(np.linspace(-90, 90.0, 13)) ax.set_xlabel("Azimuth (degrees)") ax.set_ylabel("Elevation (degrees)") # setup for 3dB contours contournum = np.ceil((plot_max - pattern_min) / 3).astype(int) levels2 = np.linspace(-contournum * 3, plot_max, contournum + 1) if pattern_min < -40: line_spec_width = 0.5 else: line_spec_width = 1 CS4 = ax.contour( az, elev, logdata, levels2, colors=("k",), linewidths=(line_spec_width,), origin=origin, ) ax.grid() if title_text != None: ax.set_title(title_text) plt.show()
def PatternPlot2D( data, az, pattern_min=-40, plot_max=0.0, logtype="amplitude", ticknum=6, line_labels=None, title_text=None, fontsize=16, ): # condition data data = np.abs(data) if data.ndim > 1: # multi line plot, condition data, as second axis should be the number of different lines. multiline = True num_lines = data.shape[1] else: multiline = False # calculate log profile if logtype == "power": logdata = 10 * np.log10(data) bar_label = "Relative Power (dB)" else: logdata = 20 * np.log10(data) if plot_max == 0.0: logdata -= np.nanmax(logdata) bar_label = "Normalised Directivity (dBi)" else: bar_label = "Directivity (dBi)" logdata[logdata <= pattern_min] = pattern_min tick_marks = np.linspace(pattern_min, plot_max, ticknum) az_rad = np.radians(az) pl.rcParams.update({"font.size": fontsize}) fig, ax = plt.subplots(subplot_kw={"projection": "polar"}) if multiline == True: for line in range(num_lines): if not (line_labels == None): ax.plot(az_rad, logdata[:, line], label=line_labels[line]) else: ax.plot(az_rad, logdata[:, line]) else: ax.plot(az_rad, logdata) ax.set_rmax(plot_max) ax.set_rticks(tick_marks, fontsize=fontsize) # Less radial ticks ax.set_rlabel_position(-22.5) # Move radial labels away from plotted line ax.grid(True) if not (line_labels == None): # legend position legend_angle = np.deg2rad(30) ax.legend( loc="lower left", bbox_to_anchor=( 0.5 + np.cos(legend_angle) / 2, 0.5 + np.sin(legend_angle) / 2, ), fontsize=fontsize, ) if not (title_text == None): ax.set_title(title_text, va="bottom", fontsize=fontsize) label_angle = np.deg2rad(280) ax.text(label_angle, plot_max * 1.2, bar_label) plt.show() # noinspection PyTypeChecker @cuda.jit(device=True) def point_directivity(Ea, Eb, az_range, el_range, interest_index): """ compute the directivity at the point of interest in the farfield pattern """ average_power = 0.0 directivity_results = cuda.local.array(shape=(3), dtype=float32) return directivity_results @cuda.jit(device=True) def EqualGainCombiningGPU(SteeringPattern, CommandIndex, weights): """ equal gain combining algorithm based on the provided steering pattern and command index """ weights = cmath.exp(-1j * np.angle(SteeringPattern[:, CommandIndex])) return weights @cuda.jit def GPUBeamformingMap(Etheta, Ephi, DirectivityMap, az_range, el_range, wavelength): """ """ az_inc, el_inc = cuda.grid(2) if az_inc < az_range.shape[0] and el_inc < el_range.shape[0]: DirectivityMap[az_inc, el_inc] = 0