//    Kepler: Bisection method for Elliptic KEPler equation

//

//    Eq.: D - e_X sin D + e_Y cos D = L

//

//    Ref: T. Fukushima, 1996, AJ, 112, 2858-2861.

//    Output Variable: bekep = D (Solution)

//    Input Variables: ex = e_X, ey = e_Y, x0 = L



public class Kepler {



    // Local Arrays:

    static int NMAX = 512;

    static double a[]= new double[NMAX+1];

    static double s[]= new double[NMAX+1];

    static double c[]= new double[NMAX+1];

    // Constants

    static int NDIM0=7, NDIM1=8, NDIM2=9, MEND=30;

    static double PIHALF=Math.PI/2, TWOPI=2*Math.PI, EPSLON=1e-15d;

    static double A3=1d/3, A6=1d/6, A12=1d/12, A30=1d/30, A90=1d/90, A620=1d/620;



    // Preparation of Auxiliary Tables

    static {                     // use static initializer because they are contants

        int ixend = 1, i = 1;

        a[1] = Math.PI;

        s[1] = 0d;

        c[1] = -1d;

        double theta = Math.PI;

        for (int n=2; n<=NDIM2; n++) {

            ixend *= 2;

            theta *= 0.5d;

            for (int ix=1; ix<=ixend; ix++) {

                i++;

                a[i] = theta*(2*ix-1); 

                s[i] = Math.sin(a[i]);

                c[i] = Math.cos(a[i]);

            }

        }

    }



    // Binary search solution of Kepler's equation (elliptic case)

    // inputs: -1 < ex < +1,  -1 < ey < +1,  (ex*ex+ey*ey)<1, 0 < x0 < 2*pi

    // output: eccentric anomaly (which is easily converted to true anomaly or x,y)

    public static double bekep (double ex, double ey, double x0) {

        // Selection of Bisection Level (Revised Nov. 11, 1996)

        int ndim;

        double e2 = ex*ex + ey*ey;

        if(e2>1d)

           return 0d;     //Hyperbolic case

        else if(e2<0.5d)

           ndim = NDIM0;

        else if(e2<0.75d)

           ndim = NDIM1;

        else

           ndim = NDIM2;



        // Reduction of 'L'

        double x;

        if(x0<ey)

            x = ey - mod(ey-x0,TWOPI) + TWOPI;

        else if(x0 > ey+TWOPI)

            x = mod(x0-ey,TWOPI) + ey;

        else

            x = x0;



        double dx = x0-x;



        // Bisection Search

        int i=1;

        double f0;

        for (int n=2; n<=ndim; n++) {

              f0 = (a[i]-(ex*s[i]-ey*c[i]))-x;

              i *= 2;

              if(f0<=0d) i++;

        }

        

        // Newton Method by using Taylor Expanded Trig. Func.

        double f2 = ex*s[i]-ey*c[i];

        double f3 = ex*c[i]+ey*s[i];

        double f1 = 1d - f3;

        f0 = (a[i]-f2)-x;

        double d = -f0/f1;

        double d2,sd,dcd,df2,df3,fx,fy=0d;

        for(int m=1; m<=MEND; m++) {

            d2 = d*d*0.5d;

            sd = d*(1d-d2*(A3-d2*(A30-d2*A620)));

            dcd= d2*(1d-d2*(A6-d2*A90));

            df2= -f2*dcd + f3*sd;

            df3= -f3*dcd - f2*sd;

            fx = f1-df3;

            fy = (f0+d)-df2;

            d  = d-fy/fx;

            if(Math.abs(fy)<EPSLON) break;

        }

        double solution;

        if(Math.abs(fy)>EPSLON) {         // No convergence

            double z = a[i]+d;

            double fz = (z - (ex*Math.sin(z) - ey*Math.cos(z))) - x;

            System.out.println("(bekep) No Conv." + ex + ey + x + fz);

            solution = 0d;

        } else {

            solution = a[i] + d;          // Converged

            if(Math.abs(dx)>0d)

               solution += dx;    // Recovery of 'L' Reduction

        }

        return solution;

    }



    private static double mod(double x, double y) {

       if(Math.abs(x)<y)

         return Math.abs(x);

       else

         return x - y*Math.floor(x/y);

    }



    public static void main(String args[]) {

        if(args.length==3) {              //enter ex,ey and L on command line to obtain D

            double ex=(new Double(args[0])).doubleValue();

            double ey=(new Double(args[1])).doubleValue();

            double L =(new Double(args[2])).doubleValue();

            double D = bekep(ex,ey,L);

            double error = (D-L)-(ex*Math.sin(D)-ey*Math.cos(D));                   

            System.out.println(" " + ex + " " + ey + " " + L/TWOPI + " " + D/TWOPI + " " + error/TWOPI);

        }

        else

            test();

    }



    public static void test() {  // Test the bekep method over a range of parameters

        int iend=8, jend=4, kend=5;

        double dw = TWOPI/iend;

        double dL = TWOPI/jend;

        double de = 1d/kend;

        System.out.println("Test Result of 'bekep'");

        System.out.println("e_X    e_Y    L/(2Pi)    D/(2Pi)   Error/(2Pi)");

        double e,w,ex,ey,L,D,error;

        for(int k=1; k<=kend; k++) {

            e = k*de;

            for(int i=1; i<=iend; i++) {

                w=(i-1)*dw;

                ex = e*Math.cos(w);

                ey = e*Math.sin(w);

                for(int j=1; j<=jend; j++) {

                    L=(j-1)*dL;

                    D = bekep(ex,ey,L);

                    error = (D-L)-(ex*Math.sin(D)-ey*Math.cos(D));

                    System.out.println(" " + ex + " " + ey + " " + L/TWOPI + " " + D/TWOPI + " " + error/TWOPI);

                }

            }

        }

    }

}