#include <cstdio>
#include <cstring>
#include <cmath>
#include <ctime>

#include "../Bamgx.h"
#include "../QuadTree.h"
#include "../SetOfE4.h"

#include "../../shared/shared.h"
#include "../../include/macros.h"
#include "../../toolkits/toolkits.h"

#undef __FUNCT__ 
#define __FUNCT__ "MatVVP2x2"

namespace bamg {

	/*Constructor*/
	/*FUNCTION MatVVP2x2::MatVVP2x2(const MetricAnIso M){{{1*/
	MatVVP2x2::MatVVP2x2(const MetricAnIso M){
		double a11=M.a11,a21=M.a21,a22=M.a22;
		const double eps = 1.e-5;
		double c11 = a11*a11, c22 = a22*a22, c21= a21*a21;
		double b=-a11-a22,c=-c21+a11*a22;
		double   delta = b*b - 4 * c ;
		double n2=(c11+c22+c21);
		if ( n2 < 1e-30) 
		 lambda1=lambda2=0,v.x=1,v.y=0;
		else if (delta < eps*n2)
		  { 
			lambda1=lambda2=-b/2, v.x=1,v.y=0;
		  }
		else 
		  {  //    ---  construction  de 2 vecteur dans (Im ( A - D(i) Id) ortogonal 
			delta = sqrt(delta);
			lambda1 = (-b-delta)/2.0,lambda2 = (-b+delta)/2.0;
			double v0 = a11-lambda1, v1 = a21,v2 = a22 - lambda1;
			double s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;

			if(s1 < s0)
			 s0=sqrt(s0),v.x=v1/s0,v.y=-v0/s0;
			else
			 s1=sqrt(s1),v.x=v2/s1,v.y=-v1/s1;
		  };
	}
	/*}}}1*/

}