diff --git a/examples/linear_algebra/R/correlation_angles_example.R b/examples/linear_algebra/R/correlation_angles_example.R
index 163ff7f2..dc96d3b9 100644
--- a/examples/linear_algebra/R/correlation_angles_example.R
+++ b/examples/linear_algebra/R/correlation_angles_example.R
@@ -1,12 +1,15 @@
 library(cmdstanr)
 fp <- file.path("./examples/linear_algebra/stan/correlation_angles_example.stan")
-mod <- cmdstan_model(fp, include_paths = "./functions/linear_algebra")
+mod <- cmdstan_model(fp, include_paths = "./functions/linear_algebra", force_recompile = T)
 
 T <- matrix(c(0.8, -0.9, -0.9,
               -0.9, 1.1, 0.3,
               -0.9, 0.4, 0.9), 3, 3)
 
-
+T <- matrix(c(1, -0.9, -0.9,
+              -0.9, 1, 0.3,
+              -0.9, 0.4, 1), 3, 3)
+chol(T)
 mod_out <- mod$optimize(
   data = list(N = 3,
               R = T)
@@ -17,12 +20,10 @@ mod_out <- mod$sample(
               R = T),
   chains = 2,
   seed = 23421,
+  
   parallel_chains = 2,
   iter_warmup = 600,
   iter_sampling = 600
 )
 
