diff --git a/functions/copula/normal_copula.stanfunctions b/functions/copula/normal_copula.stanfunctions
index 31bd81e..c631746 100644
--- a/functions/copula/normal_copula.stanfunctions
+++ b/functions/copula/normal_copula.stanfunctions
@@ -61,46 +61,6 @@ real normal_copula_vector_lpdf(vector u, vector v, real rho) {
   return a1 * x / a2 - N * a3;
 }
 
-/**
-  * Multi-Normal Cholesky copula log density
-  *
-  * \f{aligned}{
-  *    c(\mathbf{u}) &= \frac{\partial^d C}{\partial \mathbf{\Phi_1}\cdots \partial \mathbf{\Phi_d}} \\
-  *                  &= \frac{1}{\sqrt{\det \Sigma}} \exp \Bigg(-\frac{1}{2}
-  *                              \begin{pmatrix} \mathbf{\Phi}^{-1}(u_1) \\ \vdots  \\ \mathbf{\Phi}^{-1}(u_d)  \end{pmatrix}^T
-  *                              (\Sigma^{-1} - I)
-  *                              \begin{pmatrix} \mathbf{\Phi}^{-1}(u_1) \\ \vdots  \\ \mathbf{\Phi}^{-1}(u_d)  \end{pmatrix}
-  *                              \Bigg) \\
-  *                  &=  \frac{1}{\sqrt{\det \Sigma}} \exp \bigg(-\frac{1}{2}(A^T\Sigma^{-1}A - A^TA) \bigg)\f}
-  * where
-  * \f[
-  *    A = \begin{pmatrix} \mathbf{\Phi}^{-1}(u_1) \\ \vdots  \\ \mathbf{\Phi}^{-1}(u_d)  \end{pmatrix}.
-  *  \f]
-  * Note that \f$\det \Sigma  = |LL^T| = |L |^2 = \big(\prod_{i=1}^d L_{i,i}\big)^2\f$ so \f$\sqrt{\det \Sigma} = \prod_{i=1}^d L_{i,i}\f$
-  * and
-  *
-  *\f{aligned}{
-  *  c(\mathbf{u}) &= \Bigg(\prod_{i=1}^d L_{i,i}\Bigg)^{-1} \exp \bigg(-\frac{1}{2}(A^T(LL^T)^{-1}A - A^TA)\bigg) \\
-  *                &=  \Bigg(\prod_{i=1}^d L_{i,i}\Bigg)^{-1} \exp \bigg(-\frac{1}{2}\big((L^{-1}A)^TL^{-1}A - A^TA\big)\bigg)
-    \f}
-  * where \f$X = L^{-1}A\f$ can be solved with \code{.cpp} X = mdivide_left_tri_low(L, A)\endcode.
-  *
-  * @copyright Sean Pinkney, 2021
-  *
-  * @param u Matrix
-  * @param L Cholesky factor matrix
-  * @return log density
-  */
-
-real multi_normal_copula_lpdf(matrix u, matrix L) {
-  int K = rows(u);
-  int N = cols(u);
-  real inv_sqrt_det_log = sum(log(diagonal(L)));
-  matrix[K, N] x = mdivide_left_tri_low(L, u);
-  
-  return -N * inv_sqrt_det_log - 0.5 * sum(columns_dot_self(x) - columns_dot_self(u));
-}
-
 /**
   * Bivariate Normal Copula cdf
   *
@@ -133,41 +93,25 @@ real bivariate_normal_copula_cdf(vector u, real rho) {
   * Multivariate Normal Copula log density (Cholesky parameterisation)
   *
   * \f{aligned}{
-  *    c(\mathbf{u}) &= \frac{\partial^d C}{\partial \mathbf{\Phi_1}\cdots \partial \mathbf{\Phi_d}} \\
-  *                  &= \prod_{i=1}^{n} \lvert \Sigma \rvert^{-1/2} \exp \left\{-\frac{1}{2}
+  * c(\mathbf{u}) &= \frac{\partial^d C}{\partial \mathbf{\Phi_1}\cdots \partial \mathbf{\Phi_d}} \\
+  *               &= \prod_{i=1}^{n} \lvert \Sigma \rvert^{-1/2} \exp \left\{-\frac{1}{2}
   *                              \begin{pmatrix} \mathbf{\Phi}^{-1}(u_1) \\ \vdots  \\ \mathbf{\Phi}^{-1}(u_d)  \end{pmatrix}^T
   *                              (\Sigma^{-1} - I)
   *                              \begin{pmatrix} \mathbf{\Phi}^{-1}(u_1) \\ \vdots  \\ \mathbf{\Phi}^{-1}(u_d)  \end{pmatrix}
   *                              \right\} \\
   *                 &= \lvert \Sigma \rvert^{-n/2} \exp\left\{ -\frac{1}{2} \text{tr}\left( \left[ \Sigma^{-1} - I \right] \sum_{i=1}^n \Phi^{-1}(u_i) \Phi^{-1}(u_i)' \right) \right\} \hspace{0.5cm} \\
-  *                 &= \lvert L \rvert^{-n} \exp\left\{  -\frac{1}{2} \text{tr}\left( \left[ L^{-T}L^{-1}- I \right] Q'Q  \right) \\
-  *                 &= \lvert L \rvert^{-n} \exp\left\{ -\frac{1}{2} \text{tr}\left( L^{-T}L^{-1} Q'Q - Q'Q \right) \right\}  \\
-  *                 &= \lvert L \rvert^{-n} \exp\left\{ -\frac{1}{2} \text{tr}\left(  QL^{-T}L^{-1}Q' \right) + \text{tr} \left( Q'Q \right) \right\} \\
-  *                 &= \lvert L \rvert^{-n} \exp\left\{ -\frac{1}{2} \left[ \text{tr}\left(  [L^{-1}Q'][L^{-1}Q']' \right) + \text{tr} \left( Q'Q \right) \right]\right\}
+  *                 &= \lvert L \rvert^{-n} \exp\left\{  -\frac{1}{2} \text{tr}\left( \left[ L^{-T}L^{-1}- I \right] Q'Q  \right) \right\} \\
+  *                 &= \lvert L \rvert^{-n} \exp\left\{  -\frac{1}{2} \sum_{j=1}^{d}\sum_{i=1}^{n} \left( \left[ L^{-T}L^{-1}- I \right] \odot Q'Q  \right) \right\} \\
   *  \f}
   * where
   * \f[
-  *    Q_{n \times d} = \begin{pmatrix} \mathbf{\Phi}^{-1}(u_1) \\ \vdots  \\ \mathbf{\Phi}^{-1}(u_d)  \end{pmatrix}.
-  *  \f]
-  * Note that \f$\det \Sigma  = |LL^T| = |L |^2 = \big(\prod_{i=1}^d L_{i,i}\big)^2\f$ so \f$\sqrt{\det \Sigma} = \prod_{i=1}^d L_{i,i}\f$
-  * and
-  *
-  *\f{aligned}{
-  *  c(\mathbf{u}) &= \Bigg(\prod_{i=1}^d L_{i,i}\Bigg)^{-1} \exp \bigg(-\frac{1}{2}(A^T(LL^T)^{-1}A - A^TA)\bigg) \\
-  *                &=  \Bigg(\prod_{i=1}^d L_{i,i}\Bigg)^{-1} \exp \bigg(-\frac{1}{2}\big((L^{-1}A)^TL^{-1}A - A^TA\big)\bigg)
-    \f}
-  * where \f$X = L^{-1}A\f$ can be solved with \code{.cpp} X = mdivide_left_tri_low(L, A)\endcode.
+  *    Q_{n \times d} = \begin{pmatrix} \mathbf{\Phi}^{-1}(u_1) \\ \vdots  \\ \mathbf{\Phi}^{-1}(u_d)  \end{pmatrix}
+  *  \f] 
+  * and \f$ \odot \f$ is the elementwise multiplication operator.
+  * Note that \f$\det \Sigma  = |LL^T| = |L |^2 = \big(\prod_{i=1}^d L_{i,i}\big)^2\f$ so \f$\sqrt{\det \Sigma} = \prod_{i=1}^d L_{i,i}\f$.
   *
   * @copyright Ethan Alt, Sean Pinkney 2021
   *
-  * @param U Matrix where the number of rows equal observations and the columns are equal to the number of outcomes
-  * @param L Cholesky factor of the correlation matrix
-  * @return log density of the distribution
-  */
-
-/**
-  * Multivariate Normal Copula lpdf (Cholesky parameterisation)
-  *
   * @param U Matrix of outcomes from marginal calculations
   * @param L Cholesky factor of the correlation/covariance matrix
   * @return log density of distribution