From cf70603693fb76ac369e21522c2060cfc9c358d3 Mon Sep 17 00:00:00 2001 From: LyceanEM <60020395+LyceanEM@users.noreply.github.com> Date: Tue, 3 Dec 2024 20:10:22 +0000 Subject: [PATCH] Updates to all examples to ensure changes to calculate farfield are working as expected. Bug fixes to the triangle export function in the base class, and to the emrotation function within the mesh transform.. --- .../examples/02_coherently_polarised_array.py | 4 +- .../03_frequency_domain_channel_modelling.py | 104 ++++-- docs/examples/05_array_beamforming.py | 3 +- docs/examples/06_farfield_polarisation.py | 122 +++---- .../07_aperture_farfield_polarisation.py | 24 +- lyceanem/base_classes.py | 4 +- lyceanem/electromagnetics/emfunctions.py | 6 +- lyceanem/geometry/targets.py | 46 ++- lyceanem/models/frequency_domain.py | 332 ++---------------- 9 files changed, 215 insertions(+), 430 deletions(-) diff --git a/docs/examples/02_coherently_polarised_array.py b/docs/examples/02_coherently_polarised_array.py index 5579152..801d295 100644 --- a/docs/examples/02_coherently_polarised_array.py +++ b/docs/examples/02_coherently_polarised_array.py @@ -110,8 +110,8 @@ UAV_Static_Pattern = antenna_pattern( azimuth_resolution=az_res, elevation_resolution=elev_res ) -UAV_Static_Pattern.pattern[:, :, 0] = Etheta -UAV_Static_Pattern.pattern[:, :, 1] = Ephi +UAV_Static_Pattern.pattern[:, :, 0] = Etheta.reshape(elev_res,az_res) +UAV_Static_Pattern.pattern[:, :, 1] = Ephi.reshape(elev_res,az_res) UAV_Static_Pattern.display_pattern(desired_pattern='Power') diff --git a/docs/examples/03_frequency_domain_channel_modelling.py b/docs/examples/03_frequency_domain_channel_modelling.py index 6ec5a5c..df4d44a 100644 --- a/docs/examples/03_frequency_domain_channel_modelling.py +++ b/docs/examples/03_frequency_domain_channel_modelling.py @@ -21,7 +21,7 @@ # Frequency and Mesh Resolution # ------------------------------ # -freq = np.asarray(22.0e9) +freq = np.asarray(24.0e9) wavelength = 3e8 / freq mesh_resolution = 0.5 * wavelength @@ -56,14 +56,14 @@ ) 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])) +transmit_horn_structure = GF.translate_mesh(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.translate_mesh(transmitting_antenna_surface_coords,np.asarray([2.695, 0, 0])) +transmitting_antenna_surface_coords = GF.translate_mesh(transmitting_antenna_surface_coords,np.asarray([2.529, 0, 0])) # %% # Position Receiver # ------------------ @@ -71,10 +71,10 @@ 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.translate_mesh(receive_horn_structure,np.asarray([0, 1.427, 0])) +receive_horn_structure = GF.translate_mesh(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.translate_mesh(receiving_antenna_surface_coords,np.asarray([0, 1.427, 0])) +receiving_antenna_surface_coords = GF.translate_mesh(receiving_antenna_surface_coords,np.asarray([0, 1.609, 0])) @@ -123,20 +123,17 @@ # %% # Visualise the Scene Geometry # ------------------------------ -#############################################NEED TO FIX THIS with pyvista + import pyvista as pv -def structure_cells(array): - ## add collumn of 3s to beggining of each row - array = np.append(np.ones((array.shape[0], 1), dtype=np.int32) * 3, array, axis=1) - return array -pyvista_mesh = pv.PolyData(reflectorplate.points, structure_cells(reflectorplate.cells[0].data)) -pyvista_mesh2 = pv.PolyData(receive_horn_structure.points, structure_cells(receive_horn_structure.cells[0].data)) -pyvista_mesh3 = pv.PolyData(transmit_horn_structure.points, structure_cells(transmit_horn_structure.cells[0].data)) +from lyceanem.utility.mesh_functions import pyvista_to_meshio ## plot the mesh plotter = pv.Plotter() -plotter.add_mesh(pyvista_mesh, color="white", show_edges=True) -plotter.add_mesh(pyvista_mesh2, color="blue", show_edges=True) -plotter.add_mesh(pyvista_mesh3, color="red", show_edges=True) +plotter.add_mesh(pv.from_meshio(reflectorplate), color="grey") +plotter.add_mesh(pv.from_meshio(scatter_points), color="grey") +plotter.add_mesh(pv.from_meshio(receive_horn_structure), color="blue") +plotter.add_mesh(pv.from_meshio(receiving_antenna_surface_coords), color="blue") +plotter.add_mesh(pv.from_meshio(transmit_horn_structure), color="red") +plotter.add_mesh(pv.from_meshio(transmitting_antenna_surface_coords), color="red") plotter.add_axes_at_origin() plotter.show() # Specify desired Transmit Polarisation @@ -196,37 +193,70 @@ def structure_cells(array): import copy from tqdm import tqdm - +from lyceanem.electromagnetics.emfunctions import update_electric_fields, PoyntingVector for angle_inc in tqdm(range(len(angle_values))): rotation_vector = np.radians(np.asarray([0.0, 0.0, angle_values[angle_inc]])) - scatter_points_temp = GF.mesh_rotate(copy.deepcopy(scatter_points),rotation_vector) - reflectorplate_temp = GF.mesh_rotate(copy.deepcopy(reflectorplate),rotation_vector) + 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]) - # pyvista_mesh = pv.PolyData(reflectorplate_temp.points, structure_cells(reflectorplate_temp.cells[0].data)) - # pyvista_mesh2 = pv.PolyData(receive_horn_structure.points, structure_cells(receive_horn_structure.cells[0].data)) - # pyvista_mesh3 = pv.PolyData(transmit_horn_structure.points, structure_cells(transmit_horn_structure.cells[0].data)) - # pyvista_mesh4 = pv.PolyData(scatter_points_temp.points) - # ## plot the mesh + + Ex, Ey, Ez = FD.calculate_scattering( + aperture_coords=transmitting_antenna_surface_coords, + sink_coords=scatter_points_temp, + antenna_solid=blockers, + desired_E_axis=np.tile(desired_E_axis,[transmitting_antenna_surface_coords.points.shape[0],1]), + scatter_points=scatter_points_temp, + wavelength=wavelength, + scattering=0, + project_vectors=False, + beta=(2*np.pi)/wavelength + ) + scattered_field=np.array([Ex*scatter_points_temp.point_data['Area'], + Ey*scatter_points_temp.point_data['Area'], + Ez*scatter_points_temp.point_data['Area']]).transpose() + #scatter_points_temp=update_electric_fields(scatter_points_temp, + # Ex*scatter_points_temp.point_data['Area'], + # Ey*scatter_points_temp.point_data['Area'], + # Ez*scatter_points_temp.point_data['Area']) + #scatter_points_temp=PoyntingVector(scatter_points_temp) + #scatter_points_temp.point_data["Ex"]=np.abs(scatter_points_temp.point_data['Ex-Real']+1j*scatter_points_temp.point_data['Ex-Imag']) + #scatter_points_temp.point_data["Ey"]=np.abs(scatter_points_temp.point_data['Ey-Real']+1j*scatter_points_temp.point_data['Ey-Imag']) + #scatter_points_temp.point_data["Ez"]=np.abs(scatter_points_temp.point_data['Ez-Real']+1j*scatter_points_temp.point_data['Ez-Imag']) # plotter = pv.Plotter() - # plotter.add_mesh(pyvista_mesh, color="white", show_edges=True) - # plotter.add_mesh(pyvista_mesh2, color="blue", show_edges=True) - # plotter.add_mesh(pyvista_mesh3, color="red", show_edges=True) - # plotter.add_mesh(pyvista_mesh4, color="green") + # plotter.add_mesh(pv.from_meshio(reflectorplate_temp), color="grey") + # plotter.add_mesh(pv.from_meshio(scatter_points_temp), scalars="Ey",clim=[0,0.0015]) + # plotter.add_mesh(pv.from_meshio(receive_horn_structure), color="blue") + # plotter.add_mesh(pv.from_meshio(receiving_antenna_surface_coords), color="blue") + # plotter.add_mesh(pv.from_meshio(transmit_horn_structure), color="red") + # plotter.add_mesh(pv.from_meshio(transmitting_antenna_surface_coords), color="red") # plotter.add_axes_at_origin() # plotter.show() - Ex, Ey, Ez = FD.calculate_scattering( + + 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 + ) + Ex3, Ey3, Ez3 = FD.calculate_scattering( aperture_coords=transmitting_antenna_surface_coords, sink_coords=receiving_antenna_surface_coords, antenna_solid=blockers, - desired_E_axis=desired_E_axis, + desired_E_axis=np.tile(desired_E_axis,[transmitting_antenna_surface_coords.points.shape[0],1]), scatter_points=scatter_points_temp, wavelength=wavelength, - scattering=1, - project_vectors=False + scattering=0, + project_vectors=False, + beta=(2*np.pi)/wavelength ) - responsex[angle_inc] = np.sum(Ex) - responsey[angle_inc] = np.sum(Ey) - responsez[angle_inc] = np.sum(Ez) + 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"]) @@ -258,10 +288,10 @@ def structure_cells(array): ax.plot(angle_values - 45, EzdB, label="Ez") plt.xlabel("$\\theta_{N}$ (degrees)") plt.ylabel("Normalised Level (dB)") -ax.set_ylim(-60.0, 0) +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(-60, 0.0, 21)) +ax.set_yticks(np.linspace(-40, 0.0, 21)) legend = ax.legend(loc="upper right", shadow=True) plt.grid() plt.show() diff --git a/docs/examples/05_array_beamforming.py b/docs/examples/05_array_beamforming.py index ed22860..d10a5ec 100644 --- a/docs/examples/05_array_beamforming.py +++ b/docs/examples/05_array_beamforming.py @@ -78,8 +78,9 @@ az_range = np.linspace(-180, 180, az_res) el_range = np.linspace(-90, 90, elev_res) +num_elements=Etheta.shape[0] directivity_map = MaximumDirectivityMap( - Etheta, Ephi, source_coords, wavelength, az_range, el_range + Etheta.reshape(num_elements,elev_res,az_res), Ephi.reshape(num_elements,elev_res,az_res), source_coords, wavelength, az_range, el_range ) from lyceanem.electromagnetics.beamforming import PatternPlot diff --git a/docs/examples/06_farfield_polarisation.py b/docs/examples/06_farfield_polarisation.py index 4ca59b4..49c2c25 100644 --- a/docs/examples/06_farfield_polarisation.py +++ b/docs/examples/06_farfield_polarisation.py @@ -74,8 +74,8 @@ u_pattern = antenna_pattern( azimuth_resolution=az_res, elevation_resolution=elev_res ) -u_pattern.pattern[:, :, 0] = Etheta -u_pattern.pattern[:, :, 1] = Ephi +u_pattern.pattern[:, :, 0] = Etheta.reshape(elev_res,az_res) +u_pattern.pattern[:, :, 1] = Ephi.reshape(elev_res,az_res) u_pattern.display_pattern(desired_pattern='Power') # %% @@ -99,8 +99,8 @@ v_pattern = antenna_pattern( azimuth_resolution=az_res, elevation_resolution=elev_res ) -v_pattern.pattern[:, :, 0] = Etheta -v_pattern.pattern[:, :, 1] = Ephi +v_pattern.pattern[:, :, 0] = Etheta.reshape(elev_res,az_res) +v_pattern.pattern[:, :, 1] = Ephi.reshape(elev_res,az_res) v_pattern.display_pattern(desired_pattern='Power') # %% @@ -124,64 +124,66 @@ n_pattern = antenna_pattern( azimuth_resolution=az_res, elevation_resolution=elev_res ) -n_pattern.pattern[:, :, 0] = Etheta -n_pattern.pattern[:, :, 1] = Ephi +n_pattern.pattern[:, :, 0] = Etheta.reshape(elev_res,az_res) +n_pattern.pattern[:, :, 1] = Ephi.reshape(elev_res,az_res) n_pattern.display_pattern(desired_pattern='Power') # %% # The point source can then be rotated, by providing a rotation matrix, and the u,v,n directions are moved with it in a consistent way. +from scipy.spatial.transform import Rotation as R -# point_antenna.rotate_antenna(o3d.geometry.get_rotation_matrix_from_axis_angle(np.radians(np.asarray([90.0,0.0,0.0])))) - -# 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(), -# 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), -# wavelength=wavelength, -# farfield_distance=20, -# elements=False, -# project_vectors=False, -# ) -# u_pattern.pattern[:, :, 0] = Etheta -# u_pattern.pattern[:, :, 1] = Ephi -# u_pattern.display_pattern(desired_pattern='Power') - - -# 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(), -# 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), -# wavelength=wavelength, -# farfield_distance=20, -# elements=False, -# project_vectors=False, -# ) -# v_pattern.pattern[:, :, 0] = Etheta -# v_pattern.pattern[:, :, 1] = Ephi -# v_pattern.display_pattern(desired_pattern='Power') - - -# 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(), -# 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), -# wavelength=wavelength, -# farfield_distance=20, -# elements=False, -# project_vectors=False, -# ) -# n_pattern.pattern[:, :, 0] = Etheta -# n_pattern.pattern[:, :, 1] = Ephi -# n_pattern.display_pattern(desired_pattern='Power') \ No newline at end of file +r=R.from_euler('xyz', np.radians(np.asarray([90.0,0.0,0.0]))) +point_antenna.rotate_antenna(r.as_matrix()) + +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(), + 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), + wavelength=wavelength, + farfield_distance=20, + elements=False, + project_vectors=False, +) +u_pattern.pattern[:, :, 0] = Etheta.reshape(elev_res,az_res) +u_pattern.pattern[:, :, 1] = Ephi.reshape(elev_res,az_res) +u_pattern.display_pattern(desired_pattern='Power') + + +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(), + 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), + wavelength=wavelength, + farfield_distance=20, + elements=False, + project_vectors=False, +) +v_pattern.pattern[:, :, 0] = Etheta.reshape(elev_res,az_res) +v_pattern.pattern[:, :, 1] = Ephi.reshape(elev_res,az_res) +v_pattern.display_pattern(desired_pattern='Power') + + +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(), + 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), + wavelength=wavelength, + farfield_distance=20, + elements=False, + project_vectors=False, +) +n_pattern.pattern[:, :, 0] = Etheta.reshape(elev_res,az_res) +n_pattern.pattern[:, :, 1] = Ephi.reshape(elev_res,az_res) +n_pattern.display_pattern(desired_pattern='Power') \ No newline at end of file diff --git a/docs/examples/07_aperture_farfield_polarisation.py b/docs/examples/07_aperture_farfield_polarisation.py index 0340ddb..ec49544 100644 --- a/docs/examples/07_aperture_farfield_polarisation.py +++ b/docs/examples/07_aperture_farfield_polarisation.py @@ -71,8 +71,8 @@ u_pattern = antenna_pattern( azimuth_resolution=az_res, elevation_resolution=elev_res ) -u_pattern.pattern[:, :, 0] = Etheta -u_pattern.pattern[:, :, 1] = Ephi +u_pattern.pattern[:, :, 0] = Etheta.reshape(elev_res,az_res) +u_pattern.pattern[:, :, 1] = Ephi.reshape(elev_res,az_res) u_pattern.display_pattern(desired_pattern='Power') # %% @@ -97,8 +97,8 @@ v_pattern = antenna_pattern( azimuth_resolution=az_res, elevation_resolution=elev_res ) -v_pattern.pattern[:, :, 0] = Etheta -v_pattern.pattern[:, :, 1] = Ephi +v_pattern.pattern[:, :, 0] = Etheta.reshape(elev_res,az_res) +v_pattern.pattern[:, :, 1] = Ephi.reshape(elev_res,az_res) v_pattern.display_pattern(desired_pattern='Power') # %% @@ -122,8 +122,8 @@ n_pattern = antenna_pattern( azimuth_resolution=az_res, elevation_resolution=elev_res ) -n_pattern.pattern[:, :, 0] = Etheta -n_pattern.pattern[:, :, 1] = Ephi +n_pattern.pattern[:, :, 0] = Etheta.reshape(elev_res,az_res) +n_pattern.pattern[:, :, 1] = Ephi.reshape(elev_res,az_res) n_pattern.display_pattern(desired_pattern='Power') # %% @@ -148,8 +148,8 @@ project_vectors=False, beta=(2*np.pi)/wavelength ) -u_pattern.pattern[:, :, 0] = Etheta -u_pattern.pattern[:, :, 1] = Ephi +u_pattern.pattern[:, :, 0] = Etheta.reshape(elev_res,az_res) +u_pattern.pattern[:, :, 1] = Ephi.reshape(elev_res,az_res) u_pattern.display_pattern(desired_pattern='Power') @@ -167,8 +167,8 @@ project_vectors=False, beta=(2*np.pi)/wavelength ) -v_pattern.pattern[:, :, 0] = Etheta -v_pattern.pattern[:, :, 1] = Ephi +v_pattern.pattern[:, :, 0] = Etheta.reshape(elev_res,az_res) +v_pattern.pattern[:, :, 1] = Ephi.reshape(elev_res,az_res) v_pattern.display_pattern(desired_pattern='Power') @@ -186,6 +186,6 @@ project_vectors=False, beta=(2*np.pi)/wavelength ) -n_pattern.pattern[:, :, 0] = Etheta -n_pattern.pattern[:, :, 1] = Ephi +n_pattern.pattern[:, :, 0] = Etheta.reshape(elev_res,az_res) +n_pattern.pattern[:, :, 1] = Ephi.reshape(elev_res,az_res) n_pattern.display_pattern(desired_pattern='Power') \ No newline at end of file diff --git a/lyceanem/base_classes.py b/lyceanem/base_classes.py index ded33f5..8c57abc 100644 --- a/lyceanem/base_classes.py +++ b/lyceanem/base_classes.py @@ -389,8 +389,8 @@ def triangles_base_raycaster(self): """ triangles = np.empty((0), dtype=base_types.triangle_t) for item in range(len(self.solids)): - temp_object = copy.deepcopy(self.solids[item]) - temp_object = GF.mesh_transform(temp_object, self.pose, False) + #temp_object = copy.deepcopy(self.solids[item]) + temp_object = GF.mesh_transform(self.solids[item], self.pose, False) triangles = np.append(triangles, RF.convertTriangles(temp_object)) diff --git a/lyceanem/electromagnetics/emfunctions.py b/lyceanem/electromagnetics/emfunctions.py index 1f1e378..50efa70 100644 --- a/lyceanem/electromagnetics/emfunctions.py +++ b/lyceanem/electromagnetics/emfunctions.py @@ -229,9 +229,8 @@ def transform_em(field_data, r): + 1j * field_data.point_data["Ey-Imag"], field_data.point_data["Ez-Real"] + 1j * field_data.point_data["Ez-Imag"], - np.zeros((field_data.point_data["Ex-Real"].shape[0])), ] - ).transpose() + ).reshape(3,-1).transpose() rot_fields = r.apply(fields) field_data.point_data["Ex-Real"] = np.real( rot_fields[:, 0] @@ -258,9 +257,8 @@ def transform_em(field_data, r): + 1j * field_data.point_data["Ey-Imag"], field_data.point_data["Ez-Real"] + 1j * field_data.point_data["Ez-Imag"], - np.zeros((field_data.point_data["Ex-Real"].shape[0])), ] - ).transpose() + ).reshape(3,-1).transpose() rot_fields = r.apply(fields) field_data.point_data["Ex-Real"] = np.real(rot_fields[:, 0]) field_data.point_data["Ey-Real"] = np.real(rot_fields[:, 1]) diff --git a/lyceanem/geometry/targets.py b/lyceanem/geometry/targets.py index b9142e8..14251c2 100755 --- a/lyceanem/geometry/targets.py +++ b/lyceanem/geometry/targets.py @@ -116,13 +116,15 @@ def rectReflector(majorsize, minorsize, thickness): ) pv_mesh = pv_mesh.triangulate() pv_mesh.compute_normals(inplace=True, consistent_normals=False) + from ..utility.mesh_functions import pyvista_to_meshio triangles = np.reshape(np.array(pv_mesh.faces), (12, 4)) triangles = triangles[:, 1:] - mesh = meshio.Mesh(pv_mesh.points, {"triangle": triangles}) + #mesh = meshio.Mesh(pv_mesh.points, {"triangle": triangles}) + mesh=pyvista_to_meshio(pv_mesh) - mesh.point_data["Normals"] = pv_mesh.point_normals - mesh.cell_data["Normals"] = pv_mesh.cell_normals + from .geometryfunctions import compute_normals + mesh=compute_normals(mesh) red = np.zeros((8, 1), dtype=np.float32) green = np.ones((8, 1), dtype=np.float32) * 0.259 blue = np.ones((8, 1), dtype=np.float32) * 0.145 @@ -198,7 +200,8 @@ def shapeTrapezoid(x_size, y_size, length, flare_angle): mesh.point_data["Normals"] = np.asarray(pv_mesh.point_normals) mesh.cell_data["Normals"] = [np.asarray(pv_mesh.cell_normals)] - + from .geometryfunctions import compute_areas + mesh=compute_areas(mesh) red = np.zeros((8, 1), dtype=np.float32) green = np.ones((8, 1), dtype=np.float32) * 0.259 blue = np.ones((8, 1), dtype=np.float32) * 0.145 @@ -206,6 +209,7 @@ def shapeTrapezoid(x_size, y_size, length, flare_angle): mesh.point_data["red"] = red mesh.point_data["green"] = green mesh.point_data["blue"] = blue + return mesh @@ -453,6 +457,39 @@ def parabola(x): return mesh, aperture_points +def spherical_field(az_range,elev_range,outward_normals=False,field_radius=1.0): + """ + Create a spherical field of points, with normals, a triangle mesh, and areas calculated for each triangle and adjusted for each field point. + + Parameters + ---------- + az_range : array_like[float] + Azimuthal angle in degrees ``[0, 360]``. + elev_range : array_like[float] + Elevation angle in degrees ``[-90, 90]``. + outward_normals : bool, optional + If outward pointing normals are required, set as True + field_radius : float, optional + The radius of the field, default is 1.0 m + + Returns + ------- + mesh : meshio object + spherical field of points at specified azimuth and elevation angles, with meshed triangles + """ + vista_pattern = pv.grid_from_sph_coords(az_range, (90-elev_range), field_radius).extract_surface() + if outward_normals: + vista_pattern.point_data['Normals']=vista_pattern.points/(np.linalg.norm(vista_pattern.points,axis=1).reshape(-1,1)) + else: + vista_pattern.point_data['Normals'] = (vista_pattern.points / ( + np.linalg.norm(vista_pattern.points, axis=1).reshape(-1, 1)))*-1.0 + + from ..utility.mesh_functions import pyvista_to_meshio + mesh=pyvista_to_meshio(vista_pattern.triangulate()) + from ..geometry.geometryfunctions import compute_areas,theta_phi_r + mesh=theta_phi_r(compute_areas(mesh)) + + return mesh def linear_parabolic_aperture(diameter, focal_length, height, thickness, mesh_size, lip=False, lip_width=1e-3): # Define function for parabola equation (y^2 = 4*focal_length*x) @@ -676,5 +713,6 @@ def gridedReflectorPoints( ], point_data={"Normals": mesh_normals}, ) + mesh_points.point_data['Area']=np.ones((mesh_points.points.shape[0]))*(grid_resolution**2) return mesh_points diff --git a/lyceanem/models/frequency_domain.py b/lyceanem/models/frequency_domain.py index 1005cdd..5ec9fe4 100644 --- a/lyceanem/models/frequency_domain.py +++ b/lyceanem/models/frequency_domain.py @@ -142,315 +142,31 @@ def calculate_farfield( """ # create sink points for the model - azaz, elel = np.meshgrid(az_range, el_range) - theta = GF.elevationtotheta(elel) - sinks = np.zeros((len(np.ravel(azaz)), 3), dtype=np.float32) - sinks[:, 0], sinks[:, 1], sinks[:, 2] = RF.azeltocart( - np.ravel(azaz), np.ravel(elel), farfield_distance + from ..geometry.targets import spherical_field + sink_coords = spherical_field(az_range, el_range, farfield_distance) + + Ex, Ey, Ez=calculate_scattering(aperture_coords, + sink_coords, + antenna_solid, + desired_E_axis, + scatter_points=scatter_points, + wavelength=wavelength, + scattering=scattering, + elements=elements, + los=los, + mesh_resolution=mesh_resolution, + project_vectors=project_vectors, + antenna_axes=antenna_axes, + multiE=False, + alpha=alpha, + beta=beta) + # convert to etheta,ephi + etheta = ( + Ex * np.cos(sink_coords.point_data["phi_(Radians)"]) * np.cos(sink_coords.point_data["theta_(Radians)"]) + + Ey * np.sin(sink_coords.point_data["phi_(Radians)"]) * np.cos(sink_coords.point_data["theta_(Radians)"]) + - Ez * np.sin(sink_coords.point_data["theta_(Radians)"]) ) - sink_normals = np.zeros((len(np.ravel(azaz)), 3), dtype=np.float32) - origin = np.zeros((len(sinks), 3), dtype=np.float32).ravel() - lengths = np.zeros((len(np.ravel(azaz)), 1), dtype=np.float32) - sink_normals, _ = calc_dv_norm( - sinks, np.zeros((len(sinks), 3), dtype=np.float32), sink_normals, lengths - ) - sink_cloud = RF.points2pointcloud(sinks) - sink_cloud.point_data["Normals"] = sink_normals - num_sources = len(np.asarray(aperture_coords.points)) - num_sinks = len(np.asarray(sink_cloud.points)) - environment_triangles = GF.mesh_conversion(antenna_solid) - - if project_vectors: - conformal_E_vectors = EM.calculate_conformalVectors( - desired_E_axis, - np.asarray(aperture_coords.point_data["Normals"]), - antenna_axes, - ) - else: - if ( - desired_E_axis.shape[0] - == np.asarray(aperture_coords.point_data["Normals"]).shape[0] - ): - conformal_E_vectors = copy.deepcopy(desired_E_axis) - else: - conformal_E_vectors = np.repeat( - desired_E_axis.reshape(1, 3).astype(np.complex64), num_sources, axis=0 - ) - - if scattering == 0: - # only use the aperture point cloud, no scattering required. - scatter_points = meshio.Mesh(points=np.empty((0, 3)), cells=[]) - unified_model = np.append( - np.asarray(aperture_coords.points).astype(np.float32), - np.asarray(sink_cloud.points).astype(np.float32), - axis=0, - ) - unified_normals = np.append( - np.asarray(aperture_coords.point_data["Normals"]).astype(np.float32), - np.asarray(sink_cloud.point_data["Normals"]).astype(np.float32), - axis=0, - ) - unified_weights = np.ones((unified_model.shape[0], 3), dtype=np.complex64) - if source_weights is None: - unified_weights[0:num_sources, :] = ( - conformal_E_vectors # / num_sources # set total amplitude to 1 for the aperture - ) - else: - unified_weights[0:num_sources, :] = source_weights - - unified_weights[num_sources : num_sources + num_sinks, :] = ( - 1 # / num_sinks # set total amplitude to 1 for the aperture - ) - point_informationv2 = np.empty((len(unified_model)), dtype=scattering_t) - # set all sources as magnetic current sources, and permittivity and permeability as free space - point_informationv2[:]["Electric"] = True - point_informationv2[:]["permittivity"] = 8.8541878176e-12 - point_informationv2[:]["permeability"] = 1.25663706212e-6 - # set position, velocity, normal, and weight of sources - point_informationv2[0:num_sources]["px"] = np.asarray( - aperture_coords.points - ).astype(np.float32)[:, 0] - point_informationv2[0:num_sources]["py"] = np.asarray( - aperture_coords.points - ).astype(np.float32)[:, 1] - point_informationv2[0:num_sources]["pz"] = np.asarray( - aperture_coords.points - ).astype(np.float32)[:, 2] - normals = np.asarray(aperture_coords.point_data["Normals"]) - point_informationv2[0:num_sources]["nx"] = normals[:, 0] - point_informationv2[0:num_sources]["ny"] = normals[:, 1] - point_informationv2[0:num_sources]["nz"] = normals[:, 2] - # set position and velocity of sinks - point_informationv2[num_sources : (num_sources + num_sinks)]["px"] = sinks[:, 0] - point_informationv2[num_sources : (num_sources + num_sinks)]["py"] = sinks[:, 1] - point_informationv2[num_sources : (num_sources + num_sinks)]["pz"] = sinks[:, 2] - point_informationv2[num_sources : (num_sources + num_sinks)]["nx"] = ( - sink_normals[:, 0] - ) - point_informationv2[num_sources : (num_sources + num_sinks)]["ny"] = ( - sink_normals[:, 1] - ) - point_informationv2[num_sources : (num_sources + num_sinks)]["nz"] = ( - sink_normals[:, 2] - ) - - point_informationv2[:]["ex"] = unified_weights[:, 0] - point_informationv2[:]["ey"] = unified_weights[:, 1] - point_informationv2[:]["ez"] = unified_weights[:, 2] - scatter_mask = np.zeros((point_informationv2.shape[0]), dtype=np.int32) - scatter_mask[0:num_sources] = 0 - scatter_mask[(num_sources + num_sinks) :] = 0 - - else: - # create scatter points on antenna solids based upon a half wavelength square - if scatter_points is None: - scatter_points, areas = TL.source_cloud_from_shape( - antenna_solid, 1e-6, (wavelength * mesh_resolution) ** 2 / 2.0 - ) - - unified_model = np.append( - np.append( - np.asarray(aperture_coords.points).astype(np.float32), - np.asarray(sink_cloud.points).astype(np.float32), - axis=0, - ), - np.asarray(scatter_points.points).astype(np.float32), - axis=0, - ) - unified_normals = np.append( - np.append( - np.asarray(aperture_coords.point_data["Normals"]).astype(np.float32), - np.asarray(sink_cloud.point_data["Normals"]).astype(np.float32), - axis=0, - ), - np.asarray(scatter_points.point_data["Normals"]).astype(np.float32), - axis=0, - ) - unified_weights = np.ones((unified_model.shape[0], 3), dtype=np.complex64) - if source_weights is None: - unified_weights[0:num_sources, :] = ( - conformal_E_vectors # / num_sources # set total amplitude to 1 for the aperture - ) - else: - unified_weights[0:num_sources, :] = source_weights - unified_weights[num_sources : num_sources + num_sinks, :] = ( - 1 # / num_sinks # set total amplitude to 1 for the aperture - ) - unified_weights[num_sources + num_sinks :, :] = ( - scattering_weight # / len(np.asarray(scatter_points.points)) # set total amplitude to 1 for the aperture - ) - point_informationv2 = np.empty((len(unified_model)), dtype=scattering_t) - # set all sources as magnetic current sources, and permittivity and permeability as free space - point_informationv2[:]["Electric"] = True - point_informationv2[:]["permittivity"] = 8.8541878176e-12 - point_informationv2[:]["permeability"] = 1.25663706212e-6 - # set position, velocity, normal, and weight of sources - point_informationv2[0:num_sources]["px"] = np.asarray( - aperture_coords.points - ).astype(np.float32)[:, 0] - point_informationv2[0:num_sources]["py"] = np.asarray( - aperture_coords.points - ).astype(np.float32)[:, 1] - point_informationv2[0:num_sources]["pz"] = np.asarray( - aperture_coords.points - ).astype(np.float32)[:, 2] - point_informationv2[0:num_sources]["nx"] = np.asarray( - aperture_coords.point_data["Normals"] - ).astype(np.float32)[:, 0] - point_informationv2[0:num_sources]["ny"] = np.asarray( - aperture_coords.point_data["Normals"] - ).astype(np.float32)[:, 1] - point_informationv2[0:num_sources]["nz"] = np.asarray( - aperture_coords.point_data["Normals"] - ).astype(np.float32)[:, 2] - # point_informationv2[0:num_sources]['ex']=unified_weights[0:num_sources,0] - # point_informationv2[0:num_sources]['ey']=unified_weights[0:num_sources,1] - # point_informationv2[0:num_sources]['ez']=unified_weights[0:num_sources,2] - # set position and velocity of sinks - point_informationv2[num_sources : (num_sources + num_sinks)]["px"] = sinks[:, 0] - point_informationv2[num_sources : (num_sources + num_sinks)]["py"] = sinks[:, 1] - point_informationv2[num_sources : (num_sources + num_sinks)]["pz"] = sinks[:, 2] - # point_informationv2[num_sources:(num_sources+num_sinks)]['vx']=0.0 - # point_informationv2[num_sources:(num_sources+num_sinks)]['vy']=0.0 - # point_informationv2[num_sources:(num_sources+num_sinks)]['vz']=0.0 - point_informationv2[num_sources : (num_sources + num_sinks)]["nx"] = ( - sink_normals[:, 0] - ) - point_informationv2[num_sources : (num_sources + num_sinks)]["ny"] = ( - sink_normals[:, 1] - ) - point_informationv2[num_sources : (num_sources + num_sinks)]["nz"] = ( - sink_normals[:, 2] - ) - point_informationv2[(num_sources + num_sinks) :]["px"] = np.asarray( - scatter_points.points - ).astype(np.float32)[:, 0] - point_informationv2[(num_sources + num_sinks) :]["py"] = np.asarray( - scatter_points.points - ).astype(np.float32)[:, 1] - point_informationv2[(num_sources + num_sinks) :]["pz"] = np.asarray( - scatter_points.points - ).astype(np.float32)[:, 2] - scatter_normals = np.asarray(scatter_points.point_data["Normals"]) - point_informationv2[(num_sources + num_sinks) :]["nx"] = scatter_normals[:, 0] - point_informationv2[(num_sources + num_sinks) :]["ny"] = scatter_normals[:, 1] - point_informationv2[(num_sources + num_sinks) :]["nz"] = scatter_normals[:, 2] - point_informationv2[:]["ex"] = unified_weights[:, 0] - point_informationv2[:]["ey"] = unified_weights[:, 1] - point_informationv2[:]["ez"] = unified_weights[:, 2] - scatter_mask = np.zeros((point_informationv2.shape[0]), dtype=np.int32) - scatter_mask[0:num_sources] = 0 - scatter_mask[(num_sources + num_sinks) :] = scattering - - # full_index, initial_index = RF.integratedraycastersetup(num_sources, - # num_sinks, - # point_informationv2, - # RF.convertTriangles(antenna_solid), - # scatter_mask) - full_index, rays = RF.workchunkingv2( - np.asarray(aperture_coords.points).astype(np.float32), - sinks, - np.asarray(scatter_points.points).astype(np.float32), - environment_triangles, - scattering + 1, - line_of_sight=los, - ) - - if elements: - # create efiles for model - etheta = np.zeros( - (num_sources, el_range.shape[0], az_range.shape[0]), dtype=np.complex64 - ) - ephi = np.zeros( - (num_sources, el_range.shape[0], az_range.shape[0]), dtype=np.complex64 - ) - Ex = np.zeros( - (num_sources, el_range.shape[0], az_range.shape[0]), dtype=np.complex64 - ) - Ey = np.zeros( - (num_sources, el_range.shape[0], az_range.shape[0]), dtype=np.complex64 - ) - Ez = np.zeros( - (num_sources, el_range.shape[0], az_range.shape[0]), dtype=np.complex64 - ) - - for element in range(num_sources): - point_informationv2[0:num_sources]["ex"] = 0.0 - point_informationv2[0:num_sources]["ey"] = 0.0 - point_informationv2[0:num_sources]["ez"] = 0.0 - point_informationv2[element]["ex"] = conformal_E_vectors[ - element, 0 - ] # / num_sources - point_informationv2[element]["ey"] = conformal_E_vectors[ - element, 1 - ] # / num_sources - point_informationv2[element]["ez"] = conformal_E_vectors[ - element, 2 - ] # / num_sources - # unified_weights[0:num_sources, :] = 0.0 - # unified_weights[element, :] = (conformal_E_vectors[element, :] / num_sources)*v_transmit - scatter_map = EM.EMGPUFreqDomain( - num_sources, - sinks.shape[0], - full_index, - point_informationv2, - wavelength, - alpha, - beta, - ) - Ex[element, :, :] = np.dot( - np.ones((num_sources)), scatter_map[:, :, 0] - ).reshape(el_range.shape[0], az_range.shape[0]) - Ey[element, :, :] = np.dot( - np.ones((num_sources)), scatter_map[:, :, 1] - ).reshape(el_range.shape[0], az_range.shape[0]) - Ez[element, :, :] = np.dot( - np.ones((num_sources)), scatter_map[:, :, 2] - ).reshape(el_range.shape[0], az_range.shape[0]) - etheta[element, :, :] = ( - Ex[element, :, :] * np.cos(np.deg2rad(azaz)) * np.cos(np.deg2rad(theta)) - + Ey[element, :, :] - * np.sin(np.deg2rad(azaz)) - * np.cos(np.deg2rad(theta)) - - Ez[element, :, :] * np.sin(np.deg2rad(theta)) - ) - ephi[element, :, :] = -Ex[element, :, :] * np.sin(np.deg2rad(azaz)) + Ey[ - element, :, : - ] * np.cos(np.deg2rad(azaz)) - - else: - # create efiles for model - etheta = np.zeros((el_range.shape[0], az_range.shape[0]), dtype=np.complex64) - ephi = np.zeros((el_range.shape[0], az_range.shape[0]), dtype=np.complex64) - Ex = np.zeros((el_range.shape[0], az_range.shape[0]), dtype=np.complex64) - Ey = np.zeros((el_range.shape[0], az_range.shape[0]), dtype=np.complex64) - Ez = np.zeros((el_range.shape[0], az_range.shape[0]), dtype=np.complex64) - scatter_map = EM.EMGPUFreqDomain( - num_sources, - num_sinks, - full_index, - point_informationv2, - wavelength, - alpha, - beta, - ) - - Ex[:, :] = np.sum(scatter_map[:, :, 0], axis=0).reshape( - el_range.shape[0], az_range.shape[0] - ) - Ey[:, :] = np.sum(scatter_map[:, :, 1], axis=0).reshape( - el_range.shape[0], az_range.shape[0] - ) - Ez[:, :] = np.sum(scatter_map[:, :, 2], axis=0).reshape( - el_range.shape[0], az_range.shape[0] - ) - # convert to etheta,ephi - etheta = ( - Ex * np.cos(np.deg2rad(azaz)) * np.cos(np.deg2rad(theta)) - + Ey * np.sin(np.deg2rad(azaz)) * np.cos(np.deg2rad(theta)) - - Ez * np.sin(np.deg2rad(theta)) - ) - ephi = -Ex * np.sin(np.deg2rad(azaz)) + Ey * np.cos(np.deg2rad(azaz)) + ephi = -Ex * np.sin(sink_coords.point_data["phi_(Radians)"]) + Ey * np.cos(sink_coords.point_data["phi_(Radians)"]) return etheta, ephi