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


/**This class does the 3D polyhedron exchange dynamics.*/

public class PolyhedronExchange {


    /**This is a simplified version of the main routine, which
       returns the super tile one gets by doing LIM steps of
       the iteration**/

    public static Polyhedron orbitTilePartial(Vector V,int LIM) {
	double[] d=orbitGeometryPartial(V,LIM);
	Polyhedron X=reconstruct(V,d);
	return(X);
    }

    public static Polyhedron orbitTilePartialInverse(Vector V,int LIM) {
	Vector W=new Vector(2-V.x[0],2-V.x[1],2-V.x[2]);
	double[] d=orbitGeometryPartial(W,LIM);
	Polyhedron X=reconstruct(W,d);
	for(int i=0;i<X.count;++i) {
	    X.V[i]=new Vector(2-X.V[i].x[0],2-X.V[i].x[1],2-X.V[i].x[2]);
	}
	return(X);
    }


    /**This is a simplified version of the orbit geometry routine.**/

    public static double[] orbitGeometryPartial(Vector V0,int LIM) {
	double[] d=new double[18];
	for(int i=0;i<18;++i) d[i]=100;
	Vector V=new Vector(V0);
	int type=-1;
	int count=0;
	int history=-1;
	while(count<LIM) {
          if(history==-1) type=PolyhedronExchange.classify(V);
	  if(history!=-1) type=PolyhedronExchange.classify(V,history);
	  history=type;
	  d=PolyhedronExchange.updateSpectrum(d,V,type);
          V=doDynamicsPlus(V,type);
	  ++count;
	}
	return(d);
    }






    /**This routine does the same as the previous one, except that
       it returns the itinerary sequence.**/

    public static int[] itinerary(Vector V0,int LIM) {
	int[] I=new int[LIM];
	Vector V=new Vector(V0);
	int type=-1;
	int count=0;
	int history=-1;
	boolean periodic=false;
	while((count<LIM)&&(periodic==false)) {
          if(history==-1) type=PolyhedronExchange.classify(V);
	  if(history!=-1) type=PolyhedronExchange.classify(V,history);
	  history=type;
	  I[count]=type;
          V=doDynamicsPlus(V,type);  
          if(Vector.dist(V,V0)<.0000001) periodic=true;
	  ++count;
	}
	int[] I2=new int[count];
	for(int i=0;i<count;++i) I2[i]=I[i];
	return(I2);
    }


    /**This routine gets the backwards itinerary.  
       The variable type records the regions in the
       backwards partition that contain the point.**/


    public static int[] itineraryInverse(Vector V0,int LIM) {
	int[] I=new int[LIM];
	Vector V=new Vector(V0);
	int type=-1;
	int count=0;
	boolean periodic=false;
	while((count<LIM)&&(periodic==false)) {
          type=classifyInverse(V);
	  I[count]=type;
          V=doDynamicsMinus(V,type); 
          if(Vector.dist(V,V0)<.0000001) periodic=true;
	  ++count;
	}
	int[] I2=new int[count];
	for(int i=0;i<count;++i) I2[i]=I[i];
	return(I2);
    }

    /**This is the main routine. It returns the orbit
       tile that contains the vector V.  It only looks
       at the first LIM forward iterates.  There are 3 modes:

       0:  straight-up tile
       1:  the core-tile of type A. 
       2:  the core-tile of type B, iso levels
       3:  the core-tile of type B, expansion levels

       The core-tile of type A consists of the tile of points that all
       return to A in the same way.  In case the points avoid A, the
       core-tile is just the orbit tile.  Type B has a similar explanation,
       except for one complication.  For the expansion levels we need to
       use the inverse dynamics for the calculation. We avoid doing this
       by precomputing the inverse first return map, then doing the
       forward dynamics, then suitably translating the result. **/

