/*************************************************************************/
/* Program file name: Ex1_7.java                                         */
/*                                                                       */
/* © Tao Pang 2006                                                       */
/*                                                                       */
/* Last modified: February 8, 2012                                       */
/*                                                                       */
/* (1) This Java program is created for the book, "An Introduction to    */
/*     Computational Physics, 2nd Edition," written by Tao Pang and      */
/*     published by Cambridge University Press on January 19, 2006.      */
/*                                                                       */
/* (2) No warranties, express or implied, are made for this program.     */
/*                                                                       */
/*************************************************************************/

// Program for Exercise 1.7.  The distance between Halley's comet
// and the Sun is calculated with an algorithm via the 2nd-order
// Taylor expansions of the position and velocity over time.
// Using the lines commented out if only 1st-order terms are to be
// used (the Euler algorithm).

import java.lang.*;
public class Ex1_7 {
  static final int n = 200000, m = 2000;
  public static void main(String argv[]) {
    double t[] = new double [n];
    double x[] = new double [n];
    double y[] = new double [n];
    double r[] = new double [n];
    double vx[] = new double [n];
    double vy[] = new double [n];
    double gx[] = new double [n];
    double gy[] = new double [n];
    double h = 2.0/(n-1), h2 = h*h/2, k = 39.478428;

 // Initialization of the problem
    t[0] = 0;
    x[0] = 1.966843;
    y[0] = 0;
    r[0] = x[0];
    vx[0] = 0;
    vy[0] = 0.815795;
    gx[0] = -k/(r[0]*r[0]);
    gy[0] = 0;

 // Algorithm for the position and velocity with 2nd-order Taylor expansion
    for (int i=0; i<n-1; ++i) {
      t[i+1]  = h*(i+1);
//    x[i+1]  = x[i] + h*vx[i];
      x[i+1]  = x[i] + h*vx[i]
              + h2*gx[i];
//    y[i+1]  = y[i] + h*vy[i];
      y[i+1]  = y[i] + h*vy[i]
              + h2*gy[i];
      double r2 = x[i+1]*x[i+1]+y[i+1]*y[i+1];
      r[i+1] = Math.sqrt(r2);
      double r3 = r2*r[i+1];
      gx[i+1] = -k*x[i+1]/r3;
      gy[i+1] = -k*y[i+1]/r3;
//    vx[i+1] = vx[i] + h*gx[i];
      vx[i+1] = vx[i] + h*gx[i]
              + k*h2*(3*x[i]*(x[i]*vx[i]+y[i]*vy[i])/(r[i]*r[i])
              - vx[i])/(r[i]*r[i]*r[i]);
//    vy[i+1] = vy[i] + h*gy[i];
      vy[i+1] = vy[i] + h*gy[i]
              + k*h2*(3*y[i]*(y[i]*vy[i]+x[i]*vx[i])/(r[i]*r[i])
              - vy[i])/(r[i]*r[i]*r[i]);
    }
    for (int i=0; i<n-m; i+=m) {
      System.out.println(t[i] +" " + r[i]);
    }
  }
}