package ChaosDemos;
import java.awt.*;
import java.util.*;

//********************************************************************
/**                    
* Finds and characterizes periodic orbits <br>
* @version 21 April 1999
* @author Michael Cross
*/
//********************************************************************

public class periodicOrbits {

    /** Period of orbit */
    int period;
    
    /** Number of orbits */
    public int nOrbits;

    /** Increment in number of fit points */
    private int increment=100;
    
    /** Number of points to collect before fit */
    private int fitNumber;    

    /** Running store of nPeriod points */
    private double[] xSave;
    private double[] ySave;
    private double[] dataSave;

    /** index of last added point */
    private int index;

    /** iterates[period][number]
    *   to give linear fit near fixed points */
    private double[][] xIn;
    private double[][] yIn;
    private double[][] xOut;
    private double[][] yOut;
    private double[][] sig;
    private int[][] ix;

    /** Numberof iterates[period]
    *   stored for each periodic point */
    int[] nIterates;

    /** Current best estimate of periodic points */
    public double[][] fp;
    
    /** Eigenvalues */
    public double[] eig1;
    public double[] eig2;
    
    /** Eigenvectors and adjoints */
    public double[][] ev1;
    public double[][] ev2;
    public double[][] fv1;
    public double[][] fv2;
    
    /** true if found ith fixed point */
    public boolean[] foundFixedPoint; 
    
    /** true if new fixed point curve should be added to graph */
    public boolean[] addPoint;                                                 
    
    /** Tolerance for near periodic point */
    private double eps=.01;
    private double epsp=.01;

    /** true of found first point near periodic orbit */
    public boolean foundFirst=false;
    
    /** true if point near periodic orbit */
    public boolean foundOrbit=false;

    /** true if point part of working orbit */
    public boolean partOfOrbit=false;
    
