diff --git a/poke/material_data/Mathewson_Al.csv b/poke/material_data/Mathewson_Al.csv
new file mode 100644
index 0000000..977908a
--- /dev/null
+++ b/poke/material_data/Mathewson_Al.csv
@@ -0,0 +1,77 @@
+wl,n
+0.49594,0.77909
+0.50606,0.81783
+0.51660,0.85227
+0.52759,0.89107
+0.53906,0.94408
+0.55104,1.00461
+0.56356,1.06864
+0.57667,1.14503
+0.59040,1.21981
+0.60480,1.30548
+0.61992,1.40111
+0.63582,1.50675
+0.65255,1.62810
+0.67018,1.77579
+0.68880,1.94186
+0.70848,2.13059
+0.72932,2.33676
+0.75142,2.57452
+0.77490,2.87316
+0.79990,3.03948
+0.82656,2.93996
+0.85506,2.52983
+0.88560,2.13552
+0.91840,1.81758
+0.95372,1.57683
+0.99187,1.41421
+1.03320,1.31513
+1.07812,1.24965
+1.12713,1.22635
+1.18080,1.22466
+1.23984,1.25115
+1.30510,1.29739
+1.37760,1.34628
+1.45864,1.42835
+1.54980,1.50000
+1.65312,1.58999
+1.77120,1.75488
+
+wl,k
+0.49594,5.84012
+0.50606,5.93033
+0.51660,6.10134
+0.52759,6.22848
+0.53906,6.35541
+0.55104,6.51991
+0.56356,6.64394
+0.57667,6.76839
+0.59040,6.92733
+0.60480,7.08550
+0.61992,7.20855
+0.63582,7.36684
+0.65255,7.49338
+0.67018,7.65855
+0.68880,7.82757
+0.70848,8.00246
+0.72932,8.15233
+0.75142,8.25398
+0.77490,8.23135
+0.79990,8.07703
+0.82656,7.75521
+0.85506,7.60921
+0.88560,7.74987
+0.91840,8.14270
+0.95372,8.68829
+0.99187,9.19239
+1.03320,9.73291
+1.07812,10.32287
+1.12713,10.88595
+1.18080,11.55421
+1.23984,12.18874
+1.30510,12.91059
+1.37760,13.70447
+1.45864,14.52722
+1.54980,15.50000
+1.65312,16.47811
+1.77120,17.75048
diff --git a/poke/materials.py b/poke/materials.py
index 1cb1f0a..9145c80 100644
--- a/poke/materials.py
+++ b/poke/materials.py
@@ -5,7 +5,7 @@
 # Silica == Fused Silica
 # fmt: off
 avail_materials = [
-    "Al", "Ag",  # metals
+    "Al", "Ag", "Mathewson_Al",  # metals
     "HfO2", "SiO2", "Ta2O5", "TiO2", "Nb2O5",  # Oxides
     "SiN",  # Nitrides
     "MgF2", "CaF2", "LiF",  # Fluorides
@@ -59,6 +59,10 @@ def create_index_model(material, verbose=False):
         pth = get_abs_path("Ciesielski_Ag.csv")
         n_end = 333
         k_start = n_end + 3
+    elif material == "Mathewson_Al":
+        pth = get_abs_path("Mathewson_Al.csv")
+        n_end = 37
+        k_start = n_end + 3
     elif material == "HfO2":
         pth = get_abs_path("Kuhaili_HfO2.csv")
     elif material == "SiO2":
diff --git a/poke/poke_math.py b/poke/poke_math.py
index f76af2f..ad91c88 100644
--- a/poke/poke_math.py
+++ b/poke/poke_math.py
@@ -282,6 +282,31 @@ def rotation_3d(angle, axis):
     return mat
 
 
+def broadcast_kron(a,b):
+    """broadcasted kronecker product of two N,M,...,2,2 arrays. Used for jones -> mueller conversion
+    In the unbroadcasted case, this output looks like
+
+    out = [a[0,0]*b,a[0,1]*b]
+          [a[1,0]*b,a[1,1]*b]
+
+    where out is a N,M,...,4,4 array. I wrote this to work for generally shaped kronecker products where the matrix
+    is contained in the last two axes, but it's only tested for the Nx2x2 case
+
+    Parameters
+    ----------
+    a : numpy.ndarray
+        N,M,...,2,2 array used to scale b in kronecker product
+    b : numpy.ndarray
+        N,M,...,2,2 array used to form block matrices in kronecker product
+
+    Returns
+    -------
+    out
+        N,M,...,4,4 array
+    """
+
+    return np.einsum('...ik,...jl',a,b).reshape([*a.shape[:-2],int(a.shape[-2]*b.shape[-2]),int(a.shape[-1]*b.shape[-1])])
+
 "Vector Operations from Quinn Jarecki"
 
 
diff --git a/poke/polarization.py b/poke/polarization.py
index 7e9cef8..d3ea0d2 100644
--- a/poke/polarization.py
+++ b/poke/polarization.py
@@ -5,6 +5,7 @@
     vector_norm,
     vector_angle,
     rotation_3d,
+    broadcast_kron,
     vectorAngle,
     rotation3D,
 )
@@ -413,6 +414,21 @@ def jones_to_mueller(Jones):
     return M
 
 
+def jones_to_mueller_broadcast(jones):
+
+    U = np.array([[1,0,0,1],
+                  [1,0,0,-1],
+                  [0,1,1,0],
+                  [0,1j,-1j,0]])
+
+    U *= np.sqrt(1/2)
+    Uinv = np.linalg.inv(U)
+    inner = broadcast_kron(jones.conj(), jones)
+
+    M = np.real(U @ inner @ Uinv)
+    return M
+
+
 def mueller_to_jones(M):
     """Converts Mueller matrix to a relative Jones matrix. Phase aberration is relative to the Pxx component.