import java.applet.Applet;
import java.awt.*;
import java.awt.event.*;
import java.applet.*;
import java.awt.geom.*;
import java.math.*;

public class EnergyEstimator {


    /**1 means eliminate.  We add a small positive fudge factor of
          2^{-40}  to make this  test more compatible with the interval 
          arithmetic version.  Adding this
         fudge factor makes the test HARDER to pass.**/

    public static int[] eliminate(double e2,double E,BoxConfiguration X) {
	double[] e1=energyLowerBound(E,X);
	int[] test={0,(int)(e1[1])};
	if(e1[0]>e2+Math.pow(.5,40)) test[0]=1;     //pass
	return(test);
    }

    /**Here is the main energy estimate.  The second number tells which index
       makes the biggest contribution to the error term.**/

    public static double[] energyLowerBound(double E,BoxConfiguration X) {
	double e1=vertexEnergy(E,X);
	double[] e2=errorTerm(E,X);
	double t=e1-e2[0];
	double[] d={t,e2[1]};
	return(d);
    }

    public static double[] energyLowerBoundDebug(double E,BoxConfiguration X) {
	double e1=vertexEnergy(E,X);
	double[] e2=errorTerm(E,X);
	double t=e1-e2[0];
	double[] d={t,e2[1]};
	return(d);
    }

    /**This routine computes the energy at the 128 vertices of a dyadic box and
       then takes the min.**/

    public static double vertexEnergy(double E,BoxConfiguration X) {
	double e=100000;
	double test=0;
	Complex[] Z=new Complex[4];
	for(int n=0;n<128;++n) {
	    Z=X.getConfiguration(n);
	    test=Energy.energy(E,Z);
	    if(e>test) e=test;
	}
	return(e);
    }

    /**Here is the error term.  This follows the notation in the paper.**/
    /**The whole thing.  The first number is the error term.
          The second number tells which box makes the biggest contribution**/

    public static double[] errorTerm(double E, BoxConfiguration X) {
	double total=0;
	double min=0;
	double test=0;
	int index=0;

	for(int i=0;i<4;++i) {
	    test=errorTermSingleIndex(E,X,i);
	    total=total+test;
	    if(min<test) {
		index=i;
		min=test;
	    }
	}

	double[] d={total,index};
	return(d);
    }




    /**For an individual index **/

    public static double errorTermSingleIndex(double E, BoxConfiguration X,int i) {
	double total=0;
	double d=X.B[i].delta();

	    for(int j=0;j<5;++j) {
		if(j!=i) {
		    total=total+d*d*epsilon(E,X.B[i],X.B[j]);
		}
	    }
	return(total);
    }


    /**epsilon from the paper**/

    public static double epsilon(double E,Box B1,Box B2) {
	double R=Separation.psiMin(B1,B2);
	if(R<=0) return(Math.pow(2,50));
	double l1=lambda1(E,B1,B2,R);
	double l2=lambda2(E,B1,B2,R);
	double l3=Math.max(0,l1)+l2;
	return(l3);
    }




    /**lambda 1 from the paper**/

    public static double lambda1(double E,Box B1,Box B2,double R) {

	if(B1.type==1) {
	    double t1=E*(E+1)/(8*Math.pow(R,E+2));
	    double t2=E*(E+2)/(32*Math.pow(R,E));
	    double t=t1-t2;
	    return(t);
	}

	if(B1.type==0) {
                      double t1=E*(E+1)/(4*Math.pow(R,E+2));
	    double t2=E*(E+2)/(16*Math.pow(R,E));
	    double t=t1-t2;
	    return(t);
	}
	return(0);
    }

    /**lambda2 from the paper **/

    public static double lambda2(double E,Box B1,Box B2,double R) {

	if(B1.type==1) {
	    double t=E/(8*Math.pow(R,E+1));
	    return(t);
	}	    

	if(B1.type==0) {   
                      double t0=E/(7.98*Math.pow(R,E+1));
	    double x1=B1.xAbsMax();
	    double y1=B1.yAbsMax();
	    double t1=Math.sqrt(1+x1*x1);
                      double t2=Math.sqrt(1+y1*y1);
	    double t=t0*(t1+t2);
	    return(t);
	}
	return(0);
    }



}