//Jimmy Corno
//April 29, 2004
//Relaxing loop program for periodic orbits in constant flow.

class RosslerDescent{
    
    static int m = 100;      //Number of segments in the guess loop.
    static int n = 10;     //Number of points per segment
    static double[][] guessLoop = new double[m][3];
    static double[][] delta = new double[m][3];
    static double[][] L = new double[m][3];
    static double[][] Lhat = new double[m][3];
    static double[] magL = new double[m];
    static double[][][] pos = new double[m][n+1][3];
    static double[][][] vel = new double[m][n+1][3];
    static double[][] COS = new double[m][n+1];
    static double[][] magV = new double[m][n+1];
    static double[] F2 = new double[m];
    static double F2new, F2old;
    static double c = 2.5;
    static double[] Gamma = new double[m];
    static int q = 0;
    static double bunchingTolerance = 3.0;
    
    public static void main( String[] Args ) {
        
        //gamma is the modifier for the descent; smaller gamma decreases step
        //size
        double gamma = 0.5;
        
        //q is the number of steps
        
        
        //this for loop lays out the guess loop
        for(int i=0; i<m-1; i++) {
            guessLoop[i][0] = 6.2*Math.sin(1.1*i*5.6/m - 1.1) + 0.1;
            guessLoop[i][1] = 5.3*Math.sin(1.1*i*5.6/m + 3.4) -0.7;
            guessLoop[i][2] = 0.1 + 2.1*Math.exp(-(i*5.6/m-3.4)*(i*5.6/m-3.4)/0.405);
        }
        
        guessLoop[m-1][0] = (5.2*Math.sin(1.1*(m-1)*5.6/m - 1.1) + 0.1 + 5.2*Math.sin( - 1.1) + 0.1 )/2.0;
        guessLoop[m-1][1] = (4.3*Math.sin(1.1*(m-1)*5.6/m + 3.4) -0.7 + 4.3*Math.sin( 3.4) -0.7)/2.0;
        guessLoop[m-1][2] = (0.1 + 3.1*Math.exp(-((m-1)*5.6/m-3.4)*((m-1)*5.6/m-3.4)/0.405) + 0.1 + 3.1*Math.exp(-(-3.4)*(-3.4)/0.405) )/2.0;
        
        //these function calls make the necessary calculations for the guess 
        //loop
        loopSegments();
        loopPoints();
        
        //Calculating initial cost function.
        F2new=0.0;
        for(int i=0; i<m; i++ ){
            F2[i] = costFunction(i);
            F2new += F2[i];
            if(F2[i] < 0.0){
                unkink(i);
                F2new = 0.0;
                i = 0;
            }
        }
        
        for(int i=0; i<m; i++ ){
            Gamma[i] = F2[i]/F2new*gamma;
        }
        //System.out.println(F2new);
        
        //do loop relaxes guess loop to reduce cost function.
        do{
            F2old = F2new;
            
            for( int i=0; i<m; i++) {
                delta[i][0] = deltax(i);
                delta[i][1] = deltay(i);
                delta[i][2] = deltaz(i);
            }
            
            //averages loop deformations for each point
            for( int i=0; i<m; i++) {
                    guessLoop[(i+1)%m][0] -= Gamma[i]*(delta[(i+1)%m][0] 
                                            + 1.0*delta[i][0])/2.0;
                    guessLoop[(i+1)%m][1] -= Gamma[i]*(delta[(i+1)%m][1] 
                                            + 1.0*delta[i][1])/2.0;
                    guessLoop[(i+1)%m][2] -= Gamma[i]*(delta[(i+1)%m][2] 
                                            + 1.0*delta[i][2])/2.0;
            }
            
            loopSegments();
            loopPoints();
            
            F2new=0.0;
            for(int i=0; i<m; i++ ){
                F2[i] = costFunction(i);
                F2new += F2[i];
                if(F2[i] < 0.0){
                    unkink(i);
                    F2new = 0.0;
                    i = 0;
                }
            }
            
            for(int i=0; i<m; i++ ){
                Gamma[i] = F2[i]/F2new*gamma;
            }
            
            //System.out.println( F2new );
            q++;
            
        }while((F2new > 0.05) && (q<100000));
        
        //distance from the origin to a point on the loop; should be 1 for 
        //correct loop
        
        /*
        double r = Math.sqrt(guessLoop[2][0]*guessLoop[2][0] + 
                    guessLoop[2][1]*guessLoop[2][1]);
        
        System.out.println(q + "," + r);
        */
        
        //prints guess loop points to screen
        for( int i=0; i<m; i++ ) {
            System.out.print( guessLoop[i][0] + "," + guessLoop[i][1] + "," + guessLoop[i][2] );
            System.out.println( "," + F2[i]);
        }
        System.out.println(q + "\t" + F2new + "\t" + F2old );
    }
    
