public class NumericalAreaIntegrator implements IAreaIntegrator {

	double m_lfAreaEps;
	IScalarField m_scalarField;
	
	public NumericalAreaIntegrator() {
		m_lfAreaEps = 0.0001;
		m_scalarField = new ScalarFieldSumGauss();
	}
	
	public NumericalAreaIntegrator(IScalarField scalarField) {
		m_lfAreaEps = 0.0001;
		m_scalarField = scalarField.makeCopy();
	}
	
	public NumericalAreaIntegrator(IScalarField scalarField, double lfEps) {
		m_lfAreaEps = lfEps;
		m_scalarField = scalarField.makeCopy();
	}
	
	public NumericalAreaIntegrator(NumericalAreaIntegrator src) {
		m_lfAreaEps = src.m_lfAreaEps;
		m_scalarField = src.m_scalarField.makeCopy();
	}
	
	public void setAreaEps(double lfEps) {
		m_lfAreaEps = lfEps;
	}
	
	public void setScalarField(IScalarField scalarField) {
		m_scalarField = scalarField.makeCopy();
	}
	
	void test1() {
		double [] arrLfCtr = new double[2];
		double [] arrLfP1 = new double[2];
		double [] arrLfP2 = new double[2];
		for(int i = 0; i < 2; ++i) {
			arrLfCtr[i] = Math.random();
			arrLfP1[i] = Math.random();
			arrLfP2[i] = Math.random();
		}
		double lfX1 = arrLfP1[0] - arrLfCtr[0];
		double lfY1 = arrLfP1[1] - arrLfCtr[1];
		double lfX2 = arrLfP2[0] - arrLfCtr[0];
		double lfY2 = arrLfP2[1] - arrLfCtr[1];
		
		Double objLfAnalyticalArea = lfX1*lfY2 - lfY1*lfX2;
		objLfAnalyticalArea = 0.5*objLfAnalyticalArea;
		
		Double objLfNumericalArea = 
			integrateTriangle(arrLfCtr, arrLfP1, arrLfP2);
		
		
		System.out.println("NumInt Tri Test : Numerical are is " +
				objLfNumericalArea.toString() + 
				" Analytical area is " + 
				objLfAnalyticalArea.toString());
		
		System.out.println("Percent error is " +
				(new Double(200*Math.abs(objLfNumericalArea -
						objLfAnalyticalArea)/(
						Math.abs(objLfNumericalArea) +
						Math.abs(objLfAnalyticalArea))).toString()));
		
	}
	
	void test2() {
		double [] arrLfCtr = new double[2];
		double [] arrLfP1 = new double[2];
		double [] arrLfP2 = new double[2];
		for(int i = 0; i < 2; ++i) {
			arrLfCtr[i] = Math.random();
			arrLfP1[i] = 1.0+Math.random();
			arrLfP2[i] = 1.0+Math.random();
		}
		double lfRad = Math.random();
		double lfX1 = arrLfP1[0] - arrLfCtr[0];
		double lfY1 = arrLfP1[1] - arrLfCtr[1];
		double lfX2 = arrLfP2[0] - arrLfCtr[0];
		double lfY2 = arrLfP2[1] - arrLfCtr[1];
		double lfRad1 = 
			Math.sqrt(lfX1*lfX1+
					lfY1*lfY1);
		double lfRad2 =
			Math.sqrt(lfX2*lfX2+
					lfY2*lfY2);
		double lfRadProd = lfRad1*lfRad2;
		double lfCross = lfX1*lfY2 - lfY1*lfX2;
		double lfDot = lfX1*lfX2 + lfY1*lfY2;
		double lfTheta = Math.atan2(lfCross/lfRadProd,lfDot/lfRadProd);
		Double objLfAnalyticalArea = 0.5*lfTheta*lfRad*lfRad;
		
		double lfOldArea = m_lfAreaEps;
		for(int j = 0; j < 4; ++j) {
			m_lfAreaEps = 0.01*Math.pow(0.1, j);
		
			java.util.Date date1 = new java.util.Date();
			Double objLfNumericalArea = 
				this.integrateClippedTriangle(arrLfCtr, arrLfP1, 
						arrLfP2, lfRad, arrLfCtr);
			java.util.Date date2 = new java.util.Date();
			
			long lnTimeDiff = date2.getTime() - date1.getTime();

			System.out.println("Time elapsed is " + 
					(new Long(lnTimeDiff)).toString() + 
					" area eps is "+
					(new Double(m_lfAreaEps)).toString());
			System.out.println("NumInt Clipped Tri Test : Numerical are is " +
					objLfNumericalArea.toString() + 
					" Analytical area is " + 
					objLfAnalyticalArea.toString());

			System.out.println("Percent error is " +
					(new Double(200*Math.abs(objLfNumericalArea -
							objLfAnalyticalArea)/(
									Math.abs(objLfNumericalArea) +
									Math.abs(objLfAnalyticalArea))).toString())+
									"\tAngle is "+
									(new Double(lfTheta)).toString());
		}
		m_lfAreaEps = lfOldArea;
		System.out.println("Radii are " +
				(new Double(lfRad1)).toString() +"," +
				(new Double(lfRad2)).toString() + "," +
				(new Double(lfRad)).toString());
	}
	