    public periodicOrbits(int inPeriod) {

        period=inPeriod;
      
        // Create arrays
        xSave=new double[period];
        ySave=new double[period];
        dataSave=new double[2*period];

        fp=new double[period][2];
        eig1=new double[period];
        eig2=new double[period];
        ev1=new double[period][2];
        ev2=new double[period][2];
        fv1=new double[period][2];
        fv2=new double[period][2];

       
        foundFixedPoint=new boolean[period];
        addPoint=new boolean[period];
       
        xIn=new double[period][increment+1];
        yIn=new double[period][increment+1];
        xOut=new double[period][increment+1];
        yOut=new double[period][increment+1];
        sig=new double[period][increment+1];
        ix=new int[period][increment+1];
        
        index=0;
        
        nIterates=new int[period];
        for(int i=0;i<period;i++) {
            nIterates[i]=0;
            foundFixedPoint[i]=false;
            addPoint[i]=false;
        }    
        nOrbits=0;    
    }
    
//********************************************************************
/**
* Sets tolerances eps1, eps2
* @param eps1
* @param eps2
* For a period <i>p</i> orbit the <i>p</i>th iterate must be within a ditance eps1
* and intermediate points must be outside a distance eps2
*/
//********************************************************************    
    public void setTolerance(double eps1, double eps2) {
        eps=eps1;
        epsp=eps2;
    }    
    
//********************************************************************
/**
* Sets the number of points accumulated before a fit is done
* @param N: the number of points
*/
//********************************************************************    
    public void setFitNumber(int N) {
        fitNumber=N;
    }
    
//********************************************************************
/**
* Tests to see if point is part of periodic orbit and stores point for
* future test.
* @param point contains x and y coordinates as first two entries
*/
//********************************************************************    
    public int testPoint(double[] point) {
        int i;
        
        if(testPeriodicOrbit(point, xSave, ySave, index)) {
            foundOrbit=true;
            if(!foundFirst) {
                for(i=0;i<period;i++) {
                    fp[i][0]=xSave[modplus(index,i,period)];
                    fp[i][1]=ySave[modplus(index,i,period)];                    
                }
                nOrbits=1;
                foundFirst=true;
                partOfOrbit=true;
            }
            else {
                partOfOrbit=false;
                
                for(i=0;i<period;i++) {
                    if(dist(point,fp[i][0],fp[i][1]) < epsp) {
                        partOfOrbit=true;
                        if(nIterates[i]<increment) {
                            nIterates[i]++;                        
                            xIn[i][nIterates[i]]=xSave[index];
                            yIn[i][nIterates[i]]=ySave[index];
                            xOut[i][nIterates[i]]=point[0];
                            yOut[i][nIterates[i]]=point[1];
                            ix[i][nIterates[i]]=nIterates[i];
                            sig[i][nIterates[i]]=1.;
                        }
                        if(nIterates[i]>=fitNumber && !foundFixedPoint[i]) {
                            foundFixedPoint[i]=true;
                            fit(i,nIterates[i]);
                            addPoint[i]=true;    
                        }                            
                        break;
                    }                        
                }                    
                
                foundOrbit=true;    
            }
        }
        else foundOrbit=false;
        if(foundOrbit) {
            point[2]=xSave[index];
            point[3]=ySave[index];
        }
        xSave[index]=point[0];
        ySave[index]=point[1];
        dataSave[2*index]=point[0];
        dataSave[2*index+1]=point[1];
        index=modplus(index,1,period);
        if(foundOrbit)            
            if(partOfOrbit) return 1;
            else return 2;
        else return 0;
    }
    
//********************************************************************
/**
* Returns array containing periodic orbit in format for plotting
* x1,y1,x2,y2.....
*/
//********************************************************************    
    public double[] returnOrbit() {
        double[] returnData= new double[2*period];
        System.arraycopy(dataSave, 0, returnData, 0, 2*period);
        return returnData;
    }
    
//********************************************************************
/**
* Updates current best estimate of periodic orbit
* @param inFp array containing orbit in format x1,y1,x2,y2....
*/
//********************************************************************    
    public void updateFit(double[][] inFp) {
        for(int i=0;i<period;i++)
            for(int j=0;j<2;j++)
                fp[i][j]=inFp[i][j];
    }
    
//********************************************************************
/**
* Tests whether point is part of periodic orbit
* @param point x,y coordinates of point to be tested
* @param xSave array of x coordinates of periodic orbit
* @param ySave array of y coordinates of periodic orbit
*/
//********************************************************************    
    private boolean testPeriodicOrbit(double[] point, double[] xSave, double[] ySave, int index) {
        int j;
        if(dist(point,xSave[index],ySave[index]) > eps) return false;
        for(int i=1;i<period;i++) {
            j=modplus(index,i,period);
            if(dist(point,xSave[j],ySave[j]) < epsp) return false;
        }
        return true;
    }
    
//********************************************************************
/**
* Fits matrix to linearization about periodic point and calculates
* eigenvalues, eigenvectors and adjoints.
* Uses svdfitClass.
* @param ifit index of periodic orbit
* @param number of points to be fitted
*/
//********************************************************************    
    private void fit(int ifit, int n) {
        double a11,a12,a21,a22,xo,yo,det,rad,f1,f2,e1,e2;
        double[] a=new double[4];
        
        svdfitClass test=new svdfitClass(n,3);
        double chisq1=test.svdfit(ix[ifit],xOut[ifit],xIn[ifit],yIn[ifit],sig[ifit],a);        
        a11=a[1];
        a12=a[2];
        xo=a[3];
        
        double chisq2=test.svdfit(ix[ifit],yOut[ifit],xIn[ifit],yIn[ifit],sig[ifit],a);
        a21=a[1];
        a22=a[2];
        yo=a[3];
        det=(1.-a22)*(1.-a11)-a12*a21;
        fp[ifit][0]=((1.-a22)*xo+a12*yo)/det;
        fp[ifit][1]=(a21*xo+(1.-a11)*yo)/det;
        rad=Math.sqrt(Math.pow((a11+a22),2)-4.*(a11*a22-a21*a12));
        e1=0.5*((a11+a22)+rad);
        e2=0.5*((a11+a22)-rad);
        if(Math.abs(e1)>Math.abs(e2)) {
            eig1[ifit]=e1;
            eig2[ifit]=e2;
        }
        else {
            eig1[ifit]=e2;
            eig2[ifit]=e1;
        }    
            
        det=1./Math.sqrt(Math.pow(a11-eig1[ifit],2)+Math.pow(a12,2));
        ev1[ifit][0]=-a12*det;
        ev1[ifit][1]=(a11-eig1[ifit])*det;
        fv2[ifit][0]=ev1[ifit][1];
        fv2[ifit][1]=-ev1[ifit][0];
        
        det=1./Math.sqrt(Math.pow(a11-eig2[ifit],2)+Math.pow(a12,2));
        ev2[ifit][0]=-a12*det;
        ev2[ifit][1]=(a11-eig2[ifit])*det;
        fv1[ifit][0]=ev2[ifit][1];
        fv1[ifit][1]=-ev2[ifit][0];
        
    }
                                           
//********************************************************************
/**
* Calculates distance between two points
* @param p1 first point x,y
* @param x x-coordinate of second point
* @param y y-coordinate of second point
*/
//********************************************************************    
    private double dist(double[] p1, double x, double y) {
         return Math.abs(p1[0]-x)+Math.abs(p1[1]-y);
    }
    

//********************************************************************
/**
* Adds two numbers modulo third number
* @param i1 first number
* @param i2 second number
* @param i3 third number
*/
//********************************************************************    
    private int modplus(int i1, int i2, int i3) {
        int result;
        result=i1+i2;
        while(result >= i3)
            result=result-i3;
        return result;
    }                                           
}