Merge remote-tracking branch 'refs/remotes/origin/feature/o3dremove' …

…into feature/o3dremove

docs/examples/03_frequency_domain_channel_modelling.py

@@ -21,7 +21,7 @@
 # Frequency and Mesh Resolution
 # ------------------------------
 #
-freq = np.asarray(24.0e9)
+freq = np.asarray(16.0e9)
 wavelength = 3e8 / freq
 mesh_resolution = 0.5 * wavelength
 
@@ -179,7 +179,7 @@ def structure_cells(array):
 
 
 
-angle_values = np.linspace(0, 90, 91)
+angle_values = np.linspace(0, 90, 181)
 angle_increment = np.diff(angle_values)[0]
 responsex = np.zeros((len(angle_values)), dtype="complex")
 responsey = np.zeros((len(angle_values)), dtype="complex")
@@ -223,9 +223,9 @@ def structure_cells(array):
         scattering=1,
         project_vectors=False
     )
-    responsex[angle_inc] = Ex
-    responsey[angle_inc] = Ey
-    responsez[angle_inc] = Ez
+    responsex[angle_inc] = np.sum(Ex)
+    responsey[angle_inc] = np.sum(Ey)
+    responsez[angle_inc] = np.sum(Ez)
 
 
 

docs/examples/04_time_domain_channel_modelling.py

@@ -30,7 +30,7 @@
 noise_power = -80  # dbw
 mean_noise = 0
 
-model_freq = 16e9
+model_freq = 24e9
 wavelength = 3e8 / model_freq
 
 # %%

docs/examples/07_aperture_farfield_polarisation.py

@@ -37,43 +37,8 @@
 from lyceanem.geometry.targets import meshedHorn
 structure,array_points=meshedHorn(3*wavelength, 1*wavelength, 4*wavelength, 1*wavelength,np.radians(10),wavelength*0.5)
 
-def planar_aperture(u_size,v_size,wall_thickness=2e-3,depth=2e-3,mesh_size=1.0):
-    with pygmsh.occ.Geometry() as geom:
-        geom.characteristic_length_max=mesh_size
-        u_points=np.linspace(-u_size/2,u_size/2,np.ceil(u_size/mesh_size).astype(int))
-        v_points=np.linspace(-v_size/2,v_size/2,np.ceil(v_size/mesh_size).astype(int))
-        n_points=np.linspace(0,0,1)
-        u_mesh,v_mesh,n_mesh=np.meshgrid(u_points,v_points,n_points,indexing='ij')
-        aperture_points=np.array([u_mesh.ravel(),v_mesh.ravel(),n_mesh.ravel()]).transpose()
-        source_normals = np.empty(((aperture_points.shape[0], 3)), dtype=np.float32)
-        source_normals[:] = 0.0
-        source_normals[:, 2] = 1.0
-        ground_plane=geom.add_box((-(u_size+wall_thickness*2)/2.0,-(v_size+wall_thickness*2)/2.0,-depth), extents=(u_size+wall_thickness*2,v_size+wall_thickness*2,depth),mesh_size=mesh_size)
-        point_list=[]
-        for inc in range(aperture_points.shape[0]):
-            point_list.append(geom.add_point(aperture_points[inc,:]))
-            geom.in_surface(point_list[-1],ground_plane)
-
-        #aperture=geom.add_rectangle((-u_size/2.0,-v_size/2.0,0.0), u_size, v_size,mesh_size=mesh_size)
-        #ground_plane=geom.add_box((-(u_size+wall_thickness*2)/2.0,-(v_size+wall_thickness*2)/2.0,-depth), extents=(u_size+wall_thickness*2,v_size+wall_thickness*2,depth),mesh_size=mesh_size)
-
-
-        mesh = geom.generate_mesh()
-        end_point=mesh.points.shape[0]
-        mesh.point_data['Normals']=np.zeros((mesh.points.shape[0],3))
-        mesh.point_data['Normals'][mesh.points[:,2]>-depth*0.5,2]=1.0
-        mesh.point_data['Normals'][mesh.points[:,2]>-depth*0.5,2]=-1.0
-        mesh.points=np.append(mesh.points,aperture_points,axis=0)
-        mesh.point_data['Normals']=np.append(mesh.point_data['Normals'],source_normals,axis=0)
-        mesh.cell_sets['Aperture']=[None,None,None,np.linspace(end_point,aperture_points.shape[0]+end_point-1,aperture_points.shape[0]).astype(int)]
-        mesh.cells[3].data=np.append(mesh.cells[3].data,np.linspace(end_point,aperture_points.shape[0]+end_point-1,aperture_points.shape[0]).astype(int).reshape(-1,1),axis=0)
-    return mesh
-
-mesh2=planar_aperture(3*wavelength,1*wavelength,mesh_size=wavelength*0.5)
-
-
-
-horn_antenna=antenna_structures(structures(solids=[mesh2]), points(points=[array_points]))
+
+horn_antenna=antenna_structures(structures(solids=[structure]), points(points=[array_points]))
 
 
 from lyceanem.models.frequency_domain import calculate_farfield

