text: 
  /* Calculate the velocity vector of the earth
   * as it moves in its elliptical orbit.
   * The difference in velocity between this and assuming
   * a circular orbit is only about 1 part in 10**4.
   *
   * Note that this gives heliocentric, not barycentric, velocity.
   *
   * Input is Julian date.  Output left in global array vearth[].
   */
  
  double jvearth = -1.0;
  double vearth[3] = {0.0};
  
  #include "kep.h"
  
  extern struct orbit earth;
  
  int velearth( J )
  double J;
  {
  double e[3], p[3], t;
  int i;
  #if DEBUG
  double x[3], A, q;
  #endif
  
  if( J == jvearth )
  	return(0);
  
  jvearth = J;
  
  /* calculate heliocentric position of the earth
   * as of a short time ago.
   */
  t = 0.005;
  kepler( TDT-t, &earth, e, p );
  
  for( i=0; i<3; i++ )
  	vearth[i] = (rearth[i] - e[i])/t;
  
  #if DEBUG
  /* Generate display for comparison with Almanac values. */
  for( i=0; i<3; i++ )
  	{
  	q = vearth[i];
  	A += q*q;
  	x[i] = q;
  	}
  A = sqrt(A);
  precess( x, TDT, 1 );
  printf( "Vearth %.6e, X %.6f, Y %.6f, Z %.6f\n", A, x[0], x[1], x[2] );
  #endif
  return(0);
  }
  