diff --git a/lyceanem/base_classes.py b/lyceanem/base_classes.py
index 92ef62d..7023a87 100644
--- a/lyceanem/base_classes.py
+++ b/lyceanem/base_classes.py
@@ -120,7 +120,9 @@ def create_points(self, points, normals):
             """
             mesh_vertices = points.reshape(-1, 3)
             mesh_normals = normals.reshape(-1, 3)
-            new_point_cloud = meshio.Mesh(points=mesh_vertices, cells=[], point_data={"Normals": mesh_normals})
+            new_point_cloud = meshio.Mesh(
+                points=mesh_vertices, cells=[], point_data={"Normals": mesh_normals}
+            )
 
             self.add_points(new_point_cloud)
 
@@ -143,7 +145,12 @@ def rotate_points(
         """
         # warning, current commond just rotates around the origin, and until Open3D can be brought up to the
         # latest version without breaking BlueCrystal reqruiements, this will require additional code.
-        assert (rotation_vector.shape == (3,) or rotation_vector.shape == (3,1) or rotation_vector.shape == (1,3) or rotation_vector.shape == (3,3)), "Rotation vector must be a 3x1 or 3x3 array"
+        assert (
+            rotation_vector.shape == (3,)
+            or rotation_vector.shape == (3, 1)
+            or rotation_vector.shape == (1, 3)
+            or rotation_vector.shape == (3, 3)
+        ), "Rotation vector must be a 3x1 or 3x3 array"
         for item in range(len(self.points)):
             self.points[item] = GF.mesh_rotate(
                 self.points[item], rotation_vector, rotation_centre
@@ -174,52 +181,93 @@ def export_points(self, point_index=None):
         combined points
         """
         if point_index == None:
-            points=np.empty((0,3))
+            points = np.empty((0, 3))
             for item in range(len(self.points)):
                 if item == 0:
-                    points = np.append(points,np.array(self.points[item].points),axis=0)
+                    points = np.append(
+                        points, np.array(self.points[item].points), axis=0
+                    )
                 else:
                     points = np.append(points, self.points[item].points, axis=0)
             point_data = {}
             for key in self.points[0].point_data.keys():
-                if len(self.points[0].point_data[key].shape)<2:
-                    point_data[key] = np.empty((0,1))
+                if len(self.points[0].point_data[key].shape) < 2:
+                    point_data[key] = np.empty((0, 1))
                 else:
-                    point_data[key] = np.empty((0,self.points[0].point_data[key].shape[1]))
-                    
-                for item in range(len(self.points)):  
+                    point_data[key] = np.empty(
+                        (0, self.points[0].point_data[key].shape[1])
+                    )
+
+                for item in range(len(self.points)):
                     point_data_element = np.array(self.points[item].point_data[key])
-                    if len(point_data_element.shape)<2:
-                        point_data[key] = np.append(point_data[key], point_data_element.reshape(-1,1), axis=0)    
+                    if len(point_data_element.shape) < 2:
+                        point_data[key] = np.append(
+                            point_data[key], point_data_element.reshape(-1, 1), axis=0
+                        )
                     else:
-                        point_data[key] = np.append(point_data[key], point_data_element, axis=0)
-
-            combinded_points = meshio.Mesh(points, cells = [("vertex", np.array([[i, ] for i in range(points.shape[0])]))], point_data=point_data)
+                        point_data[key] = np.append(
+                            point_data[key], point_data_element, axis=0
+                        )
+
+            combinded_points = meshio.Mesh(
+                points,
+                cells=[
+                    (
+                        "vertex",
+                        np.array(
+                            [
+                                [
+                                    i,
+                                ]
+                                for i in range(points.shape[0])
+                            ]
+                        ),
+                    )
+                ],
+                point_data=point_data,
+            )
             combinded_points = GF.mesh_transform(combinded_points, self.pose, False)
             return combinded_points
-                    
-
-
 
         else:
-            points=np.empty((0,3))
+            points = np.empty((0, 3))
             for item in point_index:
                 if item == 0:
-                    points = np.append(points,np.array(self.points[item].points),axis=0)
+                    points = np.append(
+                        points, np.array(self.points[item].points), axis=0
+                    )
                 else:
                     points = np.append(points, self.points[item].points, axis=0)
             point_data = {}
             for key in self.points[point_index[0]].point_data.keys():
-                point_data[key] = np.empty((0,self.points[point_index[0]].point_data[key].shape[1]))
+                point_data[key] = np.empty(
+                    (0, self.points[point_index[0]].point_data[key].shape[1])
+                )
                 for item in point_index:
                     point_data_element = np.array(self.points[item].point_data[key])
-                    point_data[key] = np.append(point_data[key], point_data_element, axis=0)
-                
-
-            combinded_points = meshio.Mesh(points, cells = [("vertex", np.array([[i, ] for i in range(points.shape[0])]))], point_data=point_data)
+                    point_data[key] = np.append(
+                        point_data[key], point_data_element, axis=0
+                    )
+
+            combinded_points = meshio.Mesh(
+                points,
+                cells=[
+                    (
+                        "vertex",
+                        np.array(
+                            [
+                                [
+                                    i,
+                                ]
+                                for i in range(points.shape[0])
+                            ]
+                        ),
+                    )
+                ],
+                point_data=point_data,
+            )
             combinded_points = GF.mesh_transform(combinded_points, self.pose, False)
 
-        
             return combinded_points
 
 
@@ -293,7 +341,7 @@ def rotate_structures(
         Parameters
         ----------
         rotation_matrix : numpy array of appropriate shape (4,4)
-        
+
         rotation_centre : 1*3 numpy float array
             centre of rotation for the structures
 
@@ -302,7 +350,6 @@ def rotate_structures(
         None
         """
 
-
         for item in range(len(self.solids)):
             if self.solids[item] is not None:
                 self.solids[item] = GF.mesh_rotate(
@@ -324,8 +371,7 @@ def translate_structures(self, vector):
         """
         for item in range(len(self.solids)):
             if self.solids[item] is not None:
-                self.solids[item] = GF.translate_mesh(self.solids[item],vector)
-
+                self.solids[item] = GF.translate_mesh(self.solids[item], vector)
 
     def triangles_base_raycaster(self):
         """
@@ -351,7 +397,6 @@ def triangles_base_raycaster(self):
         return triangles
 
 
-
 class antenna_structures(object3d):
     """
     Dedicated class to store information on a specific antenna, including aperture points
@@ -369,11 +414,10 @@ def __init__(self, structures, points):
         super().__init__()
         self.structures = structures
         self.points = points
-        #self.antenna_xyz = o3d.geometry.TriangleMesh.create_coordinate_frame(
-           # size=1, origin=self.pose[:3, 3]
-       # )
-
+        # self.antenna_xyz = o3d.geometry.TriangleMesh.create_coordinate_frame(
+        # size=1, origin=self.pose[:3, 3]
 
+    # )
 
     def export_all_points(self, point_index=None):
         if point_index == None:
@@ -383,7 +427,6 @@ def export_all_points(self, point_index=None):
 
         point_cloud = GF.mesh_transform(point_cloud, self.pose, False)
 
-
         return point_cloud
 
     def excitation_function(
@@ -393,7 +436,7 @@ def excitation_function(
         phase_shift="none",
         wavelength=1.0,
         steering_vector=np.zeros((1, 3)),
-        transmit_power=1.0
+        transmit_power=1.0,
     ):
         # generate the local excitation function and then convert into the global coordinate frame.
         if point_index == None:
@@ -401,9 +444,9 @@ def excitation_function(
         else:
             aperture_points = self.export_all_points(point_index=point_index)
 
-        #as export all points imposes the transformation from local to global frame on the points and associated normal vectors, no rotation is required within calculate_conformalVectors
+        # as export all points imposes the transformation from local to global frame on the points and associated normal vectors, no rotation is required within calculate_conformalVectors
         aperture_weights = EM.calculate_conformalVectors(
-            desired_e_vector, aperture_points.point_data['Normals'], np.eye(3)
+            desired_e_vector, aperture_points.point_data["Normals"], np.eye(3)
         )
         if phase_shift == "wavefront":
             source_points = np.asarray(aperture_points.points)
@@ -411,23 +454,26 @@ def excitation_function(
                 source_points, steering_vector, wavelength
             )
             aperture_weights = aperture_weights * phase_weights.reshape(-1, 1)
-        if phase_shift =="coherence":
+        if phase_shift == "coherence":
             source_points = np.asarray(aperture_points.points)
-            phase_weights=BM.ArbitaryCoherenceWeights(source_points, steering_vector, wavelength
+            phase_weights = BM.ArbitaryCoherenceWeights(
+                source_points, steering_vector, wavelength
             )
             aperture_weights = aperture_weights * phase_weights.reshape(-1, 1)
 
         from .utility.math_functions import discrete_transmit_power
-        if 'Area' in aperture_points.point_data.keys():
-            areas=aperture_points.point_data['Area']
+
+        if "Area" in aperture_points.point_data.keys():
+            areas = aperture_points.point_data["Area"]
         else:
-            areas=np.zeros((aperture_points.points.shape[0]))
-            areas[:]=(wavelength*0.5)**2
+            areas = np.zeros((aperture_points.points.shape[0]))
+            areas[:] = (wavelength * 0.5) ** 2
 
-        calibrated_amplitude_density=discrete_transmit_power(aperture_weights,areas,transmit_power)
+        calibrated_amplitude_density = discrete_transmit_power(
+            aperture_weights, areas, transmit_power
+        )
         return calibrated_amplitude_density
 
-    
     def export_all_structures(self):
         objects = copy.deepcopy(self.structures.solids)
         for item in range(len(objects)):
@@ -437,11 +483,14 @@ def export_all_structures(self):
                 objects[item] = GF.mesh_transform(objects[item], self.pose, False)
 
         return objects
-    def rotate_antenna(self, rotation_vector,rotation_centre=np.zeros((3, 1), dtype=np.float32)):
-        self.structures.rotate_structures(rotation_vector,rotation_centre)
-        self.points.rotate_points(rotation_vector,rotation_centre)
 
-    def translate_antenna(self,translation_vector):
+    def rotate_antenna(
+        self, rotation_vector, rotation_centre=np.zeros((3, 1), dtype=np.float32)
+    ):
+        self.structures.rotate_structures(rotation_vector, rotation_centre)
+        self.points.rotate_points(rotation_vector, rotation_centre)
+
+    def translate_antenna(self, translation_vector):
         self.structures.translate_structures(translation_vector)
         self.points.translate_points(translation_vector)
 
@@ -458,6 +507,7 @@ def pyvista_export(self):
 
         """
         import pyvista as pv
+
         aperture_meshes = []
         structure_meshes = []
         # export and convert structures
@@ -465,7 +515,9 @@ def pyvista_export(self):
 
         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)
+            array = np.append(
+                np.ones((array.shape[0], 1), dtype=np.int32) * 3, array, axis=1
+            )
             return array
 
         for item in triangle_meshes:
@@ -481,9 +533,6 @@ def structure_cells(array):
 
         return aperture_meshes, structure_meshes
 
-  
-
-
 
 class antenna_pattern(object3d):
     """
@@ -571,16 +620,12 @@ def initilise_pattern(self):
             elev_index = np.where(elev_angles > 0)
             self.pattern[elev_index, :, 0] = 1.0
 
-
-
-
-
-
     def export_pyvista_object(self):
         """
         Return Antenna Pattern as a PyVista Structured Mesh Data Object
         """
         import pyvista as pv
+
         def cell_bounds(points, bound_position=0.5):
             """
             Calculate coordinate cell boundaries.
@@ -615,47 +660,73 @@ def cell_bounds(points, bound_position=0.5):
             return bounds
 
         self.transmute_pattern(desired_format="ExEyEz")
-        #hack for cst patterns in spatial intelligence project
-        ex=self.pattern[:,:,0]
-        ey=self.pattern[:,:,1]
-        ez=self.pattern[:,:,2]
-        az=np.linspace(0,360,self.azimuth_resolution)
-        elevation=np.linspace(-90,90,self.elevation_resolution)
+        # hack for cst patterns in spatial intelligence project
+        ex = self.pattern[:, :, 0]
+        ey = self.pattern[:, :, 1]
+        ez = self.pattern[:, :, 2]
+        az = np.linspace(0, 360, self.azimuth_resolution)
+        elevation = np.linspace(-90, 90, self.elevation_resolution)
         xx_bounds = cell_bounds(az)
-        yy_bounds = cell_bounds(90-elevation)
+        yy_bounds = cell_bounds(90 - elevation)
         self.transmute_pattern(desired_format="Etheta/Ephi")
-        et=self.pattern[:,:,0]
-        ep=self.pattern[:,:,1]
-        vista_pattern=pv.grid_from_sph_coords(xx_bounds, yy_bounds,  self.field_radius)
-        vista_pattern['Real']=np.real(np.append(np.append(ex.transpose().ravel().reshape(ex.size,1),
-                                                          ey.transpose().ravel().reshape(ex.size,1),axis=1),
-                                                ez.transpose().ravel().reshape(ex.size,1),axis=1))
-        vista_pattern['Imag']=np.imag(np.append(np.append(ex.transpose().ravel().reshape(ex.size,1),
-                                                          ey.transpose().ravel().reshape(ex.size,1),axis=1),
-                                                ez.transpose().ravel().reshape(ex.size,1),axis=1))
-        vista_pattern['E(theta) Magnitude']=np.abs(et.transpose().ravel())
-        vista_pattern['E(theta) Phase']=np.angle(et.transpose().ravel())
-        vista_pattern['E(phi) Magnitude']=np.abs(ep.transpose().ravel())
-        vista_pattern['E(phi) Phase']=np.angle(ep.transpose().ravel())
-        vista_pattern['Magnitude']=np.abs(vista_pattern['Real']+1j*vista_pattern['Imag'])
-        vista_pattern['Phase']=np.angle(vista_pattern['Real']+1j*vista_pattern['Imag'])
+        et = self.pattern[:, :, 0]
+        ep = self.pattern[:, :, 1]
+        vista_pattern = pv.grid_from_sph_coords(xx_bounds, yy_bounds, self.field_radius)
+        vista_pattern["Real"] = np.real(
+            np.append(
+                np.append(
+                    ex.transpose().ravel().reshape(ex.size, 1),
+                    ey.transpose().ravel().reshape(ex.size, 1),
+                    axis=1,
+                ),
+                ez.transpose().ravel().reshape(ex.size, 1),
+                axis=1,
+            )
+        )
+        vista_pattern["Imag"] = np.imag(
+            np.append(
+                np.append(
+                    ex.transpose().ravel().reshape(ex.size, 1),
+                    ey.transpose().ravel().reshape(ex.size, 1),
+                    axis=1,
+                ),
+                ez.transpose().ravel().reshape(ex.size, 1),
+                axis=1,
+            )
+        )
+        vista_pattern["E(theta) Magnitude"] = np.abs(et.transpose().ravel())
+        vista_pattern["E(theta) Phase"] = np.angle(et.transpose().ravel())
+        vista_pattern["E(phi) Magnitude"] = np.abs(ep.transpose().ravel())
+        vista_pattern["E(phi) Phase"] = np.angle(ep.transpose().ravel())
+        vista_pattern["Magnitude"] = np.abs(
+            vista_pattern["Real"] + 1j * vista_pattern["Imag"]
+        )
+        vista_pattern["Phase"] = np.angle(
+            vista_pattern["Real"] + 1j * vista_pattern["Imag"]
+        )
         return vista_pattern
+
     def pattern_mesh(self):
         points = self.cartesian_points()
         mesh = pv.StructuredGrid(points[:, 0], points[:, 1], points[:, 2])
         mesh.dimensions = [self.azimuth_resolution, self.elevation_resolution, 1]
         if self.arbitary_pattern_format == "Etheta/Ephi":
-            mesh.point_data['Etheta'] = self.pattern[:,:,0]
-            mesh.point_data['Ephi'] = self.pattern[:, :, 1]
+            mesh.point_data["Etheta"] = self.pattern[:, :, 0]
+            mesh.point_data["Ephi"] = self.pattern[:, :, 1]
         elif self.arbitary_pattern_format == "ExEyEz":
-            mesh.point_data['Ex'] = self.pattern[:, :, 0]
-            mesh.point_data['Ey'] = self.pattern[:, :, 1]
-            mesh.point_data['Ez'] = self.pattern[:, :, 2]
+            mesh.point_data["Ex"] = self.pattern[:, :, 0]
+            mesh.point_data["Ey"] = self.pattern[:, :, 1]
+            mesh.point_data["Ez"] = self.pattern[:, :, 2]
 
         return mesh
 
     def display_pattern(
-        self, plottype="Polar", desired_pattern="both", pattern_min=-40, plot_max=0,plotengine="matplotlib"
+        self,
+        plottype="Polar",
+        desired_pattern="both",
+        pattern_min=-40,
+        plot_max=0,
+        plotengine="matplotlib",
     ):
         """
         Displays the Antenna Pattern using :func:`lyceanem.electromagnetics.beamforming.PatternPlot`
@@ -673,18 +744,19 @@ def display_pattern(
         -------
         None
         """
-        if plotengine=="pyvista":
+        if plotengine == "pyvista":
+
             def PatternPlot(
-                    data,
-                    az,
-                    elev,
-                    pattern_min=-40,
-                    plot_max=0.0,
-                    plottype="Polar",
-                    logtype="amplitude",
-                    ticknum=6,
-                    title_text=None,
-                    shell_radius=1.0
+                data,
+                az,
+                elev,
+                pattern_min=-40,
+                plot_max=0.0,
+                plottype="Polar",
+                logtype="amplitude",
+                ticknum=6,
+                title_text=None,
+                shell_radius=1.0,
             ):
                 # points = spherical_mesh(az, elev,
                 # (data / np.nanmax(data) * shell_radius))
@@ -707,6 +779,7 @@ def PatternPlot(
                 pl.show()
 
                 return pl
+
         if self.arbitary_pattern_format == "Etheta/Ephi":
             if desired_pattern == "both":
                 BM.PatternPlot(
diff --git a/lyceanem/electromagnetics/beamforming.py b/lyceanem/electromagnetics/beamforming.py
index 0326709..a7173fe 100644
--- a/lyceanem/electromagnetics/beamforming.py
+++ b/lyceanem/electromagnetics/beamforming.py
@@ -893,7 +893,6 @@ def MaximumfieldMapDiscrete(
     return efield_map
 
 
-
 @njit(cache=True, nogil=True)
 def directivity_transform(
     Etheta,
@@ -1347,21 +1346,30 @@ def GPUBeamformingMap(Etheta, Ephi, DirectivityMap, az_range, el_range, waveleng
         DirectivityMap[az_inc, el_inc] = 0
 
 
-def create_display_mesh(field_data, field_radius=1.0, label='Poynting Vector (Magnitude)', log_type='amplitude',
-                        plot_max=None):
+def create_display_mesh(
+    field_data,
+    field_radius=1.0,
+    label="Poynting Vector (Magnitude)",
+    log_type="amplitude",
+    plot_max=None,
+):
     pattern_mesh = copy.deepcopy(field_data).clean(tolerance=1e-6)
     # pattern_mesh=pv.PolyData(normals).delaunay_2d(inplace=True)
     # pattern_mesh.compute_normals(inplace=True,flip_normals=False) # for some reason the computed normals are all inwards unless I set flip_normals to True
 
-    if log_type == 'amplitude':
+    if log_type == "amplitude":
         log_multiplier = 20.0
-    elif log_type == 'power':
+    elif log_type == "power":
         log_multiplier = 10.0
         # create Normed dB Values for labelled value
     if len(pattern_mesh.point_data[label].shape) > 1:
         # assumes complex data in iq format is to be plotted.
         logdata = log_multiplier * np.log10(
-            np.abs(pattern_mesh.point_data[label][:, 0] + 1j * pattern_mesh.point_data[label][:, 1]))
+            np.abs(
+                pattern_mesh.point_data[label][:, 0]
+                + 1j * pattern_mesh.point_data[label][:, 1]
+            )
+        )
     else:
         logdata = log_multiplier * np.log10(pattern_mesh.point_data[label])
     if np.nanmax(logdata) <= 0.0:
@@ -1379,20 +1387,37 @@ def create_display_mesh(field_data, field_radius=1.0, label='Poynting Vector (Ma
     normed_logdata[normed_logdata < pattern_min] = pattern_min
 
     norm_log = ((normed_logdata - pattern_min) / np.abs(pattern_min)) * field_radius
-    pattern_mesh.point_data['Normed Power'] = norm_log
-    pattern_mesh.point_data['Gain (dB)'] = logdata
-    pattern_mesh.points = pattern_mesh.center + pattern_mesh['Normals'] * norm_log.reshape(-1, 1)
+    pattern_mesh.point_data["Normed Power"] = norm_log
+    pattern_mesh.point_data["Gain (dB)"] = logdata
+    pattern_mesh.points = pattern_mesh.center + pattern_mesh[
+        "Normals"
+    ] * norm_log.reshape(-1, 1)
     return pattern_mesh
-def AnimatedPatterns(pattern_list, command_vectors, field_radius=1.0, plot_scalar='Poynting Vector (Magnitude)',
-                     log_type='amplitude', plot_max=None, filename="BeamformingExport.mp4"):
+
+
+def AnimatedPatterns(
+    pattern_list,
+    command_vectors,
+    field_radius=1.0,
+    plot_scalar="Poynting Vector (Magnitude)",
+    log_type="amplitude",
+    plot_max=None,
+    filename="BeamformingExport.mp4",
+):
     pattern_max = plot_max
 
     pattern_min = pattern_max - 40
     display_list = []
     for inc in range(len(pattern_list)):
         display_list.append(
-            create_display_mesh(pattern_list[inc], field_radius=field_radius, label=plot_scalar, log_type=log_type,
-                                plot_max=plot_max))
+            create_display_mesh(
+                pattern_list[inc],
+                field_radius=field_radius,
+                label=plot_scalar,
+                log_type=log_type,
+                plot_max=plot_max,
+            )
+        )
 
     time_horizontal_fraction = 0.25
     screen_size = 3008
@@ -1400,20 +1425,33 @@ def AnimatedPatterns(pattern_list, command_vectors, field_radius=1.0, plot_scala
     # Open a movie file
     pl.open_movie(filename)
     pl.ren_win.SetSize([screen_size, screen_size])
-    actor = pl.add_mesh(display_list[0], scalars='Gain (dB)', clim=[pattern_min, pattern_max])
+    actor = pl.add_mesh(
+        display_list[0], scalars="Gain (dB)", clim=[pattern_min, pattern_max]
+    )
     steering_vector = pl.add_text(
-        "Command Vector = ({:3.0f},{:3.0f})".format(command_vectors[0, 0], command_vectors[0, 1]),
-        position=(time_horizontal_fraction * screen_size, 0.95 * screen_size), font_size=50)
+        "Command Vector = ({:3.0f},{:3.0f})".format(
+            command_vectors[0, 0], command_vectors[0, 1]
+        ),
+        position=(time_horizontal_fraction * screen_size, 0.95 * screen_size),
+        font_size=50,
+    )
     pl.add_axes()
     pl.write_frame()
     from tqdm import tqdm
+
     for frame in tqdm(range(len(pattern_list))):
         pl.remove_actor(actor)
         pl.remove_actor(steering_vector)
-        actor = pl.add_mesh(display_list[frame], scalars='Gain (dB)', clim=[pattern_min, pattern_max])
+        actor = pl.add_mesh(
+            display_list[frame], scalars="Gain (dB)", clim=[pattern_min, pattern_max]
+        )
         steering_vector = pl.add_text(
-            "Command Vector = ({:3.0f},{:3.0f})".format(command_vectors[frame, 0], command_vectors[frame, 1]),
-            position=(time_horizontal_fraction * screen_size, 0.95 * screen_size), font_size=50)
+            "Command Vector = ({:3.0f},{:3.0f})".format(
+                command_vectors[frame, 0], command_vectors[frame, 1]
+            ),
+            position=(time_horizontal_fraction * screen_size, 0.95 * screen_size),
+            font_size=50,
+        )
         pl.write_frame()
 
-    pl.close()
\ No newline at end of file
+    pl.close()
diff --git a/lyceanem/electromagnetics/data/propagation_constants.py b/lyceanem/electromagnetics/data/propagation_constants.py
index 2e2a846..66f0891 100644
--- a/lyceanem/electromagnetics/data/propagation_constants.py
+++ b/lyceanem/electromagnetics/data/propagation_constants.py
@@ -1,9 +1,11 @@
 from importlib_resources import files
 import lyceanem.electromagnetics.data as data
+
+
 def oxygen_lines():
     data_lines = []
-    oxy_data=str(files(data).joinpath("Oxy.txt"))
-    with open(oxy_data, 'r') as file:
+    oxy_data = str(files(data).joinpath("Oxy.txt"))
+    with open(oxy_data, "r") as file:
         for line in file:
             if line.strip():
                 values = [float(x) for x in line.split()]
@@ -11,14 +13,15 @@ def oxygen_lines():
 
     return data_lines
 
+
 def water_vapour_lines():
 
     data_lines = []
     water_data = str(files(data).joinpath("Vapour.txt"))
-    with open(water_data, 'r') as file:
+    with open(water_data, "r") as file:
         for line in file:
             if line.strip():
                 values = [float(x) for x in line.split()]
                 data_lines.append(values[:7])
 
-    return data_lines
\ No newline at end of file
+    return data_lines
diff --git a/lyceanem/electromagnetics/emfunctions.py b/lyceanem/electromagnetics/emfunctions.py
index b025947..2381123 100644
--- a/lyceanem/electromagnetics/emfunctions.py
+++ b/lyceanem/electromagnetics/emfunctions.py
@@ -2,10 +2,8 @@
 import pyvista as pv
 
 
-
-
-def fresnel_zone(pointA,pointB,wavelength,zone=1):
-    #based on the provided points, wavelength, and zone number, calculate the fresnel zone. This is defined as an ellipsoid for which the difference in distance between the line AB (line of sight), and AP+PB (reflected wave) is a constant multiple of ($n\dfrac{\lambda}{2}$).
+def fresnel_zone(pointA, pointB, wavelength, zone=1):
+    # based on the provided points, wavelength, and zone number, calculate the fresnel zone. This is defined as an ellipsoid for which the difference in distance between the line AB (line of sight), and AP+PB (reflected wave) is a constant multiple of ($n\dfrac{\lambda}{2}$).
     foci = np.stack([pointA, pointB])
     center = np.mean(foci, axis=0)
     ellipsoid = pv.ParametricEllipsoid()
@@ -15,56 +13,112 @@ def fresnel_zone(pointA,pointB,wavelength,zone=1):
     # using binomial approximation via assuming the distance between A and B is much larger than the maximum radius of the zone
     fresnel_radius = 0.5 * ((zone * wavelength * separation) ** 0.5)
     import scipy.spatial.transform as spt
-    pose=np.eye(4)
-    pose[:3, :3] = spt.Rotation.align_vectors(major_axis / separation, [1, 0, 0])[0].as_matrix()
+
+    pose = np.eye(4)
+    pose[:3, :3] = spt.Rotation.align_vectors(major_axis / separation, [1, 0, 0])[
+        0
+    ].as_matrix()
     pose[:3, 3] = center
-    ellipsoid = pv.ParametricEllipsoid(separation * 0.5, fresnel_radius, fresnel_radius).transform(pose)
+    ellipsoid = pv.ParametricEllipsoid(
+        separation * 0.5, fresnel_radius, fresnel_radius
+    ).transform(pose)
     return ellipsoid
 
 
 def field_magnitude_phase(field_data):
     # calculate magnitude and phase representation
 
-    if all(k in field_data.point_data.keys() for k in
-           ('Ex - Real', 'Ex - Imag', 'Ey - Real', 'Ey - Imag', 'Ez - Real', 'Ez - Imag')):
+    if all(
+        k in field_data.point_data.keys()
+        for k in (
+            "Ex - Real",
+            "Ex - Imag",
+            "Ey - Real",
+            "Ey - Imag",
+            "Ez - Real",
+            "Ez - Imag",
+        )
+    ):
         # Exyz exsists in the dataset
-        Ex = field_data.point_data['Ex - Real'] + 1j * field_data.point_data['Ex - Imag']
-        Ey = field_data.point_data['Ey - Real'] + 1j * field_data.point_data['Ey - Imag']
-        Ez = field_data.point_data['Ez - Real'] + 1j * field_data.point_data['Ez - Imag']
-        field_data.point_data['Ex - Magnitude'] = np.abs(Ex)
-        field_data.point_data['Ex - Phase'] = np.angle(Ex)
-        field_data.point_data['Ey - Magnitude'] = np.abs(Ey)
-        field_data.point_data['Ey - Phase'] = np.angle(Ey)
-        field_data.point_data['Ez - Magnitude'] = np.abs(Ez)
-        field_data.point_data['Ez - Phase'] = np.angle(Ez)
-
-    if all(k in field_data.point_data.keys() for k in
-           ('E(theta) - Real', 'E(theta) - Imag', 'E(phi) - Real', 'E(phi) - Imag')):
+        Ex = (
+            field_data.point_data["Ex - Real"] + 1j * field_data.point_data["Ex - Imag"]
+        )
+        Ey = (
+            field_data.point_data["Ey - Real"] + 1j * field_data.point_data["Ey - Imag"]
+        )
+        Ez = (
+            field_data.point_data["Ez - Real"] + 1j * field_data.point_data["Ez - Imag"]
+        )
+        field_data.point_data["Ex - Magnitude"] = np.abs(Ex)
+        field_data.point_data["Ex - Phase"] = np.angle(Ex)
+        field_data.point_data["Ey - Magnitude"] = np.abs(Ey)
+        field_data.point_data["Ey - Phase"] = np.angle(Ey)
+        field_data.point_data["Ez - Magnitude"] = np.abs(Ez)
+        field_data.point_data["Ez - Phase"] = np.angle(Ez)
+
+    if all(
+        k in field_data.point_data.keys()
+        for k in (
+            "E(theta) - Real",
+            "E(theta) - Imag",
+            "E(phi) - Real",
+            "E(phi) - Imag",
+        )
+    ):
         # E(theta) and E(phi) are not present in the data
-        Etheta = field_data.point_data['E(theta) - Real'] + 1j * field_data.point_data['E(theta) - Imag']
-        Ephi = field_data.point_data['E(phi) - Real'] + 1j * field_data.point_data['E(phi) - Imag']
-        field_data.point_data['E(theta) - Magnitude'] = np.abs(Etheta)
-        field_data.point_data['E(theta) - Phase'] = np.angle(Etheta)
-        field_data.point_data['E(phi) - Magnitude'] = np.abs(Ephi)
-        field_data.point_data['E(phi) - Phase'] = np.angle(Ephi)
+        Etheta = (
+            field_data.point_data["E(theta) - Real"]
+            + 1j * field_data.point_data["E(theta) - Imag"]
+        )
+        Ephi = (
+            field_data.point_data["E(phi) - Real"]
+            + 1j * field_data.point_data["E(phi) - Imag"]
+        )
+        field_data.point_data["E(theta) - Magnitude"] = np.abs(Etheta)
+        field_data.point_data["E(theta) - Phase"] = np.angle(Etheta)
+        field_data.point_data["E(phi) - Magnitude"] = np.abs(Ephi)
+        field_data.point_data["E(phi) - Phase"] = np.angle(Ephi)
 
     return field_data
+
+
 def extract_electric_fields(field_data):
-    fields = np.array([field_data.point_data['Ex - Real'][:,0] + 1j * field_data.point_data['Ex - Imag'][:,0],
-                       field_data.point_data['Ey - Real'][:,0] + 1j * field_data.point_data['Ey - Imag'][:,0],
-                       field_data.point_data['Ez - Real'][:,0] + 1j * field_data.point_data['Ez - Imag'][:,0],
-                       ]).transpose()
+    fields = np.array(
+        [
+            field_data.point_data["Ex - Real"][:, 0]
+            + 1j * field_data.point_data["Ex - Imag"][:, 0],
+            field_data.point_data["Ey - Real"][:, 0]
+            + 1j * field_data.point_data["Ey - Imag"][:, 0],
+            field_data.point_data["Ez - Real"][:, 0]
+            + 1j * field_data.point_data["Ez - Imag"][:, 0],
+        ]
+    ).transpose()
 
     return fields
+
+
 def EthetaEphi_to_Exyz(field_data):
-    if not all(k in field_data.point_data.keys() for k in ('theta (Radians)', 'phi (Radians)')):
+    if not all(
+        k in field_data.point_data.keys() for k in ("theta (Radians)", "phi (Radians)")
+    ):
         # theta and phi don't exist in the dataset
         from lyceanem.geometry.geometryfunctions import theta_phi_r
+
         field_data = theta_phi_r(field_data)
     #    # output Ex,Ey,Ez on point data.
-    Espherical = np.array([(field_data.point_data['E(theta) - Real'] + 1j * field_data.point_data['E(theta) - Imag']).ravel(),
-                           (field_data.point_data['E(phi) - Real'] + 1j * field_data.point_data['E(phi) - Imag']).ravel(),
-                           np.zeros((field_data.points.shape[0]))]).transpose()
+    Espherical = np.array(
+        [
+            (
+                field_data.point_data["E(theta) - Real"]
+                + 1j * field_data.point_data["E(theta) - Imag"]
+            ).ravel(),
+            (
+                field_data.point_data["E(phi) - Real"]
+                + 1j * field_data.point_data["E(phi) - Imag"]
+            ).ravel(),
+            np.zeros((field_data.points.shape[0])),
+        ]
+    ).transpose()
     # local_coordinates=field_data.points-field_data.center
     # radial_distance=np.linalg.norm(local_coordinates,axis=1)
     # theta=np.arccos(local_coordinates[:,2]/radial_distance)
@@ -73,55 +127,74 @@ def EthetaEphi_to_Exyz(field_data):
     # etheta=np.zeros(field_data.points.shape[0])
     # ephi=np.zeros(field_data.points.shape[0])
     # conversion matrix from EthetaEphi to Exyz, the inverse operation is via transposing the matrix.
-    ex = Espherical[:, 2] * np.sin(field_data.point_data['theta (Radians)']) * np.cos(
-        field_data.point_data['phi (Radians)']) + Espherical[:, 0] * np.cos(
-        field_data.point_data['theta (Radians)']) * np.cos(field_data.point_data['phi (Radians)']) - Espherical[:,
-                                                                                                     1] * np.sin(
-        field_data.point_data['phi (Radians)'])
-    ey = Espherical[:, 2] * np.sin(field_data.point_data['theta (Radians)']) * np.sin(
-        field_data.point_data['phi (Radians)']) + Espherical[:, 0] * np.cos(
-        field_data.point_data['theta (Radians)']) * np.sin(field_data.point_data['phi (Radians)']) + Espherical[:,
-                                                                                                     1] * np.cos(
-        field_data.point_data['phi (Radians)'])
-    ez = Espherical[:, 2] * np.cos(field_data.point_data['theta (Radians)']) + Espherical[:, 0] * np.sin(
-        field_data.point_data['theta (Radians)'])
-
-    field_data.point_data['Ex - Real'] = np.array([np.real(ex)]).transpose()
-    field_data.point_data['Ex - Imag'] = np.array([np.imag(ex)]).transpose()
-    field_data.point_data['Ey - Real'] = np.array([np.real(ey)]).transpose()
-    field_data.point_data['Ey - Imag'] = np.array([np.imag(ey)]).transpose()
-    field_data.point_data['Ez - Real'] = np.array([np.real(ez)]).transpose()
-    field_data.point_data['Ez - Imag'] = np.array([np.imag(ez)]).transpose()
+    ex = (
+        Espherical[:, 2]
+        * np.sin(field_data.point_data["theta (Radians)"])
+        * np.cos(field_data.point_data["phi (Radians)"])
+        + Espherical[:, 0]
+        * np.cos(field_data.point_data["theta (Radians)"])
+        * np.cos(field_data.point_data["phi (Radians)"])
+        - Espherical[:, 1] * np.sin(field_data.point_data["phi (Radians)"])
+    )
+    ey = (
+        Espherical[:, 2]
+        * np.sin(field_data.point_data["theta (Radians)"])
+        * np.sin(field_data.point_data["phi (Radians)"])
+        + Espherical[:, 0]
+        * np.cos(field_data.point_data["theta (Radians)"])
+        * np.sin(field_data.point_data["phi (Radians)"])
+        + Espherical[:, 1] * np.cos(field_data.point_data["phi (Radians)"])
+    )
+    ez = Espherical[:, 2] * np.cos(
+        field_data.point_data["theta (Radians)"]
+    ) + Espherical[:, 0] * np.sin(field_data.point_data["theta (Radians)"])
+
+    field_data.point_data["Ex - Real"] = np.array([np.real(ex)]).transpose()
+    field_data.point_data["Ex - Imag"] = np.array([np.imag(ex)]).transpose()
+    field_data.point_data["Ey - Real"] = np.array([np.real(ey)]).transpose()
+    field_data.point_data["Ey - Imag"] = np.array([np.imag(ey)]).transpose()
+    field_data.point_data["Ez - Real"] = np.array([np.real(ez)]).transpose()
+    field_data.point_data["Ez - Imag"] = np.array([np.imag(ez)]).transpose()
     return field_data
 
 
 def Exyz_to_EthetaEphi(field_data):
     # this function assumes a spherical field definition, will need to write a function which works based on the poynting vector/normal vector of the point
     electric_fields = extract_electric_fields(field_data)
-    theta = field_data.point_data['theta (Radians)']
-    phi = field_data.point_data['phi (Radians)']
-    etheta = electric_fields[:, 0] * np.cos(phi) * np.cos(theta) + electric_fields[:, 1] * np.sin(phi) * np.cos(
-        theta) - electric_fields[:, 2] * np.sin(theta)
+    theta = field_data.point_data["theta (Radians)"]
+    phi = field_data.point_data["phi (Radians)"]
+    etheta = (
+        electric_fields[:, 0] * np.cos(phi) * np.cos(theta)
+        + electric_fields[:, 1] * np.sin(phi) * np.cos(theta)
+        - electric_fields[:, 2] * np.sin(theta)
+    )
     ephi = -electric_fields[:, 0] * np.sin(phi) + electric_fields[:, 1] * np.cos(phi)
-    field_data.point_data['E(theta) - Real'] = np.zeros((electric_fields.shape[0]))
-    field_data.point_data['E(phi) - Real'] = np.zeros((electric_fields.shape[0]))
-    field_data.point_data['E(theta) - Imag'] = np.zeros((electric_fields.shape[0]))
-    field_data.point_data['E(phi) - Imag'] = np.zeros((electric_fields.shape[0]))
-    field_data.point_data['E(theta) - Real'] = np.real(etheta)
-    field_data.point_data['E(theta) - Imag'] = np.imag(etheta)
-    field_data.point_data['E(phi) - Real'] = np.real(ephi)
-    field_data.point_data['E(phi) - Imag'] = np.imag(ephi)
+    field_data.point_data["E(theta) - Real"] = np.zeros((electric_fields.shape[0]))
+    field_data.point_data["E(phi) - Real"] = np.zeros((electric_fields.shape[0]))
+    field_data.point_data["E(theta) - Imag"] = np.zeros((electric_fields.shape[0]))
+    field_data.point_data["E(phi) - Imag"] = np.zeros((electric_fields.shape[0]))
+    field_data.point_data["E(theta) - Real"] = np.real(etheta)
+    field_data.point_data["E(theta) - Imag"] = np.imag(etheta)
+    field_data.point_data["E(phi) - Real"] = np.real(ephi)
+    field_data.point_data["E(phi) - Imag"] = np.imag(ephi)
     return field_data
 
 
 def field_vectors(field_data):
-    fields = np.array([field_data.point_data['Ex - Real'] + 1j * field_data.point_data['Ex - Imag'],
-                       field_data.point_data['Ey - Real'] + 1j * field_data.point_data['Ey - Imag'],
-                       field_data.point_data['Ez - Real'] + 1j * field_data.point_data['Ez - Imag'],
-                       ]).transpose()
+    fields = np.array(
+        [
+            field_data.point_data["Ex - Real"]
+            + 1j * field_data.point_data["Ex - Imag"],
+            field_data.point_data["Ey - Real"]
+            + 1j * field_data.point_data["Ey - Imag"],
+            field_data.point_data["Ez - Real"]
+            + 1j * field_data.point_data["Ez - Imag"],
+        ]
+    ).transpose()
     directions = np.abs(fields)
     return directions
 
+
 def transform_em(field_data, r):
     """
     Transform the electric current vectors into the new coordinate system defined by the rotation matrix.
@@ -135,114 +208,196 @@ def transform_em(field_data, r):
     -------
 
     """
-    if all(k in field_data.point_data.keys() for k in ("Ex - Real", "Ey - Real", "Ez - Real")):
+    if all(
+        k in field_data.point_data.keys()
+        for k in ("Ex - Real", "Ey - Real", "Ez - Real")
+    ):
         # print("Rotating Ex,Ey,Ez")
-        fields = np.array([field_data.point_data['Ex - Real'] + 1j * field_data.point_data['Ex - Imag'],
-                           field_data.point_data['Ey - Real'] + 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()
+        fields = np.array(
+            [
+                field_data.point_data["Ex - Real"]
+                + 1j * field_data.point_data["Ex - Imag"],
+                field_data.point_data["Ey - Real"]
+                + 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()
         rot_fields = r.apply(fields)
-        field_data.point_data['Ex - Real'] = np.real(
-            rot_fields[:, 0])  # np.array([np.real(rot_fields[:,0]),np.imag(rot_fields[:,0])]).transpose()
-        field_data.point_data['Ey - Real'] = np.real(rot_fields[:, 1])
-        field_data.point_data['Ez - Real'] = np.real(rot_fields[:, 2])
-        field_data.point_data['Ex - Imag'] = np.imag(
-            rot_fields[:, 0])  # np.array([np.real(rot_fields[:,0]),np.imag(rot_fields[:,0])]).transpose()
-        field_data.point_data['Ey - Imag'] = np.imag(rot_fields[:, 1])
-        field_data.point_data['Ez - Imag'] = np.imag(rot_fields[:, 2])
+        field_data.point_data["Ex - Real"] = np.real(
+            rot_fields[:, 0]
+        )  # np.array([np.real(rot_fields[:,0]),np.imag(rot_fields[:,0])]).transpose()
+        field_data.point_data["Ey - Real"] = np.real(rot_fields[:, 1])
+        field_data.point_data["Ez - Real"] = np.real(rot_fields[:, 2])
+        field_data.point_data["Ex - Imag"] = np.imag(
+            rot_fields[:, 0]
+        )  # np.array([np.real(rot_fields[:,0]),np.imag(rot_fields[:,0])]).transpose()
+        field_data.point_data["Ey - Imag"] = np.imag(rot_fields[:, 1])
+        field_data.point_data["Ez - Imag"] = np.imag(rot_fields[:, 2])
         # if all(k in field_data.point_data.keys() for k in ('E(theta)','E(phi)')):
-    elif all(k in field_data.point_data.keys() for k in ('E(theta)', 'E(phi)')):
+    elif all(k in field_data.point_data.keys() for k in ("E(theta)", "E(phi)")):
         # print("Converting Fields and Rotating Ex,Ey,Ez")
         from lyceanem.geometry.geometryfunctions import theta_phi_r
+
         field_data = theta_phi_r(field_data)
         field_data = EthetaEphi_to_Exyz(field_data)
-        fields = np.array([field_data.point_data['Ex - Real'] + 1j * field_data.point_data['Ex - Imag'],
-                           field_data.point_data['Ey - Real'] + 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()
+        fields = np.array(
+            [
+                field_data.point_data["Ex - Real"]
+                + 1j * field_data.point_data["Ex - Imag"],
+                field_data.point_data["Ey - Real"]
+                + 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()
         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])
-        field_data.point_data['Ez - Real'] = np.real(rot_fields[:, 2])
-        field_data.point_data['Ex - Imag'] = np.imag(rot_fields[:, 0])
-        field_data.point_data['Ey - Imag'] = np.imag(rot_fields[:, 1])
-        field_data.point_data['Ez - Imag'] = np.imag(rot_fields[:, 2])
+        field_data.point_data["Ex - Real"] = np.real(rot_fields[:, 0])
+        field_data.point_data["Ey - Real"] = np.real(rot_fields[:, 1])
+        field_data.point_data["Ez - Real"] = np.real(rot_fields[:, 2])
+        field_data.point_data["Ex - Imag"] = np.imag(rot_fields[:, 0])
+        field_data.point_data["Ey - Imag"] = np.imag(rot_fields[:, 1])
+        field_data.point_data["Ez - Imag"] = np.imag(rot_fields[:, 2])
     return field_data
 
-def update_electric_fields(field_data,ex,ey,ez):
-    field_data.point_data['Ex - Real']=np.zeros((field_data.points.shape[0],1))
-    field_data.point_data['Ey - Real']=np.zeros((field_data.points.shape[0],1))
-    field_data.point_data['Ez - Real']=np.zeros((field_data.points.shape[0],1))
-    field_data.point_data['Ex - Imag']=np.zeros((field_data.points.shape[0],1))
-    field_data.point_data['Ey - Imag']=np.zeros((field_data.points.shape[0],1))
-    field_data.point_data['Ez - Imag']=np.zeros((field_data.points.shape[0],1))
-    field_data.point_data['Ex - Real']=np.array([np.real(ex)]).transpose()
-    field_data.point_data['Ex - Imag']=np.array([np.imag(ex)]).transpose()
-    field_data.point_data['Ey - Real']=np.array([np.real(ey)]).transpose()
-    field_data.point_data['Ey - Imag']=np.array([np.imag(ey)]).transpose()
-    field_data.point_data['Ez - Real']=np.array([np.real(ez)]).transpose()
-    field_data.point_data['Ez - Imag']=np.array([np.imag(ez)]).transpose()
+
+def update_electric_fields(field_data, ex, ey, ez):
+    field_data.point_data["Ex - Real"] = np.zeros((field_data.points.shape[0], 1))
+    field_data.point_data["Ey - Real"] = np.zeros((field_data.points.shape[0], 1))
+    field_data.point_data["Ez - Real"] = np.zeros((field_data.points.shape[0], 1))
+    field_data.point_data["Ex - Imag"] = np.zeros((field_data.points.shape[0], 1))
+    field_data.point_data["Ey - Imag"] = np.zeros((field_data.points.shape[0], 1))
+    field_data.point_data["Ez - Imag"] = np.zeros((field_data.points.shape[0], 1))
+    field_data.point_data["Ex - Real"] = np.array([np.real(ex)]).transpose()
+    field_data.point_data["Ex - Imag"] = np.array([np.imag(ex)]).transpose()
+    field_data.point_data["Ey - Real"] = np.array([np.real(ey)]).transpose()
+    field_data.point_data["Ey - Imag"] = np.array([np.imag(ey)]).transpose()
+    field_data.point_data["Ez - Real"] = np.array([np.real(ez)]).transpose()
+    field_data.point_data["Ez - Imag"] = np.array([np.imag(ez)]).transpose()
     return field_data
 
 
 def PoyntingVector(field_data):
-    if all(k in field_data.point_data.keys() for k in (
-    'Permittivity - Real', 'Permittivity - Imag', 'Permeability - Real', 'Permeability - Imag', 'Conductivity')):
-        eta = (field_data.point['Permeability - Real'] + 1j * field_data.point['Permeability - Imag'] /
-               field_data.point['Permittivity - Real'] + 1j * field_data.point['Permittivity - Imag']) ** 0.5
+    if all(
+        k in field_data.point_data.keys()
+        for k in (
+            "Permittivity - Real",
+            "Permittivity - Imag",
+            "Permeability - Real",
+            "Permeability - Imag",
+            "Conductivity",
+        )
+    ):
+        eta = (
+            field_data.point["Permeability - Real"]
+            + 1j
+            * field_data.point["Permeability - Imag"]
+            / field_data.point["Permittivity - Real"]
+            + 1j * field_data.point["Permittivity - Imag"]
+        ) ** 0.5
     else:
         from scipy.constants import physical_constants
-        eta = np.ones((field_data.point_data['Ex - Real'].shape[0])) * \
-              physical_constants['characteristic impedance of vacuum'][0]
+
+        eta = (
+            np.ones((field_data.point_data["Ex - Real"].shape[0]))
+            * physical_constants["characteristic impedance of vacuum"][0]
+        )
 
     electric_field_vectors = np.array(
-        [field_data.point_data['Ex - Real'].ravel() + 1j * field_data.point_data['Ex - Imag'].ravel(),
-         field_data.point_data['Ey - Real'].ravel() + 1j * field_data.point_data['Ey - Imag'].ravel(),
-         field_data.point_data['Ez - Real'].ravel() + 1j * field_data.point_data['Ez - Imag'].ravel()]).transpose()
-    if all(k in field_data.point_data.keys() for k in ('Hx - Real', 'Hy - Real', 'Hz - Real')):
+        [
+            field_data.point_data["Ex - Real"].ravel()
+            + 1j * field_data.point_data["Ex - Imag"].ravel(),
+            field_data.point_data["Ey - Real"].ravel()
+            + 1j * field_data.point_data["Ey - Imag"].ravel(),
+            field_data.point_data["Ez - Real"].ravel()
+            + 1j * field_data.point_data["Ez - Imag"].ravel(),
+        ]
+    ).transpose()
+    if all(
+        k in field_data.point_data.keys()
+        for k in ("Hx - Real", "Hy - Real", "Hz - Real")
+    ):
         # magnetic field data present, so use
-        magnetic_field_vectors = np.array([field_data.point_data['Hx - Real'] + 1j * field_data.point_data['Hx - Imag'],
-                                           field_data.point_data['Hy - Real'] + 1j * field_data.point_data['Hy - Imag'],
-                                           field_data.point_data['Hz - Real'] + 1j * field_data.point_data[
-                                               'Hz - Imag']]).transpose()
+        magnetic_field_vectors = np.array(
+            [
+                field_data.point_data["Hx - Real"]
+                + 1j * field_data.point_data["Hx - Imag"],
+                field_data.point_data["Hy - Real"]
+                + 1j * field_data.point_data["Hy - Imag"],
+                field_data.point_data["Hz - Real"]
+                + 1j * field_data.point_data["Hz - Imag"],
+            ]
+        ).transpose()
         # calculate poynting vector using electric and magnetic field vectors
-        poynting_vector_complex = np.cross(electric_field_vectors, magnetic_field_vectors)
+        poynting_vector_complex = np.cross(
+            electric_field_vectors, magnetic_field_vectors
+        )
     else:
         # use normal vectors instead
         # poynting_vector_complex=field_data.point_data['Normals']*((np.linalg.norm(electric_field_vectors,axis=1)**2)/eta).reshape(-1,1)
-        poynting_vector_complex = ((np.linalg.norm(electric_field_vectors, axis=1) ** 2) / eta).reshape(-1, 1)
-    field_data.point_data['Poynting Vector (Magnitude)'] =np.zeros((field_data.points.shape[0],1))
-    field_data.point_data['Poynting Vector (Magnitude (dB))'] = np.zeros((field_data.points.shape[0], 1))
-    field_data.point_data['Poynting Vector (Magnitude)'] = np.linalg.norm(poynting_vector_complex, axis=1).reshape(-1,1)
-    field_data.point_data['Poynting Vector (Magnitude (dB))'] = 10 * np.log10(
-        np.linalg.norm(poynting_vector_complex.reshape(-1,1), axis=1))
+        poynting_vector_complex = (
+            (np.linalg.norm(electric_field_vectors, axis=1) ** 2) / eta
+        ).reshape(-1, 1)
+    field_data.point_data["Poynting Vector (Magnitude)"] = np.zeros(
+        (field_data.points.shape[0], 1)
+    )
+    field_data.point_data["Poynting Vector (Magnitude (dB))"] = np.zeros(
+        (field_data.points.shape[0], 1)
+    )
+    field_data.point_data["Poynting Vector (Magnitude)"] = np.linalg.norm(
+        poynting_vector_complex, axis=1
+    ).reshape(-1, 1)
+    field_data.point_data["Poynting Vector (Magnitude (dB))"] = 10 * np.log10(
+        np.linalg.norm(poynting_vector_complex.reshape(-1, 1), axis=1)
+    )
     return field_data
 
 
 def Directivity(field_data):
     # calculate directivity for the given pattern
 
-    if not all(k in field_data.point_data.keys() for k in ('theta (Radians)', 'phi (Radians)')):
+    if not all(
+        k in field_data.point_data.keys() for k in ("theta (Radians)", "phi (Radians)")
+    ):
         # theta and phi don't exist in the dataset
         from lyceanem.geometry.geometryfunctions import theta_phi_r
+
         field_data = theta_phi_r(field_data)
 
-    if not all(k in field_data.point_data.keys() for k in ('E(theta) - Real', 'E(phi) - Real')):
+    if not all(
+        k in field_data.point_data.keys() for k in ("E(theta) - Real", "E(phi) - Real")
+    ):
         # E(theta) and E(phi) are not present in the data
         field_data = Exyz_to_EthetaEphi()
 
-    Utheta = np.abs((field_data.point_data['E(theta) - Real'] + 1j * field_data.point_data['E(theta) - Imag']) ** 2)
-    Uphi = np.abs((field_data.point_data['E(phi) - Real'] + 1j * field_data.point_data['E(phi) - Imag']) ** 2)
+    Utheta = np.abs(
+        (
+            field_data.point_data["E(theta) - Real"]
+            + 1j * field_data.point_data["E(theta) - Imag"]
+        )
+        ** 2
+    )
+    Uphi = np.abs(
+        (
+            field_data.point_data["E(phi) - Real"]
+            + 1j * field_data.point_data["E(phi) - Imag"]
+        )
+        ** 2
+    )
     Utotal = Utheta + Uphi
-    sin_factor = np.abs(np.sin(field_data.point_data['theta (Radians)']))  # only want positive factor
+    sin_factor = np.abs(
+        np.sin(field_data.point_data["theta (Radians)"])
+    )  # only want positive factor
     power_sum = np.sum(np.abs(Utheta * sin_factor)) + np.sum(np.abs(Uphi * sin_factor))
     # need to dynamically account for the average contribution of each point, this is only true for a theta step of 1 degree, and phi step of 10 degrees
     Uav = (power_sum * (np.radians(1.0) * np.radians(10.0))) / (4 * np.pi)
     Dtheta = Utheta / Uav
     Dphi = Uphi / Uav
     Dtot = Utotal / Uav
-    field_data.point_data['D(theta)'] = Dtheta
-    field_data.point_data['D(phi)'] = Dphi
-    field_data.point_data['D(Total)'] = Dtot
+    field_data.point_data["D(theta)"] = Dtheta
+    field_data.point_data["D(phi)"] = Dphi
+    field_data.point_data["D(Total)"] = Dtot
 
-    return field_data
\ No newline at end of file
+    return field_data
diff --git a/lyceanem/electromagnetics/empropagation.py b/lyceanem/electromagnetics/empropagation.py
index 2f32eb9..d8b812e 100644
--- a/lyceanem/electromagnetics/empropagation.py
+++ b/lyceanem/electromagnetics/empropagation.py
@@ -14,7 +14,11 @@
 import lyceanem.base_types as base_types
 import lyceanem.geometry.geometryfunctions as GF
 import lyceanem.raycasting.rayfunctions as RF
-from lyceanem.electromagnetics.data.propagation_constants import water_vapour_lines, oxygen_lines
+from lyceanem.electromagnetics.data.propagation_constants import (
+    water_vapour_lines,
+    oxygen_lines,
+)
+
 
 @cuda.jit(device=True)
 def dot(ax1, ay1, az1, ax2, ay2, az2):
@@ -748,23 +752,23 @@ def scatteringkernalv3(
         scattering_matrix[sink_index, 2] = scattering_matrix[sink_index, 2] + (
             ray_component[2] * loss1
         )
+
+
 @cuda.jit(device=True)
-def clip(a,a_min,a_max):
-    if a< a_min:
-        a=a_min
-    elif a> a_max:
-        a=a_max
+def clip(a, a_min, a_max):
+    if a < a_min:
+        a = a_min
+    elif a > a_max:
+        a = a_max
 
     return a
+
+
 @cuda.jit(device=True)
-def lossy_propagation(point1,point2,alpha,beta):
+def lossy_propagation(point1, point2, alpha, beta):
     # calculate loss using improved Rayliegh-Summerfeld
     length = float(0)
-    length=calc_sep(
-        point1,
-        point2,
-        length
-    )
+    length = calc_sep(point1, point2, length)
     outgoing_dir = cuda.local.array(shape=(3), dtype=np.float32)
     calc_dv(
         point1,
@@ -775,19 +779,20 @@ def lossy_propagation(point1,point2,alpha,beta):
     normal[0] = point1["nx"]
     normal[1] = point1["ny"]
     normal[2] = point1["nz"]
-    projection_dot=dot_vec(outgoing_dir,normal)
-    front=-(1/(2*cmath.pi))
-    G=(cmath.exp(-(alpha+1j*beta)*length))/length
-    #dG=cmath.cos(angle)*(-(alpha+1j*beta)-(1/lengths))*G
-    dG=(-(alpha+1j*beta)-(1/length))*G
-    loss=front*dG*projection_dot
-    #loss = G
-
-    #test replacement with old loss funciton
-    #loss = cmath.exp(-1j * beta * lengths)
-    #loss = loss * (((2*cmath.pi)/beta) / (4 * (cmath.pi) * (lengths)))
+    projection_dot = dot_vec(outgoing_dir, normal)
+    front = -(1 / (2 * cmath.pi))
+    G = (cmath.exp(-(alpha + 1j * beta) * length)) / length
+    # dG=cmath.cos(angle)*(-(alpha+1j*beta)-(1/lengths))*G
+    dG = (-(alpha + 1j * beta) - (1 / length)) * G
+    loss = front * dG * projection_dot
+    # loss = G
+
+    # test replacement with old loss funciton
+    # loss = cmath.exp(-1j * beta * lengths)
+    # loss = loss * (((2*cmath.pi)/beta) / (4 * (cmath.pi) * (lengths)))
     return loss
 
+
 @cuda.jit
 def scatteringkernalv4(
     problem_size,
@@ -1069,7 +1074,13 @@ def polaranddistance(network_index, point_information, polar_coefficients, dista
 
 @cuda.jit
 def freqdomainkernal(
-    network_index, point_information, source_sink_index, wavelength, scattering_network, alpha, beta
+    network_index,
+    point_information,
+    source_sink_index,
+    wavelength,
+    scattering_network,
+    alpha,
+    beta,
 ):
     cu_ray_num = cuda.grid(1)  # alias for threadIdx.x + ( blockIdx.x * blockDim.x ),
     #           threadIdx.y + ( blockIdx.y * blockDim.y )
@@ -1169,27 +1180,24 @@ def freqdomainkernal(
         # print(cu_ray_num,source_sink_index[cu_ray_num,0],source_sink_index[cu_ray_num,1])
 
         # scatter_coefficient=(1/(4*cmath.pi))**(complex(scatter_index))
-        #alpha = 0.0
-        #beta = (2.0 * cmath.pi) / wavelength[0]
+        # alpha = 0.0
+        # beta = (2.0 * cmath.pi) / wavelength[0]
         loss = lossy_propagation(
-
-                point_information[network_index[cu_ray_num, 0] - 1],
-                point_information[network_index[cu_ray_num, 1] - 1],
-
-
+            point_information[network_index[cu_ray_num, 0] - 1],
+            point_information[network_index[cu_ray_num, 1] - 1],
             alpha,
-            beta
+            beta,
         )
-        for i in range(1,network_index.shape[1] - 1):
+        for i in range(1, network_index.shape[1] - 1):
             if network_index[cu_ray_num, i + 1] != 0:
 
-                loss *= lossy_propagation(point_information[network_index[cu_ray_num, i] - 1],
-                                          point_information[network_index[cu_ray_num, i + 1] - 1],
-                                          alpha,
-                                          beta
-                                          )
-            
-                
+                loss *= lossy_propagation(
+                    point_information[network_index[cu_ray_num, i] - 1],
+                    point_information[network_index[cu_ray_num, i + 1] - 1],
+                    alpha,
+                    beta,
+                )
+
         ray_component[0] *= loss
         ray_component[1] *= loss
         ray_component[2] *= loss
@@ -1275,19 +1283,20 @@ def freqdomainisokernal(
         alpha = 0.0
         beta = (2.0 * cmath.pi) / wavelength[0]
         loss = lossy_propagation(
-                point_information[network_index[cu_ray_num, 0] - 1],
-                point_information[network_index[cu_ray_num, 1] - 1],
+            point_information[network_index[cu_ray_num, 0] - 1],
+            point_information[network_index[cu_ray_num, 1] - 1],
             alpha,
-            beta
+            beta,
         )
-        for i in range(1,network_index.shape[1] - 1):
+        for i in range(1, network_index.shape[1] - 1):
             if network_index[cu_ray_num, i + 1] != 0:
 
-                loss *= lossy_propagation(point_information[network_index[cu_ray_num, i] - 1],
-                                          point_information[network_index[cu_ray_num, i + 1] - 1],
-                                          alpha,
-                                          beta
-                                          )
+                loss *= lossy_propagation(
+                    point_information[network_index[cu_ray_num, i] - 1],
+                    point_information[network_index[cu_ray_num, i + 1] - 1],
+                    alpha,
+                    beta,
+                )
         ray_component[0] *= loss
         # ray_component[1]=loss
         # ray_component[2]=loss
@@ -1327,8 +1336,8 @@ def timedomainkernal(
     arrival_time,
     wake_time,
     time_map,
-        alpha,
-        beta
+    alpha,
+    beta,
 ):
     # this kernal is planned to calculate the time domain response for a given input signal
     # for flexibility this should probably start out as smn port pairs
@@ -1440,25 +1449,24 @@ def timedomainkernal(
             i = i + 1
 
         # print(cu_ray_num,source_sink_index[cu_ray_num,0],source_sink_index[cu_ray_num,1])
-        #alpha = 0.0
-        #beta = (2.0 * cmath.pi) / wavelength[0]
+        # alpha = 0.0
+        # beta = (2.0 * cmath.pi) / wavelength[0]
         loss = lossy_propagation(
-
-                point_information[full_index[cu_ray_num, 0] - 1],
-                point_information[full_index[cu_ray_num, 1] - 1],
-
+            point_information[full_index[cu_ray_num, 0] - 1],
+            point_information[full_index[cu_ray_num, 1] - 1],
             alpha,
-            beta
+            beta,
         )
-        for i in range(1,full_index.shape[1] - 1):
+        for i in range(1, full_index.shape[1] - 1):
             if full_index[cu_ray_num, i + 1] != 0:
 
-                loss *= lossy_propagation(point_information[full_index[cu_ray_num, i] - 1],
-                                          point_information[full_index[cu_ray_num, i + 1] - 1],
-                                          alpha,
-                                          beta
-                                          )
-        loss_magnitude,loss_phase=cmath.polar(loss)
+                loss *= lossy_propagation(
+                    point_information[full_index[cu_ray_num, i] - 1],
+                    point_information[full_index[cu_ray_num, i + 1] - 1],
+                    alpha,
+                    beta,
+                )
+        loss_magnitude, loss_phase = cmath.polar(loss)
 
         ray_component[0] *= loss_magnitude
         ray_component[1] *= loss_magnitude
@@ -2169,7 +2177,9 @@ def EMGPUJointPathLengthandPolar(source_num, sink_num, full_index, point_informa
     return path_lengths, polar_coefficients
 
 
-def EMGPUFreqDomain(source_num, sink_num, full_index, point_information, wavelength, alpha, beta):
+def EMGPUFreqDomain(
+    source_num, sink_num, full_index, point_information, wavelength, alpha, beta
+):
     """
     Wrapper for the GPU EM processer
     At present, the indexing only supports processing the rays for line of sight and single or double bounces
@@ -2289,7 +2299,7 @@ def EMGPUFreqDomain(source_num, sink_num, full_index, point_information, wavelen
                 d_wavelength,
                 temp_scattering_network,
                 d_alpha,
-                d_beta
+                d_beta,
             )
             # polaranddistance(d_full_index,d_point_information,polar_c,paths)
             # cuda.profile_stop()
@@ -2358,7 +2368,7 @@ def EMGPUFreqDomain(source_num, sink_num, full_index, point_information, wavelen
             d_wavelength,
             d_scattering_network,
             d_alpha,
-            d_beta
+            d_beta,
         )
         # polaranddistance(d_full_index,d_point_information,polar_c,paths)
         # cuda.profile_stop()
@@ -2885,8 +2895,8 @@ def TimeDomainv3(
     excitation_signal,
     sampling_freq,
     num_samples,
-        alpha,
-        beta
+    alpha,
+    beta,
 ):
     """
     New wrapper to run time domain propagation on the GPU, allowing for faster simulations.
@@ -3015,7 +3025,7 @@ def TimeDomainv3(
                 d_wake_time,
                 d_temp_map,
                 d_alpha,
-                d_beta
+                d_beta,
             )
             # print(source_chunking[n],source_chunking[n+1])
 
@@ -3576,9 +3586,9 @@ def EMGPUScatteringWrapper(
                 d_wavelength,
             )
             # cuda.profile_stop()
-            temp_rays[
-                ray_chunks[n] : ray_chunks[n + 1], :
-            ] = d_scatter_matrix.copy_to_host()
+            temp_rays[ray_chunks[n] : ray_chunks[n + 1], :] = (
+                d_scatter_matrix.copy_to_host()
+            )
             # ray_components[ray_chunks[n]:ray_chunks[n+1],:]=d_scatter_matrix.copy_to_host()
             # distmap[source_chunks[n]:source_chunks[n+1],target_chunks[m]:target_chunks[m+1]] = d_distmap_chunked.copy_to_host()
             # first_ray_payload[ray_chunks[n]:ray_chunks[n+1]]=d_chunk_payload.copy_to_host()
@@ -3653,9 +3663,6 @@ def EMGPUWrapper(source_num, sink_num, full_index, point_information, wavelength
     return resultant_rays
 
 
-
-
-
 # @njit(cache=True, nogil=True)
 def vector_mapping(local_E_vector, point_normal, rotation_matrix):
     """
@@ -3680,7 +3687,7 @@ def vector_mapping(local_E_vector, point_normal, rotation_matrix):
     global_vector
 
     """
-    point_vector = np.matmul(point_normal.astype(local_E_vector.dtype),rotation_matrix)
+    point_vector = np.matmul(point_normal.astype(local_E_vector.dtype), rotation_matrix)
     local_axes = np.eye(3)
     uvn_axes = np.zeros((3, 3), dtype=local_E_vector.dtype)
     uvn_axes[2, :] = point_vector
@@ -3695,7 +3702,7 @@ def vector_mapping(local_E_vector, point_normal, rotation_matrix):
     if abs(z_orth) == 0:
         # cannot use z axis as reference, so point normal is aligned with z axis, therefore face_u should be the on the
         # antenna y_axis, therefore face_v can be used to define backwards.
-        uvn_axes[0, :] = np.cross(local_axes[0, :],point_vector) / np.linalg.norm(
+        uvn_axes[0, :] = np.cross(local_axes[0, :], point_vector) / np.linalg.norm(
             np.cross(local_axes[0, :], point_vector)
         )
 
@@ -3850,9 +3857,6 @@ def face_centric_E_vectors(sink_normals, major_axis, scatter_map):
     return new_scatter_map
 
 
-
-
-
 def importDat(fileaddress):
     datafile = pathlib.Path(fileaddress)
     # noinspection PyTypeChecker
@@ -4716,8 +4720,6 @@ def source_transform3to2(
     # return thetaphi_E_vector
 
 
-
-
 def calculate_oxygen_attenuation(frequency, pressure, temperature, oxygen_lines):
     """
     Calculate the specific attenuation due to oxygen using the ITU-R P.676-11 model.
@@ -4737,15 +4739,22 @@ def calculate_oxygen_attenuation(frequency, pressure, temperature, oxygen_lines)
 
     for line in oxygen_lines:
         f_line, a1, a2, a3, a4, a5, a6 = line
-        S = a1 * 10 ** -7 * pressure * theta ** 3 * math.exp(a2 * (1 - theta))
-        ffo = a3 * 10 ** -4 * (pressure * theta ** (0.8 - a4) + 1.1 * pressure * theta)
-        delta = (a5 + a6 * theta) * 10 ** -4 * (pressure) * theta ** 0.8
-        F = (frequency / f_line) * ((ffo - delta * (f_line - frequency)) / ((f_line - frequency) ** 2 + ffo ** 2)+(ffo - delta * (f_line + frequency)) / ((f_line + frequency) ** 2 + ffo ** 2))
+        S = a1 * 10**-7 * pressure * theta**3 * math.exp(a2 * (1 - theta))
+        ffo = a3 * 10**-4 * (pressure * theta ** (0.8 - a4) + 1.1 * pressure * theta)
+        delta = (a5 + a6 * theta) * 10**-4 * (pressure) * theta**0.8
+        F = (frequency / f_line) * (
+            (ffo - delta * (f_line - frequency)) / ((f_line - frequency) ** 2 + ffo**2)
+            + (ffo - delta * (f_line + frequency))
+            / ((f_line + frequency) ** 2 + ffo**2)
+        )
         specific_attenuation += (frequency / f_line) * S * F
 
     return specific_attenuation
 
-def calculate_water_vapor_attenuation(frequency, pressure, temperature, water_vapor_lines):
+
+def calculate_water_vapor_attenuation(
+    frequency, pressure, temperature, water_vapor_lines
+):
     """
     Calculate the specific attenuation due to water vapor using the ITU-R P.676-11 model.
 
@@ -4765,15 +4774,24 @@ def calculate_water_vapor_attenuation(frequency, pressure, temperature, water_va
 
     for line in water_vapor_lines:
         f_line, a1, a2, a3, a4, a5, a6 = line
-        S = a1 * 10 ** -1 * e * theta ** 3.5 * math.exp(a2 * (1 - theta))
-        ffo = a3 * 10 ** -4 * (pressure * theta ** a4 + a5 * e * theta ** a6)
-        F = (frequency / f_line) * ((ffo - 0 * (f_line - frequency)) / ((f_line - frequency) ** 2 + ffo ** 2) + (
-                    ffo - 0 * (f_line + frequency)) / ((f_line + frequency) ** 2 + ffo ** 2))
+        S = a1 * 10**-1 * e * theta**3.5 * math.exp(a2 * (1 - theta))
+        ffo = a3 * 10**-4 * (pressure * theta**a4 + a5 * e * theta**a6)
+        F = (frequency / f_line) * (
+            (ffo - 0 * (f_line - frequency)) / ((f_line - frequency) ** 2 + ffo**2)
+            + (ffo - 0 * (f_line + frequency)) / ((f_line + frequency) ** 2 + ffo**2)
+        )
         specific_attenuation += (frequency / f_line) * S * F
 
     return specific_attenuation
 
-def calculate_total_gaseous_attenuation(frequency, pressure, temperature, oxygen_lines=oxygen_lines(), water_vapor_lines=water_vapour_lines()):
+
+def calculate_total_gaseous_attenuation(
+    frequency,
+    pressure,
+    temperature,
+    oxygen_lines=oxygen_lines(),
+    water_vapor_lines=water_vapour_lines(),
+):
     """
     Calculate the total gaseous attenuation due to both oxygen and water vapor.
 
@@ -4788,9 +4806,15 @@ def calculate_total_gaseous_attenuation(frequency, pressure, temperature, oxygen
     float: The calculated total gaseous attenuation in Np/m.
     """
     # Calculate specific attenuation
-    oxygen_attenuation = calculate_oxygen_attenuation(frequency, pressure, temperature, oxygen_lines)
-    water_vapor_attenuation = calculate_water_vapor_attenuation(frequency, pressure, temperature, water_vapor_lines)
-    specific_attenuation = 0.1820 * frequency * (oxygen_attenuation + water_vapor_attenuation)
+    oxygen_attenuation = calculate_oxygen_attenuation(
+        frequency, pressure, temperature, oxygen_lines
+    )
+    water_vapor_attenuation = calculate_water_vapor_attenuation(
+        frequency, pressure, temperature, water_vapor_lines
+    )
+    specific_attenuation = (
+        0.1820 * frequency * (oxygen_attenuation + water_vapor_attenuation)
+    )
     specific_attenuation = specific_attenuation / (8.686 * 1000)
 
     return specific_attenuation
@@ -4819,7 +4843,8 @@ def calculate_phase_constant(frequency, temperature, pressure, water_vapor_densi
     """
     # Constants
     from scipy.constants import speed_of_light
-    #c = 3e8  # Speed of light in vacuum (m/s)
+
+    # c = 3e8  # Speed of light in vacuum (m/s)
     T0 = 273.15  # Standard temperature in Kelvin
     e_s0 = 611  # Saturation vapor pressure at T0 in Pa
     Lv = 2.5e6  # Latent heat of vaporization of water in J/kg
@@ -4838,7 +4863,7 @@ def calculate_phase_constant(frequency, temperature, pressure, water_vapor_densi
     P = pressure * 100  # Convert hPa to Pa
 
     # Refractivity N(h)
-    N = 77.6 * (P / temperature_K) + (3.73e5 * e) / (temperature_K ** 2)
+    N = 77.6 * (P / temperature_K) + (3.73e5 * e) / (temperature_K**2)
 
     # Refractive index n
     n = 1 + N * 1e-6
@@ -4849,7 +4874,9 @@ def calculate_phase_constant(frequency, temperature, pressure, water_vapor_densi
     return beta
 
 
-def calculate_atmospheric_propagation_constant(frequency, temperature, pressure, water_vapor_density):
+def calculate_atmospheric_propagation_constant(
+    frequency, temperature, pressure, water_vapor_density
+):
     """
     Calculate the propagation constant as a function of frequency (GHz), temperature (Celsius), atmospheric pressure (hectoPascals) and water vapour density (g/m^3).
 
@@ -4869,7 +4896,11 @@ def calculate_atmospheric_propagation_constant(frequency, temperature, pressure,
     propagation constant : complex
 
     """
-    alpha = calculate_attenuation_constant(frequency, temperature, pressure, water_vapor_density)
-    beta = calculate_phase_constant(frequency, temperature, pressure, water_vapor_density)
+    alpha = calculate_attenuation_constant(
+        frequency, temperature, pressure, water_vapor_density
+    )
+    beta = calculate_phase_constant(
+        frequency, temperature, pressure, water_vapor_density
+    )
     gamma = alpha + 1j * beta
-    return gamma
\ No newline at end of file
+    return gamma
diff --git a/lyceanem/geometry/geometryfunctions.py b/lyceanem/geometry/geometryfunctions.py
index b249e94..1d72a64 100644
--- a/lyceanem/geometry/geometryfunctions.py
+++ b/lyceanem/geometry/geometryfunctions.py
@@ -21,7 +21,9 @@ def tri_areas(triangle_mesh):
     Returns:
         np.ndarray: A (N,) array of the area of each triangle.
     """
-    assert triangle_mesh.cells[0].type == "triangle", "Only triangle meshes are supported."
+    assert (
+        triangle_mesh.cells[0].type == "triangle"
+    ), "Only triangle meshes are supported."
     triangle_indices = triangle_mesh.cells[0].data
     v0 = triangle_mesh.points[triangle_indices[:, 0]]
     v1 = triangle_mesh.points[triangle_indices[:, 1]]
@@ -37,17 +39,31 @@ def cell_centroids(field_data):
     """
     for inc, cell in enumerate(field_data.cells):
 
-        if cell.type == 'triangle':
+        if cell.type == "triangle":
             v0 = field_data.points[cell.data[:, 0], :]
             v1 = field_data.points[cell.data[:, 1], :]
             v2 = field_data.points[cell.data[:, 2], :]
             centroids = (1 / 3) * (v0 + v1 + v2)
-            triangle_inc=inc
-
-    centroid_cloud = meshio.Mesh(points=centroids, cells=[("vertex", np.array([[i, ] for i in range(len(centroids))]))])
+            triangle_inc = inc
+
+    centroid_cloud = meshio.Mesh(
+        points=centroids,
+        cells=[
+            (
+                "vertex",
+                np.array(
+                    [
+                        [
+                            i,
+                        ]
+                        for i in range(len(centroids))
+                    ]
+                ),
+            )
+        ],
+    )
     for key in field_data.cell_data.keys():
-        centroid_cloud.point_data[key]=field_data.cell_data[key][triangle_inc]
-    
+        centroid_cloud.point_data[key] = field_data.cell_data[key][triangle_inc]
 
     return centroid_cloud
 
@@ -79,31 +95,47 @@ def mesh_rotate(mesh, rotation, rotation_centre=np.zeros((1, 3), dtype=np.float3
     cell_data = mesh.cell_data
     point_data = mesh.point_data
 
-    if 'Normals' in mesh.point_data:
-        #rotate normals cloud
-        normals = mesh.point_data['Normals']
+    if "Normals" in mesh.point_data:
+        # rotate normals cloud
+        normals = mesh.point_data["Normals"]
         rotated_normals = r.apply(normals)
-        point_data['Normals'] = rotated_normals
-    if 'Normals' in mesh.cell_data:
-        #rotate normals cloud
-        for i, rotated_normals in enumerate(mesh.cell_data['Normals']):
+        point_data["Normals"] = rotated_normals
+    if "Normals" in mesh.cell_data:
+        # rotate normals cloud
+        for i, rotated_normals in enumerate(mesh.cell_data["Normals"]):
             rotated_normals = r.apply(rotated_normals)
-            cell_data['Normals'][i] = rotated_normals
+            cell_data["Normals"][i] = rotated_normals
 
     mesh_return = meshio.Mesh(points=rotated_points, cells=mesh.cells)
     mesh_return.point_data = point_data
     mesh_return.cell_data = cell_data
 
     # if field data is present, rotate fields
-    if all(k in mesh.point_data.keys() for k in ("Ex - Real", "Ey - Real", "Ez - Real")):
-        #rotate field data
-        fields=transform_em(copy.deepcopy(mesh),r)
-        for key in ("Ex - Real","Ex - Imag", "Ey - Real","Ey - Imag", "Ez - Real","Ez - Imag"):
-            mesh_return.point_data[key]=fields.point_data[key]
-
-    elif all(k in mesh.point_data.keys() for k in ('E(theta)', 'E(phi)')):
-        fields=transform_em(copy.deepcopy(mesh),r)
-        for key in ("Ex - Real", "Ex - Imag", "Ey - Real", "Ey - Imag", "Ez - Real", "Ez - Imag"):
+    if all(
+        k in mesh.point_data.keys() for k in ("Ex - Real", "Ey - Real", "Ez - Real")
+    ):
+        # rotate field data
+        fields = transform_em(copy.deepcopy(mesh), r)
+        for key in (
+            "Ex - Real",
+            "Ex - Imag",
+            "Ey - Real",
+            "Ey - Imag",
+            "Ez - Real",
+            "Ez - Imag",
+        ):
+            mesh_return.point_data[key] = fields.point_data[key]
+
+    elif all(k in mesh.point_data.keys() for k in ("E(theta)", "E(phi)")):
+        fields = transform_em(copy.deepcopy(mesh), r)
+        for key in (
+            "Ex - Real",
+            "Ex - Imag",
+            "Ey - Real",
+            "Ey - Imag",
+            "Ez - Real",
+            "Ez - Imag",
+        ):
             mesh_return.point_data[key] = fields.point_data[key]
 
     return mesh_return
@@ -113,20 +145,27 @@ def mesh_transform(mesh, transform_matrix, rotate_only):
     return_mesh = mesh
     if rotate_only:
         for i in range(mesh.points.shape[0]):
-            return_mesh.points[i] = np.dot(transform_matrix, np.append(mesh.points[i], 0))[:3]
-            return_mesh.point_data['Normals'][i] = np.dot(transform_matrix,
-                                                          np.append(mesh.point_data['Normals'][i], 0))[:3]
+            return_mesh.points[i] = np.dot(
+                transform_matrix, np.append(mesh.points[i], 0)
+            )[:3]
+            return_mesh.point_data["Normals"][i] = np.dot(
+                transform_matrix, np.append(mesh.point_data["Normals"][i], 0)
+            )[:3]
 
     else:
         for i in range(mesh.points.shape[0]):
-            return_mesh.points[i] = np.dot(transform_matrix, np.append(mesh.points[i], 1))[:3]
-            return_mesh.point_data['Normals'][i] = np.dot(transform_matrix,
-                                                          np.append(mesh.point_data['Normals'][i], 0))[:3]
-        if 'Normals' in mesh.cell_data:
-            for i in range(len(mesh.cell_data['Normals'])):
-                for j in range(mesh.cell_data['Normals'][i].shape[0]):
-                    return_mesh.cell_data['Normals'][i][j] = np.dot(transform_matrix,
-                                                                    np.append(mesh.cell_data['Normals'][i][j], 0))[:3]
+            return_mesh.points[i] = np.dot(
+                transform_matrix, np.append(mesh.points[i], 1)
+            )[:3]
+            return_mesh.point_data["Normals"][i] = np.dot(
+                transform_matrix, np.append(mesh.point_data["Normals"][i], 0)
+            )[:3]
+        if "Normals" in mesh.cell_data:
+            for i in range(len(mesh.cell_data["Normals"])):
+                for j in range(mesh.cell_data["Normals"][i].shape[0]):
+                    return_mesh.cell_data["Normals"][i][j] = np.dot(
+                        transform_matrix, np.append(mesh.cell_data["Normals"][i][j], 0)
+                    )[:3]
 
     return return_mesh
 
@@ -135,42 +174,70 @@ def compute_areas(field_data):
     cell_areas = []
     for inc, cell in enumerate(field_data.cells):
 
-        if cell.type == 'triangle':
+        if cell.type == "triangle":
             # Heron's Formula
-            edge1 = np.linalg.norm(field_data.points[cell.data[:, 0], :] - field_data.points[cell.data[:, 1], :],
-                                   axis=1)
-            edge2 = np.linalg.norm(field_data.points[cell.data[:, 1], :] - field_data.points[cell.data[:, 2], :],
-                                   axis=1)
-            edge3 = np.linalg.norm(field_data.points[cell.data[:, 2], :] - field_data.points[cell.data[:, 0], :],
-                                   axis=1)
+            edge1 = np.linalg.norm(
+                field_data.points[cell.data[:, 0], :]
+                - field_data.points[cell.data[:, 1], :],
+                axis=1,
+            )
+            edge2 = np.linalg.norm(
+                field_data.points[cell.data[:, 1], :]
+                - field_data.points[cell.data[:, 2], :],
+                axis=1,
+            )
+            edge3 = np.linalg.norm(
+                field_data.points[cell.data[:, 2], :]
+                - field_data.points[cell.data[:, 0], :],
+                axis=1,
+            )
             s = (edge1 + edge2 + edge3) / 2
             areas = (s * (s - edge1) * (s - edge2) * (s - edge3)) ** 0.5
             cell_areas.append(areas)
-        if cell.type == 'quad':
+        if cell.type == "quad":
             # Heron's Formula twice
-            edge1 = np.linalg.norm(field_data.points[cell.data[:, 0], :] - field_data.points[cell.data[:, 1], :],
-                                   axis=1)
-            edge2 = np.linalg.norm(field_data.points[cell.data[:, 1], :] - field_data.points[cell.data[:, 2], :],
-                                   axis=1)
-            edge3 = np.linalg.norm(field_data.points[cell.data[:, 2], :] - field_data.points[cell.data[:, 0], :],
-                                   axis=1)
-            edge4 = np.linalg.norm(field_data.points[cell.data[:, 2], :] - field_data.points[cell.data[:, 3], :],
-                                   axis=1)
-            edge5 = np.linalg.norm(field_data.points[cell.data[:, 3], :] - field_data.points[cell.data[:, 0], :],
-                                   axis=1)
+            edge1 = np.linalg.norm(
+                field_data.points[cell.data[:, 0], :]
+                - field_data.points[cell.data[:, 1], :],
+                axis=1,
+            )
+            edge2 = np.linalg.norm(
+                field_data.points[cell.data[:, 1], :]
+                - field_data.points[cell.data[:, 2], :],
+                axis=1,
+            )
+            edge3 = np.linalg.norm(
+                field_data.points[cell.data[:, 2], :]
+                - field_data.points[cell.data[:, 0], :],
+                axis=1,
+            )
+            edge4 = np.linalg.norm(
+                field_data.points[cell.data[:, 2], :]
+                - field_data.points[cell.data[:, 3], :],
+                axis=1,
+            )
+            edge5 = np.linalg.norm(
+                field_data.points[cell.data[:, 3], :]
+                - field_data.points[cell.data[:, 0], :],
+                axis=1,
+            )
 
             s1 = (edge1 + edge2 + edge3) / 2
             s2 = (edge3 + edge4 + edge5) / 2
             areas = (s1 * (s1 - edge1) * (s1 - edge2) * (s1 - edge3)) ** 0.5 + (
-                    s2 * (s2 - edge3) * (s2 - edge4) * (s2 - edge5)) ** 0.5
+                s2 * (s2 - edge3) * (s2 - edge4) * (s2 - edge5)
+            ) ** 0.5
             cell_areas.append(areas)
 
-    field_data.cell_data['Area'] = cell_areas
-    field_data.point_data['Area'] = np.zeros((field_data.points.shape[0]))
+    field_data.cell_data["Area"] = cell_areas
+    field_data.point_data["Area"] = np.zeros((field_data.points.shape[0]))
     for inc, cell in enumerate(field_data.cells):
         for point_inc in range(field_data.points.shape[0]):
-            field_data.point_data['Area'][point_inc] = np.mean(
-                field_data.cell_data['Area'][inc][np.where(field_data.cells[inc].data == point_inc)[0]])
+            field_data.point_data["Area"][point_inc] = np.mean(
+                field_data.cell_data["Area"][inc][
+                    np.where(field_data.cells[inc].data == point_inc)[0]
+                ]
+            )
 
     return field_data
 
@@ -190,52 +257,66 @@ def compute_normals(mesh):
     """
     cell_normal_list = []
     for inc, cell in enumerate(mesh.cells):
-        #print(cell.type, cell.data.shape[0])
-        if cell.type == 'vertex':
+        # print(cell.type, cell.data.shape[0])
+        if cell.type == "vertex":
             # assume outward pointing normals from centroid
-            vertex_normals = mesh.points[cell.data[:, 0], :]/np.linalg.norm(mesh.points[cell.data[:, 0], :],axis=1).reshape(-1,1)
+            vertex_normals = mesh.points[cell.data[:, 0], :] / np.linalg.norm(
+                mesh.points[cell.data[:, 0], :], axis=1
+            ).reshape(-1, 1)
             cell_normal_list.append(vertex_normals)
-        if cell.type == 'line':
+        if cell.type == "line":
             line_normals = np.zeros((cell.data.shape[0], 3))
             cell_normal_list.append(line_normals)
-        if cell.type == 'triangle':
+        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)
+            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':
+        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)
+            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
+    mesh.cell_data["Normals"] = cell_normal_list
 
-    #calculate vertex normals
+    # calculate vertex normals
     point_normals = np.empty((0, 3))
     for inc, cell in enumerate(mesh.cells):
-        if cell.type == 'triangle':
+        if cell.type == "triangle":
             for inc in range(mesh.points.shape[0]):
                 associated_cells = np.where(inc == cell.data)[0]
-                #print(associated_cells)
-                point_normals = np.append(point_normals,
-                                          np.mean(mesh.cell_data['Normals'][0][associated_cells, :], axis=0).reshape(1,
-                                                                                                                     3),
-                                          axis=0)
-        if cell.type == 'vertex':
+                # print(associated_cells)
+                point_normals = np.append(
+                    point_normals,
+                    np.mean(
+                        mesh.cell_data["Normals"][0][associated_cells, :], axis=0
+                    ).reshape(1, 3),
+                    axis=0,
+                )
+        if cell.type == "vertex":
             for inc in range(mesh.points.shape[0]):
                 associated_cells = np.where(inc == cell.data)[0]
-                #print(associated_cells)
-                point_normals = np.append(point_normals,
-                                          np.mean(mesh.cell_data['Normals'][0][associated_cells, :], axis=0).reshape(1,
-                                                                                                                     3),
-                                          axis=0)
+                # print(associated_cells)
+                point_normals = np.append(
+                    point_normals,
+                    np.mean(
+                        mesh.cell_data["Normals"][0][associated_cells, :], axis=0
+                    ).reshape(1, 3),
+                    axis=0,
+                )
 
-    mesh.point_data['Normals'] = point_normals / np.linalg.norm(point_normals, axis=1).reshape(-1, 1)
+    mesh.point_data["Normals"] = point_normals / np.linalg.norm(
+        point_normals, axis=1
+    ).reshape(-1, 1)
 
     return mesh
 
@@ -252,10 +333,15 @@ def theta_phi_r(field_data):
     -------
 
     """
-    field_data.point_data['Radial Distance (m)'] = np.linalg.norm(field_data.points - np.zeros((1, 3)), axis=1)
-    field_data.point_data['theta (Radians)'] = np.arccos(
-        field_data.points[:, 2] / field_data.point_data['Radial Distance (m)'])
-    field_data.point_data['phi (Radians)'] = np.arctan2(field_data.points[:, 1], field_data.points[:, 0])
+    field_data.point_data["Radial Distance (m)"] = np.linalg.norm(
+        field_data.points - np.zeros((1, 3)), axis=1
+    )
+    field_data.point_data["theta (Radians)"] = np.arccos(
+        field_data.points[:, 2] / field_data.point_data["Radial Distance (m)"]
+    )
+    field_data.point_data["phi (Radians)"] = np.arctan2(
+        field_data.points[:, 1], field_data.points[:, 0]
+    )
     return field_data
 
 
@@ -355,9 +441,9 @@ def calculate_align_mat(pVec_Arr):
     z_c_vec_mat = get_cross_prod_mat(z_c_vec)
 
     qTrans_Mat = (
-            np.eye(3, 3)
-            + z_c_vec_mat
-            + np.matmul(z_c_vec_mat, z_c_vec_mat) / (1 + np.dot(z_unit_Arr, pVec_Arr))
+        np.eye(3, 3)
+        + z_c_vec_mat
+        + np.matmul(z_c_vec_mat, z_c_vec_mat) / (1 + np.dot(z_unit_Arr, pVec_Arr))
     )
     qTrans_Mat *= scale
     return qTrans_Mat
diff --git a/lyceanem/geometry/parabolas.scad b/lyceanem/geometry/parabolas.scad
deleted file mode 100644
index 21c5ea5..0000000
--- a/lyceanem/geometry/parabolas.scad
+++ /dev/null
@@ -1,51 +0,0 @@
-//parabola tests
-module openscad_paraboloid (y=10, f=5, rfa=0, fc=1, detail=44){
-	// y = height of paraboloid
-	// f = focus distance 
-	// fc : 1 = center paraboloid in focus point(x=0, y=f); 0 = center paraboloid on top (x=0, y=0)
-	// rfa = radius of the focus area : 0 = point focus
-	// detail = $fn of cone
-
-	hi = (y+2*f)/sqrt(2);								// height and radius of the cone -> alpha = 45° -> sin(45°)=1/sqrt(2)
-	x =2*f*sqrt(y/f);									// x  = half size of parabola
-	
-   translate([0,0,-f*fc])								// center on focus 
-	rotate_extrude(convexity = 10,$fn=detail )		// extrude paraboild
-	translate([rfa,0,0])								// translate for fokus area	 
-	difference(){
-		union(){											// adding square for focal area
-			projection(cut = true)																			// reduce from 3D cone to 2D parabola
-				translate([0,0,f*2]) rotate([45,0,0])													// rotate cone 45° and translate for cutting
-				translate([0,0,-hi/2])cylinder(h= hi, r1=hi, r2=0, center=true, $fn=detail);   	// center cone on tip
-			translate([-(rfa+x ),0]) square ([rfa+x , y ]);											// focal area square
-		}
-		translate([-(2*rfa+x ), -1/2]) square ([rfa+x ,y +1] ); 					// cut of half at rotation center 
-	}
-}
-
-translate([0,0,0]){
-   difference(){
-openscad_paraboloid (y=0.07239,f=0.146,rfa= 0,fc=0,detail=120);
-translate([0,0,0.005])openscad_paraboloid (y=0.07239,f=0.146,rfa= 0,fc=0,detail=120);
-}
-}
-//detail=200;
-//rfa=0;
-//f=0.146;
-//y=0.07239;
-//hi = (y+2*f)/sqrt(2);	
-////hi=0.2;
-//fc=0;
-//x =2*f*sqrt(y/f);
-//   translate([0,0,-f*fc])								// center on focus 
-//	rotate_extrude(convexity = 10,$fn=detail )		// extrude paraboild
-//	translate([rfa,0,0])								// translate for fokus area	 
-//	difference(){
-//		union(){											// adding square for focal area
-//			projection(cut = true)																			// reduce from 3D cone to 2D parabola
-//				translate([0,0,f*2]) rotate([45,0,0])													// rotate cone 45° and translate for cutting
-//				translate([0,0,-hi/2])cylinder(h= hi, r1=hi, r2=0, center=true, $fn=detail);   	// center cone on tip
-//			translate([-(rfa+x ),0]) square ([rfa+x , y ]);											// focal area square
-//		}
-//		translate([-(2*rfa+x ), -1/2]) square ([rfa+x ,y +1] ); 					// cut of half at rotation center 
-//	}
diff --git a/lyceanem/geometry/targets.py b/lyceanem/geometry/targets.py
index dd35877..8d8e567 100755
--- a/lyceanem/geometry/targets.py
+++ b/lyceanem/geometry/targets.py
@@ -99,34 +99,31 @@ def converttostl():
         p.wait()
 
 
-
-
-
 def rectReflector(majorsize, minorsize, thickness):
     """
     create a primative of the right size, assuming always orientated
     with normal aligned with zenith, and major axis with x,
     adjust position so the face is centred on (0,0,0)
     """
-    print("majorsize",majorsize)
-    print("minorsize",minorsize)
-    print("thickness",thickness)
+    print("majorsize", majorsize)
+    print("minorsize", minorsize)
+    print("thickness", thickness)
 
     halfMajor = majorsize / 2.0
     halfMinor = minorsize / 2.0
-    pv_mesh = pv.Box((-halfMajor, halfMajor, -halfMinor, halfMinor, -(thickness+EPSILON),- EPSILON))
+    pv_mesh = pv.Box(
+        (-halfMajor, halfMajor, -halfMinor, halfMinor, -(thickness + EPSILON), -EPSILON)
+    )
     pv_mesh = pv_mesh.triangulate()
-    pv_mesh.compute_normals(inplace=True,consistent_normals=False)
-    triangles = np.reshape(np.array(pv_mesh.faces),(12,4))
-    triangles = triangles[:,1:]
+    pv_mesh.compute_normals(inplace=True, consistent_normals=False)
+    triangles = np.reshape(np.array(pv_mesh.faces), (12, 4))
+    triangles = triangles[:, 1:]
 
     mesh = meshio.Mesh(pv_mesh.points, {"triangle": triangles})
 
-
-
     mesh.point_data["Normals"] = pv_mesh.point_normals
     mesh.cell_data["Normals"] = pv_mesh.cell_normals
-    red = np.zeros((8, 1), dtype=np.float32) 
+    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
 
@@ -191,17 +188,18 @@ def shapeTrapezoid(x_size, y_size, length, flare_angle):
     triangle_list[11, :] = [6, 5, 7]
     mesh = meshio.Mesh(
         points=mesh_vertices,
-        cells=[("triangle", triangle_list)],)
-    #print(mesh)
+        cells=[("triangle", triangle_list)],
+    )
+    # print(mesh)
     triangle_list = np.insert(triangle_list, 0, 3, axis=1)
-    
-    pv_mesh = pv.PolyData( mesh_vertices, faces = triangle_list)
-    pv_mesh.compute_normals(inplace=True,consistent_normals=False)
+
+    pv_mesh = pv.PolyData(mesh_vertices, faces=triangle_list)
+    pv_mesh.compute_normals(inplace=True, consistent_normals=False)
 
     mesh.point_data["Normals"] = np.asarray(pv_mesh.point_normals)
     mesh.cell_data["Normals"] = [np.asarray(pv_mesh.cell_normals)]
 
-    red = np.zeros((8, 1), dtype=np.float32) 
+    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
 
@@ -211,7 +209,6 @@ def shapeTrapezoid(x_size, y_size, length, flare_angle):
     return mesh
 
 
-
 def meshedReflector(majorsize, minorsize, thickness, grid_resolution, sides="all"):
     """
 
@@ -250,10 +247,6 @@ def meshedReflector(majorsize, minorsize, thickness, grid_resolution, sides="all
     return reflector, mesh_points
 
 
-
-
-
-
 def meshedHorn(
     majorsize,
     minorsize,
@@ -291,7 +284,7 @@ def meshedHorn(
     mesh_points : meshio object
         the source points for the horn aperture
     """
-    #print("HIHIH")
+    # print("HIHIH")
     structure = shapeTrapezoid(
         majorsize + (edge_width * 2), minorsize + (edge_width * 2), length, flare_angle
     )
@@ -299,29 +292,50 @@ def meshedHorn(
         majorsize, minorsize, 1e-6, grid_resolution, sides
     )
 
-    mesh_points = GF.translate_mesh(mesh_points, [0, 0, EPSILON*2])
+    mesh_points = GF.translate_mesh(mesh_points, [0, 0, EPSILON * 2])
 
     return structure, mesh_points
 
 
-def parabolic_aperture(diameter, focal_length, thickness, mesh_size, sides='front', lip=False, lip_height=1e-3,
-                       lip_width=1e-3):
+def parabolic_aperture(
+    diameter,
+    focal_length,
+    thickness,
+    mesh_size,
+    sides="front",
+    lip=False,
+    lip_height=1e-3,
+    lip_width=1e-3,
+):
     # Define function for parabola equation (y^2 = 4*focal_length*x)
     import lyceanem.geometry.geometryfunctions as GF
     import lyceanem.utility.math_functions as math_functions
+
     def parabola(x):
-        return (1 / (4 * focal_length)) * x ** 2
+        return (1 / (4 * focal_length)) * x**2
 
     with pygmsh.occ.Geometry() as geom:
         geom.characteristic_length_max = mesh_size * 0.8
         # 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))])
+        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])
@@ -330,13 +344,26 @@ def parabola(x):
 
         # 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)
+        _, 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)
+            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)
 
         if lip:
@@ -344,7 +371,9 @@ def parabola(x):
 
             start_point = np.array([0.0, 0.0, parabola(diameter * 0.5)])
             cylinder1 = geom.add_cylinder(start_point.ravel(), axis1, diameter / 2)
-            cylinder2 = geom.add_cylinder(start_point.ravel(), axis1, diameter / 2 + lip_width)
+            cylinder2 = geom.add_cylinder(
+                start_point.ravel(), axis1, diameter / 2 + lip_width
+            )
             final = geom.boolean_difference([cylinder2], [cylinder1])
             volume_list.append(final)
 
@@ -352,10 +381,11 @@ def parabola(x):
 
         mesh_temp = geom.generate_mesh(dim=2)
     for inc, cell in enumerate(mesh_temp.cells):
-        if cell.type == 'triangle':
+        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)
     mesh = GF.compute_normals(mesh)
@@ -364,7 +394,7 @@ def parabola(x):
         (diameter / 2),
         int(np.max(np.asarray([2, np.ceil((diameter * 0.5) / (mesh_size))]))),
     )
-    z_space = (1 / (4 * focal_length)) * x_space ** 2
+    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(
         [
@@ -415,9 +445,23 @@ def parabola(x):
 
     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})
+    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
 
@@ -569,8 +613,22 @@ def gridedReflectorPoints(
         mesh_vertices = source_coords
         mesh_normals = source_normals
 
-    mesh_points = meshio.Mesh(points=mesh_vertices,
-                              cells=[("vertex", np.array([[i, ] for i in range(len(mesh_vertices))]))],
-                              point_data={"Normals": mesh_normals})
+    mesh_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_points
diff --git a/lyceanem/models/frequency_domain.py b/lyceanem/models/frequency_domain.py
index 8b1b998..44c70fb 100644
--- a/lyceanem/models/frequency_domain.py
+++ b/lyceanem/models/frequency_domain.py
@@ -17,7 +17,7 @@ def aperture_projection(
     wavelength=1.0,
     az_range=np.linspace(-180.0, 180.0, 19),
     elev_range=np.linspace(-90.0, 90.0, 19),
-    farfield_distance=2.0
+    farfield_distance=2.0,
 ):
     """
 
@@ -45,9 +45,9 @@ def aperture_projection(
         a point cloud colored according to the projected area, normalised to the total projected area of the aperture.
     """
     if environment is None:
-        blocking_triangles=GF.mesh_conversion(aperture)
+        blocking_triangles = GF.mesh_conversion(aperture)
     else:
-        blocking_triangles=GF.mesh_conversion(environment)
+        blocking_triangles = GF.mesh_conversion(environment)
 
     directivity_envelope = np.zeros(
         (elev_range.shape[0], az_range.shape[0]), dtype=np.float32
@@ -62,9 +62,9 @@ def aperture_projection(
         vertex_area=triangle_areas,
         az_range=az_range,
         elev_range=elev_range,
-        shell_range=farfield_distance
+        shell_range=farfield_distance,
     )
-    directivity_envelope[:, :] = (4 * np.pi * visible_patterns) / (wavelength ** 2)
+    directivity_envelope[:, :] = (4 * np.pi * visible_patterns) / (wavelength**2)
 
     return directivity_envelope, pcd
 
@@ -86,8 +86,8 @@ def calculate_farfield(
     los=True,
     project_vectors=False,
     antenna_axes=np.eye(3),
-        alpha=0.0,
-        beta=(np.pi*2)/1.0
+    alpha=0.0,
+    beta=(np.pi * 2) / 1.0,
 ):
     """
     Based upon the aperture coordinates and solids, predict the farfield for the antenna.
@@ -148,15 +148,20 @@ def calculate_farfield(
     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)
+    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
+            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)
+        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
@@ -164,7 +169,7 @@ def calculate_farfield(
 
     if scattering == 0:
         # only use the aperture point cloud, no scattering required.
-        scatter_points = meshio.Mesh(points= np.empty((0, 3)), cells=[])
+        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),
@@ -178,14 +183,14 @@ def calculate_farfield(
         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
+                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_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
@@ -209,15 +214,15 @@ def calculate_farfield(
         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[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]
@@ -254,16 +259,16 @@ def calculate_farfield(
         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
+                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
+        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
@@ -298,15 +303,15 @@ def calculate_farfield(
         # 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_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]
@@ -363,19 +368,25 @@ def calculate_farfield(
             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
-            )
+            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
+                num_sources,
+                sinks.shape[0],
+                full_index,
+                point_informationv2,
+                wavelength,
+                alpha,
+                beta,
             )
             Ex[element, :, :] = np.dot(
                 np.ones((num_sources)), scatter_map[:, :, 0]
@@ -405,7 +416,13 @@ def calculate_farfield(
         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
+            num_sources,
+            num_sinks,
+            full_index,
+            point_informationv2,
+            wavelength,
+            alpha,
+            beta,
         )
 
         Ex[:, :] = np.sum(scatter_map[:, :, 0], axis=0).reshape(
@@ -442,8 +459,8 @@ def calculate_scattering(
     project_vectors=False,
     antenna_axes=np.eye(3),
     multiE=False,
-        alpha=0.0,
-        beta=(np.pi*2)/1.0
+    alpha=0.0,
+    beta=(np.pi * 2) / 1.0,
 ):
     """
     calculating the scattering from the provided source coordinates, to the provided sink coordinates in the environment.
@@ -485,37 +502,51 @@ def calculate_scattering(
     num_sources = len(np.asarray(aperture_coords.points))
     num_sinks = len(np.asarray(sink_coords.points))
 
-    environment_triangles=GF.mesh_conversion(antenna_solid)
+    environment_triangles = GF.mesh_conversion(antenna_solid)
 
     if not multiE:
         if project_vectors:
             conformal_E_vectors = EM.calculate_conformalVectors(
-                desired_E_axis, np.asarray(aperture_coords.point_data["Normals"]), antenna_axes
+                desired_E_axis,
+                np.asarray(aperture_coords.point_data["Normals"]),
+                antenna_axes,
             )
         else:
-            #print("hi from here", aperture_coords.cell_data)
-            if desired_E_axis.shape[0] == np.asarray(aperture_coords.point_data["Normals"]).shape[0]:
+            # print("hi from here", aperture_coords.cell_data)
+            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
+                    desired_E_axis.reshape(1, 3).astype(np.complex64),
+                    num_sources,
+                    axis=0,
                 )
     else:
         if project_vectors:
             conformal_E_vectors = EM.calculate_conformalVectors(
-                desired_E_axis, np.asarray(aperture_coords.point_data["Normals"]), antenna_axes
+                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]:
+            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
+                    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=[])
+        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_coords.points).astype(np.float32),
@@ -528,11 +559,11 @@ def calculate_scattering(
         )
         unified_weights = np.ones((unified_model.shape[0], 3), dtype=np.complex64)
         unified_weights[0:num_sources, :] = (
-            conformal_E_vectors #/ num_sources
-        )  # set total amplitude to 1 for the aperture
+            conformal_E_vectors  # / num_sources  # set total amplitude to 1 for the aperture
+        )
         unified_weights[num_sources : num_sources + num_sinks, :] = (
-            1 #/ num_sinks
-        )  # set total amplitude to 1 for the aperture
+            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
@@ -595,7 +626,9 @@ def calculate_scattering(
             if project_vectors:
                 conformal_E_vectors = EM.calculate_conformalVectors(
                     desired_E_axis[0, :].reshape(1, 3),
-                    np.asarray(aperture_coords.point_data["Normals"]).astype(np.float32),
+                    np.asarray(aperture_coords.point_data["Normals"]).astype(
+                        np.float32
+                    ),
                 )
             else:
                 conformal_E_vectors = np.repeat(
@@ -607,7 +640,9 @@ def calculate_scattering(
             if project_vectors:
                 conformal_E_vectors = EM.calculate_conformalVectors(
                     desired_E_axis[0, :].reshape(1, 3),
-                    np.asarray(aperture_coords.point_data["Normals"]).astype(np.float32),
+                    np.asarray(aperture_coords.point_data["Normals"]).astype(
+                        np.float32
+                    ),
                 )
             else:
                 if desired_E_axis.size == 3:
@@ -639,14 +674,14 @@ def calculate_scattering(
         )
         unified_weights = np.ones((unified_model.shape[0], 3), dtype=np.complex64)
         unified_weights[0:num_sources, :] = (
-            conformal_E_vectors #/ num_sources
-        )  # set total amplitude to 1 for the aperture
+            conformal_E_vectors  # / num_sources  # set total amplitude to 1 for the aperture
+        )
         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 :, :
-        ] = 1 # / len(np.asarray(scatter_points.points))  # set total amplitude to 1 for the aperture
+            1  # / num_sinks  # set total amplitude to 1 for the aperture
+        )
+        unified_weights[num_sources + num_sinks :, :] = (
+            1  # / 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
@@ -739,34 +774,47 @@ def calculate_scattering(
     if not elements:
         # create efiles for model
         if multiE:
-            Ex = np.zeros((desired_E_axis.shape[0],num_sinks), dtype=np.complex64)
-            Ey = np.zeros((desired_E_axis.shape[0],num_sinks), dtype=np.complex64)
-            Ez = np.zeros((desired_E_axis.shape[0],num_sinks), dtype=np.complex64)
+            Ex = np.zeros((desired_E_axis.shape[0], num_sinks), dtype=np.complex64)
+            Ey = np.zeros((desired_E_axis.shape[0], num_sinks), dtype=np.complex64)
+            Ez = np.zeros((desired_E_axis.shape[0], num_sinks), dtype=np.complex64)
             for e_inc in range(desired_E_axis.shape[0]):
                 conformal_E_vectors = EM.calculate_conformalVectors(
                     desired_E_axis[e_inc, :],
-                    np.asarray(aperture_coords.point_data["Normals"]).astype(np.float32),
+                    np.asarray(aperture_coords.point_data["Normals"]).astype(
+                        np.float32
+                    ),
                 )
-                unified_weights[0:num_sources, :] = conformal_E_vectors# / num_sources
+                unified_weights[0:num_sources, :] = conformal_E_vectors  # / num_sources
                 point_informationv2[:]["ex"] = unified_weights[:, 0]
                 point_informationv2[:]["ey"] = unified_weights[:, 1]
                 point_informationv2[:]["ez"] = unified_weights[:, 2]
                 scatter_map = EM.EMGPUFreqDomain(
-                    num_sources, num_sinks, full_index, point_informationv2, wavelength, alpha, beta
+                    num_sources,
+                    num_sinks,
+                    full_index,
+                    point_informationv2,
+                    wavelength,
+                    alpha,
+                    beta,
                 )
                 Ex[e_inc] = np.dot(np.ones((num_sources)), scatter_map[:, :, 0])
                 Ey[e_inc] = np.dot(np.ones((num_sources)), scatter_map[:, :, 1])
                 Ez[e_inc] = np.dot(np.ones((num_sources)), scatter_map[:, :, 2])
         else:
             scatter_map = EM.EMGPUFreqDomain(
-                num_sources, num_sinks, full_index, point_informationv2, wavelength, alpha, beta
+                num_sources,
+                num_sinks,
+                full_index,
+                point_informationv2,
+                wavelength,
+                alpha,
+                beta,
             )
 
             Ex = np.dot(np.ones((num_sources)), scatter_map[:, :, 0])
             Ey = np.dot(np.ones((num_sources)), scatter_map[:, :, 1])
             Ez = np.dot(np.ones((num_sources)), scatter_map[:, :, 2])
 
-
         # convert to etheta,ephi
 
     else:
@@ -778,25 +826,27 @@ def calculate_scattering(
             for e_inc in range(desired_E_axis.shape[1]):
                 conformal_E_vectors = EM.calculate_conformalVectors(
                     desired_E_axis[e_inc, :],
-                    np.asarray(aperture_coords.point_data["Normals"]).astype(np.float32),
+                    np.asarray(aperture_coords.point_data["Normals"]).astype(
+                        np.float32
+                    ),
                 )
                 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
-                    )
+                    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
-                    )
+                    unified_weights[element, :] = conformal_E_vectors[
+                        element, :
+                    ]  # / num_sources
                     scatter_map = EM.EMGPUFreqDomain(
                         num_sources,
                         num_sinks,
@@ -804,7 +854,7 @@ def calculate_scattering(
                         point_informationv2,
                         wavelength,
                         alpha,
-                        beta
+                        beta,
                     )
                     Ex[element, :, e_inc] = np.dot(
                         np.ones((num_sources)), scatter_map[:, :, 0]
@@ -820,14 +870,16 @@ def calculate_scattering(
             Ey = np.zeros((num_sources, num_sinks), dtype=np.complex64)
             Ez = np.zeros((num_sources, num_sinks), dtype=np.complex64)
             scatter_map = EM.EMGPUFreqDomain(
-                num_sources, num_sinks, full_index, point_informationv2, wavelength, alpha, beta
+                num_sources,
+                num_sinks,
+                full_index,
+                point_informationv2,
+                wavelength,
+                alpha,
+                beta,
             )
             Ex = scatter_map[:, :, 0]
             Ey = scatter_map[:, :, 1]
             Ez = scatter_map[:, :, 2]
-                
 
     return Ex, Ey, Ez
-
-
-
diff --git a/lyceanem/models/time_domain.py b/lyceanem/models/time_domain.py
index 6911e4d..3627217 100644
--- a/lyceanem/models/time_domain.py
+++ b/lyceanem/models/time_domain.py
@@ -27,8 +27,8 @@ def calculate_scattering(
     mesh_resolution=0.5,
     antenna_axes=np.eye(3),
     project_vectors=True,
-        alpha=0.0,
-        beta=(np.pi*2)/1.0
+    alpha=0.0,
+    beta=(np.pi * 2) / 1.0,
 ):
     """
     Based upon the parameters given, calculate the time domain scattering for the apertures and sinks.
@@ -76,24 +76,33 @@ def calculate_scattering(
 
     num_sources = len(np.asarray(aperture_coords.points))
     num_sinks = len(np.asarray(sink_coords.points))
-    environment_triangles=GF.mesh_conversion(antenna_solid)
+    environment_triangles = GF.mesh_conversion(antenna_solid)
 
     if not multiE:
         if project_vectors:
             conformal_E_vectors = EM.calculate_conformalVectors(
-                desired_E_axis, np.asarray(aperture_coords.point_data["Normals"]), antenna_axes
+                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]:
+            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
+                    desired_E_axis.reshape(1, 3).astype(np.complex64),
+                    num_sources,
+                    axis=0,
                 )
     else:
         if project_vectors:
             conformal_E_vectors = EM.calculate_conformalVectors(
-                desired_E_axis, np.asarray(aperture_coords.point_data["Normals"]), antenna_axes
+                desired_E_axis,
+                np.asarray(aperture_coords.point_data["Normals"]),
+                antenna_axes,
             )
         else:
             if desired_E_axis.size == 3:
@@ -105,7 +114,9 @@ def calculate_scattering(
 
     if scattering == 0:
         # only use the aperture point cloud, no scattering required.
-        scatter_points = meshio.Mesh(points = np.empty((0,3), dtype=np.float32), cells=[])
+        scatter_points = meshio.Mesh(
+            points=np.empty((0, 3), dtype=np.float32), cells=[]
+        )
 
         unified_model = np.append(
             np.asarray(aperture_coords.points).astype(np.float32),
@@ -119,11 +130,11 @@ def calculate_scattering(
         )
         unified_weights = np.ones((unified_model.shape[0], 3), dtype=np.complex64)
         unified_weights[0:num_sources, :] = (
-            conformal_E_vectors #/ num_sources
-        )  # set total amplitude to 1 for the aperture
+            conformal_E_vectors  # / num_sources  # set total amplitude to 1 for the aperture
+        )
         unified_weights[num_sources : num_sources + num_sinks, :] = (
-            1 #/ num_sinks
-        )  # set total amplitude to 1 for the aperture
+            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
@@ -202,14 +213,14 @@ def calculate_scattering(
         )
         unified_weights = np.ones((unified_model.shape[0], 3), dtype=np.complex64)
         unified_weights[0:num_sources, :] = (
-            conformal_E_vectors #/ num_sources
-        )  # set total amplitude to 1 for the aperture
+            conformal_E_vectors  # / num_sources  # set total amplitude to 1 for the aperture
+        )
         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 :, :] = 1 #/ len(
-          #  np.asarray(scatter_points.points)
-        #)  # set total amplitude to 1 for the aperture
+            1  # / num_sinks  # set total amplitude to 1 for the aperture
+        )
+        unified_weights[num_sources + num_sinks :, :] = 1  # / 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
@@ -307,9 +318,11 @@ def calculate_scattering(
             for e_inc in range(desired_E_axis.shape[0]):
                 conformal_E_vectors = EM.calculate_conformalVectors(
                     desired_E_axis[e_inc, :],
-                    np.asarray(aperture_coords.point_data["Normals"]).astype(np.float32),
+                    np.asarray(aperture_coords.point_data["Normals"]).astype(
+                        np.float32
+                    ),
                 )
-                unified_weights[0:num_sources, :] = conformal_E_vectors #/ num_sources
+                unified_weights[0:num_sources, :] = conformal_E_vectors  # / num_sources
                 point_informationv2[:]["ex"] = unified_weights[:, 0]
                 point_informationv2[:]["ey"] = unified_weights[:, 1]
                 point_informationv2[:]["ez"] = unified_weights[:, 2]
@@ -325,7 +338,7 @@ def calculate_scattering(
                     sampling_freq,
                     num_samples,
                     alpha,
-                    beta
+                    beta,
                 )
                 wake_index = np.digitize(WakeTimes, time_index)
                 Ex[e_inc, wake_index:] = np.dot(
@@ -366,7 +379,7 @@ def calculate_scattering(
                 sampling_freq,
                 num_samples,
                 alpha,
-                beta
+                beta,
             )
             wake_index = np.digitize(WakeTimes, time_index)
             Ex[wake_index:] = np.dot(
@@ -405,7 +418,10 @@ def calculate_scattering(
             for e_inc in range(desired_E_axis.shape[1]):
                 conformal_E_vectors = EM.calculate_conformalVectors(
                     desired_E_axis[e_inc, :],
-                    np.asarray(aperture_coords.point_data["Normals"]).astype(np.float32),antenna_axes
+                    np.asarray(aperture_coords.point_data["Normals"]).astype(
+                        np.float32
+                    ),
+                    antenna_axes,
                 )
                 for element in range(num_sources):
                     point_informationv2[0:num_sources]["ex"] = 0.0
@@ -436,7 +452,7 @@ def calculate_scattering(
                         sampling_freq,
                         num_samples,
                         alpha,
-                        beta
+                        beta,
                     )
                     wake_index = np.digitize(WakeTimes, time_index)
                     Ex[element, wake_index:] = np.dot(
@@ -494,7 +510,7 @@ def calculate_scattering(
                     sampling_freq,
                     num_samples,
                     alpha,
-                    beta
+                    beta,
                 )
                 Ex[element, :, :] = np.dot(
                     np.ones((num_sources)),
diff --git a/lyceanem/raycasting/rayfunctions.py b/lyceanem/raycasting/rayfunctions.py
index c401676..5986ef5 100644
--- a/lyceanem/raycasting/rayfunctions.py
+++ b/lyceanem/raycasting/rayfunctions.py
@@ -293,6 +293,7 @@ def integratedRaycaster(ray_index, scattering_points, environment_local):
 
     return final_index
 
+
 # @njit
 def patterntocloud(pattern_data, shell_coords, maxarea):
     # takes the pattern_data and shell_coordinates, and creates an open3d point cloud based upon the data.
@@ -305,8 +306,10 @@ def patterntocloud(pattern_data, shell_coords, maxarea):
     point_cloud.point_data["red"] = np_colors[:, 0]
     point_cloud.point_data["green"] = np_colors[:, 1]
     point_cloud.point_data["blue"] = np_colors[:, 2]
-    
+
     return point_cloud
+
+
 def visiblespace(
     source_coords,
     source_normals,
@@ -481,7 +484,6 @@ def patterncreator(az_range, elev_range, source_coords, pointingindex, hit_index
     return visible_patterns
 
 
-
 def quickpatterncreator(
     az_range, elev_range, source_coords, angles, vertex_area, hit_index
 ):
@@ -510,9 +512,6 @@ def quickpatterncreator(
     return visible_patterns
 
 
-
-
-
 def azeltocart(az_data, el_data, radius):
     # convert from az,el and radius data to xyz
     x_data = radius * np.cos(np.deg2rad(el_data)) * np.cos(np.deg2rad(az_data))
@@ -528,11 +527,11 @@ def convertTriangles(triangle_object):
     if triangle_object == None:
         triangles = np.empty(0, dtype=base_types.triangle_t)
     else:
-        #assert triangle_object.cells[1].type == "triangle", "Not a triangle mesh"
+        # assert triangle_object.cells[1].type == "triangle", "Not a triangle mesh"
         for idx in range(len(triangle_object.cells)):
-            if (triangle_object.cells[idx].type=="triangle"):
-                triangle_cell_index=idx
-                
+            if triangle_object.cells[idx].type == "triangle":
+                triangle_cell_index = idx
+
         vertices = np.asarray(triangle_object.points)
         tri_index = np.asarray(triangle_object.cells[triangle_cell_index].data)
         triangles = np.empty(len(tri_index), dtype=base_types.triangle_t)
@@ -550,8 +549,6 @@ def convertTriangles(triangle_object):
     return triangles
 
 
-
-
 def points2pointcloud(xyz):
     """
     turns numpy array of xyz data into a meshio format point cloud
@@ -1965,9 +1962,9 @@ def launchRaycaster1Dv2(
         kernel1Dv2[grids, threads](d_chunk_payload, d_environment)
         # cuda.profile_stop()
         # distmap[source_chunks[n]:source_chunks[n+1],target_chunks[m]:target_chunks[m+1]] = d_distmap_chunked.copy_to_host()
-        first_ray_payload[
-            ray_chunks[n] : ray_chunks[n + 1]
-        ] = d_chunk_payload.copy_to_host()
+        first_ray_payload[ray_chunks[n] : ray_chunks[n + 1]] = (
+            d_chunk_payload.copy_to_host()
+        )
         # first_ray_payload[ray_chunks[n]:ray_chunks[n+1]]['intersect']=d_chunk_payload['intersect'].copy_to_host()
         # first_ray_payload[ray_chunks[n]:ray_chunks[n+1]]['dist']=d_chunk_payload['dist'].copy_to_host()
 
@@ -2121,9 +2118,9 @@ def chunkingRaycaster1Dv2(
         # distmap[source_chunks[n]:source_chunks[n+1],target_chunks[m]:target_chunks[m+1]] = d_distmap_chunked.copy_to_host()
         # second_ray_payload[ray_chunks[n]:ray_chunks[n+1]]=d_chunk_payload.copy_to_host()
         # first_ray_payload[ray_chunks[n]:ray_chunks[n+1]]=d_chunk_payload.copy_to_host()
-        second_ray_payload[
-            ray_chunks[n] : ray_chunks[n + 1]
-        ] = d_chunk_payload.copy_to_host()
+        second_ray_payload[ray_chunks[n] : ray_chunks[n + 1]] = (
+            d_chunk_payload.copy_to_host()
+        )
 
     # cuda.close()
     kernel_dt = timer() - raystart
@@ -2399,12 +2396,12 @@ def chunkingRaycaster1Dv3(
     #         # ctx = cuda.current_context()
     #         # ctx.reset()
     # else:
-        
-        # deallocate memory on gpu
-        # ctx = cuda.current_context()
-        # deallocs = ctx.deallocations
-        # deallocs.clear()
-        # print("Second Stage: Prep {:3.1f} s, Raycasting  {:3.1f} s, Path Processing {:3.1f} s".format(prep_dt,kernel_dt,mem_dt) )
+
+    # deallocate memory on gpu
+    # ctx = cuda.current_context()
+    # deallocs = ctx.deallocations
+    # deallocs.clear()
+    # print("Second Stage: Prep {:3.1f} s, Raycasting  {:3.1f} s, Path Processing {:3.1f} s".format(prep_dt,kernel_dt,mem_dt) )
     return filtered_index2, final_index2, RAYS_CAST
 
 
@@ -2474,42 +2471,42 @@ def ray_charge_core_final_vector(
         directions, norm_length = math_functions.calc_dv_norm(
             local_sources[:, -3:], targets, directions, norm_length
         )
-        temp_ray_payload[source_chunking[n] : source_chunking[n + 1]][
-            "ox"
-        ] = local_sources[:, -3]
-        temp_ray_payload[source_chunking[n] : source_chunking[n + 1]][
-            "oy"
-        ] = local_sources[:, -2]
-        temp_ray_payload[source_chunking[n] : source_chunking[n + 1]][
-            "oz"
-        ] = local_sources[:, -1]
-        temp_ray_payload[source_chunking[n] : source_chunking[n + 1]][
-            "dx"
-        ] = directions[:, 0]
-        temp_ray_payload[source_chunking[n] : source_chunking[n + 1]][
-            "dy"
-        ] = directions[:, 1]
-        temp_ray_payload[source_chunking[n] : source_chunking[n + 1]][
-            "dz"
-        ] = directions[:, 2]
+        temp_ray_payload[source_chunking[n] : source_chunking[n + 1]]["ox"] = (
+            local_sources[:, -3]
+        )
+        temp_ray_payload[source_chunking[n] : source_chunking[n + 1]]["oy"] = (
+            local_sources[:, -2]
+        )
+        temp_ray_payload[source_chunking[n] : source_chunking[n + 1]]["oz"] = (
+            local_sources[:, -1]
+        )
+        temp_ray_payload[source_chunking[n] : source_chunking[n + 1]]["dx"] = (
+            directions[:, 0]
+        )
+        temp_ray_payload[source_chunking[n] : source_chunking[n + 1]]["dy"] = (
+            directions[:, 1]
+        )
+        temp_ray_payload[source_chunking[n] : source_chunking[n + 1]]["dz"] = (
+            directions[:, 2]
+        )
         # temp_ray_payload[source_chunking[n]:source_chunking[n+1]]['tx']=targets[:,0]
         # temp_ray_payload[source_chunking[n]:source_chunking[n+1]]['ty']=targets[:,1]
         # temp_ray_payload[source_chunking[n]:source_chunking[n+1]]['tz']=targets[:,2]
-        temp_ray_payload[source_chunking[n] : source_chunking[n + 1]][
-            "dist"
-        ] = norm_length[:, 0]
+        temp_ray_payload[source_chunking[n] : source_chunking[n + 1]]["dist"] = (
+            norm_length[:, 0]
+        )
         temp_ray_payload[source_chunking[n] : source_chunking[n + 1]][
             "intersect"
         ] = False
         origins[source_chunking[n] : source_chunking[n + 1], -origin_width:] = sources[
             n, :
         ]
-        source_sink_index[
-            source_chunking[n] : source_chunking[n + 1], 0
-        ] = point_indexing[n, 0]
-        source_sink_index[
-            source_chunking[n] : source_chunking[n + 1], 1
-        ] = target_indexing.ravel()
+        source_sink_index[source_chunking[n] : source_chunking[n + 1], 0] = (
+            point_indexing[n, 0]
+        )
+        source_sink_index[source_chunking[n] : source_chunking[n + 1], 1] = (
+            target_indexing.ravel()
+        )
 
     return temp_ray_payload, origins, source_sink_index
 
@@ -2753,8 +2750,7 @@ def workchunkingv2(
     free_mem, total_mem = cuda.current_context().get_memory_info()
     max_mem = np.ceil(free_mem * 0.8).astype(np.int64)
     ray_limit = (
-        np.floor(np.floor((max_mem - environment.nbytes) / base_types.ray_t.size) )
-        
+        np.floor(np.floor((max_mem - environment.nbytes) / base_types.ray_t.size))
     ).astype(np.int64)
     # establish index boundaries
     source_index = np.arange(1, sources.shape[0] + 1).reshape(
diff --git a/lyceanem/tests/reflectordata.py b/lyceanem/tests/reflectordata.py
index fd269de..af46d8a 100644
--- a/lyceanem/tests/reflectordata.py
+++ b/lyceanem/tests/reflectordata.py
@@ -47,25 +47,29 @@ def exampleUAV(frequency):
     array = GF.mesh_rotate(array, rotation_vector1)
     array = GF.mesh_rotate(array, rotation_vector2)
 
-
-    body = GF.translate_mesh(body, np.array([0.25, 0, 0])+np.array([-0.18, 0, 0.0125]))
-    array = GF.translate_mesh(array, np.array([0.25, 0, 0])+np.array([-0.18, 0, 0.0125])
+    body = GF.translate_mesh(
+        body, np.array([0.25, 0, 0]) + np.array([-0.18, 0, 0.0125])
+    )
+    array = GF.translate_mesh(
+        array, np.array([0.25, 0, 0]) + np.array([-0.18, 0, 0.0125])
     )
 
     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)
+        array = np.append(
+            np.ones((array.shape[0], 1), dtype=np.int32) * 3, array, axis=1
+        )
         return array
+
     pyvista_array = pv.PolyData(array.points, structure_cells(array.cells[0].data))
     pyvista_body = pv.PolyData(body.points, structure_cells(body.cells[0].data))
     pyvista_array.compute_normals(inplace=True)
     pyvista_body.compute_normals(inplace=True)
 
-    array.point_data['Normals'] = pyvista_array.point_normals
-    body.point_data['Normals'] = pyvista_body.point_normals
-    array.cell_data['Normals'] = pyvista_array.cell_normals
-    body.cell_data['Normals'] = pyvista_body.cell_normals
-
+    array.point_data["Normals"] = pyvista_array.point_normals
+    body.point_data["Normals"] = pyvista_body.point_normals
+    array.cell_data["Normals"] = pyvista_array.cell_normals
+    body.cell_data["Normals"] = pyvista_body.cell_normals
 
     wavelength = 3e8 / frequency
     mesh_sep = wavelength * 0.5
@@ -195,8 +199,7 @@ def structure_cells(array):
         axis=0,
     ) + np.array([0.1025, 0, -0.025])
     source_pcd = RF.points2pointcloud(total_array)
-    source_pcd.point_data['Normals'] = total_array_normals
+    source_pcd.point_data["Normals"] = total_array_normals
     source_pcd = GF.translate_mesh(source_pcd, np.array([-0.18, 0, 0.0125]))
 
-
     return body, array, source_pcd
diff --git a/lyceanem/tests/test_base_classes.py b/lyceanem/tests/test_base_classes.py
index d0ad84d..a4064e1 100644
--- a/lyceanem/tests/test_base_classes.py
+++ b/lyceanem/tests/test_base_classes.py
@@ -1,5 +1,5 @@
 import numpy as np
-import meshio 
+import meshio
 import pytest
 from numpy.testing import assert_allclose
 from scipy.spatial.transform import Rotation as R
@@ -9,39 +9,49 @@
 
 def cube():
     # a cube centered at the origin, with side lengths of 1m (default)
-    
-
-    cube_points = np.array([[-0.5, -0.5, -0.5],
-                            [ 0.5, -0.5, -0.5],
-                            [-0.5 ,-0.5,  0.5],
-                            [ 0.5 ,-0.5 , 0.5],
-                            [-0.5,  0.5, -0.5],
-                            [ 0.5 , 0.5, -0.5],
-                            [-0.5,  0.5,  0.5],
-                            [ 0.5 , 0.5 , 0.5]])
-    cube_cells = [[4, 7, 5],
-                [4, 6, 7],
-                [0, 2, 4],
-                [2, 6, 4],
-                [0, 1, 2],
-                [1, 3, 2],
-                [1, 5, 7],
-                [1, 7, 3],
-                [2, 3, 7],
-                [2, 7, 6],
-                [0 ,4, 1],
-                [1, 4, 5]]
+
+    cube_points = np.array(
+        [
+            [-0.5, -0.5, -0.5],
+            [0.5, -0.5, -0.5],
+            [-0.5, -0.5, 0.5],
+            [0.5, -0.5, 0.5],
+            [-0.5, 0.5, -0.5],
+            [0.5, 0.5, -0.5],
+            [-0.5, 0.5, 0.5],
+            [0.5, 0.5, 0.5],
+        ]
+    )
+    cube_cells = [
+        [4, 7, 5],
+        [4, 6, 7],
+        [0, 2, 4],
+        [2, 6, 4],
+        [0, 1, 2],
+        [1, 3, 2],
+        [1, 5, 7],
+        [1, 7, 3],
+        [2, 3, 7],
+        [2, 7, 6],
+        [0, 4, 1],
+        [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)
+
+    cube = compute_normals(cube)
     return cube
 
 
 def point():
     # a single point on the +x center of the cube with consistent normal vector
 
-    pc_mesh = meshio.Mesh(points=np.array([[0.5, 0, 0]]), cells=[], point_data={"Normals": np.array([[1.0, 0, 0]])})
+    pc_mesh = meshio.Mesh(
+        points=np.array([[0.5, 0, 0]]),
+        cells=[],
+        point_data={"Normals": np.array([[1.0, 0, 0]])},
+    )
     return pc_mesh
 
 
@@ -57,7 +67,6 @@ def standard_antenna():
     return antenna()
 
 
-
 def test_excitation_function_x_u(standard_antenna):
     # test that an unrotated antenna with u (horizontal-y) polarisation gives horizontal-y polarisation
     desired_E_vector = np.array([1, 0, 0], dtype=np.complex64)
@@ -140,4 +149,4 @@ def test_excitation_function_x_mz_cp(standard_antenna):
         standard_antenna.excitation_function(desired_e_vector=desired_E_vector),
         final_vector,
         atol=1e-12,
-    )
\ No newline at end of file
+    )
diff --git a/lyceanem/tests/test_empropagation.py b/lyceanem/tests/test_empropagation.py
index 8c6de33..9fd4c6c 100644
--- a/lyceanem/tests/test_empropagation.py
+++ b/lyceanem/tests/test_empropagation.py
@@ -6,55 +6,90 @@
 
 
 def test_vector_mapping_x_u():
-    #initially test in global coordinate set
-    rotation_matrix=R.from_euler('z',0,degrees=True).as_matrix()
-    normal_vector=np.array([1.0, 0, 0])
-    desired_E_vector=np.array([1.0, 0, 0],dtype=np.complex64) #uvn axes with n being the normal vector
-    final_E_vector=np.array([0.0, 1.0, 0],dtype=np.complex64)
-    assert_allclose(vector_mapping(desired_E_vector, normal_vector, rotation_matrix), final_E_vector,
-                    atol=1e-12)
+    # initially test in global coordinate set
+    rotation_matrix = R.from_euler("z", 0, degrees=True).as_matrix()
+    normal_vector = np.array([1.0, 0, 0])
+    desired_E_vector = np.array(
+        [1.0, 0, 0], dtype=np.complex64
+    )  # uvn axes with n being the normal vector
+    final_E_vector = np.array([0.0, 1.0, 0], dtype=np.complex64)
+    assert_allclose(
+        vector_mapping(desired_E_vector, normal_vector, rotation_matrix),
+        final_E_vector,
+        atol=1e-12,
+    )
+
 
 def test_vector_mapping_x_v():
-    #initially test in global coordinate set
-    rotation_matrix=R.from_euler('z',0,degrees=True).as_matrix()
-    normal_vector=np.array([1.0, 0, 0])
-    desired_E_vector=np.array([0.0, 1.0, 0],dtype=np.complex64) #uvn axes with n being the normal vector
-    final_E_vector=np.array([0.0, 0.0, 1.0],dtype=np.complex64)
-    assert_allclose(vector_mapping(desired_E_vector, normal_vector, rotation_matrix), final_E_vector,
-                    atol=1e-12)
+    # initially test in global coordinate set
+    rotation_matrix = R.from_euler("z", 0, degrees=True).as_matrix()
+    normal_vector = np.array([1.0, 0, 0])
+    desired_E_vector = np.array(
+        [0.0, 1.0, 0], dtype=np.complex64
+    )  # uvn axes with n being the normal vector
+    final_E_vector = np.array([0.0, 0.0, 1.0], dtype=np.complex64)
+    assert_allclose(
+        vector_mapping(desired_E_vector, normal_vector, rotation_matrix),
+        final_E_vector,
+        atol=1e-12,
+    )
+
 
 def test_vector_mapping_x_n():
-    #initially test in global coordinate set
-    rotation_matrix=R.from_euler('z',0,degrees=True).as_matrix()
-    normal_vector=np.array([1.0, 0, 0])
-    desired_E_vector=np.array([0, 0, 1.0],dtype=np.complex64) #uvn axes with n being the normal vector
-    final_E_vector=np.array([1.0, 0, 0],dtype=np.complex64)
-    assert_allclose(vector_mapping(desired_E_vector, normal_vector, rotation_matrix), final_E_vector,
-                    atol=1e-12)
+    # initially test in global coordinate set
+    rotation_matrix = R.from_euler("z", 0, degrees=True).as_matrix()
+    normal_vector = np.array([1.0, 0, 0])
+    desired_E_vector = np.array(
+        [0, 0, 1.0], dtype=np.complex64
+    )  # uvn axes with n being the normal vector
+    final_E_vector = np.array([1.0, 0, 0], dtype=np.complex64)
+    assert_allclose(
+        vector_mapping(desired_E_vector, normal_vector, rotation_matrix),
+        final_E_vector,
+        atol=1e-12,
+    )
+
 
 def test_vector_mapping_x_y_u():
-    #initially test in global coordinate set
-    rotation_matrix=R.from_euler('z',-90,degrees=True).as_matrix()
-    normal_vector=np.array([1.0, 0, 0])
-    desired_E_vector=np.array([1.0, 0, 0],dtype=np.complex64) #uvn axes with n being the normal vector
-    final_E_vector=np.array([-1.0, 0.0, 0],dtype=np.complex64)
-    assert_allclose(vector_mapping(desired_E_vector, normal_vector, rotation_matrix), final_E_vector,
-                    atol=1e-12)
+    # initially test in global coordinate set
+    rotation_matrix = R.from_euler("z", -90, degrees=True).as_matrix()
+    normal_vector = np.array([1.0, 0, 0])
+    desired_E_vector = np.array(
+        [1.0, 0, 0], dtype=np.complex64
+    )  # uvn axes with n being the normal vector
+    final_E_vector = np.array([-1.0, 0.0, 0], dtype=np.complex64)
+    assert_allclose(
+        vector_mapping(desired_E_vector, normal_vector, rotation_matrix),
+        final_E_vector,
+        atol=1e-12,
+    )
+
 
 def test_vector_mapping_x_y_v():
-    #initially test in global coordinate set
-    rotation_matrix=R.from_euler('z',-90,degrees=True).as_matrix()
-    normal_vector=np.array([1.0, 0, 0])
-    desired_E_vector=np.array([0.0, 1.0, 0],dtype=np.complex64) #uvn axes with n being the normal vector
-    final_E_vector=np.array([0.0, 0.0, 1.0],dtype=np.complex64)
-    assert_allclose(vector_mapping(desired_E_vector, normal_vector, rotation_matrix), final_E_vector,
-                    atol=1e-12)
+    # initially test in global coordinate set
+    rotation_matrix = R.from_euler("z", -90, degrees=True).as_matrix()
+    normal_vector = np.array([1.0, 0, 0])
+    desired_E_vector = np.array(
+        [0.0, 1.0, 0], dtype=np.complex64
+    )  # uvn axes with n being the normal vector
+    final_E_vector = np.array([0.0, 0.0, 1.0], dtype=np.complex64)
+    assert_allclose(
+        vector_mapping(desired_E_vector, normal_vector, rotation_matrix),
+        final_E_vector,
+        atol=1e-12,
+    )
+
 
 def test_vector_mapping_x_y_n():
-    #initially test in global coordinate set
-    rotation_matrix=R.from_euler('z',-90,degrees=True).as_matrix()
-    normal_vector=np.array([1.0, 0, 0])
-    desired_E_vector=np.array([0, 0, 1.0],dtype=np.complex64) #uvn axes with n being the normal vector
-    final_E_vector=np.array([0.0, 1.0, 0],dtype=np.complex64)
-    assert_allclose(vector_mapping(desired_E_vector,normal_vector,rotation_matrix), final_E_vector,
-                    atol=1e-12)
\ No newline at end of file
+    # initially test in global coordinate set
+    rotation_matrix = R.from_euler("z", -90, degrees=True).as_matrix()
+    normal_vector = np.array([1.0, 0, 0])
+    desired_E_vector = np.array(
+        [0, 0, 1.0], dtype=np.complex64
+    )  # uvn axes with n being the normal vector
+    final_E_vector = np.array([0.0, 1.0, 0], dtype=np.complex64)
+    assert_allclose(
+        vector_mapping(desired_E_vector, normal_vector, rotation_matrix),
+        final_E_vector,
+        atol=1e-12,
+    )
diff --git a/lyceanem/tests/test_geometryfunctions.py b/lyceanem/tests/test_geometryfunctions.py
index 2cd99ce..95249d1 100644
--- a/lyceanem/tests/test_geometryfunctions.py
+++ b/lyceanem/tests/test_geometryfunctions.py
@@ -5,14 +5,13 @@
 from scipy.spatial.transform import Rotation as R
 import meshio
 from importlib_resources import files
-from  ..geometry import geometryfunctions as GF
+from ..geometry import geometryfunctions as GF
 import lyceanem.tests.data
 
 
-
 @pytest.fixture
 def standard_cube():
-    cube_location= files(lyceanem.tests.data).joinpath("cube.ply")
+    cube_location = files(lyceanem.tests.data).joinpath("cube.ply")
     cube = meshio.read(cube_location)
     return cube
 
@@ -20,67 +19,113 @@ def standard_cube():
 def are_meshes_equal(mesh1, mesh2, rtol=1e-5):
     # Check vertices
     meshes_are_equal = True
-    meshes_are_equal =  np.allclose(mesh1.points, mesh2.points, rtol=rtol)
-
+    meshes_are_equal = np.allclose(mesh1.points, mesh2.points, rtol=rtol)
 
     # Check cells
     meshes_are_equal = len(mesh1.cells) == len(mesh2.cells)
     for i in range(len(mesh1.cells)):
-        meshes_are_equal  *= np.allclose(mesh1.cells[i].data, mesh2.cells[i].data, rtol=rtol)
+        meshes_are_equal *= np.allclose(
+            mesh1.cells[i].data, mesh2.cells[i].data, rtol=rtol
+        )
     ## normals
     print("checking 'normffffals3")
 
-    if 'nx' in mesh1.point_data and 'ny' in mesh1.point_data and 'nz' in mesh1.point_data:
-        if 'nx' in mesh2.point_data and 'ny' in mesh2.point_data and 'nz' in mesh2.point_data:
-            meshes_are_equal  *= np.allclose(mesh1.point_data['nx'], mesh2.point_data['nx'], rtol=rtol)
-            meshes_are_equal  *= np.allclose(mesh1.point_data['ny'], mesh2.point_data['ny'], rtol=rtol)
-            meshes_are_equal  *= np.allclose(mesh1.point_data['nz'], mesh2.point_data['nz'], rtol=rtol)
-        elif 'Normals' in mesh2.point_data:
-            meshes_are_equal  *= np.allclose(mesh1.point_data['nx'], mesh2.point_data['Normals'][:,0], rtol=rtol)
-            meshes_are_equal  *= np.allclose(mesh1.point_data['ny'], mesh2.point_data['Normals'][:,1], rtol=rtol)
-            meshes_are_equal  *= np.allclose(mesh1.point_data['nz'], mesh2.point_data['Normals'][:,2], rtol=rtol)
+    if (
+        "nx" in mesh1.point_data
+        and "ny" in mesh1.point_data
+        and "nz" in mesh1.point_data
+    ):
+        if (
+            "nx" in mesh2.point_data
+            and "ny" in mesh2.point_data
+            and "nz" in mesh2.point_data
+        ):
+            meshes_are_equal *= np.allclose(
+                mesh1.point_data["nx"], mesh2.point_data["nx"], rtol=rtol
+            )
+            meshes_are_equal *= np.allclose(
+                mesh1.point_data["ny"], mesh2.point_data["ny"], rtol=rtol
+            )
+            meshes_are_equal *= np.allclose(
+                mesh1.point_data["nz"], mesh2.point_data["nz"], rtol=rtol
+            )
+        elif "Normals" in mesh2.point_data:
+            meshes_are_equal *= np.allclose(
+                mesh1.point_data["nx"], mesh2.point_data["Normals"][:, 0], rtol=rtol
+            )
+            meshes_are_equal *= np.allclose(
+                mesh1.point_data["ny"], mesh2.point_data["Normals"][:, 1], rtol=rtol
+            )
+            meshes_are_equal *= np.allclose(
+                mesh1.point_data["nz"], mesh2.point_data["Normals"][:, 2], rtol=rtol
+            )
         else:
-            meshes_are_equal  *= False
-    elif 'Normals' in mesh1.point_data:
-        if 'nx' and 'ny' and 'nz' in mesh2.point_data:
-            meshes_are_equal  *= np.allclose(mesh1.point_data['Normals'][:,0], mesh2.point_data['nx'], rtol=rtol)
-            meshes_are_equal  *= np.allclose(mesh1.point_data['Normals'][:,1], mesh2.point_data['ny'], rtol=rtol)
-            meshes_are_equal  *= np.allclose(mesh1.point_data['Normals'][:,2], mesh2.point_data['nz'], rtol=rtol)
+            meshes_are_equal *= False
+    elif "Normals" in mesh1.point_data:
+        if "nx" and "ny" and "nz" in mesh2.point_data:
+            meshes_are_equal *= np.allclose(
+                mesh1.point_data["Normals"][:, 0], mesh2.point_data["nx"], rtol=rtol
+            )
+            meshes_are_equal *= np.allclose(
+                mesh1.point_data["Normals"][:, 1], mesh2.point_data["ny"], rtol=rtol
+            )
+            meshes_are_equal *= np.allclose(
+                mesh1.point_data["Normals"][:, 2], mesh2.point_data["nz"], rtol=rtol
+            )
             print("checking 'normals3")
-        elif 'Normals' in mesh2.point_data:
-            meshes_are_equal  *= np.allclose(mesh1.point_data['Normals'], mesh2.point_data['Normals'], rtol=rtol)
+        elif "Normals" in mesh2.point_data:
+            meshes_are_equal *= np.allclose(
+                mesh1.point_data["Normals"], mesh2.point_data["Normals"], rtol=rtol
+            )
         else:
-            meshes_are_equal  *= False
+            meshes_are_equal *= False
     else:
-        meshes_are_equal  *= True
+        meshes_are_equal *= True
     ##normals cell data
     for i in range(len(mesh1.cells)):
-        if 'nx' and 'ny' and 'nz' in mesh1.cell_data:
-            if 'nx' and 'ny' and 'nz' in mesh2.cell_data:
-                meshes_are_equal  *= np.allclose(mesh1.cell_data['nx'], mesh2.cell_data['nx'], rtol=rtol)
-                meshes_are_equal  *= np.allclose(mesh1.cell_data['ny'], mesh2.cell_data['ny'], rtol=rtol)
-                meshes_are_equal  *= np.allclose(mesh1.cell_data['nz'], mesh2.cell_data['nz'], rtol=rtol)
-            elif 'Normals' in mesh2.cell_data:
-                meshes_are_equal  *= np.allclose(mesh1.cell_data['nx'], mesh2.cell_data['Normals'][:,0], rtol=rtol)
-                meshes_are_equal  *= np.allclose(mesh1.cell_data['ny'], mesh2.cell_data['Normals'][:,1], rtol=rtol)
-                meshes_are_equal  *= np.allclose(mesh1.cell_data['nz'], mesh2.cell_data['Normals'][:,2], rtol=rtol)
+        if "nx" and "ny" and "nz" in mesh1.cell_data:
+            if "nx" and "ny" and "nz" in mesh2.cell_data:
+                meshes_are_equal *= np.allclose(
+                    mesh1.cell_data["nx"], mesh2.cell_data["nx"], rtol=rtol
+                )
+                meshes_are_equal *= np.allclose(
+                    mesh1.cell_data["ny"], mesh2.cell_data["ny"], rtol=rtol
+                )
+                meshes_are_equal *= np.allclose(
+                    mesh1.cell_data["nz"], mesh2.cell_data["nz"], rtol=rtol
+                )
+            elif "Normals" in mesh2.cell_data:
+                meshes_are_equal *= np.allclose(
+                    mesh1.cell_data["nx"], mesh2.cell_data["Normals"][:, 0], rtol=rtol
+                )
+                meshes_are_equal *= np.allclose(
+                    mesh1.cell_data["ny"], mesh2.cell_data["Normals"][:, 1], rtol=rtol
+                )
+                meshes_are_equal *= np.allclose(
+                    mesh1.cell_data["nz"], mesh2.cell_data["Normals"][:, 2], rtol=rtol
+                )
             else:
-                meshes_are_equal  *= False
-        elif 'Normals' in mesh1.cell_data:
-            if 'nx' and 'ny' and 'nz' in mesh2.cell_data:
-                meshes_are_equal  *= np.allclose(mesh1.cell_data['Normals'][:,0], mesh2.cell_data['nx'], rtol=rtol)
-                meshes_are_equal  *= np.allclose(mesh1.cell_data['Normals'][:,1], mesh2.cell_data['ny'], rtol=rtol)
-                meshes_are_equal  *= np.allclose(mesh1.cell_data['Normals'][:,2], mesh2.cell_data['nz'], rtol=rtol)
-            elif 'Normals' in mesh2.cell_data:
-                meshes_are_equal  *= np.allclose(mesh1.cell_data['Normals'], mesh2.cell_data['Normals'], rtol=rtol)
+                meshes_are_equal *= False
+        elif "Normals" in mesh1.cell_data:
+            if "nx" and "ny" and "nz" in mesh2.cell_data:
+                meshes_are_equal *= np.allclose(
+                    mesh1.cell_data["Normals"][:, 0], mesh2.cell_data["nx"], rtol=rtol
+                )
+                meshes_are_equal *= np.allclose(
+                    mesh1.cell_data["Normals"][:, 1], mesh2.cell_data["ny"], rtol=rtol
+                )
+                meshes_are_equal *= np.allclose(
+                    mesh1.cell_data["Normals"][:, 2], mesh2.cell_data["nz"], rtol=rtol
+                )
+            elif "Normals" in mesh2.cell_data:
+                meshes_are_equal *= np.allclose(
+                    mesh1.cell_data["Normals"], mesh2.cell_data["Normals"], rtol=rtol
+                )
             else:
-                meshes_are_equal  *= False
+                meshes_are_equal *= False
         else:
-            meshes_are_equal  *= True
+            meshes_are_equal *= True
     return meshes_are_equal
-        
-
-        
 
 
 def test_tri_areas_2():
@@ -89,87 +134,120 @@ def test_tri_areas_2():
     result = GF.tri_areas(input)
     assert np.allclose(result, refrence)
 
+
 def test_tri_centroids():
-    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'))
+    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(files(lyceanem.tests.data).joinpath('rotated_recieve.ply'))
-    input = meshio.read(files(lyceanem.tests.data).joinpath('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])
- 
 
     result = GF.mesh_rotate(input, rotation_vector1, centre)
-    
+
     assert are_meshes_equal(result, refrence)
 
+
 def test_tri_areas(standard_cube):
-    #define a cube with side length 1m, and test that the triangle areas are correct
-    result=np.full((12),0.5)
-    assert np.all(GF.tri_areas(standard_cube)==result)
+    # define a cube with side length 1m, and test that the triangle areas are correct
+    result = np.full((12), 0.5)
+    assert np.all(GF.tri_areas(standard_cube) == result)
+
 
 def test_axes_from_normal_x_x():
     boresight = np.zeros((3), dtype=np.float32)
     boresight[0] = 1
-    rotation_matrix = R.from_euler('z',0,degrees=True).as_matrix()
-    assert_allclose(GF.axes_from_normal(boresight,boresight_along='x'),rotation_matrix,atol=1e-12)
+    rotation_matrix = R.from_euler("z", 0, degrees=True).as_matrix()
+    assert_allclose(
+        GF.axes_from_normal(boresight, boresight_along="x"), rotation_matrix, atol=1e-12
+    )
+
 
 def test_axes_from_normal_x_y():
     boresight = np.zeros((3), dtype=np.float32)
     boresight[1] = 1
-    rotation_matrix = R.from_euler('z',90,degrees=True).as_matrix()
-    assert_allclose(GF.axes_from_normal(boresight,boresight_along='x'),rotation_matrix,atol=1e-12)
+    rotation_matrix = R.from_euler("z", 90, degrees=True).as_matrix()
+    assert_allclose(
+        GF.axes_from_normal(boresight, boresight_along="x"), rotation_matrix, atol=1e-12
+    )
+
 
 def test_axes_from_normal_x_z():
     boresight = np.zeros((3), dtype=np.float32)
     boresight[2] = 1
-    rotation_matrix = R.from_euler('y',-90,degrees=True).as_matrix()
-    assert_allclose(GF.axes_from_normal(boresight,boresight_along='x'),rotation_matrix,atol=1e-12)
+    rotation_matrix = R.from_euler("y", -90, degrees=True).as_matrix()
+    assert_allclose(
+        GF.axes_from_normal(boresight, boresight_along="x"), rotation_matrix, atol=1e-12
+    )
+
 
 def test_axes_from_normal_x_nz():
     boresight = np.zeros((3), dtype=np.float32)
     boresight[2] = -1
-    rotation_matrix = R.from_euler('y',90,degrees=True).as_matrix()
-    assert_allclose(GF.axes_from_normal(boresight,boresight_along='x'),rotation_matrix,atol=1e-12)
+    rotation_matrix = R.from_euler("y", 90, degrees=True).as_matrix()
+    assert_allclose(
+        GF.axes_from_normal(boresight, boresight_along="x"), rotation_matrix, atol=1e-12
+    )
 
 
 def test_axes_from_normal_y_x():
     boresight = np.zeros((3), dtype=np.float32)
     boresight[0] = 1
-    rotation_matrix = R.from_euler('z',-90,degrees=True).as_matrix()
-    assert_allclose(GF.axes_from_normal(boresight,boresight_along='y'),rotation_matrix,atol=1e-12)
+    rotation_matrix = R.from_euler("z", -90, degrees=True).as_matrix()
+    assert_allclose(
+        GF.axes_from_normal(boresight, boresight_along="y"), rotation_matrix, atol=1e-12
+    )
+
 
 def test_axes_from_normal_y_y():
     boresight = np.zeros((3), dtype=np.float32)
     boresight[1] = 1
-    rotation_matrix = R.from_euler('z',0,degrees=True).as_matrix()
-    assert_allclose(GF.axes_from_normal(boresight,boresight_along='y'),rotation_matrix,atol=1e-12)
+    rotation_matrix = R.from_euler("z", 0, degrees=True).as_matrix()
+    assert_allclose(
+        GF.axes_from_normal(boresight, boresight_along="y"), rotation_matrix, atol=1e-12
+    )
+
 
 def test_axes_from_normal_y_z():
     boresight = np.zeros((3), dtype=np.float32)
     boresight[2] = 1
-    rotation_matrix = R.from_euler('x',90,degrees=True).as_matrix()
-    assert_allclose(GF.axes_from_normal(boresight,boresight_along='y'),rotation_matrix,atol=1e-12)
+    rotation_matrix = R.from_euler("x", 90, degrees=True).as_matrix()
+    assert_allclose(
+        GF.axes_from_normal(boresight, boresight_along="y"), rotation_matrix, atol=1e-12
+    )
+
 
 def test_axes_from_normal_z_x():
     boresight = np.zeros((3), dtype=np.float32)
     boresight[0] = 1
-    rotation_matrix = R.from_euler('y', 90, degrees=True).as_matrix()
-    assert_allclose(GF.axes_from_normal(boresight,boresight_along='z'),rotation_matrix,atol=1e-12)
+    rotation_matrix = R.from_euler("y", 90, degrees=True).as_matrix()
+    assert_allclose(
+        GF.axes_from_normal(boresight, boresight_along="z"), rotation_matrix, atol=1e-12
+    )
 
 
 def test_axes_from_normal_z_y():
     boresight = np.zeros((3), dtype=np.float32)
     boresight[1] = 1
-    rotation_matrix = R.from_euler('x',-90,degrees=True).as_matrix()
-    assert_allclose(GF.axes_from_normal(boresight,boresight_along='z'),rotation_matrix,atol=1e-12)
+    rotation_matrix = R.from_euler("x", -90, degrees=True).as_matrix()
+    assert_allclose(
+        GF.axes_from_normal(boresight, boresight_along="z"), rotation_matrix, atol=1e-12
+    )
+
 
 def test_axes_from_normal_z_z():
     boresight = np.zeros((3), dtype=np.float32)
     boresight[2] = 1
-    rotation_matrix = R.from_euler('x',0,degrees=True).as_matrix()
-    assert_allclose(GF.axes_from_normal(boresight,boresight_along='z'),rotation_matrix,atol=1e-12)
+    rotation_matrix = R.from_euler("x", 0, degrees=True).as_matrix()
+    assert_allclose(
+        GF.axes_from_normal(boresight, boresight_along="z"), rotation_matrix, atol=1e-12
+    )
diff --git a/lyceanem/tests/test_ray_functions.py b/lyceanem/tests/test_ray_functions.py
index fc67be5..bb278d0 100644
--- a/lyceanem/tests/test_ray_functions.py
+++ b/lyceanem/tests/test_ray_functions.py
@@ -1,4 +1,4 @@
-from ..raycasting import rayfunctions as RF 
+from ..raycasting import rayfunctions as RF
 import numpy as np
 import pytest
 import lyceanem.base_types as base_types
@@ -6,8 +6,11 @@
 import lyceanem.tests.data
 import meshio
 
+
 def test_convertTriangles():
-    reference = np.load(files(lyceanem.tests.data).joinpath("triangle_type_array_ref.npy"))
+    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):
@@ -20,6 +23,3 @@ def test_convertTriangles():
         assert np.allclose(triangles[i]["v2x"], reference[i]["v2x"])
         assert np.allclose(triangles[i]["v2y"], reference[i]["v2y"])
         assert np.allclose(triangles[i]["v2z"], reference[i]["v2z"])
-
-
-
diff --git a/lyceanem/utility/math_functions.py b/lyceanem/utility/math_functions.py
index 68119d6..1bdf2fb 100644
--- a/lyceanem/utility/math_functions.py
+++ b/lyceanem/utility/math_functions.py
@@ -1,4 +1,3 @@
-
 import math
 
 import numpy as np
@@ -7,7 +6,7 @@
 
 @njit
 def cart2pol(x, y):
-    rho = np.sqrt(x ** 2 + y ** 2)
+    rho = np.sqrt(x**2 + y**2)
     phi = np.arctan2(y, x)
     return (rho, phi)
 
@@ -18,17 +17,19 @@ def pol2cart(rho, phi):
     y = rho * np.sin(phi)
     return (x, y)
 
-def polar2cart(theta,phi,r):
-    x = r*np.sin(phi) * np.cos(theta)
-    y = r*np.sin(phi) * np.sin(theta)
-    z = r*np.cos(phi)
-    return x,y,z
 
-def cart2polar(x,y,z):
-    r=np.sqrt(x ** 2 + y ** 2 + z ** 2)
-    theta=np.arctan2(y, x)
-    phi =np.arctan2(np.sqrt(x ** 2 + y ** 2), z)
-    return theta,phi,r
+def polar2cart(theta, phi, r):
+    x = r * np.sin(phi) * np.cos(theta)
+    y = r * np.sin(phi) * np.sin(theta)
+    z = r * np.cos(phi)
+    return x, y, z
+
+
+def cart2polar(x, y, z):
+    r = np.sqrt(x**2 + y**2 + z**2)
+    theta = np.arctan2(y, x)
+    phi = np.arctan2(np.sqrt(x**2 + y**2), z)
+    return theta, phi, r
 
 
 @njit
@@ -66,7 +67,7 @@ def calc_normals(T):
     diry = e1z * e2x - e1x * e2z
     dirz = e1x * e2y - e1y * e2x
 
-    normconst = math.sqrt(dirx ** 2 + diry ** 2 + dirz ** 2)
+    normconst = math.sqrt(dirx**2 + diry**2 + dirz**2)
 
     T["normx"] = dirx / normconst
     T["normy"] = diry / normconst
@@ -74,7 +75,10 @@ def calc_normals(T):
 
     return T
 
-def discrete_transmit_power(weights,element_area,transmit_power=100.0,impedance=np.pi*120.0):
+
+def discrete_transmit_power(
+    weights, element_area, transmit_power=100.0, impedance=np.pi * 120.0
+):
     """
     Calculate the transmitting aperture amplitude density required for a given transmit power in watts.
     Parameters
@@ -88,16 +92,25 @@ def discrete_transmit_power(weights,element_area,transmit_power=100.0,impedance=
     -------
 
     """
-    #Calculate the Power at each mesh point from the Intensity at each point
-    power_at_point=np.linalg.norm(weights,axis=1).reshape(-1,1)*element_area.reshape(-1,1) # calculate power
-    #integrate power over aperture and normalise to desired power for whole aperture
-    power_normalisation=transmit_power/np.sum(power_at_point)
-    transmit_power_density=(power_at_point*power_normalisation)/element_area.reshape(-1,1)
-    #calculate amplitude (V/m)
-    transmit_amplitude=((transmit_power_density*impedance)**0.5)*(weights/np.linalg.norm(weights,axis=1).reshape(-1,1))
-    #transmit_excitation=transmit_amplitude_density.reshape(-1,1)*element_area.reshape(-1,1)
+    # Calculate the Power at each mesh point from the Intensity at each point
+    power_at_point = np.linalg.norm(weights, axis=1).reshape(
+        -1, 1
+    ) * element_area.reshape(
+        -1, 1
+    )  # calculate power
+    # integrate power over aperture and normalise to desired power for whole aperture
+    power_normalisation = transmit_power / np.sum(power_at_point)
+    transmit_power_density = (
+        power_at_point * power_normalisation
+    ) / element_area.reshape(-1, 1)
+    # calculate amplitude (V/m)
+    transmit_amplitude = ((transmit_power_density * impedance) ** 0.5) * (
+        weights / np.linalg.norm(weights, axis=1).reshape(-1, 1)
+    )
+    # transmit_excitation=transmit_amplitude_density.reshape(-1,1)*element_area.reshape(-1,1)
     return transmit_amplitude
 
+
 @guvectorize(
     [(float32[:], float32[:], float32[:], float32)],
     "(n),(n)->(n),()",
@@ -107,7 +120,7 @@ def fast_calc_dv(source, target, dv, normconst):
     dirx = target[0] - source[0]
     diry = target[1] - source[1]
     dirz = target[2] - source[2]
-    normconst = math.sqrt(dirx ** 2 + diry ** 2 + dirz ** 2)
+    normconst = math.sqrt(dirx**2 + diry**2 + dirz**2)
     dv = np.array([dirx, diry, dirz]) / normconst
 
 
@@ -116,7 +129,7 @@ def calc_dv(source, target):
     dirx = target[0] - source[0]
     diry = target[1] - source[1]
     dirz = target[2] - source[2]
-    normconst = np.sqrt(dirx ** 2 + diry ** 2 + dirz ** 2)
+    normconst = np.sqrt(dirx**2 + diry**2 + dirz**2)
     dv = np.array([dirx, diry, dirz]) / normconst
     return dv[0], dv[1], dv[2], normconst
 
@@ -131,6 +144,11 @@ def calc_dv_norm(source, target, direction, length):
     direction = (target - source) / length
     return direction, length
 
-def FindAlignedVector(command_vector,directions):
-    index_vector=np.argmax(np.einsum('nm,nm->n',np.tile(command_vector,(directions.shape[0],1)),directions))
-    return index_vector
\ No newline at end of file
+
+def FindAlignedVector(command_vector, directions):
+    index_vector = np.argmax(
+        np.einsum(
+            "nm,nm->n", np.tile(command_vector, (directions.shape[0], 1)), directions
+        )
+    )
+    return index_vector