lyceanem/electromagnetics/empropagation.py

@@ -1266,10 +1266,27 @@ def freqdomainisokernal(
         # print(cu_ray_num,source_sink_index[cu_ray_num,0],source_sink_index[cu_ray_num,1])
         wave_vector = (2.0 * cmath.pi) / wavelength
         # scatter_coefficient=(1/(4*cmath.pi))**(complex(scatter_index))
-        loss = cmath.exp(-lengths * wave_vector * 1j) * (
-            wavelength / (4 * cmath.pi * lengths)
+        alpha = 0.0
+        beta = (2.0 * cmath.pi) / wavelength[0]
+        loss = lossy_propagation(
+            calc_sep(
+                point_information[network_index[cu_ray_num, 0] - 1],
+                point_information[network_index[cu_ray_num, 1] - 1],
+                float(0)
+            ),
+            alpha,
+            beta
         )
-        ray_component[0] = loss
+        for i in range(1,network_index.shape[1] - 1):
+            if network_index[cu_ray_num, i + 1] != 0:
+
+                loss *= lossy_propagation(calc_sep(point_information[network_index[cu_ray_num, i] - 1],
+                                                   point_information[network_index[cu_ray_num, i + 1] - 1],
+                                                   float(0)),
+                                          alpha,
+                                          beta
+                                          )
+        ray_component[0] *= loss
         # ray_component[1]=loss
         # ray_component[2]=loss
         # print(ray_component[0].real,ray_component[1].real,ray_component[2].real)
@@ -1419,9 +1436,26 @@ def timedomainkernal(
             i = i + 1
 
         # print(cu_ray_num,source_sink_index[cu_ray_num,0],source_sink_index[cu_ray_num,1])
-        wave_vector = (2.0 * cmath.pi) / wavelength[0]
-        # scatter_coefficient=(1/(4*cmath.pi))**(complex(scatter_index))
-        loss = wavelength[0] / (4 * cmath.pi * lengths)
+        alpha = 0.0
+        beta = (2.0 * cmath.pi) / wavelength[0]
+        loss = lossy_propagation(
+            calc_sep(
+                point_information[full_index[cu_ray_num, 0] - 1],
+                point_information[full_index[cu_ray_num, 1] - 1],
+                float(0)
+            ),
+            alpha,
+            beta
+        )
+        for i in range(1,full_index.shape[1] - 1):
+            if full_index[cu_ray_num, i + 1] != 0:
+
+                loss *= lossy_propagation(calc_sep(point_information[full_index[cu_ray_num, i] - 1],
+                                                   point_information[full_index[cu_ray_num, i + 1] - 1],
+                                                   float(0)),
+                                          alpha,
+                                          beta
+                                          )
 
         ray_component[0] *= loss
         ray_component[1] *= loss

lyceanem/geometry/geometryfunctions.py

@@ -92,6 +92,44 @@ def mesh_transform(mesh, transform_matrix, rotate_only):
 
     return return_mesh
 
+def compute_normals(mesh):
+    """
+    Computes the Cell Normals for meshio mesh objects, this does not currently calculate the point normals, but this will be done soon.
+
+    Parameters
+    ----------
+    mesh
+
+    Returns
+    -------
+
+    """
+    cell_normal_list = []
+    for inc, cell in enumerate(mesh.cells):
+        print(cell.type, cell.data.shape[0])
+        if cell.type == 'vertex':
+            vertex_normals = np.zeros((cell.data.shape[0], 3))
+            cell_normal_list.append(vertex_normals)
+        if cell.type == 'line':
+            line_normals = np.zeros((cell.data.shape[0], 3))
+            cell_normal_list.append(line_normals)
+        if cell.type == 'triangle':
+            # print(inc)
+            edge1 = mesh.points[cell.data[:, 0], :] - mesh.points[cell.data[:, 1], :]
+            edge2 = mesh.points[cell.data[:, 0], :] - mesh.points[cell.data[:, 2], :]
+            tri_cell_normals = np.cross(edge1, edge2)
+            tri_cell_normals *= (1 / np.linalg.norm(tri_cell_normals, axis=1)).reshape(-1, 1)
+            cell_normal_list.append(tri_cell_normals)
+        if cell.type == 'tetra':
+            # print(inc)
+            edge1 = mesh.points[cell.data[:, 0], :] - mesh.points[cell.data[:, 1], :]
+            edge2 = mesh.points[cell.data[:, 0], :] - mesh.points[cell.data[:, 2], :]
+            tetra_cell_normals = np.cross(edge1, edge2)
+            tetra_cell_normals *= (1 / np.linalg.norm(tetra_cell_normals, axis=1)).reshape(-1, 1)
+            cell_normal_list.append(tetra_cell_normals)
+
+    mesh.cell_data['Normals'] = cell_normal_list
+    return mesh
 
 def mesh_conversion(conversion_object):
     """

lyceanem/geometry/targets.py

@@ -5,8 +5,7 @@
 import numpy as np
 import meshio
 import pyvista as pv
-import scipy.stats
-import solid as sd
+import pygmsh
 from importlib_resources import files
 from scipy.spatial.transform import Rotation as R
 
@@ -305,6 +304,114 @@ def meshedHorn(
     return structure, mesh_points
 
 
+def parabolic_aperture(diameter, focal_length, thickness, mesh_size, sides='front'):
+    # Define function for parabola equation (y^2 = 4*focal_length*x)
+    def parabola(x):
+        return (1 / (4 * focal_length)) * x ** 2
+
+    with pygmsh.occ.Geometry() as geom:
+        geom.characteristic_length_max = mesh_size
+        # Define points
+        cp1 = geom.add_point([0, 0, 0])  # Center point
+        cp2 = geom.add_point([diameter * 0.5 * (1 / 6), 0, parabola(diameter * 0.5 * (1 / 6))])
+        cp3 = geom.add_point([diameter * 0.5 * (2 / 6), 0, parabola(diameter * 0.5 * (2 / 6))])
+        cp4 = geom.add_point([diameter * 0.5 * (3 / 6), 0, parabola(diameter * 0.5 * (3 / 6))])
+        cp5 = geom.add_point([diameter * 0.5 * (4 / 6), 0, parabola(diameter * 0.5 * (4 / 6))])
+        cp6 = geom.add_point([diameter * 0.5 * (5 / 6), 0, parabola(diameter * 0.5 * (5 / 6))])
+        cp7 = geom.add_point([diameter * 0.5 * (6 / 6), 0, parabola(diameter * 0.5 * (6 / 6))])
+
+        # Define top line based on points
+        line = geom.add_bspline([cp1, cp2, cp3, cp4, cp5, cp6, cp7])
+
+        _, surface, _ = geom.extrude(line, translation_axis=[0.0, 0.0, -thickness])
+
+        # Revolve line to create revolution surface
+        volume_list = []
+        _, b, _ = geom.revolve(surface, rotation_axis=[0.0, 0.0, 1.0], point_on_axis=[0.0, 0.0, 0.0],
+                               angle=0.25 * np.pi)
+        volume_list.append(b)
+        for inc in range(7):
+
+            geom.rotate(surface, point=[0.0, 0.0, 0.0], angle=(1 / 4) * np.pi, axis=[0.0, 0.0, 1.0])
+            _, b2, _ = geom.revolve(surface, rotation_axis=[0.0, 0.0, 1.0], point_on_axis=[0.0, 0.0, 0.0],
+                                    angle=0.25 * np.pi)
+            volume_list.append(b2)
+
+        full_reflector = geom.boolean_union(volume_list)
+
+        mesh_temp = geom.generate_mesh(dim=2)
+    for inc, cell in enumerate(mesh_temp.cells):
+        if cell.type == 'triangle':
+            triangle_index = inc
+
+    import meshio
+    triangle_cells = [("triangle", mesh_temp.cells[triangle_index].data)]
+    mesh = meshio.Mesh(mesh_temp.points, triangle_cells)
+
+    x_space = np.linspace(
+        mesh_size,
+        (diameter / 2),
+        int(np.max(np.asarray([2, np.ceil((diameter * 0.5) / (mesh_size))]))),
+    )
+    z_space = (1 / (4 * focal_length)) * x_space ** 2
+    c_space = np.ceil((2 * np.pi * x_space) / mesh_size).astype(int)
+    normal_gradiant_vector = np.array(
+        [
+            np.ones((len(x_space))),
+            np.zeros((len(x_space))),
+            -1 / (1 / (2 * focal_length) * x_space),
+        ]
+    )
+    source = np.array([x_space, np.zeros((len(x_space))), z_space]).transpose()
+    target = (normal_gradiant_vector * -x_space).transpose()
+    base_directions = np.zeros((x_space.shape[0], 3), dtype=np.float32)
+    norm_length = np.zeros((x_space.shape[0], 1), dtype=np.float32)
+    base_directions, norm_length = math_functions.calc_dv_norm(
+        source, target, base_directions, norm_length
+    )
+    source_coords = np.empty(((1, 3)), dtype=np.float32)
+    source_coords[0, :] = 0
+    source_normals = np.empty(((1, 3)), dtype=np.float32)
+    source_normals[0, :] = 0
+    source_normals[0, 2] = 1
+    for r_index in range(x_space.shape[0]):
+        source_coords = np.append(
+            source_coords,
+            np.array(
+                [
+                    x_space[r_index]
+                    * np.cos(np.linspace(0, (2 * np.pi), c_space[r_index])[0:-1]),
+                    x_space[r_index]
+                    * np.sin(np.linspace(0, (2 * np.pi), c_space[r_index])[0:-1]),
+                    z_space[r_index] * np.ones(c_space[r_index] - 1),
+                ]
+            ).transpose(),
+            axis=0,
+        )
+        source_normals = np.append(
+            source_normals,
+            np.array(
+                [
+                    base_directions[r_index, 0]
+                    * np.cos(np.linspace(0, (2 * np.pi), c_space[r_index])[0:-1]),
+                    base_directions[r_index, 0]
+                    * np.sin(np.linspace(0, (2 * np.pi), c_space[r_index])[0:-1]),
+                    base_directions[r_index, 2] * np.ones(c_space[r_index] - 1),
+                ]
+            ).transpose(),
+            axis=0,
+        )
+
+    mesh_vertices = source_coords + np.array([0, 0, 1e-6])
+    mesh_normals = source_normals
+    aperture_points = meshio.Mesh(points=mesh_vertices,
+                                  cells=[("vertex", np.array([[i, ] for i in range(len(mesh_vertices))]))],
+                                  point_data={'Normals': mesh_normals})
+
+    return mesh, aperture_points
+
+
+
 
 
 

lyceanem/tests/test_base_classes.py

@@ -33,7 +33,8 @@ def cube():
                 [1, 4, 5]]
     ## put into meshio mesh
     cube = meshio.Mesh(points=cube_points, cells=[("triangle", cube_cells)])
-
+    from ..geometry.geometryfunctions import compute_normals
+    cube=compute_normals(cube)
     return cube
 
 

lyceanem/tests/test_geometryfunctions.py

@@ -1,17 +1,19 @@
 import numpy as np
-
+import lyceanem
 import pytest
 from numpy.testing import assert_allclose
 from scipy.spatial.transform import Rotation as R
 import meshio
-
+from importlib_resources import files
 from  ..geometry import geometryfunctions as GF
+import lyceanem.tests.data
 
 
 
 @pytest.fixture
 def standard_cube():
-    cube = meshio.read('data/cube.ply')
+    cube_location= files(lyceanem.tests.data).joinpath("cube.ply")
+    cube = meshio.read(cube_location)
     return cube
 
 
@@ -82,21 +84,21 @@ def are_meshes_equal(mesh1, mesh2, rtol=1e-5):
 
 
 def test_tri_areas_2():
-    refrence = np.load('data/areas.npy')
-    input = meshio.read('data/receive_horn.ply')
+    refrence = np.load(files(lyceanem.tests.data).joinpath("areas.npy"))
+    input = meshio.read(files(lyceanem.tests.data).joinpath("receive_horn.ply"))
     result = GF.tri_areas(input)
     assert np.allclose(result, refrence)
 
 def test_tri_centroids():
-    refrence = np.load('data/centroid_numpy.npy')
-    refrenece_mesh = meshio.read('data/centroid_mesh.ply')
-    input = meshio.read('data/receive_horn.ply')
+    refrence = np.load(files(lyceanem.tests.data).joinpath('centroid_numpy.npy'))
+    refrenece_mesh = meshio.read(files(lyceanem.tests.data).joinpath('centroid_mesh.ply'))
+    input = meshio.read(files(lyceanem.tests.data).joinpath('receive_horn.ply'))
     result_numpy, result_mesh = GF.tri_centroids(input)
     assert np.allclose(result_numpy, refrence)
     assert are_meshes_equal(result_mesh, refrenece_mesh)
 def test_mesh_rotate():
-    refrence = meshio.read('data/rotated_recieve.ply')
-    input = meshio.read('data/receive_horn.ply')
+    refrence = meshio.read(files(lyceanem.tests.data).joinpath('rotated_recieve.ply'))
+    input = meshio.read(files(lyceanem.tests.data).joinpath('receive_horn.ply'))
     rotation_vector1 = np.radians(np.asarray([90.0, 0.0, 0.0]))
     centre = np.array([1.0, 1.0, 1.0])
 

lyceanem/tests/test_ray_functions.py

@@ -2,12 +2,13 @@
 import numpy as np
 import pytest
 import lyceanem.base_types as base_types
-
+from importlib_resources import files
+import lyceanem.tests.data
 import meshio
 
 def test_convertTriangles():
-    reference = np.load("data/triangle_type_array_ref.npy")
-    mesh = meshio.read("data/receive_horn.ply")
+    reference = np.load(files(lyceanem.tests.data).joinpath("triangle_type_array_ref.npy"))
+    mesh = meshio.read(files(lyceanem.tests.data).joinpath("receive_horn.ply"))
     triangles = RF.convertTriangles(mesh)
     for i in range(triangles.size):
         assert np.allclose(triangles[i]["v0x"], reference[i]["v0x"])