   public static Polyhedron orbitTile(int mode,double bot,double top,Vector V,int LIM) {
       if(mode<3) return(orbitTileBasic(mode,bot,top,V,LIM));  
        int[] listB=VerifySupport.itineraryAB(V,false,false,true);
	GoldenVector GW=VerifySupport.getTotalMoveInverse(listB);  //vector is scaled up by 2.
	Vector W=GW.toVector();
	W=W.scale(.5);  //scale back to correct size
	Vector RV=Vector.plus(V,W);
	RV=TorusMap.fundamentalDomain(RV);
	Polyhedron P=orbitTileBasic(2,bot,top,RV,LIM);
	Vector CM=P.getCenter();
	boolean high=false;
	if(CM.x[1]>1) high=true;
	for(int i=0;i<P.count;++i) {
	   P.V[i]=TorusMap.fundamentalDomain(Vector.minus(P.V[i],W));
	   if((high==true)&&(P.V[i].x[1]<.000001)) P.V[i].x[1]=P.V[i].x[1]+2;
	   if((high==false)&&(P.V[i].x[1]>1.999999)) P.V[i].x[1]=P.V[i].x[1]-2;
	}
	return(P);
   }

    /**Here is the basic routine. It works for the modes 0,1,2.  The
       mode-3 case requires the correction mentioned above.*/

    public static Polyhedron orbitTileBasic(int mode,double bot,double top,Vector V,int LIM) {
	double[] d=orbitGeometry(mode,V,LIM);
	if(d==null) return(null);
	d[0]=Math.min(d[0],top-V.x[2]);  //chop the top 
	d[9]=Math.min(d[9],V.x[2]-bot);  //chop the bottom
	Polyhedron X=reconstruct(V,d);
	X.orbit=(int)(d[18]);
	return(X);
    }



    /**This routine looks at the orbit and records how far each
       point in the orbit lies from each of the faces of the
       polyhedron containing it.  For each direction, the minimum
       over the orbit is computed.  These 18 numbers are then
       used to compute the orbit tile.  There are 3 options:
 
       mode = 0:  straight-up orbit
       mode = 1:  track the orbit until it lies in renorm set A
       mode = 2:  track the orbit until it lies in renorm set B
    **/

    public static double[] orbitGeometry(int mode,Vector V0,int LIM) {
	double[] d=new double[18];
	for(int i=0;i<18;++i) d[i]=100;
	Vector V=new Vector(V0);
	int history=-1;
	int type=-1;
	int count=0;
	int periodic=0;
	boolean test=false;
	boolean inside=false;
	while((count<LIM)&&(test==false)) {

	    /**the chopped polyhedra**/
	    if(mode==1) {
		int[] q=DataRenorm.indexA(V);
		if(q!=null) {
		    d=updateSpectrumA(d,V,q);
		    if(count>0) test=true;
		}
	    }

	    /**the chopped polyhedra**/
	    if(mode==2) {
	            int[] q=DataRenorm.indexB(V);
		    if(q!=null) {
		    d=updateSpectrumB(d,V,q);
                    if(count>0) test=true;
		}
	    }

	    if(test==false) {
	      if(history==-1) type=PolyhedronExchange.classify(V);
	      if(history!=-1) type=PolyhedronExchange.classify(V,history);
	      history=type;
	      d=PolyhedronExchange.updateSpectrum(d,V,type);
              V=doDynamicsPlus(V,type);
	      if(Vector.dist(V,V0)<.0000001) {
                 test=true;
	         periodic=1;
	      }
	      ++count;
	    }
	}  

	if(test==false) return(null);
	double[] e=new double[20];
	for(int i=0;i<18;++i) e[i]=d[i];
	e[18]=count; //the number of iterates
	return(e);
    }










    /**The vector is supposed to lie inside the qth
       polyhedron.  The list of 18 numbers is the
       dot product with the outer normals of the polyhedron. 
       Most entries will be 0.**/