	public void unitTest() {
		IScalarField fldBackup = m_scalarField;
		m_scalarField = new ScalarFieldUtils.LinearScalarField(1.0);
		
		int nTestIter = 0;
		for(nTestIter = 0; nTestIter < 10; ++nTestIter) {
			test1();
		}
		
		for(nTestIter = 0; nTestIter < 10; ++nTestIter) {
			test2();
		}

		m_scalarField = fldBackup;
		System.out.println("Error : Haven't written all unit tests");
	}
	
	public double integrateTriangle(double[] arrLfCenter, double[] arrLfPt1,
			double[] arrLfPt2) {
		double v1x = (arrLfPt1[0]-arrLfCenter[0]);
		double v1y = (arrLfPt1[1]-arrLfCenter[1]);
		double v2x = (arrLfPt2[0]-arrLfCenter[0]);
		double v2y = (arrLfPt2[1]-arrLfCenter[1]);
		double lfArea = v1x*v2y-v1y*v2x;
		if(Math.abs(lfArea) < m_lfAreaEps) {
			return lfArea*(m_scalarField.getVal(arrLfCenter)+
					m_scalarField.getVal(arrLfPt1) +
					m_scalarField.getVal(arrLfPt2))/6.0;
		}
		
		// Okay so recursive integration technique is not hte most efficient
		// will replace later
		double arrLfMidC1[] = new double[2];
		arrLfMidC1[0] = 0.5*(arrLfCenter[0]+arrLfPt1[0]);
		arrLfMidC1[1] = 0.5*(arrLfCenter[1]+arrLfPt1[1]);
		double arrLfMidC2[] = new double[2];
		arrLfMidC2[0] = 0.5*(arrLfCenter[0]+arrLfPt2[0]);
		arrLfMidC2[1] = 0.5*(arrLfCenter[1]+arrLfPt2[1]);
		double arrLfMid12[] = new double[2];
		arrLfMid12[0] = 0.5*(arrLfPt1[0]+arrLfPt2[0]);
		arrLfMid12[1] = 0.5*(arrLfPt1[1]+arrLfPt2[1]);
		
		//               1
		//       c1
		// c             12
		//       c2
		//               2
		// c-1-2 = 
		// c-c1-c2 + c1-1-12 + c2-12-2 + c2-c1-12
		return integrateTriangle(arrLfCenter,arrLfMidC1,arrLfMidC2) +
			integrateTriangle(arrLfMidC1,arrLfPt1,arrLfMid12) +
			integrateTriangle(arrLfMidC2,arrLfMid12,arrLfPt2) +
			integrateTriangle(arrLfMidC2,arrLfMidC1,arrLfMid12);
	}

