///////////////////////////////////////////////////////////////////////////
//                                                                       //
// Program file name: Comet.java                                         //
//                                                                       //
// © Tao Pang 2006                                                       //
//                                                                       //
// Last modified: January 18, 2006                                       //
//                                                                       //
// (1) This Java program is part of 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.     //
//                                                                       //
///////////////////////////////////////////////////////////////////////////

// An example to study the time-dependent position and
// velocity of Halley's comet via the Verlet algorithm.

import java.lang.*;
public class Comet {
  static final int n = 20000, m = 200;
  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;

 // Verlet algorithm for the position and velocity
    for (int i=0; i<n-1; ++i) {
      t[i+1]  = h*(i+1);
      x[i+1]  = x[i]+h*vx[i]+h2*gx[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+1]+gx[i])/2;
      vy[i+1]= vy[i]+h*(gy[i+1]+gy[i])/2;
    }
    for (int i=0; i<n-m; i+=m) {
      System.out.println(t[i]);
      System.out.println(r[i]);
      System.out.println();
    }
  }
}