    public static double[] updateSpectrum(double[] d1,Vector W,int q) {	
        GoldenPolyhedron P=DataPartition.getGoldenPolyhedron(q);
	int[] N=DataPartitionRaw.polyNormal(q);
	int[] A=DataPartitionRaw.faceAnchor(q);
	double[] d2=new double[18];
	for(int i=0;i<18;++i) d2[i]=d1[i];
	for(int i=0;i<A.length;++i) {
            GoldenVector gn=DataPartition.getNormal(N[i]);
	    Vector n=gn.toVector();
	    Vector p0=P.V[A[i]].toVector();
	    double test=Vector.dot(n,Vector.minus(p0,W));
	    if(d2[N[i]]>test) d2[N[i]]=test;
	}
	return(d2);
    }

    /**This does the same thing but incorporates the
       polyhedra from the A renorm set.*/

    public static double[] updateSpectrumA(double[] d1,Vector V,int[] q) {
	Vector W=DataRenorm.AtoB(q[0],q[1],q[2],V);
        GoldenPolyhedron P=DataRenorm.getGoldenB(q[1],q[2]);
	int[] N=DataRenormRaw.polyNormal(q[1],q[2]);
	int[] A=DataRenormRaw.faceAnchor(q[1],q[2]);
	double[] d2=new double[18];
	for(int i=0;i<18;++i) d2[i]=d1[i];
	for(int i=0;i<A.length;++i) {
            GoldenVector gn=DataPartition.getNormal(N[i]);
	    Vector n=gn.toVector();
	    Vector p0=P.V[A[i]].toVector();
	    double test=Vector.dot(n,Vector.minus(p0,W));
	    if(q[0]%2==0) test=test*GoldenRatio.phi(-3);
	    if(d2[N[i]]>test) d2[N[i]]=test;
	}
	return(d2);
    }




    /**This does the same thing but incorporates the
       polyhedra from the B renorm set.*/

    public static double[] updateSpectrumB(double[] d1,Vector W,int[] q) {	
        GoldenPolyhedron P=DataRenorm.getGoldenB(q[0],q[1]);
	int[] N=DataRenormRaw.polyNormal(q[0],q[1]);
	int[] A=DataRenormRaw.faceAnchor(q[0],q[1]);
	double[] d2=new double[18];
	for(int i=0;i<18;++i) d2[i]=d1[i];
	for(int i=0;i<A.length;++i) {
            GoldenVector gn=DataPartition.getNormal(N[i]);
	    Vector n=gn.toVector();
	    Vector p0=P.V[A[i]].toVector();
	    double test=Vector.dot(n,Vector.minus(p0,W));
	    if(d2[N[i]]>test) d2[N[i]]=test;
	}
	return(d2);
    }




    /**DYNAMICS ROUTINES**/

    public static Vector doDynamicsPlus(Vector V) {
	int type=classify(V);
	return(doDynamicsPlus(V,type));
    }

    public static Vector doDynamicsMinus(Vector V) {
	int type=classifyInverse(V);
	return(doDynamicsMinus(V,type));
    }

    public static Vector doDynamicsPlus(Vector V,int type) {
	Complex z=DataPartition.getMove3D(type);
	Vector W=Vector.plus(V,new Vector(z.x,z.y,0));
	W=TorusMap.fundamentalDomain(W);
	return(W);
    }

    public static Vector doDynamicsMinus(Vector V,int type) {
	Complex z=DataPartition.getMove3D(type);
	Vector W=Vector.minus(V,new Vector(z.x,z.y,0));
	W=TorusMap.fundamentalDomain(W);
	return(W);
    }
	/**CLASSIFICATION ROUTINES**/

    /**This tells which of the 64 regions contains
       the point W.*/

    public static int classify(Vector W) {
	for(int i=0;i<64;++i) {
	    if(testInside(W,i)==true) return(i);
	}
	return(-1);
    }

    /**This routine speeds up the preceding one.  If
       we know the history of the map, i.e. the previous
       iterate, then we can predict the possible
       indices which will contain the vector.**/