	public double integrateClippedTriangle(double[] arrLfCenter, double[] arrLfPt1,
			double[] arrLfPt2, double lfRadius, 
			double [] arrLfCircleCenter) {
		double v1x = (arrLfPt1[0]-arrLfCenter[0]);
		double v1y = (arrLfPt1[1]-arrLfCenter[1]);
		double v2x = (arrLfPt2[0]-arrLfCenter[0]);
		double v2y = (arrLfPt2[1]-arrLfCenter[1]);
		double lfArea = v1x*v2y-v1y*v2x;

		double lfTmpX = arrLfCenter[0] - arrLfCircleCenter[0];
		double lfTmpY = arrLfCenter[1] - arrLfCircleCenter[1];
		double lfDistCSqrd = lfTmpX*lfTmpX + lfTmpY*lfTmpY;
			
		lfTmpX = arrLfPt1[0] - arrLfCircleCenter[0];
		lfTmpY = arrLfPt1[1] - arrLfCircleCenter[1];
		double lfDist1Sqrd = lfTmpX*lfTmpX + lfTmpY*lfTmpY;
		
		lfTmpX = arrLfPt2[0] - arrLfCircleCenter[0];
		lfTmpY = arrLfPt2[1] - arrLfCircleCenter[1];
		double lfDist2Sqrd = lfTmpX*lfTmpX + lfTmpY*lfTmpY;
		
		double lfRadSqrd = lfRadius*lfRadius;
		boolean bCin = (lfDistCSqrd < lfRadSqrd);
		boolean b1in = (lfDist1Sqrd < lfRadSqrd);
		boolean b2in = (lfDist2Sqrd < lfRadSqrd);
		
		if(!bCin && !b1in && !b2in) {
			lfTmpX = arrLfPt1[0] - arrLfCenter[0];
			lfTmpY = arrLfPt1[1] - arrLfCenter[1];
			double lfTriRad = lfTmpX*lfTmpX + lfTmpY*lfTmpY;
			lfTmpX = arrLfPt2[0] - arrLfCenter[0];
			lfTmpY = arrLfPt2[1] - arrLfCenter[1];
			lfTriRad = Math.max(lfTriRad,lfTmpX*lfTmpX + lfTmpY*lfTmpY);
			lfTmpX = arrLfPt2[0] - arrLfPt1[0];
			lfTmpY = arrLfPt2[1] - arrLfPt1[1];
			lfTriRad = Math.max(lfTriRad,lfTmpX*lfTmpX + lfTmpY*lfTmpY);
			lfTriRad = 0.5*Math.sqrt(lfTriRad);
			
			if(Math.sqrt(lfDistCSqrd) > lfTriRad + lfRadius &&
					Math.sqrt(lfDist1Sqrd) > lfTriRad + lfRadius &&
					Math.sqrt(lfDist2Sqrd) > lfTriRad + lfRadius) {
				return 0.0;
			}
		}
		
		
//		if(!bCin && !b1in && !b2in) {
//			return 0.0;
//		}
		
		// really crude hack, fix later. --------------------v
		if(Math.abs(lfArea) < m_lfAreaEps*Math.min(lfRadSqrd,1.0)) {
			double lfSum = 0.0;
			if(bCin) {
				lfSum += m_scalarField.getVal(arrLfCenter);
			}
			if(b1in) {
				lfSum += m_scalarField.getVal(arrLfPt1);
			}
			if(b2in) {
				lfSum += m_scalarField.getVal(arrLfPt2);
			}
			
			return lfArea*lfSum/6.0;
		}
		
		// Okay so recursive integration technique is not hte most efficient
		// will replace later
		double arrLfMidC1[] = new double[2];
		arrLfMidC1[0] = 0.5*(arrLfCenter[0]+arrLfPt1[0]);
		arrLfMidC1[1] = 0.5*(arrLfCenter[1]+arrLfPt1[1]);
		double arrLfMidC2[] = new double[2];
		arrLfMidC2[0] = 0.5*(arrLfCenter[0]+arrLfPt2[0]);
		arrLfMidC2[1] = 0.5*(arrLfCenter[1]+arrLfPt2[1]);
		double arrLfMid12[] = new double[2];
		arrLfMid12[0] = 0.5*(arrLfPt1[0]+arrLfPt2[0]);
		arrLfMid12[1] = 0.5*(arrLfPt1[1]+arrLfPt2[1]);
		
		//               1
		//       c1
		// c             12
		//       c2
		//               2
		// c-1-2 = 
		// c-c1-c2 + c1-1-12 + c2-12-2 + c2-c1-12
		return integrateClippedTriangle(arrLfCenter,arrLfMidC1,
				arrLfMidC2, lfRadius, arrLfCircleCenter) +
			integrateClippedTriangle(arrLfMidC1,arrLfPt1,
					arrLfMid12, lfRadius, arrLfCircleCenter) +
			integrateClippedTriangle(arrLfMidC2,arrLfMid12,
					arrLfPt2, lfRadius, arrLfCircleCenter) +
			integrateClippedTriangle(arrLfMidC2,arrLfMidC1,
					arrLfMid12, lfRadius, arrLfCircleCenter);
	}
	
	
	// override?
	public double integrateWedge(double[] arrLfCenter, double[] arrLfPt1,
			double lfAngle) {
		double lfX = arrLfPt1[0] - arrLfCenter[0];
		double lfY = arrLfPt1[1] - arrLfCenter[1];
		double lfRad = Math.sqrt(lfX*lfX+lfY*lfY);
		double lfSum = 0.0;
		double lfCurrAngle = 0.0;
		double lfStep = Math.sqrt(m_lfAreaEps)/lfRad;
		double lfSinS = Math.sin(lfStep);
		double lfCosS = Math.cos(lfStep);
		double [] arrLfLastPt = new double[2];
		double [] arrLfNextPt = new double[2];
		double [] tmp = null;
		arrLfLastPt[0] = arrLfCenter[0] + lfX;
		arrLfLastPt[1] = arrLfCenter[1] + lfY;
		for(; lfCurrAngle + lfStep < lfAngle ; lfCurrAngle += lfStep) {
			double lfXNew = lfX*lfCosS - lfY*lfSinS;
			double lfYNew = lfX*lfSinS + lfY*lfCosS;
			lfX = lfXNew;
			lfY = lfYNew;
			arrLfNextPt[0] = arrLfCenter[0] + lfX;
			arrLfNextPt[1] = arrLfCenter[1] + lfY;
			
			lfSum += integrateTriangle(arrLfCenter, 
					arrLfLastPt, arrLfNextPt);
			tmp = arrLfLastPt; arrLfLastPt = arrLfNextPt;
			arrLfNextPt = tmp;
		}
		lfStep = lfAngle - lfCurrAngle;
		lfSinS = Math.sin(lfStep);
		lfCosS = Math.cos(lfStep);
		
		double lfXNew = lfX*lfCosS - lfY*lfSinS;
		double lfYNew = lfX*lfSinS + lfY*lfCosS;
		lfX = lfXNew;
		lfY = lfYNew;
		arrLfNextPt[0] = arrLfCenter[0] + lfX;
		arrLfNextPt[1] = arrLfCenter[1] + lfY;
		
		lfSum += integrateTriangle(arrLfCenter, 
				arrLfLastPt, arrLfNextPt);
	
		return lfSum;
	}

