diff --git a/include/pmmc.h b/include/pmmc.h
index 269561a3..10d9772b 100644
--- a/include/pmmc.h
+++ b/include/pmmc.h
@@ -3511,7 +3511,7 @@ inline double pmmc_CubeSurfaceInterpValue(DoubleArray &CubeValues, DTMutableList
 }
 //--------------------------------------------------------------------------------------------------------
 inline double pmmc_CubeCurveInterpValue(DoubleArray &CubeValues, DoubleArray &CurveValues, 
-										DTMutableList<Point> &Points, int npts)
+										DTMutableList<Point> &Points, int i, int j, int k, int npts)
 {
 	int p;
 	Point A,B;
@@ -3535,7 +3535,10 @@ inline double pmmc_CubeCurveInterpValue(DoubleArray &CubeValues, DoubleArray &Cu
 
 	for (p=0; p<npts; p++){
 		A = Points(p);
-		CurveValues(p) = a + b*A.x + c*A.y+d*A.z + e*A.x*A.y + f*A.x*A.z + g*A.y*A.z + h*A.x*A.y*A.z;
+		x = A.x-1.0*i;
+		y = A.y-1.0*j;
+		z = A.z-1.0*k;
+		CurveValues(p) = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
 	}
 	
 	integral = 0.0;
@@ -3543,9 +3546,11 @@ inline double pmmc_CubeCurveInterpValue(DoubleArray &CubeValues, DoubleArray &Cu
 		// Extract the line segment
 		A = Points(p);
 		B = Points(p+1);
+		vA = CurveValues(p);
+		vB = CurveValues(p+1);
 		// Compute the length of the segment
 		length = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
-		integral += 0.5*length*(CurveValues(p) + CurveValues(p+1));
+		integral += 0.5*length*(vA + vB);
 	}
 	return integral;
 }
@@ -3558,6 +3563,7 @@ inline double pmmc_CubeContactAngle(DoubleArray &CubeValues, DoubleArray &CurveV
 	int p;
 	Point A,B;
 	double vA,vB;
+	double x,y,z;
 	double s,s1,s2,s3,temp;
 	double a,b,c,d,e,f,g,h;
 	double integral;
@@ -3600,9 +3606,12 @@ inline double pmmc_CubeContactAngle(DoubleArray &CubeValues, DoubleArray &CurveV
 	g = CubeValues(0,1,1)-a-c-d;
 	h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
 
-	for (int p=0; p<npts; p++){
+	for (p=0; p<npts; p++){
 		A = Points(p);
-		CurveValues(p) = a + b*A.x + c*A.y+d*A.z + e*A.x*A.y + f*A.x*A.z + g*A.y*A.z + h*A.x*A.y*A.z;
+		x = A.x-1.0*i;
+		y = A.y-1.0*j;
+		z = A.z-1.0*k;
+		CurveValues(p) = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
 	}
 	
 	integral = 0.0;
@@ -3610,10 +3619,187 @@ inline double pmmc_CubeContactAngle(DoubleArray &CubeValues, DoubleArray &CurveV
 		// Extract the line segment
 		A = Points(p);
 		B = Points(p+1);
+		vA = CurveValues(p);
+		vB = CurveValues(p+1);
 		// Compute the length of the segment
 		length = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
-		integral += 0.5*length*(CurveValues(p) + CurveValues(p+1));
+		integral += 0.5*length*(vA + vB);
 	}
-	
 	return integral;
-}//--------------------------------------------------------------------------------------------------------
+}
+//--------------------------------------------------------------------------------------------------------
+inline void pmmc_CubeSurfaceInterpVector(DoubleArray &Vec_x, DoubleArray &Vec_y, DoubleArray &Vec_z,
+									DoubleArray &CubeValues, DTMutableList<Point> &Points, IntArray &Triangles,
+									  DoubleArray &SurfaceVector, DoubleArray &VecAvg, int i, int j, int k, int npts, int ntris)
+{
+	Point A,B,C;
+	int p;
+	double vA,vB,vC;
+	double x,y,z;
+	double s,s1,s2,s3,temp;
+	double a,b,c,d,e,f,g,h;
+	double integral;
+
+	// ................x component .............................
+	// Copy the curvature values for the cube
+	CubeValues(0,0,0) = Vec_x(i,j,k);
+	CubeValues(1,0,0) = Vec_x(i+1,j,k);
+	CubeValues(0,1,0) = Vec_x(i,j+1,k);
+	CubeValues(1,1,0) = Vec_x(i+1,j+1,k);
+	CubeValues(0,0,1) = Vec_x(i,j,k+1);
+	CubeValues(1,0,1) = Vec_x(i+1,j,k+1);
+	CubeValues(0,1,1) = Vec_x(i,j+1,k+1);
+	CubeValues(1,1,1) = Vec_x(i+1,j+1,k+1);
+
+	// trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
+	a = CubeValues(0,0,0);
+	b = CubeValues(1,0,0)-a;
+	c = CubeValues(0,1,0)-a;
+	d = CubeValues(0,0,1)-a;
+	e = CubeValues(1,1,0)-a-b-c;
+	f = CubeValues(1,0,1)-a-b-d;
+	g = CubeValues(0,1,1)-a-c-d;
+	h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
+
+	for (p=0; p<npts; p++){
+		A = Points(p);
+		x = A.x-1.0*i;
+		y = A.y-1.0*j;
+		z = A.z-1.0*k;
+		SurfaceVector(p) = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
+	}
+
+	// ................y component .............................
+	// Copy the curvature values for the cube
+	CubeValues(0,0,0) = Vec_y(i,j,k);
+	CubeValues(1,0,0) = Vec_y(i+1,j,k);
+	CubeValues(0,1,0) = Vec_y(i,j+1,k);
+	CubeValues(1,1,0) = Vec_y(i+1,j+1,k);
+	CubeValues(0,0,1) = Vec_y(i,j,k+1);
+	CubeValues(1,0,1) = Vec_y(i+1,j,k+1);
+	CubeValues(0,1,1) = Vec_y(i,j+1,k+1);
+	CubeValues(1,1,1) = Vec_y(i+1,j+1,k+1);
+
+	// trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
+	a = CubeValues(0,0,0);
+	b = CubeValues(1,0,0)-a;
+	c = CubeValues(0,1,0)-a;
+	d = CubeValues(0,0,1)-a;
+	e = CubeValues(1,1,0)-a-b-c;
+	f = CubeValues(1,0,1)-a-b-d;
+	g = CubeValues(0,1,1)-a-c-d;
+	h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
+
+	for (p=0; p<npts; p++){
+		A = Points(p);
+		x = A.x-1.0*i;
+		y = A.y-1.0*j;
+		z = A.z-1.0*k;
+		SurfaceVector(npts+p) = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
+	}
+
+	// ................z component .............................
+	// Copy the curvature values for the cube
+	CubeValues(0,0,0) = Vec_z(i,j,k);
+	CubeValues(1,0,0) = Vec_z(i+1,j,k);
+	CubeValues(0,1,0) = Vec_z(i,j+1,k);
+	CubeValues(1,1,0) = Vec_z(i+1,j+1,k);
+	CubeValues(0,0,1) = Vec_z(i,j,k+1);
+	CubeValues(1,0,1) = Vec_z(i+1,j,k+1);
+	CubeValues(0,1,1) = Vec_z(i,j+1,k+1);
+	CubeValues(1,1,1) = Vec_z(i+1,j+1,k+1);
+
+	// trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
+	a = CubeValues(0,0,0);
+	b = CubeValues(1,0,0)-a;
+	c = CubeValues(0,1,0)-a;
+	d = CubeValues(0,0,1)-a;
+	e = CubeValues(1,1,0)-a-b-c;
+	f = CubeValues(1,0,1)-a-b-d;
+	g = CubeValues(0,1,1)-a-c-d;
+	h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
+
+	for (p=0; p<npts; p++){
+		A = Points(p);
+		x = A.x-1.0*i;
+		y = A.y-1.0*j;
+		z = A.z-1.0*k;
+		SurfaceVector(2*npts+p) = a + b*x + c*y + d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
+	}
+	//.............................................................................
+
+	for (int r=0; r<ntris; r++){
+		A = Points(Triangles(0,r));
+		B = Points(Triangles(1,r));
+		C = Points(Triangles(2,r));
+		// Compute length of sides (assume dx=dy=dz)
+		s1 = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
+		s2 = sqrt((A.x-C.x)*(A.x-C.x)+(A.y-C.y)*(A.y-C.y)+(A.z-C.z)*(A.z-C.z));
+		s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
+		s = 0.5*(s1+s2+s3);
+		temp = s*(s-s1)*(s-s2)*(s-s3);
+		if (temp > 0.0){
+			// Increment the averaged values
+			// x component
+			vA = SurfaceVector(Triangles(0,r));
+			vB = SurfaceVector(Triangles(1,r));
+			vC = SurfaceVector(Triangles(2,r));
+			VecAvg(0) += sqrt(temp)*0.33333333333333333*(vA+vB+vC);
+			// y component
+			vA = SurfaceVector(npts+Triangles(0,r));
+			vB = SurfaceVector(npts+Triangles(1,r));
+			vC = SurfaceVector(npts+Triangles(2,r));
+			VecAvg(1) += sqrt(temp)*0.33333333333333333*(vA+vB+vC);
+			// z component
+			vA = SurfaceVector(2*npts+Triangles(0,r));
+			vB = SurfaceVector(2*npts+Triangles(1,r));
+			vC = SurfaceVector(2*npts+Triangles(2,r));
+			VecAvg(2) += sqrt(temp)*0.33333333333333333*(vA+vB+vC);
+		}
+	}
+}
+//--------------------------------------------------------------------------------------------------------
+inline double pmmc_CubeCurveInterpVector(DoubleArray &CubeValues, DoubleArray &CurveValues,
+										DTMutableList<Point> &Points, int i, int j, int k, int npts)
+{
+	int p;
+	Point A,B;
+	double vA,vB;
+	double x,y,z;
+	double s,s1,s2,s3,temp;
+	double a,b,c,d,e,f,g,h;
+	double integral;
+	double length;
+
+	// trilinear coefficients: f(x,y,z) = a+bx+cy+dz+exy+fxz+gyz+hxyz
+	// Evaluate the coefficients
+	a = CubeValues(0,0,0);
+	b = CubeValues(1,0,0)-a;
+	c = CubeValues(0,1,0)-a;
+	d = CubeValues(0,0,1)-a;
+	e = CubeValues(1,1,0)-a-b-c;
+	f = CubeValues(1,0,1)-a-b-d;
+	g = CubeValues(0,1,1)-a-c-d;
+	h = CubeValues(1,1,1)-a-b-c-d-e-f-g;
+
+	for (p=0; p<npts; p++){
+		A = Points(p);
+		x = A.x-1.0*i;
+		y = A.y-1.0*j;
+		z = A.z-1.0*k;
+		CurveValues(p) = a + b*x + c*y+d*z + e*x*y + f*x*z + g*y*z + h*x*y*z;
+	}
+
+	integral = 0.0;
+	for (p=0; p < npts-1; p++){
+		// Extract the line segment
+		A = Points(p);
+		B = Points(p+1);
+		vA = CurveValues(p);
+		vB = CurveValues(p+1);
+		// Compute the length of the segment
+		length = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
+		integral += 0.5*length*(vA + vB);
+	}
+	return integral;
+}
diff --git a/tests/TestContactAngle.cpp b/tests/TestContactAngle.cpp
new file mode 100644
index 00000000..b04d8c42
--- /dev/null
+++ b/tests/TestContactAngle.cpp
@@ -0,0 +1,283 @@
+#include <iostream>
+#include <math.h>
+#include "pmmc.h"
+//#include "PointList.h"
+//#include "Array.h"
+
+#define RADIUS 15
+#define HEIGHT 15.5
+#define N 60
+#define PI 3.14159
+
+int main (int argc, char *argv[])
+{
+	
+	//  printf("Radius = %s \n,"RADIUS);  
+	int SIZE = N*N*N;
+	int Nx,Ny,Nz;
+	Nx = Ny = Nz = N;
+	int i,j,k,p,q,r;
+	
+//	double *Solid; // cylinder
+//	double *Phase; // region of the cylinder	
+//	Solid = new double [SIZE];
+//	Phase = new double [SIZE];
+	DoubleArray SignDist(Nx,Ny,Nz);
+	DoubleArray Phase(Nx,Ny,Nz);
+	double fluid_isovalue = 0.0;
+	double solid_isovalue = 0.0;
+
+	
+	/* ****************************************************************
+	 VARIABLES FOR THE PMMC ALGORITHM
+	 ****************************************************************** */
+	//...........................................................................
+	// Averaging variables
+	//...........................................................................
+	double awn,ans,aws,lwns,nwp_volume;
+	double As;
+	double dEs,dAwn,dAns;			 // Global surface energy (calculated by rank=0)
+	double awn_global,ans_global,aws_global,lwns_global,nwp_volume_global;	
+	double As_global;
+	//	bool add=1;			// Set to false if any corners contain nw-phase ( F > fluid_isovalue)
+	int cube[8][3] = {{0,0,0},{1,0,0},{0,1,0},{1,1,0},{0,0,1},{1,0,1},{0,1,1},{1,1,1}};  // cube corners
+	//	int count_in=0,count_out=0;
+	//	int nodx,nody,nodz;
+	// initialize lists for vertices for surfaces, common line
+	DTMutableList<Point> nw_pts(20);
+	DTMutableList<Point> ns_pts(20);
+	DTMutableList<Point> ws_pts(20);
+	DTMutableList<Point> nws_pts(20);
+	// initialize triangle lists for surfaces
+	IntArray nw_tris(3,20);
+	IntArray ns_tris(3,20);
+	IntArray ws_tris(3,20);
+	// initialize list for line segments
+	IntArray nws_seg(2,20);
+	
+	DTMutableList<Point> tmp(20);
+	//	IntArray store;
+	
+	int n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0;
+	int n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;
+	
+	double s,s1,s2,s3;		// Triangle sides (lengths)
+	Point A,B,C,P;
+	//	double area;
+	
+	// Initialize arrays for local solid surface
+	DTMutableList<Point> local_sol_pts(20);
+	int n_local_sol_pts = 0;
+	IntArray local_sol_tris(3,18);
+	int n_local_sol_tris;
+	DoubleArray values(20);
+	DTMutableList<Point> local_nws_pts(20);
+	int n_local_nws_pts;
+	
+	int c;
+	//...........................................................................
+	int ncubes = (Nx-2)*(Ny-2)*(Nz-2);	// Exclude the "upper" halo
+	IntArray cubeList(3,ncubes);
+	pmmc_CubeListFromMesh(cubeList, ncubes, Nx, Ny, Nz);
+	//...........................................................................
+	double Cx,Cy,Cz;
+	double dist1,dist2;
+	Cx = Cy = Cz = N*0.51;
+	for (k=0; k<N; k++){
+		for (j=0; j<N; j++){
+			for (i=0; i<N; i++){
+				dist1 = sqrt((i-Cx)*(i-Cx)+(j-Cy)*(j-Cy)) - RADIUS;
+				dist2 = fabs(Cz-k)-HEIGHT;
+				//	printf("distances = %f, %f \n",dist1,dist2);
+				//Solid.data[k*Nx*Ny+j*Nx+i] = dist1;
+				//Phase[k*Nx*Ny+j*Nx+i] = dist2;
+				SignDist(i,j,k) = -dist1;
+				Phase(i,j,k) = dist2;
+			}
+		}   
+	}
+
+	FILE *STRIS;
+	STRIS = fopen("solid-triangles.out","w");
+	
+	FILE *WN_TRIS;
+	WN_TRIS = fopen("wn-tris.out","w");
+	
+	FILE *NS_TRIS;
+	NS_TRIS = fopen("ns-tris.out","w");
+	
+	FILE *WS_TRIS;
+	WS_TRIS = fopen("ws-tris.out","w");
+	
+	FILE *WNS_PTS;
+	WNS_PTS = fopen("wns-pts.out","w");
+	
+	// End of the loop to set the values
+	awn = aws = ans = lwns = 0.0;
+	nwp_volume = 0.0;
+	As = 0.0;
+	
+	for (c=0;c<ncubes;c++){
+		// Get cube from the list
+		i = cubeList(0,c);
+		j = cubeList(1,c);
+		k = cubeList(2,c);
+		
+		for (p=0;p<8;p++){
+			if ( Phase(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 
+				&&  SignDist(i+cube[p][0],j+cube[p][1],k+cube[p][2]) > 0 ){
+				nwp_volume += 0.125;
+			}
+		}
+		
+		// Run PMMC
+		n_local_sol_tris = 0;
+		n_local_sol_pts = 0;
+		n_local_nws_pts = 0;
+		
+		n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0;
+		n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;
+
+		// Construct the interfaces and common curve
+		pmmc_ConstructLocalCube(SignDist, Phase, solid_isovalue, fluid_isovalue,
+				nw_pts, nw_tris, values, ns_pts, ns_tris, ws_pts, ws_tris,
+				local_nws_pts, nws_pts, nws_seg, local_sol_pts, local_sol_tris,
+				n_local_sol_tris, n_local_sol_pts, n_nw_pts, n_nw_tris,
+				n_ws_pts, n_ws_tris, n_ns_tris, n_ns_pts, n_local_nws_pts, n_nws_pts, n_nws_seg,
+				i, j, k, Nx, Ny, Nz);
+
+		
+		//*******************************************************************
+		// Compute the Interfacial Areas, Common Line length for blob p
+		// nw surface
+		double temp;
+		for (r=0;r<n_nw_tris;r++){
+			A = nw_pts(nw_tris(0,r));
+			B = nw_pts(nw_tris(1,r));
+			C = nw_pts(nw_tris(2,r));
+			// Compute length of sides (assume dx=dy=dz)
+			s1 = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
+			s2 = sqrt((A.x-C.x)*(A.x-C.x)+(A.y-C.y)*(A.y-C.y)+(A.z-C.z)*(A.z-C.z));
+			s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
+			s = 0.5*(s1+s2+s3);
+			temp = s*(s-s1)*(s-s2)*(s-s3);
+			if (temp > 0.0) awn += sqrt(temp);
+			
+		}
+		for (r=0;r<n_ns_tris;r++){
+			A = ns_pts(ns_tris(0,r));
+			B = ns_pts(ns_tris(1,r));
+			C = ns_pts(ns_tris(2,r));
+			// Compute length of sides (assume dx=dy=dz)
+			s1 = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
+			s2 = sqrt((A.x-C.x)*(A.x-C.x)+(A.y-C.y)*(A.y-C.y)+(A.z-C.z)*(A.z-C.z));
+			s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
+			s = 0.5*(s1+s2+s3);
+			//ans=ans+sqrt(s*(s-s1)*(s-s2)*(s-s3));
+			temp = s*(s-s1)*(s-s2)*(s-s3);
+			if (temp > 0.0) ans += sqrt(temp);
+		}
+		for (r=0;r<n_ws_tris;r++){
+			A = ws_pts(ws_tris(0,r));
+			B = ws_pts(ws_tris(1,r));
+			C = ws_pts(ws_tris(2,r));
+			// Compute length of sides (assume dx=dy=dz)
+			s1 = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
+			s2 = sqrt((A.x-C.x)*(A.x-C.x)+(A.y-C.y)*(A.y-C.y)+(A.z-C.z)*(A.z-C.z));
+			s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
+			s = 0.5*(s1+s2+s3);
+			//aws=aws+sqrt(s*(s-s1)*(s-s2)*(s-s3));
+			temp = s*(s-s1)*(s-s2)*(s-s3);
+			if (temp > 0.0) aws += sqrt(temp);
+		}
+		for (r=0;r<n_local_sol_tris;r++){
+			A = local_sol_pts(local_sol_tris(0,r));
+			B = local_sol_pts(local_sol_tris(1,r));
+			C = local_sol_pts(local_sol_tris(2,r));
+			// Compute length of sides (assume dx=dy=dz)
+			s1 = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
+			s2 = sqrt((A.x-C.x)*(A.x-C.x)+(A.y-C.y)*(A.y-C.y)+(A.z-C.z)*(A.z-C.z));
+			s3 = sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
+			s = 0.5*(s1+s2+s3);
+			//aws=aws+sqrt(s*(s-s1)*(s-s2)*(s-s3));
+			temp = s*(s-s1)*(s-s2)*(s-s3);
+			if (temp > 0.0) As += sqrt(temp);
+		}
+		for (p=0; p < n_local_nws_pts-1; p++){
+			// Extract the line segment
+			A = local_nws_pts(p);
+			B = local_nws_pts(p+1);
+			// Compute the length of the segment
+			s = sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
+			// Add the length to the common line 
+			lwns += s;
+		}
+		//.......................................................................................
+		// Write the triangle lists to text file
+		for (r=0;r<n_nw_tris;r++){
+			A = nw_pts(nw_tris(0,r));
+			B = nw_pts(nw_tris(1,r));
+			C = nw_pts(nw_tris(2,r));
+			fprintf(WN_TRIS,"%f %f %f %f %f %f %f %f %f \n",A.x,A.y,A.z,B.x,B.y,B.z,C.x,C.y,C.z);
+		}		
+		for (r=0;r<n_ws_tris;r++){
+			A = ws_pts(ws_tris(0,r));
+			B = ws_pts(ws_tris(1,r));
+			C = ws_pts(ws_tris(2,r));
+			fprintf(WS_TRIS,"%f %f %f %f %f %f %f %f %f \n",A.x,A.y,A.z,B.x,B.y,B.z,C.x,C.y,C.z);
+		}
+		for (r=0;r<n_ns_tris;r++){
+			A = ns_pts(ns_tris(0,r));
+			B = ns_pts(ns_tris(1,r));
+			C = ns_pts(ns_tris(2,r));
+			fprintf(NS_TRIS,"%f %f %f %f %f %f %f %f %f \n",A.x,A.y,A.z,B.x,B.y,B.z,C.x,C.y,C.z);
+		}
+		for (p=0; p < n_nws_pts; p++){
+			P = nws_pts(p);
+			fprintf(WNS_PTS,"%f %f %f \n",P.x, P.y, P.z);
+		}
+		
+		//.......................................................................................
+		
+		for (r=0;r<n_local_sol_tris;r++){
+			A = local_sol_pts(local_sol_tris(0,r));
+			B = local_sol_pts(local_sol_tris(1,r));
+			C = local_sol_pts(local_sol_tris(2,r));
+			fprintf(STRIS,"%f %f %f %f %f %f %f %f %f \n",A.x,A.y,A.z,B.x,B.y,B.z,C.x,C.y,C.z);
+		}
+		//*******************************************************************
+		// Reset the triangle counts to zero
+		n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0;
+		n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;
+		//	n_ns_tris_beg = 0;//n_ns_tris;
+		//	n_ws_tris_beg = 0;//n_ws_tris;
+		//	n_nws_seg_beg = n_nws_seg;
+		//*******************************************************************
+	}
+	fclose(WN_TRIS);
+	fclose(NS_TRIS);
+	fclose(WS_TRIS);
+	fclose(WNS_PTS);
+	fclose(STRIS);
+
+	printf("-------------------------------- \n");
+	printf("NWP volume = %f \n", nwp_volume);
+	printf("Area wn = %f, Analytical = %f \n", awn,2*PI*RADIUS*RADIUS);
+	printf("Area ns = %f, Analytical = %f \n", ans, 2*PI*RADIUS*(N-2)-4*PI*RADIUS*HEIGHT);
+	printf("Area ws = %f, Analytical = %f \n", aws, 4*PI*RADIUS*HEIGHT);
+	printf("Area s = %f, Analytical = %f \n", As, 2*PI*RADIUS*(N-2));
+	printf("Length wns = %f, Analytical = %f \n", lwns, 4*PI*RADIUS);
+	printf("-------------------------------- \n");	
+	//.........................................................................
+	
+/*	FILE *PHASE;
+	PHASE = fopen("Phase.in","wb");
+	fwrite(Phase,8,SIZE,PHASE);
+	fclose(PHASE);
+	
+	FILE *SOLID;
+	SOLID = fopen("Distance.in","wb");
+	fwrite(Solid,8,SIZE,SOLID);
+	fclose(SOLID);
+*/
+}
diff --git a/tests/TestCylinderAreas.cpp b/tests/TestCylinderAreas.cpp
index b04d8c42..bfe31ff9 100644
--- a/tests/TestCylinderAreas.cpp
+++ b/tests/TestCylinderAreas.cpp
@@ -145,7 +145,6 @@ int main (int argc, char *argv[])
 				n_local_sol_tris, n_local_sol_pts, n_nw_pts, n_nw_tris,
 				n_ws_pts, n_ws_tris, n_ns_tris, n_ns_pts, n_local_nws_pts, n_nws_pts, n_nws_seg,
 				i, j, k, Nx, Ny, Nz);
-
 		
 		//*******************************************************************
 		// Compute the Interfacial Areas, Common Line length for blob p
diff --git a/tests/TestSphereCurvature.cpp b/tests/TestSphereCurvature.cpp
index 219d0dd3..a30a9aea 100644
--- a/tests/TestSphereCurvature.cpp
+++ b/tests/TestSphereCurvature.cpp
@@ -141,7 +141,7 @@ int main(int argc, char **argv)
 
 	printf("Mean Curvature Average =  %f, Analytical = %f \n", wn_curvature_sum/wn_area_sum, 2.0/rad[0]/101 );
 	
-	FILE *CURVATURE;
+/*	FILE *CURVATURE;
 	CURVATURE = fopen("Curvature.dat","wb");
 	fwrite(MeanCurvature.data,8,N,CURVATURE);
 	fclose(CURVATURE);
@@ -150,5 +150,5 @@ int main(int argc, char **argv)
 	DISTANCE = fopen("SignDist.dat","wb");
 	fwrite(Phase.data,8,N,DISTANCE);
 	fclose(DISTANCE);
-	
+*/
 }
diff --git a/tests/TestSphereOrientation.cpp b/tests/TestSphereOrientation.cpp
new file mode 100644
index 00000000..e0bc8cff
--- /dev/null
+++ b/tests/TestSphereOrientation.cpp
@@ -0,0 +1,145 @@
+#include <iostream>
+#include <math.h>
+#include "pmmc.h"
+#include "Domain.h"
+
+using namespace std;
+
+/*
+ *  Compare the measured and analytical curvature for a sphere
+ *
+ */
+
+int main(int argc, char **argv)
+{
+	int i,j,k;
+	int Nx,Ny,Nz,N;
+	double Lx,Ly,Lz;
+	double fluid_isovalue=0.0;
+	double solid_isovalue=0.0;
+
+	Lx = Ly = Lz = 1.0;
+	Nx = Ny = Nz = 102;
+	N = Nx*Ny*Nz;
+
+	//...........................................................................
+	// Set up the cube list
+	//...........................................................................
+	int ncubes = (Nx-2)*(Ny-2)*(Nz-2);	// Exclude the "upper" halo
+	IntArray cubeList(3,ncubes);
+	pmmc_CubeListFromMesh(cubeList, ncubes, Nx, Ny, Nz);
+	//...........................................................................
+
+	//****************************************************************************************
+	// Create the structures needed to carry out the PMMC algorithm
+	//****************************************************************************************
+	DTMutableList<Point> nw_pts(20);
+	DTMutableList<Point> ns_pts(20);
+	DTMutableList<Point> ws_pts(20);
+	DTMutableList<Point> nws_pts(20);
+	// initialize triangle lists for surfaces
+	IntArray nw_tris(3,20);
+	IntArray ns_tris(3,20);
+	IntArray ws_tris(3,20);
+	// initialize list for line segments
+	IntArray nws_seg(2,20);
+
+	DTMutableList<Point> tmp(20);
+
+	DoubleArray wn_curvature(20);
+
+	int n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0;
+	int n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;
+
+	// Initialize arrays for local solid surface
+	DTMutableList<Point> local_sol_pts(20);
+	int n_local_sol_pts = 0;
+	IntArray local_sol_tris(3,18);
+	int n_local_sol_tris;
+	DoubleArray values(20);
+	DTMutableList<Point> local_nws_pts(20);
+	int n_local_nws_pts;
+	//****************************************************************************************
+
+	int nspheres=1;
+	// Set up the signed distance function for a sphere
+	double *cx,*cy,*cz,*rad;
+	cx = new double[nspheres];
+	cy = new double[nspheres];
+	cz = new double[nspheres];
+	rad = new double[nspheres];
+
+	// Sphere in the middle of the domain
+	cx[0] = cy[0] = cz[0] = 0.5;
+	rad[0] = 0.3;
+
+	DoubleArray SignDist(Nx,Ny,Nz);
+	DoubleArray Phase(Nx,Ny,Nz);
+	DoubleArray Phase_x(Nx,Ny,Nz);
+	DoubleArray Phase_y(Nx,Ny,Nz);
+	DoubleArray Phase_z(Nx,Ny,Nz);
+	DoubleArray CubeValues(2,2,2);
+
+	// Compute the signed distance function
+	SignedDistance(Phase.data,nspheres,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,0,0,0,1,1,1);
+
+	for (k=0; k<Nz; k++){
+		for (j=0; j<Ny; j++){
+			for (i=0; i<Nx; i++){
+				SignDist(i,j,k) = 100.0;
+				Phase(i,j,k) = sqrt((1.0*i-0.5*Nx)*(1.0*i-0.5*Nx)+(1.0*j-0.5*Ny)*(1.0*j-0.5*Ny)+(1.0*k-0.5*Nz)*(1.0*k-0.5*Nz))-0.3*Nx;
+			}
+		}
+	}
+	SignedDistance(SignDist.data,0,cx,cy,cz,rad,Lx,Ly,Lz,Nx,Ny,Nz,0,0,0,1,1,1);
+
+	double Gxx_sum,Gyy_sum,Gzz_sum,Gxy_sum,Gxz_sum,Gyz_sum;
+	double wn_curvature_sum = 0.0;
+	double wn_area_sum = 0.0;
+
+	pmmc_MeshGradient(Phase, Phase_x, Phase_y, Phase_z, Nx, Ny, Nz);
+
+	for (int c=0;c<ncubes;c++){
+
+		// Get cube from the list
+		i = cubeList(0,c);
+		j = cubeList(1,c);
+		k = cubeList(2,c);
+
+		// Run PMMC
+		n_local_sol_tris = 0;
+		n_local_sol_pts = 0;
+		n_local_nws_pts = 0;
+
+		n_nw_pts=0,n_ns_pts=0,n_ws_pts=0,n_nws_pts=0, map=0;
+		n_nw_tris=0, n_ns_tris=0, n_ws_tris=0, n_nws_seg=0;
+
+		// Construct the interfaces and common curve
+		pmmc_ConstructLocalCube(SignDist, Phase, solid_isovalue, fluid_isovalue,
+				nw_pts, nw_tris, values, ns_pts, ns_tris, ws_pts, ws_tris,
+				local_nws_pts, nws_pts, nws_seg, local_sol_pts, local_sol_tris,
+				n_local_sol_tris, n_local_sol_pts, n_nw_pts, n_nw_tris,
+				n_ws_pts, n_ws_tris, n_ns_tris, n_ns_pts, n_local_nws_pts, n_nws_pts, n_nws_seg,
+				i, j, k, Nx, Ny, Nz);
+
+		// Copy the curvature values for the cube
+		CubeValues(0,0,0) = Phase_x(i,j,k);
+		CubeValues(1,0,0) = Phase_x(i+1,j,k);
+		CubeValues(0,1,0) = Phase_x(i,j+1,k);
+		CubeValues(1,1,0) = Phase_x(i+1,j+1,k);
+		CubeValues(0,0,1) = Phase_x(i,j,k+1);
+		CubeValues(1,0,1) = Phase_x(i+1,j,k+1);
+		CubeValues(0,1,1) = Phase_x(i,j+1,k+1);
+		CubeValues(1,1,1) = Phase_x(i+1,j+1,k+1);
+
+		// Interpolate the curvature onto the surface
+//		wn_curvature_sum += pmmc_CubeSurfaceInterpValue(CubeValues, nw_pts, nw_tris,
+//									wn_curvature, i, j, k, n_nw_pts, n_nw_tris);
+
+		wn_area_sum += pmmc_CubeSurfaceArea(nw_pts, nw_tris, n_nw_tris);
+
+	}
+
+	printf("Gxx =  %f, Analytical = 1/3 \n", wn_curvature_sum/wn_area_sum);
+	
+}