diff --git a/lyceanem/electromagnetics/emfunctions.py b/lyceanem/electromagnetics/emfunctions.py
index e15d2fb..a8afa9a 100644
--- a/lyceanem/electromagnetics/emfunctions.py
+++ b/lyceanem/electromagnetics/emfunctions.py
@@ -50,9 +50,9 @@ def field_magnitude_phase(field_data):
 
     return field_data
 def extract_electric_fields(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'],
+    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
@@ -93,9 +93,9 @@ def EthetaEphi_to_Exyz(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()
-    theta = field_data.point_Data['theta (Radians)']
-    phi = field_data.point_Data['phi (Radians)']
+    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)
     ephi = -electric_fields[:, 0] * np.sin(phi) + electric_fields[:, 1] * np.cos(phi)
@@ -169,7 +169,7 @@ def PoyntingVector(field_data):
 
 def Directivity(field_data):
     # calculate directivity for the given pattern
-    field_area = field_data.area
+
     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
         field_data = theta_phi_r(field_data)
diff --git a/lyceanem/geometry/geometryfunctions.py b/lyceanem/geometry/geometryfunctions.py
index 88e51d9..9656680 100644
--- a/lyceanem/geometry/geometryfunctions.py
+++ b/lyceanem/geometry/geometryfunctions.py
@@ -107,6 +107,49 @@ def mesh_transform(mesh, transform_matrix, rotate_only):
 
     return return_mesh
 
+
+def compute_areas(field_data):
+    cell_areas = []
+    for inc, cell in enumerate(field_data.cells):
+
+        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)
+            s = (edge1 + edge2 + edge3) / 2
+            areas = (s * (s - edge1) * (s - edge2) * (s - edge3)) ** 0.5
+            cell_areas.append(areas)
+        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)
+
+            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
+            cell_areas.append(areas)
+
+    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]])
+
+    return field_data
+
 def compute_normals(mesh):
     """
     Computes the Cell Normals for meshio mesh objects, the point normals will also be calculated for all points which are connected to a triangle, isolated points will be propulated with a normal vector of nan.