class Perihelion {
    // Mass, distance, and time are in centimeters
    static long r = 10000000;   // Geodesic step = 10000000 cm
    static long mass = 300000;  // Mass of star = 300000 cm (Sun)
    static long x0 = 100000000; // Planet initial position: x = 100000000 cm,
    static long y0 = 0;         // y = 0
    static double vx0 = 0;      // Initial velocity, vx = 0, vy = 0.02 c, 
    static double vy0 = 0.02;   // used to compute the second point

  static double d (double at, double ax, double ay, 
                   double bt, double bx, double by ) {
   double r,gtt,gxx,gxy,gyy,ut,ux,uy;
   r = Math.sqrt(ax*ax+ay*ay);
   gtt = 1 - mass/r; 
   gxx = - (ax * ax / (r * (r - mass)) + ay * ay / (r * r));
   gxy = - (mass * ax * ay / (r * r * (r - mass)));
   gyy = - (ax * ax / (r * r) + ay * ay / (r * (r - mass)));
   ut = bt-at;
   ux = bx-ax;
   uy = by-ay;
   return Math.sqrt(ut*gtt*ut+ux*gxx*ux+ux*gxy*uy + uy*gxy*ux + uy*gyy*uy);}

  static double l (long at, long ax, long ay, 
                   long bt, long bx, long by,
                   long ct, long cx, long cy) {
    return d(at,ax,ay,bt,bx,by) + d(bt,bx,by,ct,cx,cy);
  }

  static double w (long at, long ax, long ay, 
                   long bt, long bx, long by,
                   long ct, long cx, long cy) {
    double rt,rx,ry;
    rt = (l(at,ax,ay,bt+1,bx,by,ct,cx,cy) - l(at,ax,ay,bt-1,bx,by,ct,cx,cy))/2;
    rx = (l(at,ax,ay,bt,bx+1,by,ct,cx,cy) - l(at,ax,ay,bt,bx-1,by,ct,cx,cy))/2;
    ry = (l(at,ax,ay,bt,bx,by+1,ct,cx,cy) - l(at,ax,ay,bt,bx,by-1,ct,cx,cy))/2;
    return rt*rt + rx*rx + ry*ry;
  }

  static long[] bestneighbor (long at, long ax, long ay,
                              long bt, long bx, long by, 
                              long ct, long cx, long cy) {
    long [] t;
    long x, y;
    double best,min;
    t = new long [2];
    t[0] = cx;
    t[1] = cy;
    best = w(at,ax,ay,bt,bx,by,ct,cx,cy);
    for (x = cx-1; x <= cx+1; x++) {
      for (y = cy-1; y <= cy+1; y++) {
        min = w(at,ax,ay,bt,bx,by,ct,x,y);
        if (min < best) {
          best = min;
          t[0] = x; t[1] = y;}}}
    return t;
  }

  static long[] nextpoint (long at, long ax, long ay,
                           long bt, long bx, long by, 
                           long ct, long cx, long cy) {
    boolean found;
    long [] t;
    found = false;
    t = bestneighbor(at,ax,ay,bt,bx,by,ct,cx,cy);
    while (!found) {
      t = bestneighbor(at,ax,ay,bt,bx,by,ct,cx,cy);
      found = (t[0] == cx) && (t[1] == cy);
      cx = t[0]; 
      cy = t[1];}
    return t;
  }

  public static void main (String [] args) {
    long at, ax, ay, bt, bx, by, ct, cx, cy, dt, dx, dy;
    double alpha, alphaold, l, v, psi, aph, per, min, shift;
    long [] t;

    Isn.initDrawing("Gravitation",10,10,400,400);
    Isn.paintCircle(200,200,4,255,0,0);

    //initial data
    bt = 0;
    bx = x0;
    by = y0;

    ct = r;
    cx = x0 + Math.round(vx0 * r); 
    cy = y0 + Math.round(vy0 * r); 

    //fakepastpoint, used only for the computation of the tentative next point
    at = 2*bt - ct;
    ax = 2*bx - cx;
    ay = 2*by - cy;

    //things to be computed
    alphaold = 0.0;
    aph = 1.0;
    per = 1.0;

    while (true) {

      // compute next point
      dt = ct + r;
      dx = ax - 3 * bx + 3 * cx; 
      dy = ay - 3 * by + 3 * cy;
      t = nextpoint(bt,bx,by,ct,cx,cy,dt,dx,dy);
      dx = t[0];
      dy = t[1];

      //move forward
      at = bt;
      ax = bx;
      ay = by;

      bt = ct;
      bx = cx;
      by = cy;

      ct = dt;
      cx = dx;
      cy = dy;

      // draw the point
      Isn.drawPixel(200 + cx / 1e6,200 - cy / 1e6,0,0,0);

      //detect Aphaelion, Perihelion, etc.
      if ((bx*bx + by*by - (cx*cx+cy*cy) > 0) && 
          (bx*bx + by*by - (ax*ax + ay*ay) > 0)) {
        Isn.paintCircle(200+bx/1e6,200-by/1e6,2,0,0,255);
        System.out.println("Aphelion");
        System.out.println("t = " + bt + " cm" + " x = " + bx + " cm"
                           + " y = " + by + " cm");
        alpha = Math.atan2(by,bx) / 3.14159 * 180;
        l = Math.sqrt(bx*bx+by*by);
        v = Math.sqrt((bx-ax)*(bx-ax)+(by-ay)*(by-ay))/(bt-at);
        System.out.println("angle = " + alpha + " deg");
        System.out.println("distance = " + l + " cm");
        System.out.println("velocity = " + v + " c");
        aph = l;
        psi = 1.5 * 180 * mass * (1/per + 1/aph);
        System.out.println("theoretical shift = " + psi + " deg");
        shift = alpha - alphaold;
        if (shift < 0) shift = shift + 360;
        System.out.println("observed shift = " + shift + " deg");
        System.out.println("---------------------------------------------");
        alphaold = alpha;
      }

      if ((bx*bx+ by*by - (cx*cx+cy*cy) < 0) && 
          (bx*bx + by*by - (ax*ax + ay*ay) < 0))   {
        alpha = Math.atan2(by,bx) / 3.14159 * 180;
        l = Math.sqrt(bx*bx+by*by);
        v = Math.sqrt((bx-ax)*(bx-ax)+(by-ay)*(by-ay))/(bt-at);
        System.out.println("Perihelion");
        System.out.println("t = " + bt + " cm" + " x = " + bx
                           + " cm" + " y = " + by + " cm");
        System.out.println("angle = " + alpha + " deg");
        System.out.println("distance = " + l + " cm");
        System.out.println("velocity = " + v + " c");
        per = l;
        System.out.println("---------------------------------------------");
      }
    }
  }
}