Keplerian Orbital Coordinates

The qbasic example program for creating motion files uses the following subroutine to generate position data for an ellipsoidal orbit. Time t is expressed in revolutions. a and b are the semimajor and semiminor axes of an ellipse in the XZ plane. x and z receive the position coordinates. Based on a routine for solving Kepler's equation in Peter Duffett-Smith, Astronomy with Your Personal Computer, 2nd ed., Cambridge University Press, 1990, ISBN 0-521-38995-X.

SUB kepler (t AS SINGLE, a AS SINGLE, b AS SINGLE, x AS SINGLE, z AS SINGLE)

   CONST pi2 = 6.283185

   DIM ec AS SINGLE           'eccentricity
   DIM am AS SINGLE           'mean anomaly
   DIM ae AS SINGLE           'eccentric anomaly
   DIM at AS SINGLE           'true anomaly
   DIM m  AS SINGLE
   DIM d  AS SINGLE

   ec = SQR(a * a - b * b) / a

   am = pi2 * t
   m = am - pi2 * INT(am / pi2)
   ae = m

   DO
      d = ae - (ec * SIN(ae)) - m
      IF ABS(d) < .000001 THEN EXIT DO
      d = d / (1 - (ec * COS(ae)))
      ae = ae - d
   LOOP

   at = 2 * ATN(SQR((1 + ec) / (1 - ec)) * TAN(ae / 2))
   x = a * COS(at) - (ec * a)
   z = b * SIN(at)

END SUB