    public static int classify(Vector W,int history) {
	int[] A=DataPartitionRaw.successor(history);
	for(int i=0;i<A.length;++i) {
	    if(testInside(W,A[i])==true) return(A[i]);
	}
	return(-1);
    }


    public static int classifyInverse(Vector W) {
	for(int i=0;i<64;++i) {
	    if(testInsideInverse(W,i)==true) return(i);
	}
	return(-1);
    }

    public static boolean testInsideInverse(Vector W,int q) {
	Complex z=DataPartition.getMove3D(q);
	Vector X=Vector.minus(W,new Vector(z.x,z.y,0));
	X=TorusMap.fundamentalDomain(X);
	return(testInside(X,q));
    }

    /**tests if a vectors is inside the qth polyhedron*/

    public static boolean testInside(Vector W,int q) {
	GoldenPolyhedron P=DataPartition.getGoldenPolyhedron(q);
	int[] N=DataPartitionRaw.polyNormal(q);
	int[] A=DataPartitionRaw.faceAnchor(q);
	for(int i=0;i<A.length;++i) {
	    GoldenVector gn=DataPartition.getNormal(N[i]);
	    Vector n=gn.toVector();
	    Vector p0=P.V[A[i]].toVector();
	    double test=Vector.dot(n,Vector.minus(p0,W));
	    if(test>-.0000000001) test=0;
	    if(test<0) return(false);
	}
	return(true);
    }





    /**RECONSTRUCTION**/

    /**This routine reconstructs the polyhedron from the
       list of faces.**/


    public static Polyhedron reconstruct(Vector V,double[] d) {

	Vector[] LIST=new Vector[100];
	int total=0;
	for(int i=0;i<18;++i) {
	    for(int j=i+1;j<18;++j) {
		for(int k=j+1;k<18;++k) {
		    int[] n={i,j,k};
		    double[] dd={d[i],d[j],d[k]};
		    if((d[i]<100)&&(d[j]<100)&&(d[k]<100)) {
			Vector W=intersect(n,dd);
			boolean test=pass(d,W);
			if(test==true) {
			  W=Vector.plus(W,V);
			  if(Lists.match(W,LIST,total)==false) {
			    LIST[total]=new Vector(W);
			    ++total;
			  }
			}
		    }
		}
	    }
	}	
        Polyhedron X=new Polyhedron();
	X.count=total;
	for(int i=0;i<total;++i) X.V[i]=LIST[i];
	return(X);
    }

    /*Here n is a list of 3 indices, for outward normal vectors,
      and d is a list of 3 positive values.  Let Pj be the plane
      of vectors whose dot product with nj is dj.  Then the point returned
      is P0 intersect P1 intersect P2.  We first check to make sure that
      the normals are in general position.*/

    public static Vector intersect(int[] n,double[] d) {
	GoldenVector[] GN=new GoldenVector[3];
	Vector[] N=new Vector[3];
	for(int i=0;i<3;++i) {
            GN[i]=DataPartition.getNormal(n[i]);
	    N[i]=GN[i].toVector();
	}
	GoldenReal test=GoldenVector.tripleProduct(GN[0],GN[1],GN[2]);
	if(test.isZero()==true) return(null);
	Matrix M=new Matrix(N);
	M=M.inverse();
	Vector W=new Vector(d[0],d[1],d[2]);
	Vector X=Matrix.act(M,W);
	return(X);
    }

    /**Checks that the vector does not lie outside the
       polyhedron**/

    public static boolean pass(double[] d,Vector V) {
	if(V==null) return(false);
	for(int i=0;i<18;++i) {
	    Vector W=DataPartition.getNormalVector(i);
	    if(d[i]<100) {
		double d1=d[i];
		double d2=Vector.dot(V,W);
		double test=d1-d2;
		if(test<-.000000000001) return(false);
	    }
	}
	return(true);
    }

}