diff --git a/contrib/anisotropic_plagioclase/plagioclase_data.py b/contrib/anisotropic_plagioclase/plagioclase_data.py
index 4b2c84b3..c5246ebb 100644
--- a/contrib/anisotropic_plagioclase/plagioclase_data.py
+++ b/contrib/anisotropic_plagioclase/plagioclase_data.py
@@ -80,6 +80,10 @@ def get_data():
     b = np.moveaxis(b, 1, 0)
     b_err = np.moveaxis(b_err, 1, 0)
 
+    # Brown reports the sum in Voigt form (Section 3.1)
+    b[:, 3:] = b[:, 3:] / 2
+    b_err[:, 3:] = b_err[:, 3:] / 2
+
     data["beta"] = {
         "P": d[:, 0] * 0.0 + 1.0e5,
         "T": d[:, 0] * 0.0 + 298.15,
diff --git a/contrib/anisotropic_plagioclase/plagioclase_model.py b/contrib/anisotropic_plagioclase/plagioclase_model.py
index 82b5385d..adb1875f 100644
--- a/contrib/anisotropic_plagioclase/plagioclase_model.py
+++ b/contrib/anisotropic_plagioclase/plagioclase_model.py
@@ -5,6 +5,7 @@
 import numpy as np
 from copy import deepcopy
 from burnman.utils.unitcell import cell_parameters_to_vectors
+from tabulate import tabulate
 
 
 def make_scalar_model(args):
@@ -14,6 +15,11 @@ def make_scalar_model(args):
     ab_scalar = SLB_2022.albite()
     an_scalar = SLB_2022.anorthite()
 
+    ab_scalar.params["V_0"] = ab_scalar.params["V_0"] + dV_ab
+    an_scalar.params["V_0"] = an_scalar.params["V_0"] + dV_an
+    ab_scalar.params["K_0"] = ab_scalar.params["K_0"] + dK_ab
+    an_scalar.params["K_0"] = an_scalar.params["K_0"] + dK_an
+
     aban_linear = SLB_2022.anorthite()
     aban_linear.params["V_0"] = 0.5 * (
         ab_scalar.params["V_0"] + an_scalar.params["V_0"]
@@ -21,17 +27,34 @@ def make_scalar_model(args):
     aban_linear.params["K_0"] = 0.5 * (
         ab_scalar.params["K_0"] + an_scalar.params["K_0"]
     )
+    aban_linear.params["Kprime_0"] = 0.5 * (
+        ab_scalar.params["Kprime_0"] + an_scalar.params["Kprime_0"]
+    )
 
     aban_scalar = deepcopy(aban_linear)
 
-    ab_scalar.params["V_0"] = ab_scalar.params["V_0"] + dV_ab
-    an_scalar.params["V_0"] = an_scalar.params["V_0"] + dV_an
     aban_scalar.params["V_0"] = aban_scalar.params["V_0"] + dV_aban
-
-    ab_scalar.params["K_0"] = ab_scalar.params["K_0"] + dK_ab
-    an_scalar.params["K_0"] = an_scalar.params["K_0"] + dK_an
     aban_scalar.params["K_0"] = aban_scalar.params["K_0"] + dK_aban
 
+    table = [
+        ["", "ab", "an", "an$_{50}$ (1)", "an$_{50}$ (2)"],
+        [
+            "V_0 (cm$^3$/mol)",
+            ab_scalar.params["V_0"] * 1.0e6,
+            an_scalar.params["V_0"] * 1.0e6,
+            aban_linear.params["V_0"] * 1.0e6,
+            aban_scalar.params["V_0"] * 1.0e6,
+        ],
+        [
+            "K_0 (GPa)",
+            ab_scalar.params["K_0"] / 1.0e9,
+            an_scalar.params["K_0"] / 1.0e9,
+            aban_linear.params["K_0"] / 1.0e9,
+            aban_scalar.params["K_0"] / 1.0e9,
+        ],
+    ]
+    print(tabulate(table, headers="firstrow", tablefmt="latex_raw", floatfmt=".6e"))
+
     class plagioclase_scalar(Solution):
         def __init__(self, molar_fractions=None):
             self.name = "plagioclase (NCAS)"
@@ -77,6 +100,16 @@ def make_anisotropic_model(scalar_args, cell_args, elastic_args):
     f = np.cbrt(an_scalar.params["V_0"] / np.linalg.det(M))
     an_cell_parameters[:3] = an_cell_parameters[:3] * f
 
+    table = [
+        ["", "$a$", "$b$", "$c$", "$\\alpha$", "$\\beta$", "$\\gamma$"],
+        ["ab"],
+        ["an"],
+    ]
+    table[1].extend(ab_cell_parameters)
+    table[2].extend(an_cell_parameters)
+
+    print(tabulate(table, headers="firstrow", tablefmt="latex_raw", floatfmt=".6e"))
+
     an_a = np.zeros((6, 6))
     ab_a = np.zeros((6, 6))
 
@@ -136,6 +169,12 @@ def make_anisotropic_model(scalar_args, cell_args, elastic_args):
     an_a[0, 0] = 1.0 - np.sum(an_a[:3, :3])
     ab_a[0, 0] = 1.0 - np.sum(ab_a[:3, :3])
 
+    table = ab_a
+    print(tabulate(table, tablefmt="latex_raw", floatfmt=".5f"))
+
+    table = an_a
+    print(tabulate(table, tablefmt="latex_raw", floatfmt=".5f"))
+
     # print(an_a)
     # print(ab_a)
     # exit()
diff --git a/contrib/anisotropic_plagioclase/plagioclase_model_plots.py b/contrib/anisotropic_plagioclase/plagioclase_model_plots.py
index e6b372b6..a858a0c0 100644
--- a/contrib/anisotropic_plagioclase/plagioclase_model_plots.py
+++ b/contrib/anisotropic_plagioclase/plagioclase_model_plots.py
@@ -3,6 +3,9 @@
 from plagioclase_parameters import scalar_args, cell_args, elastic_args
 from plagioclase_model import make_anisotropic_model
 from plagioclase_data import get_data
+from burnman.minerals.SLB_2022 import plagioclase
+
+ss_SLB = plagioclase()
 
 ss = make_anisotropic_model(scalar_args, cell_args, elastic_args)
 
@@ -26,10 +29,20 @@
     molar_fractions,
 )
 
+prps_SLB = ss_SLB.evaluate(
+    ["molar_volume", "isothermal_bulk_modulus_reuss"],
+    pressures,
+    temperatures,
+    molar_fractions,
+)
+
+
 labels = ["V", "KTR", "CT", "CN", "ST", "betaT", "cell"]
 prps = {labels[i]: prps[i] for i in range(7)}
 prps["psi"] = np.einsum("ijk, i->ijk", prps["ST"], prps["KTR"])
 
+prps_SLB = {labels[i]: prps_SLB[i] for i in range(2)}
+
 d = get_data()
 
 
@@ -40,7 +53,15 @@
 
 mask2 = p_ans < 0.5
 
-ln = ax[0].plot(p_ans[mask2], prps["V"][mask2] * 1.0e6)
+ln = ax[0].plot(p_ans[mask2], prps["V"][mask2] * 1.0e6, label="C$\\bar{1}$, this study")
+ln = ax[0].plot(
+    p_ans,
+    prps_SLB["V"] * 1.0e6,
+    linestyle="--",
+    color=ln[0].get_color(),
+    label="SLB2022",
+)
+
 ax[0].plot(
     p_ans[~mask2], prps["V"][~mask2] * 1.0e6, color=ln[0].get_color(), linestyle=":"
 )
@@ -72,7 +93,8 @@
 )
 
 