    //vx and vy and vz are methods to calculate the flow velocity
    public static double vx( double x, double y, double z ){
        double VX = -y - z;
        return(VX);
    }
    
    public static double vy( double x, double y, double z ){
        double VY = x + 0.2*y;
        return(VY);
    }
    
        public static double vz( double x, double y, double z ){
        double VZ = 0.2 + z*(x - c);
        return(VZ);
    }
    
    //loopSegments picks out the loop segments from the guess loop
    public static void loopSegments() {
        for(int i=0; i<m; i++) {

            L[i][0] = guessLoop[(i+1)%m][0] - guessLoop[i][0];
            L[i][1] = guessLoop[(i+1)%m][1] - guessLoop[i][1];
            L[i][2] = guessLoop[(i+1)%m][2] - guessLoop[i][2];
            magL[i] = Math.sqrt( L[i][0]*L[i][0] + L[i][1]*L[i][1] 
                      + L[i][2]*L[i][2] );
            Lhat[i][0] = L[i][0] / magL[i];
            Lhat[i][1] = L[i][1] / magL[i];
            Lhat[i][2] = L[i][2] / magL[i];
            
            if(magL[i]/magL[(i+m-1)%m] > bunchingTolerance ){
                unkink(i);
            }
            
        }
    }
    
    //same as loopSegments, but this one is for specific segments.
        public static void loopSegments( int i) {
            int k = (i+m-1)%m;

            L[i][0] = guessLoop[(i+1)%m][0] - guessLoop[i][0];
            L[i][1] = guessLoop[(i+1)%m][1] - guessLoop[i][1];
            L[i][2] = guessLoop[(i+1)%m][2] - guessLoop[i][2];
            magL[i] = Math.sqrt( L[i][0]*L[i][0] + L[i][1]*L[i][1] 
                      + L[i][2]*L[i][2] );
            Lhat[i][0] = L[i][0] / magL[i];
            Lhat[i][1] = L[i][1] / magL[i];
            Lhat[i][2] = L[i][2] / magL[i];
            
            L[k][0] = guessLoop[(k+1)%m][0] - guessLoop[k][0];
            L[k][1] = guessLoop[(k+1)%m][1] - guessLoop[k][1];
            L[k][2] = guessLoop[(k+1)%m][2] - guessLoop[k][2];
            magL[k] = Math.sqrt( L[k][0]*L[k][0] + L[k][1]*L[k][1] 
                      + L[k][2]*L[k][2] );
            Lhat[k][0] = L[k][0] / magL[k];
            Lhat[k][1] = L[k][1] / magL[k];
            Lhat[k][2] = L[k][2] / magL[k];
    }
    
    //loopPoints determines the position, flow velocity, and cosine between
    //flow velocity and loop velocity for the integration points.
    public static void loopPoints() {
        for( int i=0; i<m; i++ ) {
            for( int j=0; j<=n; j++ ) {
                pos[i][j][0] = guessLoop[i][0] + j*L[i][0] / n;
                pos[i][j][1] = guessLoop[i][1] + j*L[i][1] / n;
                pos[i][j][2] = guessLoop[i][2] + j*L[i][2] / n;
                vel[i][j][0] = vx( pos[i][j][0], pos[i][j][1], pos[i][j][2] );
                vel[i][j][1] = vy( pos[i][j][0], pos[i][j][1], pos[i][j][2] );
                vel[i][j][2] = vz( pos[i][j][0], pos[i][j][1], pos[i][j][2] );
                magV[i][j] = Math.sqrt( vel[i][j][0]*vel[i][j][0] + 
                                        vel[i][j][1]*vel[i][j][1] + 
                                        vel[i][j][2]*vel[i][j][2] );
                COS[i][j] = (L[i][0]*vel[i][j][0] + L[i][1]*vel[i][j][1]
                            + L[i][2]*vel[i][j][2])/(magL[i] * magV[i][j]);
                
/*                if(COS[i][j] <= 0.0 ) {
                    System.out.println( COS[i][j] + "," + pos[i][j][0] + ","
                                            + pos[i][j][1]);
                }
 */
            }
        }

    }
    
