www.pudn.com > Image_segment.rar > LevelSetMethods.cpp

//Copyright (c) 2004-2005, Baris Sumengen 
//All rights reserved. 
// 
// CIMPL Matrix Performance Library 
// 
//Redistribution and use in source and binary 
//forms, with or without modification, are 
//permitted provided that the following 
//conditions are met: 
// 
//    * No commercial use is allowed.  
//    This software can only be used 
//    for non-commercial purposes. This  
//    distribution is mainly intended for 
//    academic research and teaching. 
//    * Redistributions of source code must 
//    retain the above copyright notice, this 
//    list of conditions and the following 
//    disclaimer. 
//    * Redistributions of binary form must 
//    mention the above copyright notice, this 
//    list of conditions and the following 
//    disclaimer in a clearly visible part  
//    in associated product manual,  
//    readme, and web site of the redistributed  
//    software. 
//    * Redistributions in binary form must 
//    reproduce the above copyright notice, 
//    this list of conditions and the 
//    following disclaimer in the 
//    documentation and/or other materials 
//    provided with the distribution. 
//    * The name of Baris Sumengen may not be 
//    used to endorse or promote products 
//    derived from this software without 
//    specific prior written permission. 
// 
//THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT 
//HOLDERS AND CONTRIBUTORS "AS IS" AND ANY 
//EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT 
//NOT LIMITED TO, THE IMPLIED WARRANTIES OF 
//MERCHANTABILITY AND FITNESS FOR A PARTICULAR 
//PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE 
//CONTRIBUTORS BE LIABLE FOR ANY 
//DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 
//EXEMPLARY, OR CONSEQUENTIAL DAMAGES 
//(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT 
//OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
//DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 
//HOWEVER CAUSED AND ON ANY THEORY OF 
//LIABILITY, WHETHER IN CONTRACT, STRICT 
//LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR 
//OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
//OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
//POSSIBILITY OF SUCH DAMAGE. 
 
 
 
#include "./LevelSetMethods.h" 
 
 
namespace LevelSetMethods 
{ 
 
 
	Matrix ExtractCurve(Matrix phi) 
	{ 
		Matrix curve(phi.Rows(), phi.Columns(), 0); 
 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int i=4;i<=endR-4;i++) 
		{ 
			for(int j=4;j<=endC-4;j++) 
			{ 
				if(phi.ElemNC(i+1,j)*phi.ElemNC(i,j) <= 0) 
				{ 
					curve.ElemNC(i+1,j) = 1; 
					curve.ElemNC(i,j) = 1; 
				} 
				if(phi.ElemNC(i-1,j)*phi.ElemNC(i,j) <= 0) 
				{ 
					curve.ElemNC(i-1,j) = 1; 
					curve.ElemNC(i,j) = 1; 
				} 
				if(phi.ElemNC(i,j+1)*phi.ElemNC(i,j) <= 0) 
				{ 
					curve.ElemNC(i,j+1) = 1; 
					curve.ElemNC(i,j) = 1; 
				} 
				if(phi.ElemNC(i,j-1)*phi.ElemNC(i,j) <= 0) 
				{ 
					curve.ElemNC(i,j-1) = 1; 
					curve.ElemNC(i,j) = 1; 
				} 
			} 
		} 
 