-mod_out$summary("R_hat")
-mat <- matrix(mod_out$summary("R_hat")$mean, 3, 3)
-mat
-chol(mat)
+mod_out$summary()
diff --git a/examples/linear_algebra/R/correlation_angles_example2.R b/examples/linear_algebra/R/correlation_angles_example2.R
new file mode 100644
index 00000000..50655c1b
--- /dev/null
+++ b/examples/linear_algebra/R/correlation_angles_example2.R
@@ -0,0 +1,120 @@
+library(cmdstanr)
+fp <- file.path("./examples/linear_algebra/stan/correlation_angles_constrain_example.stan")
+mod <- cmdstan_model(fp, include_paths = "./functions/linear_algebra", force_recompile = T)
+
+library(rstan)
+expose_stan_functions(fp)
+
+library(StanHeaders)
+stanFunction("multiply_lower_tri_self_transpose", x = L)
+
+
+N <- 3
+test_mat <- matrix(c(0.8, -0.9, -0.9,
+              -0.9, 1.1, 0.3,
+              -0.9, 0.4, 0.9), 3, 3)
+
+T <- matrix(c(1, -0.9, -0.9,
+              -0.9, 1, 0.3,
+              -0.9, 0.4, 1), 3, 3)
+
+
+test_mat <- matrix(c(1.0000,   0,         0,         0,   -0.9360,
+                     0,    1.0000,   -0.5500,   -0.3645,   -0.5300,
+                     0,   -0.5500,    1.0000,   0,    0.0875,
+                     0,   -0.3645,    0,    1.0000,    0.4557,
+                     -0.9360,   -0.5300,    0.0875,    0.4557,    1.0000),
+                   byrow = T, 5, 5)
+
+where_zero <- (test_mat[upper.tri(test_mat)] == 0) * 1
+sum(where_zero) * 1
+
+fit <- stan(file = fp, 
+            data = list(N = nrow(test_mat),
+                        R = test_mat,
+                        is_symmetric = 1,
+                        Z = sum(where_zero),
+                        K = ( nrow(test_mat) * (nrow(test_mat) - 1 )) / 2,
+                        where_zero = where_zero ),
+            chains = 1)
+
+K <- 5
+theta_test <- build_angle_mat(where_zero, runif(6, 0, pi), K)
+L = zero_constrain(theta_test, K)
+tcrossprod(L)
+
+angle_raw <-  runif(6)
+N <- length(where_zero);
+mat <- matrix(0, N, N)
+count <- 1;
+raw_count <- 1;
+
+for (k in 2:K){
+  mat[k, k] = 0;
+  for (j in 1:k - 1) {
+    if (where_zero[raw_count] != 1) {
+      mat[k, j] = angle_raw[count];
+      count <= count +  1;
+    }
+    raw_count = raw_count + 1;
+  }
+}
+
+count <- 0
+
+for (k in 2:K){
+  for (j in 1:(k - 1)) {
+    count <- count + 1
+    print(count)
+  }}
+
+
+
+mod_out <- mod$sample(
+  data = list(N = nrow(test_mat),
+              R = test_mat,
+              is_symmetric = 1,
+              Z = sum(where_zero),
+              K = ( nrow(test_mat) * (nrow(test_mat) - 1 )) / 2,
+              where_zero = where_zero ),
+  chains = 2,
+  #init = 0,
+  #seed = 23421,
+  #adapt_delta = 0.8,
+  #max_treedepth = 10,
+  parallel_chains = 4,
+  iter_warmup = 400,
+  iter_sampling = 400
+)
+
+N <- nrow(test_mat)
+
+round(matrix(mod_out$summary("R_out")$mean,nrow(test_mat), nrow(test_mat)), 3)
+
+
+chol(test_mat)
+
+mod_out <- mod$optimize(
+  data = list(N = nrow(test_mat),
+                      R = test_mat,
+                      is_symmetric = 1,
+                      Z = sum(where_zero),
+                      K = ( nrow(test_mat) * (nrow(test_mat) - 1 )) / 2,
+                      where_zero = where_zero ))
+
+round(matrix(mod_out$summary("R_out")$estimate,nrow(test_mat), nrow(test_mat)), 3)
+
+mod_out <- mod$sample(
+  data = list(N = nrow(test_mat),
+              R = test_mat,
+              is_symmetric = 1),
+  chains = 2,
+  seed = 23421,
+  adapt_delta = 0.9,
+  max_treedepth = 10,
+  parallel_chains = 2,
+  iter_warmup = 200,
+  iter_sampling = 200
+)
+mod_out$summary("R_out")
+round(matrix(mod_out$summary("R_out")$mean, N, N), 3)
diff --git a/examples/linear_algebra/stan/correlation_angles_constrain.stan b/examples/linear_algebra/stan/correlation_angles_constrain.stan
new file mode 100644
index 00000000..ad3d76b1
--- /dev/null
+++ b/examples/linear_algebra/stan/correlation_angles_constrain.stan
@@ -0,0 +1,100 @@
+functions {
+#include correlation_angles.stan
+#include triangular.stan
+ vector lower_elements(matrix M){
+    int n = rows(M);
+    int k = n * (n - 1) / 2;
+    int counter = 1;
+    vector[k] out;
+
+    for (i in 2:n){
+      for (j in 1:(i - 1)) {
+        out[counter] = M[i, j];
+        counter += 1;
+      }
+    }
+    return out;
+  }
+
+   vector lower_elements_constrain(matrix M, int A){
+    int n = rows(M);
+    int counter = 1;
+    vector[A] out;
+
+    for (i in 2:n){
+      for (j in 1:(i - 1)) {
+        if(M[i, j] > 0){
+         out[counter] = M[i, j];
+         counter += 1;
+        }
+      }
+    }
+    return out;
+  }
+
+ matrix build_angle_mat (vector where_zero, vector angle_raw, int K) {
+  int N = num_elements(where_zero);
+  matrix[K, K] mat;
+  int count = 1;
+  int raw_count = 0;
+
+  mat[1, 1] = 0.0;
+  for (k in 2:K){
+    mat[k, k] = 0.0;
+    for (j in 1:k - 1) {
+      raw_count += 1;
+      if (where_zero[raw_count] != 1) {
+         mat[k, j] = angle_raw[count];
+         count += 1;
+      }
+    }
+  }
+
+  return mat;
+}
+
+ // set the upper triangle to 0
+ // only looking at strictly lower tri part
+ vector sparse_cholesky_lp (vector angle_raw, int[] csr_rows, int[] csr_cols, int Z, int N){
+  vector[Z + N] sparse_chol; // Z + N is the number of non-zero values in lower tri plus
+                             // the diagonal since the angles do not include the diagonal
+  int R = size(csr_rows);
+  int C = size(csr_cols);
+  int S[R, C] = append_array(csr_rows, csr_cols);
+  int count = 1;
+
+  sparse_chol[count] = 1;
+
+  // traversing in col-major order
+  // S[2, 1]           skips S[2, 2]
+  // S[3, 1], S[3, 2]  skips S[3, 3]
+  // etc
+  for (r in S) {
+    int this_rows_column_num = 0; 
+    for (c in r) {
+      count += 1;
+      this_rows_column_num += 1;
+
+      if(this_rows_column_num == 1)
+        sparse_chol[count] = cos(angle_raw[count]); // constrain first column
+
+
+    }
+
+    if (i > 2) {
+      for (j in 2:(i - 1)) {
+        real prod_sines = prod(sin(angle[i, 1:j - 1]));
+        real cos_theta;
+        if (is_nan(angle_raw[i, j]) == 1) {
+          cos_theta = -(dot_product(inv_mat[j, 1:j - 1], inv_mat[i, 1:j - 1]) ) / ( inv_mat[j, j] * prod_sines );
+            if ( cos_theta < -1 || cos_theta > 1 ) reject("cos_theta is ", cos_theta, " and must be in [-1, 1]"); // cos_theta = 0;
+           // else if( cos_theta > 1) cos_theta = 1;
+            angle[i, j] = acos( cos_theta );
+          }
+        inv_mat[i, j] = cos(angle[i, j]) * prod_sines;
+      }
+    }
+    inv_mat[i, i] = prod(sin(angle[i, 1:(i - 1)]));
+  }
+  return inv_mat;
+}
\ No newline at end of file
diff --git a/examples/linear_algebra/stan/correlation_angles_constrain.stanfunctions b/examples/linear_algebra/stan/correlation_angles_constrain.stanfunctions
new file mode 100644
index 00000000..ad3d76b1
--- /dev/null
+++ b/examples/linear_algebra/stan/correlation_angles_constrain.stanfunctions
@@ -0,0 +1,100 @@
+functions {
+#include correlation_angles.stan
+#include triangular.stan
+ vector lower_elements(matrix M){
+    int n = rows(M);
+    int k = n * (n - 1) / 2;
+    int counter = 1;
+    vector[k] out;
+
+    for (i in 2:n){
+      for (j in 1:(i - 1)) {
+        out[counter] = M[i, j];
+        counter += 1;
+      }
+    }
+    return out;
+  }
+
+   vector lower_elements_constrain(matrix M, int A){
+    int n = rows(M);
+    int counter = 1;
+    vector[A] out;
+
+    for (i in 2:n){
+      for (j in 1:(i - 1)) {
+        if(M[i, j] > 0){
+         out[counter] = M[i, j];
+         counter += 1;
+        }
+      }
+    }
+    return out;
+  }
+
+ matrix build_angle_mat (vector where_zero, vector angle_raw, int K) {
+  int N = num_elements(where_zero);
+  matrix[K, K] mat;
+  int count = 1;
+  int raw_count = 0;
+
+  mat[1, 1] = 0.0;
+  for (k in 2:K){
+    mat[k, k] = 0.0;
+    for (j in 1:k - 1) {
+      raw_count += 1;
+      if (where_zero[raw_count] != 1) {
+         mat[k, j] = angle_raw[count];
+         count += 1;
+      }
+    }
+  }
+
+  return mat;
+}
+
+ // set the upper triangle to 0
+ // only looking at strictly lower tri part
+ vector sparse_cholesky_lp (vector angle_raw, int[] csr_rows, int[] csr_cols, int Z, int N){
+  vector[Z + N] sparse_chol; // Z + N is the number of non-zero values in lower tri plus
+                             // the diagonal since the angles do not include the diagonal
+  int R = size(csr_rows);
+  int C = size(csr_cols);
+  int S[R, C] = append_array(csr_rows, csr_cols);
+  int count = 1;
+
+  sparse_chol[count] = 1;
+
+  // traversing in col-major order
+  // S[2, 1]           skips S[2, 2]
+  // S[3, 1], S[3, 2]  skips S[3, 3]
+  // etc
+  for (r in S) {
+    int this_rows_column_num = 0; 
+    for (c in r) {
+      count += 1;
+      this_rows_column_num += 1;
+
+      if(this_rows_column_num == 1)
+        sparse_chol[count] = cos(angle_raw[count]); // constrain first column
+
+
+    }
+
+    if (i > 2) {
+      for (j in 2:(i - 1)) {
+        real prod_sines = prod(sin(angle[i, 1:j - 1]));
+        real cos_theta;
+        if (is_nan(angle_raw[i, j]) == 1) {
+          cos_theta = -(dot_product(inv_mat[j, 1:j - 1], inv_mat[i, 1:j - 1]) ) / ( inv_mat[j, j] * prod_sines );
+            if ( cos_theta < -1 || cos_theta > 1 ) reject("cos_theta is ", cos_theta, " and must be in [-1, 1]"); // cos_theta = 0;
+           // else if( cos_theta > 1) cos_theta = 1;
+            angle[i, j] = acos( cos_theta );
+          }
+        inv_mat[i, j] = cos(angle[i, j]) * prod_sines;
+      }
+    }
+    inv_mat[i, i] = prod(sin(angle[i, 1:(i - 1)]));
+  }
+  return inv_mat;
+}
\ No newline at end of file
diff --git a/examples/linear_algebra/stan/correlation_angles_constrain_example.stan b/examples/linear_algebra/stan/correlation_angles_constrain_example.stan
new file mode 100644
index 00000000..e9627e1d
--- /dev/null
+++ b/examples/linear_algebra/stan/correlation_angles_constrain_example.stan
@@ -0,0 +1,88 @@
+functions {
+ #include correlation_angles.stan
+ #include triangular.stan
+ vector lower_elements(matrix M, int tri_size){
+    int n = rows(M);
+    int counter = 1;
+    vector[tri_size] out;
+
+    for (i in 2:n){
+      for (j in 1:(i - 1)) {
+        out[counter] = M[i, j];
+        counter += 1;
+      }
+    }
+    return out;
+  }
+
+   vector lower_elements_constrain(matrix M, int A){
+    int n = rows(M);
+    int counter = 1;
+    vector[A] out;
+
+    for (i in 2:n){
+      for (j in 1:(i - 1)) {
+        if(M[i, j] > 0){
+         out[counter] = M[i, j];
+         counter += 1;
+        }
+      }
+    }
+    return out;
+  }
+
+ matrix build_angle_mat (vector where_zero, vector angle_raw, int K) {
+  int N = num_elements(where_zero);
+  matrix[K, K] mat;
+  int count = 1;
+  int raw_count = 0;
+
+  mat[1, 1] = 0.0;
+  for (k in 2:K){
+    mat[k, k] = 0.0;
+    for (j in 1:k - 1) {
+      raw_count += 1;
+      if (where_zero[raw_count] != 1) {
+         mat[k, j] = angle_raw[count];
+         count += 1;
+      }
+    }
+  }
+
+  return mat;
+}
+
+ matrix zero_constrain (matrix angle_raw, int N){
+  matrix[N, N] inv_mat;
+  matrix[N, N] angle = angle_raw;
+
+  inv_mat[1, 1] = 1;
+
+  for (i in 2:N) {
+    // constrain first column
+    // if C = BB^T then for the first column
+    // c_{i, 1} = b_{i, 1} = 0 so 
+    // if c_{i, 1} = 0 then cos(\theta_{i, 1}) = 0
+    // acos(0) = \theta_{i, 1} = pi / 2
+     if (is_nan(angle_raw[i, 1]) == 1) {
+       inv_mat[i, 1] = 0;
+       angle[i, 1] = pi() / 2;
+     } else inv_mat[i, 1] = cos(angle[i, 1]);
+
+    if (i > 2) {
+      for (j in 2:(i - 1)) {
+        real prod_sines = prod(sin(angle[i, 1:j - 1]));
+        real cos_theta;
+        if (is_nan(angle_raw[i, j]) == 1) {
+          cos_theta = -(dot_product(inv_mat[j, 1:j - 1], inv_mat[i, 1:j - 1]) ) / ( inv_mat[j, j] * prod_sines );
+            if ( cos_theta < -1 || cos_theta > 1 ) reject("cos_theta is ", cos_theta, " and must be in [-1, 1]"); // cos_theta = 0;
+           // else if( cos_theta > 1) cos_theta = 1;
+            angle[i, j] = acos( cos_theta );
+          }
+        inv_mat[i, j] = cos(angle[i, j]) * prod_sines;
+      }
+    }
+    inv_mat[i, i] = prod(sin(angle[i, 1:(i - 1)]));
+  }
+  return inv_mat;
+ }
\ No newline at end of file
diff --git a/examples/linear_algebra/stan/correlation_angles_constrain_working.stan b/examples/linear_algebra/stan/correlation_angles_constrain_working.stan
new file mode 100644
index 00000000..0a9e9acc
--- /dev/null
+++ b/examples/linear_algebra/stan/correlation_angles_constrain_working.stan
@@ -0,0 +1,100 @@
+functions {
+#include correlation_angles.stan
+#include triangular.stan
+ vector lower_elements(matrix M){
+    int n = rows(M);
+    int k = n * (n - 1) / 2;
+    int counter = 1;
+    vector[k] out;
+
+    for (i in 2:n){
+      for (j in 1:(i - 1)) {
+        out[counter] = M[i, j];
+        counter += 1;
+      }
+    }
+    return out;
+  }
+
+   vector lower_elements_constrain(matrix M, int A){
+    int n = rows(M);
+    int counter = 1;
+    vector[A] out;
+
+    for (i in 2:n){
+      for (j in 1:(i - 1)) {
+        if(M[i, j] > 0){
+         out[counter] = M[i, j];
+         counter += 1;
+        }
+      }
+    }
+    return out;
+  }
+
+ matrix build_angle_mat (vector where_zero, vector angle_raw, int K) {
+  int N = num_elements(where_zero);
+  matrix[K, K] mat;
+  int count = 1;
+  int raw_count = 0;
+
+  mat[1, 1] = 0.0;
+  for (k in 2:K){
+    mat[k, k] = 0.0;
+    for (j in 1:k - 1) {
+      raw_count += 1;
+      if (where_zero[raw_count] != 1) {
+         mat[k, j] = angle_raw[count];
+         count += 1;
+      }
+    }
+  }
+
+  return mat;
+}
+
+ // set the upper triangle to 0
+ // only looking at strictly lower tri part
+ vector sparse_cholesky_lp (vector angle_raw, int[] csr_rows, int[] csr_cols, int Z, int N){
+  vector[Z + N] sparse_chol; // Z + N is the number of non-zero values in lower tri plus
+                             // the diagonal since the angles do not include the diagonal
+  int R = size(csr_rows);
+  int C = size(csr_cols);
+  int S[R, C] = append_array(csr_rows, csr_cols);
+  int count = 1;
+
+  sparse_chol[count] = 1;
+
+  // traversing in col-major order
+  // S[2, 1]           skips S[2, 2]
+  // S[3, 1], S[3, 2]  skips S[3, 3]
+  // etc
+  for (r in S) {
+    int this_rows_column_num = 0; 
+    for (c in r) {
+      count += 1;
+      this_rows_column_num += 1;
+
+      if(this_rows_column_num == 1)
+        sparse_chol[count] = cos(angle_raw[count]); // constrain first column
+
+
+    }
+
+    if (i > 2) {
+      for (j in 2:(i - 1)) {
+        real prod_sines = prod(sin(angle[i, 1:j - 1]));
+        real cos_theta;
+        if (is_nan(angle_raw[i, j]) == 1) {
+          cos_theta = -(dot_product(inv_mat[j, 1:j - 1], inv_mat[i, 1:j - 1]) ) / ( inv_mat[j, j] * prod_sines );
+            if ( cos_theta < -1 || cos_theta > 1 ) reject("cos_theta is ", cos_theta, " and must be in [-1, 1]"); // cos_theta = 0;
+           // else if( cos_theta > 1) cos_theta = 1;
+            angle[i, j] = acos( cos_theta );
+          }
+        inv_mat[i, j] = cos(angle[i, j]) * prod_sines;
+      }
+    }
+    inv_mat[i, i] = prod(sin(angle[i, 1:(i - 1)]));
+  }
+  return inv_mat;
+ }
\ No newline at end of file
diff --git a/examples/linear_algebra/stan/correlation_angles_example.stan b/examples/linear_algebra/stan/correlation_angles_example.stan
index cab8d43b..b69ffc46 100644
--- a/examples/linear_algebra/stan/correlation_angles_example.stan
+++ b/examples/linear_algebra/stan/correlation_angles_example.stan
@@ -5,20 +5,35 @@ functions {
 data {
   int N;
   matrix[N, N] R;
+  int<lower=0, upper=1> is_symmetric;
 }
 transformed data {
   int K = N * (N - 1) / 2;
+  row_vector[N] m;
+  
+  for (n in 1:N)
+    m[n] = N - n + 1;
 }
 parameters {
-  vector<lower = -pi(), upper = pi()>[K] theta;
-}
+  vector<lower = 0, upper = pi()>[K - (N - 2)] theta_free;
+} 
 transformed parameters {
-  matrix[N, N] R_hat = multiply_lower_tri_self_transpose(angle2chol(angle_vec2angle_mat(theta, N)));
+  vector[K] theta = append_row(rep_vector(0.5 * pi(), N - 2), theta_free);
 }
 model {
-  vector[N] R_vec = to_vector(R);
-  vector[N] R_hat_vec = to_vector(R_hat);
-  
-  squared_distance(R_vec, R_hat_vec) ~ chi_square(1);
-  target += sum(log(fabs(2 * (R_vec - R_hat_vec))));
+   matrix[N, N] L = angle2chol(angle_vec2angle_mat(theta, N));
+   matrix[N, N] R_hat = multiply_lower_tri_self_transpose(L);
+   
+  // log_absd_jacs 
+  // sin(theta) is always > 0 since theta in (0, pi)
+  // because of this the diagonal of L is always > 0 also
+   target += 2 * sum(log(sin(theta)));            // angle2chol
+   target += N * log(2) + m * diagonal(log(L));   // cholesky tcrossprod
+   
+   if (is_symmetric == 1)
+    lower_elements(R, K) ~ normal(lower_elements(R_hat, K), 0.001);
+    else to_vector(R) ~ normal(to_vector(R_hat), 0.001);
+}
+generated quantities {
+   matrix[N, N] R_out = multiply_lower_tri_self_transpose(angle2chol(angle_vec2angle_mat(theta, N)));
 }
diff --git a/examples/linear_algebra/stan/correlation_angles_first_col_constrain_example.stan b/examples/linear_algebra/stan/correlation_angles_first_col_constrain_example.stan
new file mode 100644
index 00000000..91a67aeb
--- /dev/null
+++ b/examples/linear_algebra/stan/correlation_angles_first_col_constrain_example.stan
@@ -0,0 +1,39 @@
+functions {
+ #include correlation_angles.stan
+ #include triangular.stan
+}
+data {
+  int N;
+  matrix[N, N] R;
+  int<lower=0, upper=1> is_symmetric;
+}
+transformed data {
+  int K = N * (N - 1) / 2;
+  row_vector[N] m;
+  
+  for (n in 1:N)
+    m[n] = N - n + 1;
+}
+parameters {
+  vector<lower = 0, upper = pi()>[K - (N - 2)] theta_free;
+} 
+transformed parameters {
+  vector[K] theta = append_row(rep_vector(0.5 * pi(), N - 2), theta_free);
+}
+model {
+   matrix[N, N] L = angle2chol(angle_vec2angle_mat(theta, N));
+   matrix[N, N] R_hat = multiply_lower_tri_self_transpose(L);
+   
+  // log_absd_jacs 
+  // sin(theta) is always > 0 since theta in (0, pi)
+  // because of this the diagonal of L is always > 0 also
+   target += 2 * sum(log(sin(theta)));            // angle2chol
+   target += N * log(2) + m * diagonal(log(L));   // cholesky tcrossprod
+   
+   if (is_symmetric == 1)
+    lower_elements(R, K) ~ normal(lower_elements(R_hat, K), 0.001);
+    else to_vector(R) ~ normal(to_vector(R_hat), 0.001);
+}
+generated quantities {
+   matrix[N, N] R_out = multiply_lower_tri_self_transpose(angle2chol(angle_vec2angle_mat(theta, N)));
+}
\ No newline at end of file