From: kaufman@Xenon.Stanford.EDU (Marc T. Kaufman)
Subject: Re: NO NUMBERS! (or even SHORT INTS!)
Reply-To: kaufman@CS.Stanford.EDU
Organization: CS Department, Stanford University, California, USA
Date:  5 Oct 92 16:10:34 GMT

markh@csd4.csd.uwm.edu (Mark) writes:

>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.

You're right, somewhat.  However, we still have to get the 85 magic
observations into the program so that we can compute the Fourier series...

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

...

>#define GM (4.0E14)
>void Gravity(Vector r) { /* g = D - GM r/|r|^3 */
>   double r2, s, a; Vector n;
           ^^ clearely this should be r_squared (use of '2' is forbidden)
>   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);
                           ^^^ tsk, tsk.  What if you were computing this
in another universe with different laws of physics?

>   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);
         ^^^^^ time to take ANSI to task for using a digit
                         ^^^^^^^^ surely this should have been
#define radians_to_degrees 180.0/PI
>   tV = atan2(v.Y, v.X)*180.0/PI; mV = sqrt(v.X*v.X + v.Y*v.Y);
>   printf("%7.3f   ", t);
             ^^^ using string concatenation, we can defer the specification
of format widths by placing the 7.3 in a #define
>   printf("%10.3f %5.2f", mR - R, tR);
>   tV -= (tR - 90.0); Dir = 1;
                ^^^^         ^ You should use an enum (North, East, South West)
>   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) {
             ^^^     ^^^ Don't you remember, it was the size of ints that
started this thread...
>   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);
           ^^^ in the old days, we read these values in from a data file
rather than embedding them in the program.  That way you don't have to
recompile just to get a new value.
>}
-- 
Marc Kaufman (kaufman@CS.Stanford.EDU)