		return curve; 
	} 
 
 
	Matrix ExtractCurve(Matrix phi) 
	{ 
		Matrix curve(phi.Rows(), phi.Columns(), 0); 
 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int i=4;i<=endR-4;i++) 
		{ 
			for(int j=4;j<=endC-4;j++) 
			{ 
				if(phi.ElemNC(i+1,j)*phi.ElemNC(i,j) <= 0) 
				{ 
					curve.ElemNC(i+1,j) = 1; 
					curve.ElemNC(i,j) = 1; 
				} 
				if(phi.ElemNC(i-1,j)*phi.ElemNC(i,j) <= 0) 
				{ 
					curve.ElemNC(i-1,j) = 1; 
					curve.ElemNC(i,j) = 1; 
				} 
				if(phi.ElemNC(i,j+1)*phi.ElemNC(i,j) <= 0) 
				{ 
					curve.ElemNC(i,j+1) = 1; 
					curve.ElemNC(i,j) = 1; 
				} 
				if(phi.ElemNC(i,j-1)*phi.ElemNC(i,j) <= 0) 
				{ 
					curve.ElemNC(i,j-1) = 1; 
					curve.ElemNC(i,j) = 1; 
				} 
			} 
		} 
 
		return curve; 
	} 
 
 
 
 
	Matrix& ExtendConst2D(Matrix& phi, int size) 
	{ 
		if(size < 1 || size >3) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Extension size should be between 1 and 3!"); 
		} 
 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		// Extrapolate first layer 
		for(int i=3;i<=endR-3;i++) 
		{ 
			phi.ElemNC(i,2) = phi.ElemNC(i,3); 
			phi.ElemNC(i,endC-2) = phi.ElemNC(i,endC-3); 
		} 
		for(int j=3;j<=endC-3;j++) 
		{ 
			phi.ElemNC(2,j) = phi.ElemNC(3,j); 
			phi.ElemNC(endR-2,j) = phi.ElemNC(endR-3,j); 
		} 
		// May be useful for kappa  
		phi.ElemNC(2,2) = phi.ElemNC(3,3); 
		phi.ElemNC(2,endC-2) = phi.ElemNC(3,endC-3); 
		phi.ElemNC(endR-2,2) = phi.ElemNC(endR-3,3); 
		phi.ElemNC(endR-2,endC-2) = phi.ElemNC(endR-3,endC-3); 
 
		// Extend second layer 
		if(size>=2) 
		{ 
			for(int i=2;i<=endR-2;i++) 
			{ 
				phi.ElemNC(i,1) = phi.ElemNC(i,2); 
				phi.ElemNC(i,endC-1) = phi.ElemNC(i,endC-2); 
			} 
			for(int j=2;j<=endC-2;j++) 
			{ 
				phi.ElemNC(1,j) = phi.ElemNC(2,j); 
				phi.ElemNC(endR-1,j) = phi.ElemNC(endR-2,j); 
			} 
		} 
 
 
		// Extend third layer 
		if(size>=3) 
		{ 
			for(int i=1;i<=endR-1;i++) 
			{ 
				phi.ElemNC(i,0) = phi.ElemNC(i,1); 
				phi.ElemNC(i,endC-0) = phi.ElemNC(i,endC-1); 
			} 
			for(int j=1;j<=endC-1;j++) 
			{ 
				phi.ElemNC(0,j) = phi.ElemNC(1,j); 
				phi.ElemNC(endR-0,j) = phi.ElemNC(endR-1,j); 
			} 
		} 
 
 
		return phi; 
	} 
 
	Matrix& Extend2D(Matrix& phi, int size) 
	{ 
		if(size < 1 || size >3) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Extension size should be between 1 and 3!"); 
		} 
 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		// Extrapolate first layer 
		for(int i=3;i<=endR-3;i++) 
		{ 
			phi.ElemNC(i,2) = 2*phi.ElemNC(i,3) - phi.ElemNC(i,4); 
			phi.ElemNC(i,endC-2) = 2*phi.ElemNC(i,endC-3) - phi.ElemNC(i,endC-4); 
		} 
		for(int j=3;j<=endC-3;j++) 
		{ 
			phi.ElemNC(2,j) = 2*phi.ElemNC(3,j) - phi.ElemNC(4,j); 
			phi.ElemNC(endR-2,j) = 2*phi.ElemNC(endR-3,j) - phi.ElemNC(endR-4,j); 
		} 
		// May be useful for kappa  
		phi.ElemNC(2,2) = 2*phi.ElemNC(3,3) - phi.ElemNC(4,4); 
		phi.ElemNC(2,endC-2) = 2*phi.ElemNC(3,endC-3) - phi.ElemNC(4,endC-4); 
		phi.ElemNC(endR-2,2) = 2*phi.ElemNC(endR-3,3) - phi.ElemNC(endR-4,4); 
		phi.ElemNC(endR-2,endC-2) = 2*phi.ElemNC(endR-3,endC-3) - phi.ElemNC(endR-4,endC-4); 
 
 
		// Extrapolate second layer 
		if(size>=2) 
		{ 
			for(int i=2;i<=endR-2;i++) 
			{ 
				phi.ElemNC(i,1) = 2*phi.ElemNC(i,2) - phi.ElemNC(i,3); 
				phi.ElemNC(i,endC-1) = 2*phi.ElemNC(i,endC-2) - phi.ElemNC(i,endC-3); 
			} 
			for(int j=2;j<=endC-2;j++) 
			{ 
				phi.ElemNC(1,j) = 2*phi.ElemNC(2,j) - phi.ElemNC(3,j); 
				phi.ElemNC(endR-1,j) = 2*phi.ElemNC(endR-2,j) - phi.ElemNC(endR-3,j); 
			} 
		} 
 
 
		// Extrapolate third layer 
		if(size>=3) 
		{ 
			for(int i=1;i<=endR-1;i++) 
			{ 
				phi.ElemNC(i,0) = 2*phi.ElemNC(i,1) - phi.ElemNC(i,2); 
				phi.ElemNC(i,endC-0) = 2*phi.ElemNC(i,endC-1) - phi.ElemNC(i,endC-2); 
			} 
			for(int j=1;j<=endC-1;j++) 
			{ 
				phi.ElemNC(0,j) = 2*phi.ElemNC(1,j) - phi.ElemNC(2,j); 
				phi.ElemNC(endR-0,j) = 2*phi.ElemNC(endR-1,j) - phi.ElemNC(endR-2,j); 
			} 
		} 
 
		return phi; 
	} 
 
 
	Matrix& Extend2D(Matrix& phi, int size) 
	{ 
		if(size < 1 || size >3) 
		{ 
			cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl; 
			Utility::RunTimeError("Extension size should be between 1 and 3!"); 
		} 
 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		// Extrapolate first layer 
		for(int i=3;i<=endR-3;i++) 
		{ 
			phi.ElemNC(i,2) = 2*phi.ElemNC(i,3) - phi.ElemNC(i,4); 
			phi.ElemNC(i,endC-2) = 2*phi.ElemNC(i,endC-3) - phi.ElemNC(i,endC-4); 
		} 
		for(int j=3;j<=endC-3;j++) 
		{ 
			phi.ElemNC(2,j) = 2*phi.ElemNC(3,j) - phi.ElemNC(4,j); 
			phi.ElemNC(endR-2,j) = 2*phi.ElemNC(endR-3,j) - phi.ElemNC(endR-4,j); 
		} 
		// May be useful for kappa  
		phi.ElemNC(2,2) = 2*phi.ElemNC(3,3) - phi.ElemNC(4,4); 
		phi.ElemNC(2,endC-2) = 2*phi.ElemNC(3,endC-3) - phi.ElemNC(4,endC-4); 
		phi.ElemNC(endR-2,2) = 2*phi.ElemNC(endR-3,3) - phi.ElemNC(endR-4,4); 
		phi.ElemNC(endR-2,endC-2) = 2*phi.ElemNC(endR-3,endC-3) - phi.ElemNC(endR-4,endC-4); 
 
 
		// Extrapolate second layer 
		if(size>=2) 
		{ 
			for(int i=2;i<=endR-2;i++) 
			{ 
				phi.ElemNC(i,1) = 2*phi.ElemNC(i,2) - phi.ElemNC(i,3); 
				phi.ElemNC(i,endC-1) = 2*phi.ElemNC(i,endC-2) - phi.ElemNC(i,endC-3); 
			} 
			for(int j=2;j<=endC-2;j++) 
			{ 
				phi.ElemNC(1,j) = 2*phi.ElemNC(2,j) - phi.ElemNC(3,j); 
				phi.ElemNC(endR-1,j) = 2*phi.ElemNC(endR-2,j) - phi.ElemNC(endR-3,j); 
			} 
		} 
 
 
		// Extrapolate third layer 
		if(size>=3) 
		{ 
			for(int i=1;i<=endR-1;i++) 
			{ 
				phi.ElemNC(i,0) = 2*phi.ElemNC(i,1) - phi.ElemNC(i,2); 
				phi.ElemNC(i,endC-0) = 2*phi.ElemNC(i,endC-1) - phi.ElemNC(i,endC-2); 
			} 
			for(int j=1;j<=endC-1;j++) 
			{ 
				phi.ElemNC(0,j) = 2*phi.ElemNC(1,j) - phi.ElemNC(2,j); 
				phi.ElemNC(endR-0,j) = 2*phi.ElemNC(endR-1,j) - phi.ElemNC(endR-2,j); 
			} 
		} 
 
		return phi; 
	} 
 
 
 
 
 
 
 
 
 
 
	void Evolve2DNormalENO1(Matrix& phi, float dx, float dy, Matrix& Vn, Matrix& delta, Matrix& H1Abs, Matrix& H2Abs) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				float phi_x_m = DerENO1Minus(phi.ElemNC(i,j-1), phi.ElemNC(i,j), dx); 
				float phi_x_p = DerENO1Plus(phi.ElemNC(i,j), phi.ElemNC(i,j+1), dx); 
				float phi_x = SelectDerNormal(phi_x_m, phi_x_p, Vn.ElemNC(i,j)); 
 
				float phi_y_m = DerENO1Minus(phi.ElemNC(i-1,j), phi.ElemNC(i,j), dy); 
				float phi_y_p = DerENO1Plus(phi.ElemNC(i,j), phi.ElemNC(i+1,j), dy); 
				float phi_y = SelectDerNormal(phi_y_m, phi_y_p, Vn.ElemNC(i,j)); 
 
				float absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y); 
				H1Abs.ElemNC(i,j) = (float)fabs( Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				H2Abs.ElemNC(i,j) = (float)fabs( Vn.ElemNC(i,j)*phi_y*phi_y/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				delta.ElemNC(i,j) = Vn.ElemNC(i,j)*absGradPhi; 
			} 
		} 
	} 
 
	void Evolve2DNormalENO2(Matrix& phi, float dx, float dy, Matrix& Vn, Matrix& delta, Matrix& H1Abs, Matrix& H2Abs) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				float phi_x_m = DerENO2Minus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), dx); 
				float phi_x_p = DerENO2Plus(phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx); 
				float phi_x = SelectDerNormal(phi_x_m, phi_x_p, Vn.ElemNC(i,j)); 
 
				float phi_y_m = DerENO2Minus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), dy); 
				float phi_y_p = DerENO2Plus(phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy); 
				float phi_y = SelectDerNormal(phi_y_m, phi_y_p, Vn.ElemNC(i,j)); 
 
				float absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y); 
				H1Abs.ElemNC(i,j) = (float)fabs( Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				H2Abs.ElemNC(i,j) = (float)fabs( Vn.ElemNC(i,j)*phi_y*phi_y/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				delta.ElemNC(i,j) = Vn.ElemNC(i,j)*absGradPhi; 
			} 
		} 
	} 
 
	void Evolve2DNormalENO3(Matrix& phi, float dx, float dy, Matrix& Vn, Matrix& delta, Matrix& H1Abs, Matrix& H2Abs) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				float phi_x_m = DerENO3Minus(phi.ElemNC(i,j-3), phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx); 
				float phi_x_p = DerENO3Plus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), phi.ElemNC(i,j+3), dx); 
				float phi_x = SelectDerNormal(phi_x_m, phi_x_p, Vn.ElemNC(i,j)); 
 
				float phi_y_m = DerENO3Minus(phi.ElemNC(i-3,j), phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy); 
				float phi_y_p = DerENO3Plus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), phi.ElemNC(i+3,j), dy); 
				float phi_y = SelectDerNormal(phi_y_m, phi_y_p, Vn.ElemNC(i,j)); 
 
				float absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y); 
				H1Abs.ElemNC(i,j) = (float)fabs( Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				H2Abs.ElemNC(i,j) = (float)fabs( Vn.ElemNC(i,j)*phi_y*phi_y/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				delta.ElemNC(i,j) = Vn.ElemNC(i,j)*absGradPhi; 
			} 
		} 
	} 
 
	void Evolve2DNormalWENO(Matrix& phi, float dx, float dy, Matrix& Vn, Matrix& delta, Matrix& H1Abs, Matrix& H2Abs) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				float phi_x_m = DerWENOMinus(phi.ElemNC(i,j-3), phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx); 
				float phi_x_p = DerWENOPlus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), phi.ElemNC(i,j+3), dx); 
				float phi_x = SelectDerNormal(phi_x_m, phi_x_p, Vn.ElemNC(i,j)); 
 
				float phi_y_m = DerWENOMinus(phi.ElemNC(i-3,j), phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy); 
				float phi_y_p = DerWENOPlus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), phi.ElemNC(i+3,j), dy); 
				float phi_y = SelectDerNormal(phi_y_m, phi_y_p, Vn.ElemNC(i,j)); 
 
				float absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y); 
				H1Abs.ElemNC(i,j) = (float)fabs( Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				H2Abs.ElemNC(i,j) = (float)fabs( Vn.ElemNC(i,j)*phi_y*phi_y/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				delta.ElemNC(i,j) = Vn.ElemNC(i,j)*absGradPhi; 
			} 
		} 
	} 
 
 
 
 
 
 
 
	void Evolve2DNormalENO1(Matrix& phi, double dx, double dy, Matrix& Vn, Matrix& delta, Matrix& H1Abs, Matrix& H2Abs) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				double phi_x_m = DerENO1Minus(phi.ElemNC(i,j-1), phi.ElemNC(i,j), dx); 
				double phi_x_p = DerENO1Plus(phi.ElemNC(i,j), phi.ElemNC(i,j+1), dx); 
				double phi_x = SelectDerNormal(phi_x_m, phi_x_p, Vn.ElemNC(i,j)); 
 
				double phi_y_m = DerENO1Minus(phi.ElemNC(i-1,j), phi.ElemNC(i,j), dy); 
				double phi_y_p = DerENO1Plus(phi.ElemNC(i,j), phi.ElemNC(i+1,j), dy); 
				double phi_y = SelectDerNormal(phi_y_m, phi_y_p, Vn.ElemNC(i,j)); 
 
				double absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y); 
				H1Abs.ElemNC(i,j) = fabs( Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				H2Abs.ElemNC(i,j) = fabs( Vn.ElemNC(i,j)*phi_y*phi_y/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				delta.ElemNC(i,j) = Vn.ElemNC(i,j)*absGradPhi; 
			} 
		} 
	} 
 
	void Evolve2DNormalENO2(Matrix& phi, double dx, double dy, Matrix& Vn, Matrix& delta, Matrix& H1Abs, Matrix& H2Abs) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				double phi_x_m = DerENO2Minus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), dx); 
				double phi_x_p = DerENO2Plus(phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx); 
				double phi_x = SelectDerNormal(phi_x_m, phi_x_p, Vn.ElemNC(i,j)); 
 
				double phi_y_m = DerENO2Minus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), dy); 
				double phi_y_p = DerENO2Plus(phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy); 
				double phi_y = SelectDerNormal(phi_y_m, phi_y_p, Vn.ElemNC(i,j)); 
 
				double absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y); 
				H1Abs.ElemNC(i,j) = fabs( Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				H2Abs.ElemNC(i,j) = fabs( Vn.ElemNC(i,j)*phi_y*phi_y/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				delta.ElemNC(i,j) = Vn.ElemNC(i,j)*absGradPhi; 
			} 
		} 
	} 
 
	void Evolve2DNormalENO3(Matrix& phi, double dx, double dy, Matrix& Vn, Matrix& delta, Matrix& H1Abs, Matrix& H2Abs) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				double phi_x_m = DerENO3Minus(phi.ElemNC(i,j-3), phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx); 
				double phi_x_p = DerENO3Plus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), phi.ElemNC(i,j+3), dx); 
				double phi_x = SelectDerNormal(phi_x_m, phi_x_p, Vn.ElemNC(i,j)); 
 
				double phi_y_m = DerENO3Minus(phi.ElemNC(i-3,j), phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy); 
				double phi_y_p = DerENO3Plus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), phi.ElemNC(i+3,j), dy); 
				double phi_y = SelectDerNormal(phi_y_m, phi_y_p, Vn.ElemNC(i,j)); 
 
				double absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y); 
				H1Abs.ElemNC(i,j) = fabs( Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				H2Abs.ElemNC(i,j) = fabs( Vn.ElemNC(i,j)*phi_y*phi_y/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				delta.ElemNC(i,j) = Vn.ElemNC(i,j)*absGradPhi; 
			} 
		} 
	} 
 
	void Evolve2DNormalWENO(Matrix& phi, double dx, double dy, Matrix& Vn, Matrix& delta, Matrix& H1Abs, Matrix& H2Abs) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				double phi_x_m = DerWENOMinus(phi.ElemNC(i,j-3), phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx); 
				double phi_x_p = DerWENOPlus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), phi.ElemNC(i,j+3), dx); 
				double phi_x = SelectDerNormal(phi_x_m, phi_x_p, Vn.ElemNC(i,j)); 
 
				double phi_y_m = DerWENOMinus(phi.ElemNC(i-3,j), phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy); 
				double phi_y_p = DerWENOPlus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), phi.ElemNC(i+3,j), dy); 
				double phi_y = SelectDerNormal(phi_y_m, phi_y_p, Vn.ElemNC(i,j)); 
 
				double absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y); 
				H1Abs.ElemNC(i,j) = fabs( Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				H2Abs.ElemNC(i,j) = fabs( Vn.ElemNC(i,j)*phi_y*phi_y/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				delta.ElemNC(i,j) = Vn.ElemNC(i,j)*absGradPhi; 
			} 
		} 
	} 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
	void Evolve2DVectorENO1(Matrix& phi, float dx, float dy, Matrix& u, Matrix& v, Matrix& delta) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				float phi_x, phi_y; 
				if(u.ElemNC(i,j)>0) 
				{ 
					phi_x = DerENO1Minus(phi.ElemNC(i,j-1), phi.ElemNC(i,j), dx); 
				} 
				else if(u.ElemNC(i,j)<0) 
				{ 
					phi_x = DerENO1Plus(phi.ElemNC(i,j), phi.ElemNC(i,j+1), dx); 
				} 
				else 
				{ 
					phi_x = 0; 
				} 
 
				if(v.ElemNC(i,j)>0) 
				{ 
					phi_y = DerENO1Minus(phi.ElemNC(i-1,j), phi.ElemNC(i,j), dy); 
				} 
				else if(v.ElemNC(i,j)<0) 
				{ 
					phi_y = DerENO1Plus(phi.ElemNC(i,j), phi.ElemNC(i+1,j), dy); 
				} 
				else 
				{ 
					phi_y = 0; 
				} 
 
				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y; 
			} 
		} 
	} 
 
 
	void Evolve2DVectorENO2(Matrix& phi, float dx, float dy, Matrix& u, Matrix& v, Matrix& delta) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				float phi_x, phi_y; 
				if(u.ElemNC(i,j)>0) 
				{ 
					phi_x = DerENO2Minus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), dx); 
				} 
				else if(u.ElemNC(i,j)<0) 
				{ 
					phi_x = DerENO2Plus(phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx); 
				} 
				else 
				{ 
					phi_x = 0; 
				} 
 
				if(v.ElemNC(i,j)>0) 
				{ 
					phi_y = DerENO2Minus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), dy); 
				} 
				else if(v.ElemNC(i,j)<0) 
				{ 
					phi_y = DerENO2Plus(phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy); 
				} 
				else 
				{ 
					phi_y = 0; 
				} 
 
				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y; 
			} 
		} 
	} 
 
 
	void Evolve2DVectorENO3(Matrix& phi, float dx, float dy, Matrix& u, Matrix& v, Matrix& delta) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				float phi_x, phi_y; 
				if(u.ElemNC(i,j)>0) 
				{ 
					phi_x = DerENO3Minus(phi.ElemNC(i,j-3), phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx); 
				} 
				else if(u.ElemNC(i,j)<0) 
				{ 
					phi_x = DerENO3Plus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), phi.ElemNC(i,j+3), dx); 
				} 
				else 
				{ 
					phi_x = 0; 
				} 
 
 
				if(v.ElemNC(i,j)>0) 
				{ 
					phi_y = DerENO3Minus(phi.ElemNC(i-3,j), phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy); 
				} 
				else if(v.ElemNC(i,j)<0) 
				{ 
					phi_y = DerENO3Plus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), phi.ElemNC(i+3,j), dy); 
				} 
				else 
				{ 
					phi_y = 0; 
				} 
 
				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y; 
			} 
		} 
	} 
 
	void Evolve2DVectorWENO(Matrix& phi, float dx, float dy, Matrix& u, Matrix& v, Matrix& delta) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				float phi_x, phi_y; 
				if(u.ElemNC(i,j)>0) 
				{ 
					phi_x = DerWENOMinus(phi.ElemNC(i,j-3), phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx); 
				} 
				else if(u.ElemNC(i,j)<0) 
				{ 
					phi_x = DerWENOPlus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), phi.ElemNC(i,j+3), dx); 
				} 
				else 
				{ 
					phi_x = 0; 
				} 
 
 
				if(v.ElemNC(i,j)>0) 
				{ 
					phi_y = DerWENOMinus(phi.ElemNC(i-3,j), phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy); 
				} 
				else if(v.ElemNC(i,j)<0) 
				{ 
					phi_y = DerWENOPlus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), phi.ElemNC(i+3,j), dy); 
				} 
				else 
				{ 
					phi_y = 0; 
				} 
 
				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y; 
			} 
		} 
	} 
 
 
 
 
 
 
	void Evolve2DVectorENO1(Matrix& phi, double dx, double dy, Matrix& u, Matrix& v, Matrix& delta) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				double phi_x, phi_y; 
				if(u.ElemNC(i,j)>0) 
				{ 
					phi_x = DerENO1Minus(phi.ElemNC(i,j-1), phi.ElemNC(i,j), dx); 
				} 
				else if(u.ElemNC(i,j)<0) 
				{ 
					phi_x = DerENO1Plus(phi.ElemNC(i,j), phi.ElemNC(i,j+1), dx); 
				} 
				else 
				{ 
					phi_x = 0; 
				} 
 
				if(v.ElemNC(i,j)>0) 
				{ 
					phi_y = DerENO1Minus(phi.ElemNC(i-1,j), phi.ElemNC(i,j), dy); 
				} 
				else if(v.ElemNC(i,j)<0) 
				{ 
					phi_y = DerENO1Plus(phi.ElemNC(i,j), phi.ElemNC(i+1,j), dy); 
				} 
				else 
				{ 
					phi_y = 0; 
				} 
 
				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y; 
			} 
		} 
	} 
 
 
	void Evolve2DVectorENO2(Matrix& phi, double dx, double dy, Matrix& u, Matrix& v, Matrix& delta) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				double phi_x, phi_y; 
				if(u.ElemNC(i,j)>0) 
				{ 
					phi_x = DerENO2Minus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), dx); 
				} 
				else if(u.ElemNC(i,j)<0) 
				{ 
					phi_x = DerENO2Plus(phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx); 
				} 
				else 
				{ 
					phi_x = 0; 
				} 
 
				if(v.ElemNC(i,j)>0) 
				{ 
					phi_y = DerENO2Minus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), dy); 
				} 
				else if(v.ElemNC(i,j)<0) 
				{ 
					phi_y = DerENO2Plus(phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy); 
				} 
				else 
				{ 
					phi_y = 0; 
				} 
 
				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y; 
			} 
		} 
	} 
 
 
	void Evolve2DVectorENO3(Matrix& phi, double dx, double dy, Matrix& u, Matrix& v, Matrix& delta) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				double phi_x, phi_y; 
				if(u.ElemNC(i,j)>0) 
				{ 
					phi_x = DerENO3Minus(phi.ElemNC(i,j-3), phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx); 
				} 
				else if(u.ElemNC(i,j)<0) 
				{ 
					phi_x = DerENO3Plus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), phi.ElemNC(i,j+3), dx); 
				} 
				else 
				{ 
					phi_x = 0; 
				} 
 
 
				if(v.ElemNC(i,j)>0) 
				{ 
					phi_y = DerENO3Minus(phi.ElemNC(i-3,j), phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy); 
				} 
				else if(v.ElemNC(i,j)<0) 
				{ 
					phi_y = DerENO3Plus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), phi.ElemNC(i+3,j), dy); 
				} 
				else 
				{ 
					phi_y = 0; 
				} 
 
				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y; 
			} 
		} 
	} 
 
	void Evolve2DVectorWENO(Matrix& phi, double dx, double dy, Matrix& u, Matrix& v, Matrix& delta) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				double phi_x, phi_y; 
				if(u.ElemNC(i,j)>0) 
				{ 
					phi_x = DerWENOMinus(phi.ElemNC(i,j-3), phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx); 
				} 
				else if(u.ElemNC(i,j)<0) 
				{ 
					phi_x = DerWENOPlus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), phi.ElemNC(i,j+3), dx); 
				} 
				else 
				{ 
					phi_x = 0; 
				} 
 
 
				if(v.ElemNC(i,j)>0) 
				{ 
					phi_y = DerWENOMinus(phi.ElemNC(i-3,j), phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy); 
				} 
				else if(v.ElemNC(i,j)<0) 
				{ 
					phi_y = DerWENOPlus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), phi.ElemNC(i+3,j), dy); 
				} 
				else 
				{ 
					phi_y = 0; 
				} 
 
				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y; 
			} 
		} 
	} 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
	void Evolve2DNormalVectorENO1SignedDistance(Matrix& phi, float dx, float dy, Matrix& u, Matrix& v, Matrix& Vn, Matrix& delta, Matrix& H1Abs, Matrix& H2Abs) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				float phi_x_m = DerENO1Minus(phi.ElemNC(i,j-1), phi.ElemNC(i,j), dx); 
				float phi_x_p = DerENO1Plus(phi.ElemNC(i,j), phi.ElemNC(i,j+1), dx); 
				float phi_x = SelectDerNormalVectorSD(phi_x_m, phi_x_p, u.ElemNC(i,j), Vn.ElemNC(i,j)); 
 
				float phi_y_m = DerENO1Minus(phi.ElemNC(i-1,j), phi.ElemNC(i,j), dy); 
				float phi_y_p = DerENO1Plus(phi.ElemNC(i,j), phi.ElemNC(i+1,j), dy); 
				float phi_y = SelectDerNormalVectorSD(phi_y_m, phi_y_p, v.ElemNC(i,j), Vn.ElemNC(i,j)); 
 
				float absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y); 
				H1Abs.ElemNC(i,j) = (float)fabs( u.ElemNC(i,j) + Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				H2Abs.ElemNC(i,j) = (float)fabs( v.ElemNC(i,j) + Vn.ElemNC(i,j)*phi_y*phi_y/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y + Vn.ElemNC(i,j)*absGradPhi; 
			} 
		} 
	} 
 
	void Evolve2DNormalVectorENO2SignedDistance(Matrix& phi, float dx, float dy, Matrix& u, Matrix& v, Matrix& Vn, Matrix& delta, Matrix& H1Abs, Matrix& H2Abs) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				float phi_x_m = DerENO2Minus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), dx); 
				float phi_x_p = DerENO2Plus(phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx); 
				float phi_x = SelectDerNormalVectorSD(phi_x_m, phi_x_p, u.ElemNC(i,j), Vn.ElemNC(i,j)); 
 
				float phi_y_m = DerENO2Minus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), dy); 
				float phi_y_p = DerENO2Plus(phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy); 
				float phi_y = SelectDerNormalVectorSD(phi_y_m, phi_y_p, v.ElemNC(i,j), Vn.ElemNC(i,j)); 
 
				float absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y); 
				H1Abs.ElemNC(i,j) = (float)fabs( u.ElemNC(i,j) + Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				H2Abs.ElemNC(i,j) = (float)fabs( v.ElemNC(i,j) + Vn.ElemNC(i,j)*phi_y*phi_y/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y + Vn.ElemNC(i,j)*absGradPhi; 
			} 
		} 
	} 
 
 
	void Evolve2DNormalVectorENO3SignedDistance(Matrix& phi, float dx, float dy, Matrix& u, Matrix& v, Matrix& Vn, Matrix& delta, Matrix& H1Abs, Matrix& H2Abs) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				float phi_x_m = DerENO3Minus(phi.ElemNC(i,j-3), phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx); 
				float phi_x_p = DerENO3Plus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), phi.ElemNC(i,j+3), dx); 
				float phi_x = SelectDerNormalVectorSD(phi_x_m, phi_x_p, u.ElemNC(i,j), Vn.ElemNC(i,j)); 
 
				float phi_y_m = DerENO3Minus(phi.ElemNC(i-3,j), phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy); 
				float phi_y_p = DerENO3Plus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), phi.ElemNC(i+3,j), dy); 
				float phi_y = SelectDerNormalVectorSD(phi_y_m, phi_y_p, v.ElemNC(i,j), Vn.ElemNC(i,j)); 
 
				float absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y); 
				H1Abs.ElemNC(i,j) = (float)fabs( u.ElemNC(i,j) + Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				H2Abs.ElemNC(i,j) = (float)fabs( v.ElemNC(i,j) + Vn.ElemNC(i,j)*phi_y*phi_y/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y + Vn.ElemNC(i,j)*absGradPhi; 
			} 
		} 
	} 
 
 
	void Evolve2DNormalVectorWENOSignedDistance(Matrix& phi, float dx, float dy, Matrix& u, Matrix& v, Matrix& Vn, Matrix& delta, Matrix& H1Abs, Matrix& H2Abs) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				float phi_x_m = DerWENOMinus(phi.ElemNC(i,j-3), phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx); 
				float phi_x_p = DerWENOPlus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), phi.ElemNC(i,j+3), dx); 
				float phi_x = SelectDerNormalVectorSD(phi_x_m, phi_x_p, u.ElemNC(i,j), Vn.ElemNC(i,j)); 
 
				float phi_y_m = DerWENOMinus(phi.ElemNC(i-3,j), phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy); 
				float phi_y_p = DerWENOPlus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), phi.ElemNC(i+3,j), dy); 
				float phi_y = SelectDerNormalVectorSD(phi_y_m, phi_y_p, v.ElemNC(i,j), Vn.ElemNC(i,j)); 
 
				float absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y); 
				H1Abs.ElemNC(i,j) = (float)fabs( u.ElemNC(i,j) + Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				H2Abs.ElemNC(i,j) = (float)fabs( v.ElemNC(i,j) + Vn.ElemNC(i,j)*phi_y*phi_y/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y + Vn.ElemNC(i,j)*absGradPhi; 
			} 
		} 
	} 
 
 
 
 
 
	void Evolve2DNormalVectorENO1SignedDistance(Matrix& phi, double dx, double dy, Matrix& u, Matrix& v, Matrix& Vn, Matrix& delta, Matrix& H1Abs, Matrix& H2Abs) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				double phi_x_m = DerENO1Minus(phi.ElemNC(i,j-1), phi.ElemNC(i,j), dx); 
				double phi_x_p = DerENO1Plus(phi.ElemNC(i,j), phi.ElemNC(i,j+1), dx); 
				double phi_x = SelectDerNormalVectorSD(phi_x_m, phi_x_p, u.ElemNC(i,j), Vn.ElemNC(i,j)); 
 
				double phi_y_m = DerENO1Minus(phi.ElemNC(i-1,j), phi.ElemNC(i,j), dy); 
				double phi_y_p = DerENO1Plus(phi.ElemNC(i,j), phi.ElemNC(i+1,j), dy); 
				double phi_y = SelectDerNormalVectorSD(phi_y_m, phi_y_p, v.ElemNC(i,j), Vn.ElemNC(i,j)); 
 
				double absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y); 
				H1Abs.ElemNC(i,j) = fabs( u.ElemNC(i,j) + Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				H2Abs.ElemNC(i,j) = fabs( v.ElemNC(i,j) + Vn.ElemNC(i,j)*phi_y*phi_y/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y + Vn.ElemNC(i,j)*absGradPhi; 
			} 
		} 
	} 
 
	void Evolve2DNormalVectorENO2SignedDistance(Matrix& phi, double dx, double dy, Matrix& u, Matrix& v, Matrix& Vn, Matrix& delta, Matrix& H1Abs, Matrix& H2Abs) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				double phi_x_m = DerENO2Minus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), dx); 
				double phi_x_p = DerENO2Plus(phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx); 
				double phi_x = SelectDerNormalVectorSD(phi_x_m, phi_x_p, u.ElemNC(i,j), Vn.ElemNC(i,j)); 
 
				double phi_y_m = DerENO2Minus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), dy); 
				double phi_y_p = DerENO2Plus(phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy); 
				double phi_y = SelectDerNormalVectorSD(phi_y_m, phi_y_p, v.ElemNC(i,j), Vn.ElemNC(i,j)); 
 
				double absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y); 
				H1Abs.ElemNC(i,j) = fabs( u.ElemNC(i,j) + Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				H2Abs.ElemNC(i,j) = fabs( v.ElemNC(i,j) + Vn.ElemNC(i,j)*phi_y*phi_y/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y + Vn.ElemNC(i,j)*absGradPhi; 
			} 
		} 
	} 
 
 
	void Evolve2DNormalVectorENO3SignedDistance(Matrix& phi, double dx, double dy, Matrix& u, Matrix& v, Matrix& Vn, Matrix& delta, Matrix& H1Abs, Matrix& H2Abs) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				double phi_x_m = DerENO3Minus(phi.ElemNC(i,j-3), phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx); 
				double phi_x_p = DerENO3Plus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), phi.ElemNC(i,j+3), dx); 
				double phi_x = SelectDerNormalVectorSD(phi_x_m, phi_x_p, u.ElemNC(i,j), Vn.ElemNC(i,j)); 
 
				double phi_y_m = DerENO3Minus(phi.ElemNC(i-3,j), phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy); 
				double phi_y_p = DerENO3Plus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), phi.ElemNC(i+3,j), dy); 
				double phi_y = SelectDerNormalVectorSD(phi_y_m, phi_y_p, v.ElemNC(i,j), Vn.ElemNC(i,j)); 
 
				double absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y); 
				H1Abs.ElemNC(i,j) = fabs( u.ElemNC(i,j) + Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				H2Abs.ElemNC(i,j) = fabs( v.ElemNC(i,j) + Vn.ElemNC(i,j)*phi_y*phi_y/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y + Vn.ElemNC(i,j)*absGradPhi; 
			} 
		} 
	} 
 
 
	void Evolve2DNormalVectorWENOSignedDistance(Matrix& phi, double dx, double dy, Matrix& u, Matrix& v, Matrix& Vn, Matrix& delta, Matrix& H1Abs, Matrix& H2Abs) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				double phi_x_m = DerWENOMinus(phi.ElemNC(i,j-3), phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), dx); 
				double phi_x_p = DerWENOPlus(phi.ElemNC(i,j-2), phi.ElemNC(i,j-1), phi.ElemNC(i,j), phi.ElemNC(i,j+1), phi.ElemNC(i,j+2), phi.ElemNC(i,j+3), dx); 
				double phi_x = SelectDerNormalVectorSD(phi_x_m, phi_x_p, u.ElemNC(i,j), Vn.ElemNC(i,j)); 
 
				double phi_y_m = DerWENOMinus(phi.ElemNC(i-3,j), phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), dy); 
				double phi_y_p = DerWENOPlus(phi.ElemNC(i-2,j), phi.ElemNC(i-1,j), phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+2,j), phi.ElemNC(i+3,j), dy); 
				double phi_y = SelectDerNormalVectorSD(phi_y_m, phi_y_p, v.ElemNC(i,j), Vn.ElemNC(i,j)); 
 
				double absGradPhi = sqrt(phi_x*phi_x+phi_y*phi_y); 
				H1Abs.ElemNC(i,j) = (double)fabs( u.ElemNC(i,j) + Vn.ElemNC(i,j)*phi_x*phi_x/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				H2Abs.ElemNC(i,j) = (double)fabs( v.ElemNC(i,j) + Vn.ElemNC(i,j)*phi_y*phi_y/(absGradPhi+(absGradPhi>1e-19?0:1)) ); 
				delta.ElemNC(i,j) = u.ElemNC(i,j)*phi_x + v.ElemNC(i,j)*phi_y + Vn.ElemNC(i,j)*absGradPhi; 
			} 
		} 
	} 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
	float Getdt2DNormal(float alpha, float dx, float dy, Matrix& H1Abs, Matrix& H2Abs) 
	{ 
		int endR = H1Abs.Rows()-1; 
		int endC = H1Abs.Columns()-1; 
		float maxs = 0.0; 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				float temp = H1Abs.ElemNC(i,j)/dx + H2Abs.ElemNC(i,j)/dy; 
				maxs = maxs>=temp?maxs:temp; 
			} 
		} 
		return (float)(alpha/(maxs+(maxs>1e-19?0:1))); 
	} 
 
	float Getdt2DNormalVector(float alpha, float dx, float dy, Matrix& H1Abs, Matrix& H2Abs) 
	{ 
		int endR = H1Abs.Rows()-1; 
		int endC = H1Abs.Columns()-1; 
		float maxs = 0.0; 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				float temp = H1Abs.ElemNC(i,j)/dx + H2Abs.ElemNC(i,j)/dy; 
				maxs = maxs>=temp?maxs:temp; 
			} 
		} 
		return (float)(alpha/(maxs+(maxs>1e-19?0:1))); 
	} 
 
 
	float Getdt2DVector(float alpha, float dx, float dy, Matrix& u, Matrix& v) 
	{ 
		int endR = u.Rows()-1; 
		int endC = u.Columns()-1; 
		float maxs = 0.0; 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				float temp = fabs(u.ElemNC(i,j))/dx + fabs(v.ElemNC(i,j))/dy; 
				maxs = maxs>=temp?maxs:temp; 
			} 
		} 
		return (float)(alpha/(maxs+(maxs>1e-19?0:1))); 
	} 
 
 
 
	float Getdt2DKappa(float alpha, float dx2, float dy2, Matrix& b) 
	{ 
		int endR = b.Rows()-1; 
		int endC = b.Columns()-1; 
		float maxs = 0.0; 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				float temp = (2*b.ElemNC(i,j))/dx2 + (2*b.ElemNC(i,j))/dy2; 
				maxs = maxs>=temp?maxs:temp; 
			} 
		} 
		return (float)(alpha/(maxs+(maxs>1e-19?0:1))); 
	} 
 
	float Getdt2DVectorKappa(float alpha, float dx, float dy, float dx2, float dy2, Matrix& u, Matrix& v, Matrix& b) 
	{ 
		int endR = b.Rows()-1; 
		int endC = b.Columns()-1; 
		float maxs = 0.0; 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				float temp = fabs(u.ElemNC(i,j))/dx + fabs(v.ElemNC(i,j))/dy +  
					(2*b.ElemNC(i,j))/dx2 + (2*b.ElemNC(i,j))/dy2; 
				maxs = maxs>=temp?maxs:temp; 
			} 
		} 
		return (float)(alpha/(maxs+(maxs>1e-19?0:1))); 
	} 
 
	float Getdt2DNormalKappa(float alpha, float dx, float dy, float dx2, float dy2, Matrix& H1Abs, Matrix& H2Abs, Matrix& b) 
	{ 
		int endR = b.Rows()-1; 
		int endC = b.Columns()-1; 
		float maxs = 0.0; 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				float temp = H1Abs.ElemNC(i,j)/dx + H2Abs.ElemNC(i,j)/dy +  
					(2*b.ElemNC(i,j))/dx2 + (2*b.ElemNC(i,j))/dy2; 
				maxs = maxs>=temp?maxs:temp; 
			} 
		} 
		return (float)(alpha/(maxs+(maxs>1e-19?0:1))); 
	} 
 
	float Getdt2DNormalVectorKappa(float alpha, float dx, float dy, float dx2, float dy2, Matrix& H1Abs, Matrix& H2Abs, Matrix& b) 
	{ 
		int endR = b.Rows()-1; 
		int endC = b.Columns()-1; 
		float maxs = 0.0; 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				float temp = H1Abs.ElemNC(i,j)/dx + H2Abs.ElemNC(i,j)/dy +  
					(2*b.ElemNC(i,j))/dx2 + (2*b.ElemNC(i,j))/dy2; 
				maxs = maxs>=temp?maxs:temp; 
			} 
		} 
		return (float)(alpha/(maxs+(maxs>1e-19?0:1))); 
	} 
 
 
 
 
 
 
 
 
	double Getdt2DNormal(double alpha, double dx, double dy, Matrix& H1Abs, Matrix& H2Abs) 
	{ 
		int endR = H1Abs.Rows()-1; 
		int endC = H1Abs.Columns()-1; 
		double maxs = 0.0; 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				double temp = H1Abs.ElemNC(i,j)/dx + H2Abs.ElemNC(i,j)/dy; 
				maxs = maxs>=temp?maxs:temp; 
			} 
		} 
		return alpha/(maxs+(maxs>1e-19?0:1)); 
	} 
 
	double Getdt2DNormalVector(double alpha, double dx, double dy, Matrix& H1Abs, Matrix& H2Abs) 
	{ 
		int endR = H1Abs.Rows()-1; 
		int endC = H1Abs.Columns()-1; 
		double maxs = 0.0; 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				double temp = H1Abs.ElemNC(i,j)/dx + H2Abs.ElemNC(i,j)/dy; 
				maxs = maxs>=temp?maxs:temp; 
			} 
		} 
		return alpha/(maxs+(maxs>1e-19?0:1)); 
	} 
 
 
	double Getdt2DVector(double alpha, double dx, double dy, Matrix& u, Matrix& v) 
	{ 
		int endR = u.Rows()-1; 
		int endC = u.Columns()-1; 
		double maxs = 0.0; 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				double temp = fabs(u.ElemNC(i,j))/dx + fabs(v.ElemNC(i,j))/dy; 
				maxs = maxs>=temp?maxs:temp; 
			} 
		} 
		return alpha/(maxs+(maxs>1e-19?0:1)); 
	} 
 
 
 
	double Getdt2DKappa(double alpha, double dx2, double dy2, Matrix& b) 
	{ 
		int endR = b.Rows()-1; 
		int endC = b.Columns()-1; 
		double maxs = 0.0; 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				double temp = (2*b.ElemNC(i,j))/dx2 + (2*b.ElemNC(i,j))/dy2; 
				maxs = maxs>=temp?maxs:temp; 
			} 
		} 
		return alpha/(maxs+(maxs>1e-19?0:1)); 
	} 
 
	double Getdt2DVectorKappa(double alpha, double dx, double dy, double dx2, double dy2, Matrix& u, Matrix& v, Matrix& b) 
	{ 
		int endR = b.Rows()-1; 
		int endC = b.Columns()-1; 
		double maxs = 0.0; 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				double temp = fabs(u.ElemNC(i,j))/dx + fabs(v.ElemNC(i,j))/dy +  
					(2*b.ElemNC(i,j))/dx2 + (2*b.ElemNC(i,j))/dy2; 
				maxs = maxs>=temp?maxs:temp; 
			} 
		} 
		return alpha/(maxs+(maxs>1e-19?0:1)); 
	} 
 
	double Getdt2DNormalKappa(double alpha, double dx, double dy, double dx2, double dy2, Matrix& H1Abs, Matrix& H2Abs, Matrix& b) 
	{ 
		int endR = b.Rows()-1; 
		int endC = b.Columns()-1; 
		double maxs = 0.0; 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				double temp = H1Abs.ElemNC(i,j)/dx + H2Abs.ElemNC(i,j)/dy +  
					(2*b.ElemNC(i,j))/dx2 + (2*b.ElemNC(i,j))/dy2; 
				maxs = maxs>=temp?maxs:temp; 
			} 
		} 
		return alpha/(maxs+(maxs>1e-19?0:1)); 
	} 
 
	double Getdt2DNormalVectorKappa(double alpha, double dx, double dy, double dx2, double dy2, Matrix& H1Abs, Matrix& H2Abs, Matrix& b) 
	{ 
		int endR = b.Rows()-1; 
		int endC = b.Columns()-1; 
		double maxs = 0.0; 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				double temp = H1Abs.ElemNC(i,j)/dx + H2Abs.ElemNC(i,j)/dy +  
					(2*b.ElemNC(i,j))/dx2 + (2*b.ElemNC(i,j))/dy2; 
				maxs = maxs>=temp?maxs:temp; 
			} 
		} 
		return alpha/(maxs+(maxs>1e-19?0:1)); 
	} 
 
 
 
 
 
 
 
 
 
 
 
 
 
	void Evolve2DKappa(Matrix& phi, float dx, float dy, float dx2, float dy2, Matrix& b, Matrix& delta) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				float kappaAbsPhi = KappaAbsPhi2D(phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+1,j+1),  
					phi.ElemNC(i,j+1), phi.ElemNC(i-1,j+1), phi.ElemNC(i-1,j), phi.ElemNC(i-1,j-1),  
					phi.ElemNC(i,j-1), phi.ElemNC(i+1,j-1), dx, dy, dx2, dy2); 
				delta.ElemNC(i,j) = b.ElemNC(i,j)*kappaAbsPhi; 
			} 
		} 
	} 
 
 
 
 
	inline float KappaAbsPhi2D(float phi_i_j, float phi_ip1_j, float phi_ip1_jp1,  
		float phi_i_jp1, float phi_im1_jp1, float phi_im1_j, float phi_im1_jm1,  
		float phi_i_jm1, float phi_ip1_jm1, float dx, float dy, float dx2, float dy2) 
	{ 
		float phi_x = (phi_i_jp1 - phi_i_jm1)/(2*dx); 
		float phi_y = (phi_ip1_j - phi_im1_j)/(2*dy); 
		float phi_xx = (phi_i_jp1 - 2*phi_i_j + phi_i_jm1)/dx2; 
		float phi_yy = (phi_ip1_j - 2*phi_i_j + phi_im1_j)/dy2; 
		float phi_xy = (phi_ip1_jp1 + phi_im1_jm1 - phi_ip1_jm1 - phi_im1_jp1)/(4*dx*dy); 
		float absGradPhiSq = (phi_x*phi_x + phi_y*phi_y); 
		float kappaAbsPhi = (phi_xx*phi_y*phi_y - 2*phi_y*phi_x*phi_xy  
			+ phi_yy*phi_x*phi_x)/(absGradPhiSq + (float)(absGradPhiSq>1e-19?0:1)); 
		return kappaAbsPhi; 
 
	} 
 
 
 
 
	void Evolve2DKappa(Matrix& phi, double dx, double dy, double dx2, double dy2, Matrix& b, Matrix& delta) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				double kappaAbsPhi = KappaAbsPhi2D(phi.ElemNC(i,j), phi.ElemNC(i+1,j), phi.ElemNC(i+1,j+1),  
					phi.ElemNC(i,j+1), phi.ElemNC(i-1,j+1), phi.ElemNC(i-1,j), phi.ElemNC(i-1,j-1),  
					phi.ElemNC(i,j-1), phi.ElemNC(i+1,j-1), dx, dy, dx2, dy2); 
				delta.ElemNC(i,j) = b.ElemNC(i,j)*kappaAbsPhi; 
			} 
		} 
	} 
 
 
 
 
	inline double KappaAbsPhi2D(double phi_i_j, double phi_ip1_j, double phi_ip1_jp1,  
		double phi_i_jp1, double phi_im1_jp1, double phi_im1_j, double phi_im1_jm1,  
		double phi_i_jm1, double phi_ip1_jm1, double dx, double dy, double dx2, double dy2) 
	{ 
		double phi_x = (phi_i_jp1 - phi_i_jm1)/(2*dx); 
		double phi_y = (phi_ip1_j - phi_im1_j)/(2*dy); 
		double phi_xx = (phi_i_jp1 - 2*phi_i_j + phi_i_jm1)/dx2; 
		double phi_yy = (phi_ip1_j - 2*phi_i_j + phi_im1_j)/dy2; 
		double phi_xy = (phi_ip1_jp1 + phi_im1_jm1 - phi_ip1_jm1 - phi_im1_jp1)/(4*dx*dy); 
		double absGradPhiSq = (phi_x*phi_x + phi_y*phi_y); 
		double kappaAbsPhi = (phi_xx*phi_y*phi_y - 2*phi_y*phi_x*phi_xy  
			+ phi_yy*phi_x*phi_x)/(absGradPhiSq + (double)(absGradPhiSq>1e-19?0:1)); 
		return kappaAbsPhi; 
 
	} 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
	void ReinitSignedDistanceENO1(Matrix& phi, float dx, float dy, float alpha, int iterations) 
	{ 
		Matrix S_phi_0 = phi*InvSqrtI(phi*phi + dx*dx); 
		Matrix delta(phi.Rows(), phi.Columns()); 
		Matrix H1Abs(phi.Rows(), phi.Columns()); 
		Matrix H2Abs(phi.Rows(), phi.Columns()); 
		for(int k=0;k& phi, float dx, float dy, float alpha, int iterations) 
	{ 
		Matrix S_phi_0 = phi*InvSqrtI(phi*phi + dx*dx); 
		Matrix delta(phi.Rows(), phi.Columns()); 
		Matrix H1Abs(phi.Rows(), phi.Columns()); 
		Matrix H2Abs(phi.Rows(), phi.Columns()); 
		for(int k=0;k& phi, float dx, float dy, float alpha, int iterations) 
	{ 
		Matrix S_phi_0 = phi*InvSqrtI(phi*phi + dx*dx); 
		Matrix delta(phi.Rows(), phi.Columns()); 
		Matrix H1Abs(phi.Rows(), phi.Columns()); 
		Matrix H2Abs(phi.Rows(), phi.Columns()); 
		for(int k=0;k& phi, float dx, float dy, float alpha, int iterations) 
	{ 
		Matrix S_phi_0 = phi*InvSqrtI(phi*phi + dx*dx); 
		Matrix delta(phi.Rows(), phi.Columns()); 
		Matrix H1Abs(phi.Rows(), phi.Columns()); 
		Matrix H2Abs(phi.Rows(), phi.Columns()); 
		for(int k=0;k& phi, double dx, double dy, double alpha, int iterations) 
	{ 
		Matrix S_phi_0 = phi*InvSqrtI(phi*phi + dx*dx); 
		Matrix delta(phi.Rows(), phi.Columns()); 
		Matrix H1Abs(phi.Rows(), phi.Columns()); 
		Matrix H2Abs(phi.Rows(), phi.Columns()); 
		for(int k=0;k& phi, double dx, double dy, double alpha, int iterations) 
	{ 
		Matrix S_phi_0 = phi*InvSqrtI(phi*phi + dx*dx); 
		Matrix delta(phi.Rows(), phi.Columns()); 
		Matrix H1Abs(phi.Rows(), phi.Columns()); 
		Matrix H2Abs(phi.Rows(), phi.Columns()); 
		for(int k=0;k& phi, double dx, double dy, double alpha, int iterations) 
	{ 
		Matrix S_phi_0 = phi*InvSqrtI(phi*phi + dx*dx); 
		Matrix delta(phi.Rows(), phi.Columns()); 
		Matrix H1Abs(phi.Rows(), phi.Columns()); 
		Matrix H2Abs(phi.Rows(), phi.Columns()); 
		for(int k=0;k& phi, double dx, double dy, double alpha, int iterations) 
	{ 
		Matrix S_phi_0 = phi*InvSqrtI(phi*phi + dx*dx); 
		Matrix delta(phi.Rows(), phi.Columns()); 
		Matrix H1Abs(phi.Rows(), phi.Columns()); 
		Matrix H2Abs(phi.Rows(), phi.Columns()); 
		for(int k=0;k& phi, Matrix& delta, float dt) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				phi.ElemNC(i,j) -= dt*delta.ElemNC(i,j); 
			} 
		} 
	} 
	void AddPhi2D(Matrix& phi, Matrix& delta, float dt) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				phi.ElemNC(i,j) += dt*delta.ElemNC(i,j); 
			} 
		} 
	} 
 
 
	void SubtractPhi2D(Matrix& phi, Matrix& delta, double dt) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				phi.ElemNC(i,j) -= dt*delta.ElemNC(i,j); 
			} 
		} 
	} 
	void AddPhi2D(Matrix& phi, Matrix& delta, double dt) 
	{ 
		int endR = phi.Rows()-1; 
		int endC = phi.Columns()-1; 
		for(int j=3;j<=endC-3;j++) 
		{ 
			for(int i=3;i<=endR-3;i++) 
			{ 
				phi.ElemNC(i,j) += dt*delta.ElemNC(i,j); 
			} 
		} 
	} 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
	inline float SelectDerNormal(float phi_x_m, float phi_x_p, float Vn) 
	{ 
		float phi_x; 
		float VnPhi_m = Vn*phi_x_m; 
		float VnPhi_p = Vn*phi_x_p; 
		if(VnPhi_m <= 0 && VnPhi_p <= 0) 
		{ 
			phi_x = phi_x_p; 
		} 
		else if(VnPhi_m >= 0 && VnPhi_p >= 0) 
		{ 
			phi_x = phi_x_m; 
		} 
		else if(VnPhi_m <= 0 && VnPhi_p >= 0) 
		{ 
			phi_x = 0; 
		} 
		else if(VnPhi_m >= 0 && VnPhi_p <= 0) 
		{ 
			if(fabs(VnPhi_p) >= fabs(VnPhi_m)) 
			{ 
				phi_x = phi_x_p; 
			} 
			else 
			{ 
				phi_x = phi_x_m; 
			} 
		} 
		return phi_x; 
	} 
 
 
	inline float SelectDerNormalVectorSD(float phi_x_m, float phi_x_p, float u, float Vn) 
	{ 
		float phi_x; 
 
		float H1_m = u + Vn*phi_x_m; 
		float H1_p = u + Vn*phi_x_p; 
 
		if(H1_m <= 0 && H1_p <= 0) 
		{ 
			phi_x = phi_x_p; 
		} 
		else if(H1_m >= 0 && H1_p >= 0) 
		{ 
			phi_x = phi_x_m; 
		} 
		else if(H1_m <= 0 && H1_p >= 0) 
		{ 
			phi_x = -u/Vn; 
		} 
		else if(H1_m >= 0 && H1_p <= 0) 
		{ 
			if(fabs(H1_p) >= fabs(H1_m)) 
			{ 
				phi_x = phi_x_p; 
			} 
			else 
			{ 
				phi_x = phi_x_m; 
			} 
		} 
		return phi_x; 
	} 
 
	inline double SelectDerNormal(double phi_x_m, double phi_x_p, double Vn) 
	{ 
		double phi_x; 
		double VnPhi_m = Vn*phi_x_m; 
		double VnPhi_p = Vn*phi_x_p; 
		if(VnPhi_m <= 0 && VnPhi_p <= 0) 
		{ 
			phi_x = phi_x_p; 
		} 
		else if(VnPhi_m >= 0 && VnPhi_p >= 0) 
		{ 
			phi_x = phi_x_m; 
		} 
		else if(VnPhi_m <= 0 && VnPhi_p >= 0) 
		{ 
			phi_x = 0; 
		} 
		else if(VnPhi_m >= 0 && VnPhi_p <= 0) 
		{ 
			if(fabs(VnPhi_p) >= fabs(VnPhi_m)) 
			{ 
				phi_x = phi_x_p; 
			} 
			else 
			{ 
				phi_x = phi_x_m; 
			} 
		} 
		return phi_x; 
	} 
 
 
	inline double SelectDerNormalVectorSD(double phi_x_m, double phi_x_p, double u, double Vn) 
	{ 
		double phi_x; 
 
		double H1_m = u + Vn*phi_x_m; 
		double H1_p = u + Vn*phi_x_p; 
 
		if(H1_m <= 0 && H1_p <= 0) 
		{ 
			phi_x = phi_x_p; 
		} 
		else if(H1_m >= 0 && H1_p >= 0) 
		{ 
			phi_x = phi_x_m; 
		} 
		else if(H1_m <= 0 && H1_p >= 0) 
		{ 
			phi_x = -u/Vn; 
		} 
		else if(H1_m >= 0 && H1_p <= 0) 
		{ 
			if(fabs(H1_p) >= fabs(H1_m)) 
			{ 
				phi_x = phi_x_p; 
			} 
			else 
			{ 
				phi_x = phi_x_m; 
			} 
		} 
		return phi_x; 
	} 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
	// for (+) derivative v1 = phi_i, for (-) derivative v2 = phi_i 
	inline float DerENO1Minus(float v1, float v2, float dx) 
	{ 
		return (v2-v1)/dx; 
	} 
	inline float DerENO1Plus(float v1, float v2, float dx) 
	{ 
		return (v2-v1)/dx; 
	} 
 
	// v3 = phi_i 
	inline float DerENO2Minus(float v1, float v2, float v3, float v4, float dx) 
	{ 
		float D1_1 = (v2-v1); 
		float D1_2 = (v3-v2); 
		float D1_3 = (v4-v3); 
 
		float D2_1 = 0.5f*(D1_2-D1_1); 
		float D2_2 = 0.5f*(D1_3-D1_2); 
 
		float absD2_1 = fabs(D2_1); 
		float absD2_2 = fabs(D2_2); 
		float Q2p; 
		if(absD2_1 <= absD2_2) 
		{ 
			Q2p = D2_1; 
		} 
		else 
		{ 
			Q2p = D2_2; 
		} 
 
		return (D1_2 + Q2p)/dx; 
	} 
 
	// v2 = phi_i 
	inline float DerENO2Plus(float v1, float v2, float v3, float v4, float dx) 
	{ 
		float D1_1 = (v2-v1); 
		float D1_2 = (v3-v2); 
		float D1_3 = (v4-v3); 
 
		float D2_1 = 0.5f*(D1_2-D1_1); 
		float D2_2 = 0.5f*(D1_3-D1_2); 
 
		float absD2_1 = fabs(D2_1); 
		float absD2_2 = fabs(D2_2); 
		float Q2p; 
		if(absD2_1 <= absD2_2) 
		{ 
			Q2p = -D2_1; 
		} 
		else 
		{ 
			Q2p = -D2_2; 
		} 
 
		return (D1_2 + Q2p)/dx; 
	} 
 
	// v4 = phi_i 
	inline float DerENO3Minus(float v1, float v2, float v3, float v4, float v5, float v6, float dx) 
	{ 
		float D1_1 = (v2-v1); 
		float D1_2 = (v3-v2); 
		float D1_3 = (v4-v3); 
		float D1_4 = (v5-v4); 
		float D1_5 = (v6-v5); 
 
		float D2_1 = 0.5f*(D1_2-D1_1); 
		float D2_2 = 0.5f*(D1_3-D1_2); 
		float D2_3 = 0.5f*(D1_4-D1_3); 
		float D2_4 = 0.5f*(D1_5-D1_4); 
 
		float D3_1 = (D2_2-D2_1)/3; 
		float D3_2 = (D2_3-D2_2)/3; 
		float D3_3 = (D2_4-D2_3)/3; 
 
		float absD2_2 = fabs(D2_2); 
		float absD2_3 = fabs(D2_3); 
		float Q2p; 
		float Q3p; 
		if(absD2_2 <= absD2_3) 
		{ 
			Q2p = D2_2; 
			float cstar; 
			float absD3_1 = fabs(D3_1); 
			float absD3_2 = fabs(D3_2); 
			if(absD3_1 <= absD3_2) 
			{ 
				cstar = D3_1; 
			} 
			else 
			{ 
				cstar = D3_2; 
			} 
			Q3p = cstar*2; 
		} 
		else 
		{ 
			Q2p = D2_3; 
			float cstar; 
			float absD3_2 = fabs(D3_2); 
			float absD3_3 = fabs(D3_3); 
			if(absD3_2 <= absD3_3) 
			{ 
				cstar = D3_2; 
			} 
			else 
			{ 
				cstar = D3_3; 
			} 
			Q3p = -cstar; 
		} 
 
		return (D1_3 + Q2p + Q3p)/dx; 
	} 
 
	// v3 = phi_i 
	inline float DerENO3Plus(float v1, float v2, float v3, float v4, float v5, float v6, float dx) 
	{ 
		float D1_1 = (v2-v1); 
		float D1_2 = (v3-v2); 
		float D1_3 = (v4-v3); 
		float D1_4 = (v5-v4); 
		float D1_5 = (v6-v5); 
 
		float D2_1 = 0.5f*(D1_2-D1_1); 
		float D2_2 = 0.5f*(D1_3-D1_2); 
		float D2_3 = 0.5f*(D1_4-D1_3); 
		float D2_4 = 0.5f*(D1_5-D1_4); 
 
		float D3_1 = (D2_2-D2_1)/3; 
		float D3_2 = (D2_3-D2_2)/3; 
		float D3_3 = (D2_4-D2_3)/3; 
 
		float absD2_2 = fabs(D2_2); 
		float absD2_3 = fabs(D2_3); 
		float Q2p; 
		float Q3p; 
		if(absD2_2 <= absD2_3) 
		{ 
			Q2p = -D2_2; 
			float cstar; 
			float absD3_1 = fabs(D3_1); 
			float absD3_2 = fabs(D3_2); 
			if(absD3_1 <= absD3_2) 
			{ 
				cstar = D3_1; 
			} 
			else 
			{ 
				cstar = D3_2; 
			} 
			Q3p = -cstar; 
		} 
		else 
		{ 
			Q2p = -D2_3; 
			float cstar; 
			float absD3_2 = fabs(D3_2); 
			float absD3_3 = fabs(D3_3); 
			if(absD3_2 <= absD3_3) 
			{ 
				cstar = D3_2; 
			} 
			else 
			{ 
				cstar = D3_3; 
			} 
			Q3p = cstar*2; 
		} 
 
		return (D1_3 + Q2p + Q3p)/dx; 
	} 
 
 
	// vv4 = phi_i 
	inline float DerWENOMinus(float vv1, float vv2, float vv3, float vv4, float vv5, float vv6, float dx) 
	{ 
		float v1 = (vv2-vv1)/dx; 
		float v2 = (vv3-vv2)/dx; 
		float v3 = (vv4-vv3)/dx; 
		float v4 = (vv5-vv4)/dx; 
		float v5 = (vv6-vv5)/dx; 
 
		float data_x_1 = (2*v1 - 7*v2 + 11*v3)/6; 
		float data_x_2 = (-v2 + 5*v3 + 2*v4)/6; 
		float data_x_3 = (2*v3 + 5*v4 - v5)/6; 
 
		float dummy1 = v1-2*v2+v3; 
		float dummy2 = v1-4*v2+3*v3; 
		float dummy3 = v2-2*v3+v4; 
		float dummy4 = v2-v4; 
		float dummy5 = v3-2*v4+v5; 
		float dummy6 = 3*v3-4*v4+v5; 
		float dummy_c = 13/12; 
		double S1 = dummy_c*dummy1*dummy1 + 0.25*dummy2*dummy2; 
		double S2 = dummy_c*dummy3*dummy3 + 0.25*dummy4*dummy4; 
		double S3 = dummy_c*dummy5*dummy5 + 0.25*dummy6*dummy6; 
 
		float m = v1*v1; 
		float temp = v2*v2; 
		m = m>=temp?m:temp; 
		temp = v3*v3; 
		m = m>=temp?m:temp; 
		temp = v4*v4; 
		m = m>=temp?m:temp; 
		temp = v5*v5; 
		m = m>=temp?m:temp; 
		double epsil = m*1e-6 + 1e-199; 
 
		double alpha1 = 0.1/((S1+epsil)*(S1+epsil)); 
		double alpha2 = 0.6/((S2+epsil)*(S2+epsil)); 
		double alpha3 = 0.3/((S3+epsil)*(S3+epsil)); 
		return (float)((alpha1*data_x_1 + alpha2*data_x_2 + alpha3*data_x_3)/(alpha1+alpha2+alpha3)); 
	} 
 
	// v3 = phi_i 
	inline float DerWENOPlus(float vv1, float vv2, float vv3, float vv4, float vv5, float vv6, float dx) 
	{ 
 
		float v5 = (vv2-vv1)/dx; 
		float v4 = (vv3-vv2)/dx; 
		float v3 = (vv4-vv3)/dx; 
		float v2 = (vv5-vv4)/dx; 
		float v1 = (vv6-vv5)/dx; 
 
		float data_x_1 = (2*v1 - 7*v2 + 11*v3)/6; 
		float data_x_2 = (-v2 + 5*v3 + 2*v4)/6; 
		float data_x_3 = (2*v3 + 5*v4 - v5)/6; 
 
		float dummy1 = v1-2*v2+v3; 
		float dummy2 = v1-4*v2+3*v3; 
		float dummy3 = v2-2*v3+v4; 
		float dummy4 = v2-v4; 
		float dummy5 = v3-2*v4+v5; 
		float dummy6 = 3*v3-4*v4+v5; 
		float dummy_c = 13/12; 
		double S1 = dummy_c*dummy1*dummy1 + 0.25*dummy2*dummy2; 
		double S2 = dummy_c*dummy3*dummy3 + 0.25*dummy4*dummy4; 
		double S3 = dummy_c*dummy5*dummy5 + 0.25*dummy6*dummy6; 
 
		float m = v1*v1; 
		float temp = v2*v2; 
		m = m>=temp?m:temp; 
		temp = v3*v3; 
		m = m>=temp?m:temp; 
		temp = v4*v4; 
		m = m>=temp?m:temp; 
		temp = v5*v5; 
		m = m>=temp?m:temp; 
		double epsil = m*1e-6 + 1e-199; 
 
		double alpha1 = 0.1/((S1+epsil)*(S1+epsil)); 
		double alpha2 = 0.6/((S2+epsil)*(S2+epsil)); 
		double alpha3 = 0.3/((S3+epsil)*(S3+epsil)); 
		return (float)((alpha1*data_x_1 + alpha2*data_x_2 + alpha3*data_x_3)/(alpha1+alpha2+alpha3)); 
	} 
 
 
 
 
 
 
	// for (+) derivative v1 = phi_i, for (-) derivative v2 = phi_i 
	inline double DerENO1Minus(double v1, double v2, double dx) 
	{ 
		return (v2-v1)/dx; 
	} 
	inline double DerENO1Plus(double v1, double v2, double dx) 
	{ 
		return (v2-v1)/dx; 
	} 
 
	// v3 = phi_i 
	inline double DerENO2Minus(double v1, double v2, double v3, double v4, double dx) 
	{ 
		double D1_1 = (v2-v1); 
		double D1_2 = (v3-v2); 
		double D1_3 = (v4-v3); 
 
		double D2_1 = 0.5f*(D1_2-D1_1); 
		double D2_2 = 0.5f*(D1_3-D1_2); 
 
		double absD2_1 = fabs(D2_1); 
		double absD2_2 = fabs(D2_2); 
		double Q2p; 
		if(absD2_1 <= absD2_2) 
		{ 
			Q2p = D2_1; 
		} 
		else 
		{ 
			Q2p = D2_2; 
		} 
 
		return (D1_2 + Q2p)/dx; 
	} 
 
	// v2 = phi_i 
	inline double DerENO2Plus(double v1, double v2, double v3, double v4, double dx) 
	{ 
		double D1_1 = (v2-v1); 
		double D1_2 = (v3-v2); 
		double D1_3 = (v4-v3); 
 
		double D2_1 = 0.5f*(D1_2-D1_1); 
		double D2_2 = 0.5f*(D1_3-D1_2); 
 
		double absD2_1 = fabs(D2_1); 
		double absD2_2 = fabs(D2_2); 
		double Q2p; 
		if(absD2_1 <= absD2_2) 
		{ 
			Q2p = -D2_1; 
		} 
		else 
		{ 
			Q2p = -D2_2; 
		} 
 
		return (D1_2 + Q2p)/dx; 
	} 
 
	// v4 = phi_i 
	inline double DerENO3Minus(double v1, double v2, double v3, double v4, double v5, double v6, double dx) 
	{ 
		double D1_1 = (v2-v1); 
		double D1_2 = (v3-v2); 
		double D1_3 = (v4-v3); 
		double D1_4 = (v5-v4); 
		double D1_5 = (v6-v5); 
 
		double D2_1 = 0.5f*(D1_2-D1_1); 
		double D2_2 = 0.5f*(D1_3-D1_2); 
		double D2_3 = 0.5f*(D1_4-D1_3); 
		double D2_4 = 0.5f*(D1_5-D1_4); 
 
		double D3_1 = (D2_2-D2_1)/3; 
		double D3_2 = (D2_3-D2_2)/3; 
		double D3_3 = (D2_4-D2_3)/3; 
 
		double absD2_2 = fabs(D2_2); 
		double absD2_3 = fabs(D2_3); 
		double Q2p; 
		double Q3p; 
		if(absD2_2 <= absD2_3) 
		{ 
			Q2p = D2_2; 
			double cstar; 
			double absD3_1 = fabs(D3_1); 
			double absD3_2 = fabs(D3_2); 
			if(absD3_1 <= absD3_2) 
			{ 
				cstar = D3_1; 
			} 
			else 
			{ 
				cstar = D3_2; 
			} 
			Q3p = cstar*2; 
		} 
		else 
		{ 
			Q2p = D2_3; 
			double cstar; 
			double absD3_2 = fabs(D3_2); 
			double absD3_3 = fabs(D3_3); 
			if(absD3_2 <= absD3_3) 
			{ 
				cstar = D3_2; 
			} 
			else 
			{ 
				cstar = D3_3; 
			} 
			Q3p = -cstar; 
		} 
 
		return (D1_3 + Q2p + Q3p)/dx; 
	} 
 
	// v3 = phi_i 
	inline double DerENO3Plus(double v1, double v2, double v3, double v4, double v5, double v6, double dx) 
	{ 
		double D1_1 = (v2-v1); 
		double D1_2 = (v3-v2); 
		double D1_3 = (v4-v3); 
		double D1_4 = (v5-v4); 
		double D1_5 = (v6-v5); 
 
		double D2_1 = 0.5f*(D1_2-D1_1); 
		double D2_2 = 0.5f*(D1_3-D1_2); 
		double D2_3 = 0.5f*(D1_4-D1_3); 
		double D2_4 = 0.5f*(D1_5-D1_4); 
 
		double D3_1 = (D2_2-D2_1)/3; 
		double D3_2 = (D2_3-D2_2)/3; 
		double D3_3 = (D2_4-D2_3)/3; 
 
		double absD2_2 = fabs(D2_2); 
		double absD2_3 = fabs(D2_3); 
		double Q2p; 
		double Q3p; 
		if(absD2_2 <= absD2_3) 
		{ 
			Q2p = -D2_2; 
			double cstar; 
			double absD3_1 = fabs(D3_1); 
			double absD3_2 = fabs(D3_2); 
			if(absD3_1 <= absD3_2) 
			{ 
				cstar = D3_1; 
			} 
			else 
			{ 
				cstar = D3_2; 
			} 
			Q3p = -cstar; 
		} 
		else 
		{ 
			Q2p = -D2_3; 
			double cstar; 
			double absD3_2 = fabs(D3_2); 
			double absD3_3 = fabs(D3_3); 
			if(absD3_2 <= absD3_3) 
			{ 
				cstar = D3_2; 
			} 
			else 
			{ 
				cstar = D3_3; 
			} 
			Q3p = cstar*2; 
		} 
 
		return (D1_3 + Q2p + Q3p)/dx; 
	} 
 
 
	// vv4 = phi_i 
	inline double DerWENOMinus(double vv1, double vv2, double vv3, double vv4, double vv5, double vv6, double dx) 
	{ 
		double v1 = (vv2-vv1)/dx; 
		double v2 = (vv3-vv2)/dx; 
		double v3 = (vv4-vv3)/dx; 
		double v4 = (vv5-vv4)/dx; 
		double v5 = (vv6-vv5)/dx; 
 
		double data_x_1 = (2*v1 - 7*v2 + 11*v3)/6; 
		double data_x_2 = (-v2 + 5*v3 + 2*v4)/6; 
		double data_x_3 = (2*v3 + 5*v4 - v5)/6; 
 
		double dummy1 = v1-2*v2+v3; 
		double dummy2 = v1-4*v2+3*v3; 
		double dummy3 = v2-2*v3+v4; 
		double dummy4 = v2-v4; 
		double dummy5 = v3-2*v4+v5; 
		double dummy6 = 3*v3-4*v4+v5; 
		double dummy_c = 13/12; 
		double S1 = dummy_c*dummy1*dummy1 + 0.25*dummy2*dummy2; 
		double S2 = dummy_c*dummy3*dummy3 + 0.25*dummy4*dummy4; 
		double S3 = dummy_c*dummy5*dummy5 + 0.25*dummy6*dummy6; 
 
		double m = v1*v1; 
		double temp = v2*v2; 
		m = m>=temp?m:temp; 
		temp = v3*v3; 
		m = m>=temp?m:temp; 
		temp = v4*v4; 
		m = m>=temp?m:temp; 
		temp = v5*v5; 
		m = m>=temp?m:temp; 
		double epsil = m*1e-6 + 1e-199; 
 
		double alpha1 = 0.1/((S1+epsil)*(S1+epsil)); 
		double alpha2 = 0.6/((S2+epsil)*(S2+epsil)); 
		double alpha3 = 0.3/((S3+epsil)*(S3+epsil)); 
		return (double)((alpha1*data_x_1 + alpha2*data_x_2 + alpha3*data_x_3)/(alpha1+alpha2+alpha3)); 
	} 
 
	// v3 = phi_i 
	inline double DerWENOPlus(double vv1, double vv2, double vv3, double vv4, double vv5, double vv6, double dx) 
	{ 
 
		double v5 = (vv2-vv1)/dx; 
		double v4 = (vv3-vv2)/dx; 
		double v3 = (vv4-vv3)/dx; 
		double v2 = (vv5-vv4)/dx; 
		double v1 = (vv6-vv5)/dx; 
 
		double data_x_1 = (2*v1 - 7*v2 + 11*v3)/6; 
		double data_x_2 = (-v2 + 5*v3 + 2*v4)/6; 
		double data_x_3 = (2*v3 + 5*v4 - v5)/6; 
 
		double dummy1 = v1-2*v2+v3; 
		double dummy2 = v1-4*v2+3*v3; 
		double dummy3 = v2-2*v3+v4; 
		double dummy4 = v2-v4; 
		double dummy5 = v3-2*v4+v5; 
		double dummy6 = 3*v3-4*v4+v5; 
		double dummy_c = 13/12; 
		double S1 = dummy_c*dummy1*dummy1 + 0.25*dummy2*dummy2; 
		double S2 = dummy_c*dummy3*dummy3 + 0.25*dummy4*dummy4; 
		double S3 = dummy_c*dummy5*dummy5 + 0.25*dummy6*dummy6; 
 
		double m = v1*v1; 
		double temp = v2*v2; 
		m = m>=temp?m:temp; 
		temp = v3*v3; 
		m = m>=temp?m:temp; 
		temp = v4*v4; 
		m = m>=temp?m:temp; 
		temp = v5*v5; 
		m = m>=temp?m:temp; 
		double epsil = m*1e-6 + 1e-199; 
 
		double alpha1 = 0.1/((S1+epsil)*(S1+epsil)); 
		double alpha2 = 0.6/((S2+epsil)*(S2+epsil)); 
		double alpha3 = 0.3/((S3+epsil)*(S3+epsil)); 
		return (double)((alpha1*data_x_1 + alpha2*data_x_2 + alpha3*data_x_3)/(alpha1+alpha2+alpha3)); 
	} 
 
 
 
 
 
 
 
 
};