diff --git a/src/lscm.cpp b/src/lscm.cpp
index d6e1b6c..1d11d9c 100644
--- a/src/lscm.cpp
+++ b/src/lscm.cpp
@@ -1,10 +1,57 @@
 #include "lscm.h"
+#include "vector_area_matrix.h"
+#include "igl/repdiag.h"
+#include "igl/cotmatrix.h"
+#include "igl/massmatrix.h"
+#include "igl/eigs.h"
+#include <Eigen/SVD>
 
 void lscm(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
   Eigen::MatrixXd & U)
-{
-  // Replace with your code
-  U = V.leftCols(2);
+{ 
+  // Construct matrix Q = (L 0; 0 L) - A where L is the laplacian and A is the vector area
+  Eigen::SparseMatrix<double> L;
+  igl::cotmatrix(V, F, L);
+  Eigen::SparseMatrix<double> L_d;
+  igl::repdiag(L, 2, L_d);
+  Eigen::SparseMatrix<double> A;
+  vector_area_matrix(F, A);
+  Eigen::SparseMatrix<double> Q = L_d  - 2.0 * A;
+
+  // Construct the matrix B = (M 0; 0 M) where M is the mass matrix
+  Eigen::SparseMatrix<double> M;
+  igl::massmatrix(V, F, igl::MassMatrixType::MASSMATRIX_TYPE_BARYCENTRIC, M);
+  Eigen::SparseMatrix<double> B;
+  igl::repdiag(M, 2, B);
+
+  // Find parameterization coordinates by a eigenvector of Qv = lambda * B * v with v non-zero eigenvalue
+  Eigen::VectorXd eigen_vals;
+  Eigen::MatrixXd eigen_vecs;
+  igl::eigs(Q, B, 5, igl::EigsType::EIGS_TYPE_SM, eigen_vecs, eigen_vals);
+  
+  int idx = 1; // the second eigenvector calculated should be the Fielder vector
+   double tol = 1e-8; // Pick an eigenvector sufficiently away from 0
+  // hope that an eigen_vec with non-zero eigen_val has been found among the calculated evals
+  // otherwise just take the largest evec with the largest evec. 
+  while (abs(eigen_vals[idx]) <= tol) {
+	idx++;
+  }
+  if (idx >= 5) {
+	idx = 4;
+  }
+
+  // Reshape eigenvector into U
+  int n_v = F.maxCoeff() + 1;
+  U.resize(n_v, 2);
+  U.col(0) = eigen_vecs.col(idx).head(n_v);
+  U.col(1) = eigen_vecs.col(idx).tail(n_v);
+  
+  // Do SVD on coordinates U to find a rotation on the coordinates U
+  Eigen::Matrix2d covar = U.transpose() * U;
+  Eigen::JacobiSVD<Eigen::Matrix2d> svd(covar, Eigen::ComputeFullU);
+  for (int i = 0; i < U.rows(); i++) {
+	U.row(i) = U.row(i) * svd.matrixU(); 
+  }
 }
diff --git a/src/tutte.cpp b/src/tutte.cpp
index 6f6109b..371beb5 100644
--- a/src/tutte.cpp
+++ b/src/tutte.cpp
@@ -1,11 +1,55 @@
 #include "tutte.h"
+#include "igl/boundary_loop.h"
+#include "igl/map_vertices_to_circle.h"
+#include "igl/edges.h"
+#include "igl/min_quad_with_fixed.h"
 
 void tutte(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
   Eigen::MatrixXd & U)
 {
-  // Replace with your code
-  U = V.leftCols(2);
+  // Determine the longest boundary loop and map its vertices onto the unit circle
+  Eigen::VectorXi bdry;
+  igl::boundary_loop(F, bdry);
+  Eigen::MatrixXd bdry_pos;
+  igl::map_vertices_to_circle(V, bdry, bdry_pos);
+
+  // Calculate the weighted Laplacian matrix
+  Eigen::MatrixXi E;
+  igl::edges(F, E);
+  Eigen::SparseMatrix<double> L;
+  typedef Eigen::Triplet<double> T;
+  std::vector<T> triples;
+  triples.reserve(4 * E.rows());
+  for (int i = 0; i < E.rows(); i++) {
+	double weight = 1.0 / ((V.row(E(i,0)) - V.row(E(i,1))).norm());
+        triples.push_back(T(E(i,0), E(i,1), weight));
+        triples.push_back(T(E(i,1), E(i,0), weight));
+	triples.push_back(T(E(i,0), E(i,0),-weight));
+	triples.push_back(T(E(i,1), E(i,1),-weight));
+  }  	  
+  int n_v = F.maxCoeff() + 1;
+  L.resize(n_v, n_v);
+  L.setFromTriplets(triples.begin(), triples.end());
+
+  // Now solve the weighted optimization problem  U_x^T L U_x + U_y^T L U_y subject to boundary constraints calculated earlier
+  igl::min_quad_with_fixed_data<double> data;
+
+  // No linear constraints or linear term in optimization objective;
+  Eigen::SparseMatrix<double> Z;
+  Z.setZero();
+  Eigen::VectorXd B;
+  B.setZero();
+  Eigen::VectorXd z(n_v);
+  z.setZero();
+
+  igl::min_quad_with_fixed_precompute(L, bdry, Z, false, data); // no linear equality constraints
+  U.resize(n_v, 2);
+  for (int i = 0; i <= 1; i++) {
+	Eigen::VectorXd cur_col;
+	igl::min_quad_with_fixed_solve(data, z, bdry_pos.col(i), B, cur_col); // no linear term in objective or linear inequality contraints
+  	U.col(i) = cur_col;
+   }
 }
 
diff --git a/src/vector_area_matrix.cpp b/src/vector_area_matrix.cpp
index 1d60ef8..2195b99 100644
--- a/src/vector_area_matrix.cpp
+++ b/src/vector_area_matrix.cpp
@@ -1,11 +1,38 @@
 #include "vector_area_matrix.h"
+#include "igl/boundary_loop.h"
 
 void vector_area_matrix(
   const Eigen::MatrixXi & F,
   Eigen::SparseMatrix<double>& A)
 {
-  // Replace with your code
   int V_size = F.maxCoeff()+1;
-  A.resize(V_size*2,V_size*2);
+  Eigen::SparseMatrix<double> A_u;
+  A_u.resize(V_size*2,V_size*2);
+
+  // Determine the longest boundary loop
+  std::vector<std::vector<int>> bdry_loops;
+  igl::boundary_loop(F, bdry_loops);
+  
+  // Each edge (i, j) in the loop contributes 1/2 (x_i y_j - x_j y_i) to the vector area
+  typedef Eigen::Triplet<double> T;
+  std::vector<T> triples;
+  triples.reserve(bdry_loops[0].size() * 4);
+  for (int l = 0; l < bdry_loops.size() ; l++) {
+  std::vector<int> bdry = bdry_loops[l];
+  for (int i = 0; i < bdry.size(); i++) {
+	if (i == bdry.size() - 1) {
+		triples.push_back(T(bdry[i], bdry[0] + V_size, 0.5));
+   		triples.push_back(T(bdry[i] + V_size, bdry[0], -0.5));
+	} 
+	else {
+		triples.push_back(T(bdry[i], bdry[i+1] + V_size, 0.5));
+   		triples.push_back(T(bdry[i] + V_size, bdry[i+1], -0.5));
+	}
+  }
+}
+  A_u.setFromTriplets(triples.begin(), triples.end());
+  Eigen::SparseMatrix<double> B;
+  B = A_u.transpose();
+  A = (A_u + B) / 2;
 }