    //same as loopPoints, but this one is for an individual segment instead of 
    //the whole loop
    public static void loopPoints( int i) {
        int k = (i+m-1)%m;
        
        for( int j=0; j<=n; j++ ) {
            pos[k][j][0] = guessLoop[k][0] + j*L[k][0] / n;
            pos[k][j][1] = guessLoop[k][1] + j*L[k][1] / n;
            pos[k][j][2] = guessLoop[k][2] + j*L[k][2] / n;
            vel[k][j][0] = vx( pos[k][j][0], pos[k][j][1], pos[k][j][2] );
            vel[k][j][1] = vy( pos[k][j][0], pos[k][j][1], pos[k][j][2] );
            vel[k][j][2] = vz( pos[k][j][0], pos[k][j][1], pos[k][j][2] );
            magV[k][j] = Math.sqrt( vel[k][j][0]*vel[k][j][0] + 
                                        vel[k][j][1]*vel[k][j][1] + 
                                        vel[k][j][2]*vel[k][j][2] );
            COS[k][j] = (L[k][0]*vel[k][j][0] + L[k][1]*vel[k][j][1]
                            + L[k][2]*vel[k][j][2])/(magL[k] * magV[k][j]);
                
            pos[i][j][0] = guessLoop[i][0] + j*L[i][0] / n;
            pos[i][j][1] = guessLoop[i][1] + j*L[i][1] / n;
            pos[i][j][2] = guessLoop[i][2] + j*L[i][2] / n;
            vel[i][j][0] = vx( pos[i][j][0], pos[i][j][1], pos[i][j][2] );
            vel[i][j][1] = vy( pos[i][j][0], pos[i][j][1], pos[i][j][2] );
            vel[i][j][2] = vz( pos[i][j][0], pos[i][j][1], pos[i][j][2] );
            magV[i][j] = Math.sqrt( vel[i][j][0]*vel[i][j][0] + 
                                        vel[i][j][1]*vel[i][j][1] + 
                                        vel[i][j][2]*vel[i][j][2] );
            COS[i][j] = (L[i][0]*vel[i][j][0] + L[i][1]*vel[i][j][1]
                            + L[i][2]*vel[i][j][2])/(magL[i] * magV[i][j]);
            }

    }
    
    //costFunction calls vps(V Perp Squared) to perform the cost function 
    //integration
    public static double costFunction( int i ) {
        double F = c0*vps(i,0) + c1*vps(i,1) + c2*vps(i,2) + 
                    c3*vps(i,3) + c3*vps(i,n-3) + c2*vps(i,n-2) + 
                    c1*vps(i,n-1) + c0*vps(i,n);
        for( int j=4; j<n-3; j++ ){
            F += vps(i,j);
        }
        return F*magL[i]/n;
    }
    
    //(V Perp Squared)/magnitude of loop velocity (so that integration is over 
    //loop segment length
    public static double vps( int i, int j ){
        double f = magV[i][j] * ( 1.0/COS[i][j] - COS[i][j] );

        return f;
    }
    
    //deltax and deltay call vpx and vpy (x and y components of v perp, 
    //respectively) to integrate the descent values
    public static double deltax( int i ) {
        double F = c0*vpx(i,0) + c1*vpx(i,1) + c2*vpx(i,2) + 
                    c3*vpx(i,3) + c3*vpx(i,n-3) + c2*vpx(i,n-2) + 
                    c1*vpx(i,n-1) + c0*vpx(i,n);
        for( int j=4; j<n-3; j++ ){
            F += vpx(i,j);
        }
        return F*magL[i]/n;
    }
    
    public static double deltay( int i ) {
        double F = c0*vpy(i,0) + c1*vpy(i,1) + c2*vpy(i,2) + 
                    c3*vpy(i,3) + c3*vpy(i,n-3) + c2*vpy(i,n-2) + 
                    c1*vpy(i,n-1) + c0*vpy(i,n);
        for( int j=4; j<n-3; j++ ){
            F += vpy(i,j);
        }
        return F*magL[i]/n;
    }
    
        public static double deltaz( int i ) {
        double F = c0*vpz(i,0) + c1*vpz(i,1) + c2*vpz(i,2) + 
                    c3*vpz(i,3) + c3*vpz(i,n-3) + c2*vpz(i,n-2) + 
                    c1*vpz(i,n-1) + c0*vpz(i,n);
        for( int j=4; j<n-3; j++ ){
            F += vpz(i,j);
        }
        return F*magL[i]/n;
    }
    
    //vpx is the x component of v perp/magnitude of loop velocity
    public static double vpx( int i, int j ){
        double f = Lhat[i][0] - vel[i][j][0]/magV[i][j]*COS[i][j];
        return f;
    }
    
    public static double vpy( int i, int j ){
        double f = Lhat[i][1] - vel[i][j][1]/magV[i][j]*COS[i][j];
        return f;
    }
    
    public static double vpz( int i, int j ){
        double f = Lhat[i][2] - vel[i][j][2]/magV[i][j]*COS[i][j];
        return f;
    }
    
    //unkink removes kinks or spaces out bunched up loop points
    public static void unkink( int i ){
        int k;
        if(i==0){
            k=m-1;
        }
        else{
            k=(i-1)%m;
        }
        for( int j=0; j<3; j++) {
            guessLoop[i][j] = (guessLoop[(i+1)%m][j] + guessLoop[k][j])/2.0;
        }
        loopSegments(i);
        loopPoints(i);
        //System.out.println("Unkinked." + q);
    }
    

    
    //These are integration constants for my extended Simpson's rule integration
    private static double c0 = 17.0/48.0;
    private static double c1 = 59.0/48.0;
    private static double c2 = 43.0/48.0;
    private static double c3 = 49.0/48.0;
    
}