/*************************************************************************/
/* Program file name: comet.c                                            */
/*                                                                       */
/* © Tao Pang 2006                                                       */
/*                                                                       */
/* Last modified: July 7, 2006                                           */
/*                                                                       */
/* (1) This C 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.     */
/*                                                                       */
/*************************************************************************/

#include <stdio.h>
#include <math.h>
#define n 20000
#define m 200

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

main() {
  int i;
  double h = 2.0/(n-1), h2 = h*h/2, k = 39.478428, r2, r3;
  double t[n], x[n], y[n],r[n], vx[n], vy[n], gx[n], gy[n];

// 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 (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];
    r2 = x[i+1]*x[i+1]+y[i+1]*y[i+1];
    r[i+1] = sqrt(r2);
    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 (i=0; i<n-m; i+=m)
    printf("%16.8lf %16.8lf\n", t[i], r[i]);
}