From: markh@csd4.csd.uwm.edu (Mark)
Newsgroups: comp.programming,comp.lang.c,comp.lang.pascal,comp.lang.objective-c
Subject: Re: NO NUMBERS! (or even SHORT INTS!)
Date: 5 Oct 1992 00:43:53 GMT
Organization: Computing Services Division, University of Wisconsin - Milwaukee
NNTP-Posting-Host: 129.89.7.4

In article <mcdonald.355@aries.scs.uiuc.edu> mcdonald@aries.scs.uiuc.edu (J. D. McDonald) writes:
>OK look at this, which, like the nutation of the moon, is a better example:
...

(Calculation intensive routine)
>#include <math.h>
>#include "astro.h"
>/* ....jupiter */
>p_jupiter (t, s, dml, ds, dm, da)
...

There's the problem.  You could end up using the same routine for every planet
and all the comets, if only the orbital parameters and masses had been
extracted out.

It's a mistake precomputing all the 85 magic numbers in the calculations of
this routine -- that's what the computer is for.  Likewise, it's a mistake to
use a source the precomputes those numbers -- those sources are meant for
hand computation.  If the author made an error -- highly probable, or you made
an error (or the scanner), the result is a hidden bug.

So you lose (1) reliability, (2) reuseability, and (3) readibility.

There's no reason not to step back to the perturbation theory itself and
write the program directly on the basis of that, integrating the necessary
Fourier Series calculations right into the software itself.  You don't need
to precompute Fourier coefficients any more.

The thing that authors prior to the computer age tend to do that's really
bothersome is they actualy expand out all the terms in full until they get
numbers by plugging in the orbital constants.  Don't "simplify" those equations
by plugging in the constants, keep them out and let the software do the
calculation directly from the equations instead.

The result will likely be even more clean and concise than the tables
themselves, and on top of that, easily verifiable.

One approach around the calculation bottleneck is to use complex
polynomials -- a language like C++ is ideally suited for this, since you can
reuse the operators +, -, *, /, etc in the software itself after defining your
own classes, such as for vectors, complex numbers or polynomials or otherwise.

Then you won't need to break down the vector equations into components
anymore.  You can directly program those nice looking pristine vector
equations into your software and never have to look at components again.

This is something I wrote that does short-term orbital extrapolation.
Apparently this is using 4th order Runge-Kutta extrapolation or maybe even a
predictor-corrector method (I don't have my book handy), and is printing out
the time, radius and velocity vector (in geocentric polar form) once each
second.  The D vector appears to be a perturbation simulating a rocket thrust,
which would mean this is giving you a second-by-second snapshot of a launch.

The Add() vector routine could have been written v += t*w (combining vector
addition and scalar multiplication), if this were in C++.  The New() routine
would have been written using C++'s memory allocation operator, new.

#include <stdio.h>
#include <math.h>
#include <float.h>
typedef struct { double X, Y; } Vector;
 
Vector g, v, r, D; double t;

/* v += t*w */ 
#define Add(v, t, w) ((v).X += (t)*(w).X, (v).Y += (t)*(w).Y)
#define New(v, x, y) ((v).X = (x), (v).Y = (y))

#define GM (4.0E14)
void Gravity(Vector r) { /* g = D - GM r/|r|^3 */
   double r2, s, a; Vector n;
   r2 = r.X*r.X + r.Y*r.Y; a = GM/r2;
   s = sqrt(r2); New(n, r.X/s, r.Y/s);
   g = D; Add(g, -a, n);
}
 
void Update(double dt) { /* t, r(t), v(t) -> t + dt, r(t + dt), v(t + dt) */
   Vector v0, vC, r0, rC;
   r0 = r; v0 = v;
   Gravity(r0); Add(v, dt/6.0, g);
   vC = v0; Add(vC, dt/2.0, g); Add(r, dt/3.0, vC);
   rC = r0; Add(rC, dt/2.0, vC); Gravity(rC); Add(v, dt/3.0, g);
   vC = v0; Add(vC, dt/2.0, g); Add(r, dt/3.0, vC);
   rC = r0; Add(rC, dt/2.0, v0); Gravity(rC); Add(v, dt/3.0, g);
   vC = v0; Add(vC, dt/2.0, g); Add(r, dt/3.0, vC);
   rC = r0; Add(rC, dt, vC); Gravity(rC); Add(v, dt/6.0, g);
   t += dt;
}
 
#define R 6378000.0
#define PI 3.1415926
void Stats(Vector r, Vector v, double t) {
   double tR, tV, mR, mV; char Dir;
   tR = atan2(r.Y, r.X)*180.0/PI; mR = sqrt(r.X*r.X + r.Y*r.Y);
   tV = atan2(v.Y, v.X)*180.0/PI; mV = sqrt(v.X*v.X + v.Y*v.Y);
   printf("%7.3f   ", t);
   printf("%10.3f %5.2f", mR - R, tR);
   tV -= (tR - 90.0); Dir = 1;
   if (tV < -180.0) tV += 360.0; else if (tV > 180.0) tV -= 360.0;
   if (tV < -90.0) Dir = 0, tV = -180.0 - tV;
   else if (tV > 90.0) Dir = 0, tV = 180.0 - tV;
   printf("   %10.3f %5.2f%c\n", mV, tV, Dir? 'F': 'B');
}

void Thrust(int T0, int T1, double X, double Y) {
   int I;
   New(D, X, Y);
   for (I = T0; I < T1; I++) Update(1.0), Stats(r, v, t);
}

void main(void) {
   int I;
   t = 0.0;
   New(r, 0.0, R), New(v, 500.0, 0.0);
   Stats(r, v, t);
   Thrust(  0,   150,  0.0, 20.0);
   Thrust(150,   335, 40.0,  0.0);
   Thrust(335, 10000,  0.0,  0.0);
}