/*************************    Program 7.1    ****************************/
/*                                                                      */
/************************************************************************/
/* Please Note:                                                         */
/*                                                                      */
/* (1) This computer program is written by Tao Pang in conjunction with */
/*     his book, "An Introduction to Computational Physics," published  */
/*     by Cambridge University Press in 1997.                           */
/*                                                                      */
/* (2) No warranties, express or implied, are made for this program.    */
/*                                                                      */
/************************************************************************/

#include <stdio.h>
#include <math.h>
#define N0 20001
#define NP 10000
#define N1   200

main()
/* Program for the time-dependent position of Halley's comet
   with the velocity version of the Verlet algorithm.
   Copyright (c) Tao Pang 1997. */
{
int i,n,m;
double h,kappa,h2,r2,r3;
double x[N0],y[N0],vx[N0],vy[N0],gx[N0],gy[N0],t[N0],r[N0];

/* Initialization of the problem */

n = N0;
m = N1;
h = 1.0/NP;
h2 = h*h/2;
kappa = 39.478428;
t[0] = 0;
x[0] = 1.966843;
y[0] = 0;
r[0] = x[0];
vx[0] = 0;
vy[0] = 0.815795;
gx[0] = -kappa*x[0]/pow(r[0],3);
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] = -kappa*x[i+1]/r3;
  gy[i+1] = -kappa*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; i = i+m)
  {
  printf("%16.8lf %16.8lf\n", t[i],r[i]);
  }
}