//    1, NORDER =  1, precision 1, Zienkiewicz #1.
//    2, NORDER =  3, precision 2, Strang and Fix formula #1.
//    3, NORDER =  3, precision 2, Strang and Fix formula #2, Zienkiewicz #2.
//    4, NORDER =  4, precision 3, Strang and Fix formula #3, Zienkiewicz #3.
//    5, NORDER =  6, precision 3, Strang and Fix formula #4.
//    6, NORDER =  6, precision 3, Stroud formula T2:3-1.
//    7, NORDER =  6, precision 4, Strang and Fix formula #5.
//    8, NORDER =  7, precision 4, Strang and Fix formula #6.
//    9, NORDER =  7, precision 5, Strang and Fix formula #7,
//       Stroud formula T2:5-1, Zienkiewicz #4, Schwarz Table 2.2.
//   10, NORDER =  9, precision 6, Strang and Fix formula #8.
//   11, NORDER = 12, precision 6, Strang and Fix formula #9.
//   12, NORDER = 13, precision 7, Strang and Fix formula #10.
//   13, NORDER =  7, precision ?.
//   14, NORDER = 16, precision 7, conical product Gauss, Stroud formula T2:7-1.
//   15, NORDER = 64, precision 15, triangular product Gauss rule.
//   16, NORDER = 19, precision 8, from CUBTRI, ACM TOMS #584.
//   17, NORDER = 19, precision 9, from TRIEX, Lyness and Jespersen.
//   18, NORDER = 28, precision 11, from TRIEX, Lyness and Jespersen.
//   19, NORDER = 37, precision 13, from ACM TOMS #706.