-ax[1].plot(p_ans[mask2], prps["KTR"][mask2] / 1.0e9)
+ax[1].plot(p_ans[mask2], prps["KTR"][mask2] / 1.0e9, color=ln[0].get_color())
+ax[1].plot(p_ans, prps_SLB["KTR"] / 1.0e9, linestyle="--", color=ln[0].get_color())
 ax[1].plot(
     p_ans[~mask2], prps["KTR"][~mask2] / 1.0e9, color=ln[0].get_color(), linestyle=":"
 )
@@ -108,9 +130,11 @@
 for i in range(2):
     ax[i].set_xlabel("$p_{an}$")
 
-ax[0].set_ylabel("$V$ (m$^3$/mol)")
+ax[0].set_ylabel("$V$ (cm$^3$/mol)")
 ax[1].set_ylabel("$K_{\\text{TR}}$ (GPa)")
 
+ax[0].legend()
+
 fig.set_tight_layout(True)
 fig.savefig("plag_V_KT.pdf")
 
@@ -189,7 +213,7 @@
 for i in range(7):
     ax[i].legend()
     ax[i].set_xlabel("$p_{an}$")
-    ax[i].set_ylabel("modulus (GPa)")
+    ax[i].set_ylabel("$C_{Nij}$ (GPa)")
 
 fig.set_tight_layout(True)
 fig.savefig("plag_stiffnesses.pdf")
@@ -198,7 +222,7 @@
 fig = plt.figure(figsize=(10, 8))
 ax = [fig.add_subplot(3, 3, i) for i in range(1, 8)]
 for irow, (i, j) in enumerate(inds):
-    name = f"$d\\Psi_{{{i+1}{j+1}}} / df$"
+    name = f"$S_{{N{i+1}{j+1}}} / \\beta_{{NR}}$"
     axi = int((irow - (irow) % 3) / 3)
     ln = ax[axi].plot(p_ans[mask2], prps["psi"][mask2, i, j])
     ax[axi].plot(
@@ -220,7 +244,7 @@
 for i in range(7):
     ax[i].legend()
     ax[i].set_xlabel("$p_{an}$")
-    ax[i].set_ylabel("$d\\Psi / df$")
+    ax[i].set_ylabel("$S_{Nij} / \\beta_{NR}$")
 
 fig.set_tight_layout(True)
 fig.savefig("plag_psi.pdf")
@@ -329,7 +353,7 @@
         alpha=0.5,
     )
 
-    ax[axi].set_ylabel(f"$\\beta_{{{axi}}}$ (m)")
+    ax[axi].set_ylabel(f"$\\beta_{{T{axi+1}}}$ (GPa$^{{-1}}$)")
     ax[axi].set_xlabel("$p_{an}$")
 
 fig.set_tight_layout(True)