	public double integrateWedgeOffCenter(double[] arrLfCenter,
			double[] arrLfIntegrateCtr, double[] arrLfPt1, double lfAngle) {
		double lfX = arrLfPt1[0] - arrLfCenter[0];
		double lfY = arrLfPt1[1] - arrLfCenter[1];
		double lfRad = Math.sqrt(lfX*lfX+lfY*lfY);
		double lfSum = 0.0;
		double lfCurrAngle = 0.0;
		double lfStep = Math.sqrt(m_lfAreaEps)/lfRad;
		double lfSinS = Math.sin(lfStep);
		double lfCosS = Math.cos(lfStep);
		double [] arrLfLastPt = new double[2];
		double [] arrLfNextPt = new double[2];
		double [] tmp = null;
		arrLfLastPt[0] = arrLfCenter[0] + lfX;
		arrLfLastPt[1] = arrLfCenter[1] + lfY;
		for(; lfCurrAngle + lfStep < lfAngle ; lfCurrAngle += lfStep) {
			double lfXNew = lfX*lfCosS - lfY*lfSinS;
			double lfYNew = lfX*lfSinS + lfY*lfCosS;
			lfX = lfXNew;
			lfY = lfYNew;
			arrLfNextPt[0] = arrLfCenter[0] + lfX;
			arrLfNextPt[1] = arrLfCenter[1] + lfY;
			
			lfSum += integrateTriangle(arrLfIntegrateCtr, 
					arrLfLastPt, arrLfNextPt);
			tmp = arrLfLastPt; arrLfLastPt = arrLfNextPt;
			arrLfNextPt = tmp;
		}
		lfStep = lfAngle - lfCurrAngle;
		lfSinS = Math.sin(lfStep);
		lfCosS = Math.cos(lfStep);
		
		double lfXNew = lfX*lfCosS - lfY*lfSinS;
		double lfYNew = lfX*lfSinS + lfY*lfCosS;
		lfX = lfXNew;
		lfY = lfYNew;
		arrLfNextPt[0] = arrLfCenter[0] + lfX;
		arrLfNextPt[1] = arrLfCenter[1] + lfY;
		
		lfSum += integrateTriangle(arrLfIntegrateCtr, 
				arrLfLastPt, arrLfNextPt);
	
		return lfSum;	}

}