// Mutual cofficient of potential between two triangular patches
// with uniform charge
//
// This version is fully numerical and uses Strang and Fix, formula #3,
// 4 points, 3th order, for numerical quadrature (from Stroud)
//
// vertexes1 is an array of 3 3D point coordinates,
//           corresponding to the vertexes of the first triangle
//
// vertexes2 is an array of 3 3D point coordinates,
//           corresponding to the vertexes of the second triangle
//
double CPotential::Mutual_3thOrd_FullNum(C3DVector_float vertexes1[3],
						  C3DVector_float vertexes2[3], bool divideByArea)
{
	C3DVector t1side1, t1side2, t1side3, t1normal, t1x, t1y, t1z;
	C3DVector t2side1, t2side2, t2side3, t2normal, t2x, t2y, t2z;
	C2DVector t1local2D[3], t2local2D[3];
	C3DVector t1p3d[POTENTIAL_PTS_3TH_ORDER], t2p3d[POTENTIAL_PTS_3TH_ORDER];
	C3DVector t1transform[2], t2transform[2];
	double mod_t1side1, mod_t2side1, mod_t1normal, mod_t2normal;
	double t1jacobian, t2jacobian, pot, result, point_x, point_y;
	int i, j, k;

	// numerical 2D quadrature over a triangle
	// uses Stroud 4 points, 3th order
	//
	// initialize matrices of integration points and weights
	// (in homogeneous coordinates)
	static const double points[POTENTIAL_PTS_3TH_ORDER][3] = {
		{0.33333333333333333333333333333333, 0.33333333333333333333333333333333, 1.0},
		{0.6,                                0.2,                                1.0},
		{0.2,                                0.6,                                1.0},
		{0.2,                                0.2,                                1.0}};

	static const double weights2[POTENTIAL_W2_3TH_ORDER] =
		{ 0.31640625, -0.29296875,                         -0.29296875,                         -0.29296875,
		 -0.29296875,  0.27126736111111111111111111111111,  0.27126736111111111111111111111111,  0.27126736111111111111111111111111,
		 -0.29296875,  0.27126736111111111111111111111111,  0.27126736111111111111111111111111,  0.27126736111111111111111111111111,
		 -0.29296875,  0.27126736111111111111111111111111,  0.27126736111111111111111111111111,  0.27126736111111111111111111111111};

	//
	// calculate local coordinate frame for triangle 1
	//

	// compute side vectors of triangle
	t1side1 = vertexes1[1] - vertexes1[0];
	t1side2 = vertexes1[2] - vertexes1[1];
	t1side3 = vertexes1[0] - vertexes1[2];

	// compute unit vector normal to triangle plane (uses vector product)
	t1normal = CrossProd(t1side1, t1side2);
	mod_t1normal = Mod(t1normal);
	t1z = t1normal / mod_t1normal;

	// compute local x, y axes versors; x is along side1,
	// y is perp. to side1 in counter-clockwise direction
	mod_t1side1 = Mod(t1side1);
	t1x = t1side1 / mod_t1side1;
	// note that, since x and z are unit vectors, y is too
	t1y = CrossProd(t1x, t1z);

	//
	// calculate local coordinate frame for triangle 2
	//

	// compute side vectors of triangle
	t2side1 = vertexes2[1] - vertexes2[0];
	t2side2 = vertexes2[2] - vertexes2[1];
	t2side3 = vertexes2[0] - vertexes2[2];

	// compute unit vector normal to triangle plane (uses vector product)
	t2normal = CrossProd(t2side1, t2side2);
	mod_t2normal = Mod(t2normal);
	t2z = t2normal / mod_t2normal;

	// compute local x, y axes versors; x is along side1,
	// y is perp. to side1 in counter-clockwise direction
	mod_t2side1 = Mod(t2side1);
	t2x = t2side1 / mod_t2side1;
	// note that, since x and z are unit vectors, y is too
	t2y = CrossProd(t2x, t2z);


	// calculates the coordinates of the triangles in local 2D frame;
	// the origin is located on first vertex.
	// Remark: it is assumed that C2DVertex coordinates defualt to zero

	t1local2D[1].x = mod_t1side1;
	t1local2D[2].x = -DotProd(t1side3, t1x);
	t1local2D[2].y = -DotProd(t1side3, t1y);

	t2local2D[1].x = mod_t2side1;
	t2local2D[2].x = -DotProd(t2side3, t2x);
	t2local2D[2].y = -DotProd(t2side3, t2y);

	// calculates the coordinate transformation needed
	// to map the first triangle to a unit triangle
	// Basically, this is the following operation
	//
	//	static const double unittranmtx[9] = { 0, 1, 0,
	//	                                      1, 0, 0,
	//							             -1,-1, 1};
	// 	transform = local2D * unittranmtx;
	//
	// where unittranmtx is the inverse of [0,1,0;1,0,0;1,1,1],
	// the matrix representing the unit triangle
	//
	// Explanation: we are in 2D, but we use homogeneous coordinates.
	// The mapping is a linear transformation, i.e. to map a point (xt, yt, 1)
	// from a unit triangle to a point (x0, y0, 1) on a general triangle we use:
	// x0 = a*xt + b*yt + c
	// y0 = d*xt + e*yt + f
	// Therefore to map the corners of a unit triangle to the corners of a general
	// triangle we have:
	// [ x0 x1 x2 ]    [ a b c ]   [ 0 1 0 ]
	// [ y0 y1 y2 ]  = [ d e f ] * [ 1 0 0 ]
	// [  1  1  1 ]    [ 0 0 1 ]   [ 1 1 1 ]
	// where the last matrix is the matrix representing the vertex of the unit
	// triangle in homogeneous coordinates
	// We would like to know the matrix Trans = [a b c; d e f; 0 0 1] to be able to map
	// integration points in the unit triangle to integration points in our triangle.
	// The matrix is therefore Trans = MyTri * inv(UnitTri)

	// calculate transform for triangle 1 and 2
	// (only significative values are calculated, so last line of transform mtx
	// is trivial, moreover t1local2D[0] is (0.0, 0.0) and t1local2D[1].y is 0.0)

	t1transform[0].x = t1local2D[1].x - t1local2D[2].x;
	t1transform[0].y = -t1local2D[2].x;
	t1transform[0].z = t1local2D[2].x;
	t1transform[1].x = -t1local2D[2].y;
	t1transform[1].y = -t1local2D[2].y;
	t1transform[1].z = t1local2D[2].y;

	t2transform[0].x = t2local2D[1].x - t2local2D[2].x;
	t2transform[0].y = -t2local2D[2].x;
	t2transform[0].z = t2local2D[2].x;
	t2transform[1].x = -t2local2D[2].y;
	t2transform[1].y = -t2local2D[2].y;
	t2transform[1].z = t2local2D[2].y;


	// this is the jacobian of the transforms (a determinant)

	t1jacobian = t1transform[0].x*t1transform[1].y-
				t1transform[1].x*t1transform[0].y;

	t2jacobian = t2transform[0].x*t2transform[1].y-
				t2transform[1].x*t2transform[0].y;


	// calculate the transformed 3D points for triangle 1 and 2
	for(i=0; i<POTENTIAL_PTS_3TH_ORDER; i++) {

		// triangle 1

		// calculates the transformed 2D numerical integration point
		point_x = t1transform[0].x * points[i][0] +
			 	  t1transform[0].y * points[i][1] +
			 	  t1transform[0].z * points[i][2];
		point_y = t1transform[1].x * points[i][0] +
				  t1transform[1].y * points[i][1] +
				  t1transform[1].z * points[i][2];

		// calculate 3D point location
		t1p3d[i] = vertexes1[0] + t1x*point_x + t1y*point_y;

		// triangle 2

		// calculates the transformed 2D numerical integration point
		point_x = t2transform[0].x * points[i][0] +
			 	  t2transform[0].y * points[i][1] +
			 	  t2transform[0].z * points[i][2];
		point_y = t2transform[1].x * points[i][0] +
				  t2transform[1].y * points[i][1] +
				  t2transform[1].z * points[i][2];

		// calculate 3D point location
		t2p3d[i] = vertexes2[0] + t2x*point_x + t2y*point_y;
	}

	// calculate the integral
	result = 0;
	for(i=0, k=0; i<POTENTIAL_PTS_3TH_ORDER; i++) {
		for(j=0; j<POTENTIAL_PTS_3TH_ORDER; j++) {
			// calculate the kernel value between given integration points
			pot = 1/Mod(t1p3d[i] - t2p3d[j]);
			// sum up the result
			result = result + weights2[k]*pot;
			// increment weight counter
			k++;
		}
	}

	// do not forget the jacobian!
	// note that we should multiply by the unit triangle area, i.e. 1/2
	// (could be included in the weights; it is a side effect of how, in Stroud,
	// the transformed coordinates are calculated); however we do NOT perform
	// this operation here, since we should then divide by the two triangles area,
	// which include a divide by 2 operation, so they compensate
	result = result * t1jacobian * t2jacobian;

	// The modulus of the triangle normal is 2 * area of triangle
	if( divideByArea == true )
		result /= mod_t1normal*mod_t2normal;

	// return result
	return( result / FOUR_PI_TIMES_E0 );
}