Volume Number: 14 (1998)
Issue Number: 5
Column Tag: Beginner's Corner

Scientific Programming in PowerPlant

by Michael McLaughlin, McLean, VA

PowerPlant provides a simple application shell with which we can do serious number crunching, and still have a user-friendly, Macintosh interface

"I can do everything I need in C."

Does that sound familiar? I hear it a lot, not just with reference to C but also to Fortran and, once, to Basic! Scientists and engineers, in particular, often tend to view a visit to the compiler much as they would a trip to the dentist -- a rueful necessity, best done as quickly and painlessly as possible. The result, unfortunately, is seldom a nice, friendly application with a pleasing graphical user interface (GUI) but, rather, HelloWorld plus footnotes, usually accompanied by lots of global variables, predefined arrays, and nary a pointer in sight.

The standard excuse is that they are in a hurry and, anyhow, C++ is useful only for software that manipulates objects (whatever they are) and is totally inappropriate for the sort of mathematical calisthenics that cause physicists, astronomers, and others to chew up hours of clock cycles in search of a single number.

There are several reasons why computations can last for hours. They might involve a large number of mini-problems, or they might address a single problem comprised of many steps. In some cases, the solution could require iteration to convergence, and this might take a while. The example described below contains all of these elements. However, its size is small enough for demonstration without compromising the intended objectives.

This article has three objectives: to demonstrate that C++ not only is appropriate but, indeed, advantageous when it comes to writing mathematical software; to show that PowerPlant is an easy way to develop Macintosh applications rapidly; finally, to present a simple template, generally suitable for adaptation to CPU-intensive applications.

We begin with the Solar System at midnight (GMT) on the Fourth of July, 1976. We are given the Cartesian coordinates for all of the planets and the Sun, along with some additional data to be discussed below. The problem is to use Newton's Law of Gravitation and his Second Law of motion to predict the configuration of the Solar System 365 days later, keeping the system center-of-mass as a constant origin.

Given Newton's Laws, this task appears straightforward. It consists of a set of 27 second-order differential equations (three coordinates per planet), each of which has the following form:

where q(n) is a coordinate of body[n], t is time, m(j) is the mass of body[j], r(ij) is the radial distance between body[i] and body[j], and G is the gravitational constant. Time is measured in days and distance in astronomical units (AU). (1 AU is, approximately, the average distance between the Earth and the Sun.) This formulation ignores several important factors, namely, Special and General Relativity, nonsphericity of the planets, and a handful of asteroids that perturb the orbit of Mars significantly, but these all are secondary effects, complicated though they are.

Solving this problem requires integrating all 27 differential equations over successive time intervals, h, until the target time is reached. Here, we shall use Numerov's method, in

where R is the error for this coordinate for interval k. Since F(q)(k+1) appears on the right-hand side of this equation, even though q(k+1) is obviously unknown, iteration to convergence will be needed.

With every iteration, the center-of-mass constraint determines the coordinates of the Sun. For each coordinate, we compute

Numerov's method is not meant to achieve world-class accuracy. In fact, its most redeeming quality is its brevity. This technique requires two very accurate starting configurations as well as values for the masses of the 10 bodies (as is customary, Earth and Moon are combined). These values are listed in the data file (init.dat) accompanying the source code package available from ftp://ftp.mactech.com/src/mactech/volume14_1998/14.05.sit.

And now, it's time to start coding.

What follows assumes familiarity with CodeWarrior and Constructor (or, at least, ResEdit). Knowledge of C++ is not assumed although a passing acquaintance with its syntax will be helpful. This example was compiled using CW Pro 2. Earlier versions should require few, if any, changes.

Start by creating a new project, named SolarSys.µ, from the Basic PowerPlant PPC (or 68K) stationery. This ensures that it gets the requisite PP files and libraries. It also means that almost all of the project headers will be precompiled. Change the names of the project and target to SolarSys. Select Run to compile the PP files and to ascertain that everything works.

The Main Routine

This default app does not do exactly what we want, so we shall have to modify it a little bit. Open the file CPPStarterApp.h and replace CPPStarterApp with SolarSysApp everywhere. Then, after the #include, add the line

#include "Flags.h"

Just before the protected block, add the line

Uint32   Flags;

Save As SolarSysApp.h and you're done. Now, for the source file. In CPPStarterApp.cp, replace CPPStarterApp with SolarSysApp everywhere, then add

#include "CSetup.h"

Delete the sample ResIDT.

We are going to make our own custom dialog box, which must be "registered." Therefore, change the constructor to read

RegisterAllPPClasses();
// Register my ppob classes
URegistrar::RegisterClass(CSetup::class_ID,
   (ClassCreatorFunc) CSetup::CreateSetupStream);

Modify cmd_New, in ObeyCommand (shown below), to get our own dialog box (that we haven't created yet ) and to initialize the Flags. Don't worry that Setup is only a local variable; PowerPlant maintains internal links between the dialog box and its SuperCommander, the application.

case cmd_New:
   LDialogBox   *Setup;
   Setup = (LDialogBox*)
      LWindow::CreateWindow(rRidL_Setup, this);
   Flags = WORKING;
   break;

Finally, modify cmd_New in FindCommandStatus as follows:

case cmd_New:
   outEnabled = (Flags & WORKING) ? false : true;
   break;

Believe it or not, that's all there is. Just Save As SolarSysApp.cp. Don't try to compile yet; some pieces still are missing.

Custom Stuff

Classical procedural languages, like C, are built around verbs. "Do this. Then, do that. If this happens, go there." and so on. Object-oriented languages, such as C++, are built around nouns. They talk about "things" and the functionality of those things. When written properly, a C++ program solves a problem using the same words that are found in the statement of the problem. In this case, we are going to invent four new nouns (classes) that are especially suitable to our specific problem. They are CVector, CMatrix, CSetup, and CSolver, the "C" indicating a class that we devise, in contrast to those, such as LWindow and URegistrar, which come with system libraries and utilities, respectively. Any specific instance of a class that we might declare is called an "object." Were this article a lot longer, we could elaborate our design by defining CPlanet, CForce, etc.

This discussion is limited to the main features of these four classes. All of the code is included in the code listings, except for SolarSysApp.cp, which is nearly unchanged from the stationery. As you examine the various functions, you may think of enhancements and/or ways to make them better. A lot of obvious features have been deliberately left out to encourage you to experiment.

CVector and CMatrix

By far, the most useful classes in this project are CVector and CMatrix. The words, "Vector" and "Matrix" do not appear in the statement of the problem but they are so ubiquitous in mathematical computing that they deserve to be separate classes. If you write a lot of such programs in C and this is your initial entry into the world of C++, then it's likely that these two classes alone are worth the price of admission. You may never have to calloc() again!

CVector and CMatrix constitute the module SimpleArrays. They are just that -- a one-frill implementation of one and two-dimensional arrays, the one frill being runtime array-bounds checking. Much more elaborate implementations can be found at several websites. A good starting point is the Object-Oriented Numerics Page at http://monet.uwaterloo.ca/oon/.

These two classes are template classes so they can be used to declare arrays of anything, even themselves! Using them is trivial. For instance,

CVector<double>   factor(Sun, Pluto);

declares a vector of doubles with indices from Sun to Pluto (enumerated elsewhere) while

CMatrix<double>   initial_state;

declares a matrix of doubles with rows and columns having the default indices 1 to 10. Of course, all arrays are resizable, such as

initial_state.resize(Sun, Pluto, X, Z);

Arrays of higher dimension are also easy to construct. Here are two different ways to declare a three-dimensional matrix, the second making copies of two existing matrices.

CVector< CMatrix<double> > state(-1, 1);
CMatrix<double>   trial_state[2] = 
      {state[current], state[next]};

There are several examples of these classes scattered throughout this project and their utility should be apparent.

One caution: SimpleArrays.cp must be #included in one of the other source files, not added to the project. This is a consequence of the template code. More often, template code is simply merged into one huge, typically unreadable, header file instead of being placed in its own source file.

The remaining source files, for CSolver, CSolver_I_O, and CSetup must all be added to the project.

CSolver

Even if you go to the trouble of making your own event loop, chances are that it will not suffice for a problem like this one. In CPU-intensive tasks, typically your code launches into a for or while loop containing the essential calculations and your computer is then frozen out until that loop finishes. This could take hours or, even, days.

To alleviate this problem and create an application that behaves reasonably, we introduce the class CSolver. This class inherits functionality from both LPeriodical and LBroadcaster.

Here is the beginning of the CSolver definition:

// CSolver is an "idler" LPeriodical which broadcasts status messages.
class   CSolver : public LPeriodical, public LBroadcaster {
public:
// These members are generic to all solvers:
                      CSolver(LListener *setup);
   virtual          ~CSolver() {};
            void      SpendTime(const EventRecord &inEvt);
            void      doInput();
            void      doOutput(int mode = 0);
private:
            void      Initialize();   // resize arrays, etc.
            void      oneCycle();     // called by SpendTime()
            Uint32   *flags;

LPeriodical functionality allows the MainEventLoop() built into LApplication to perform as you would expect. Among other things, this means letting you Quit whenever you want to, not just when your computation is finished. The trick is to make Solver (the CSolver object) an "idler" periodical by calling its method StartIdling(). Then, its virtual function SpendTime() will not be called until the application's event queue is empty. (Here, there is also a READY flag for extra protection.) SpendTime(), in turn, will trigger one cycle of the computation. Therefore, all we have to do is to break up our task into little pieces, short enough to be completed before returning to MainEventLoop(). This is usually easy to do since, as noted above, long computations are almost always formed from a large number of smaller computations.

Solver terminates each cycle by broadcasting one of three status messages. These indicate that a phase of the computation has finished, or that the integrator diverged, or that there is more work to be done.

The I/O functions are generic components but their details are task-specific. These two methods could be a lot fancier, for instance, with standard file dialogs instead of assuming a default folder. Even so, note how much simpler and cleaner they look compared to their usual C counterparts. This code includes the large, non-PP header <fstream>. In general, there could be several I/O functions requiring this header. Combining them in a separate module is one way to keep the size of a project to a minimum.

The rest of CSolver is also specific to our task. Most of this will become clearer once all of the project's components are in place.

CSetup

We are now ready for the final piece. This is the Setup/Status dialog, a customized dialog box which will function as our GUI and as the "Commander/Listener" for Solver. In a sense, Solver belongs to Setup just as Setup belongs to SolarSysApp. It is, therefore, to Setup that Solver sends its messages. This partitioning of responsibility may not sound like much but it is, in fact, crucial to the simplification of complex code and is one of the key characteristics of object-oriented programming.

We shall first describe the main features of CSetup and, then, say a few words about building the dialog using Constructor.

In this project, the principal differences between CSetup and the default LDialogBox, from which it inherits, are found in the methods FinishCreateSelf() and ListenToMessage(). To FinishCreateSelf(), we add some code to tell the dialog box that it has more than just the OK and Cancel buttons. We also grab a pointer to the application's Flags variable. This pointer is later passed to Solver. Thus, there is only a single Flags variable, each bit of which is a separate flag. (See Flags.h.)

ListenToMessage() does all the real work. Whenever you click in Setup, you send a message. As Figure 1 shows, there are buttons and checkboxes to click and editfields in which to type. (Your Appearance may differ.) The editfield labeled "Steps to Go" doubles as a progress indicator during the integration.

Figure 1. Setup/Status Dialog.

In addition to messages received from the dialog box, Setup listens for one of the three messages sent by Solver and takes the appropriate action. Altogether, five messages have been defined and/or customized. (See CSetup.cp. )

In the ListenToMessage() function, DONE1 and PHASE2 refer to the first and second phases of the integration. As the figure indicates, you can enter any timestep, delta_t, in the topmost editfield. However, the datafile read in at the start provides only a very tiny step, dt = 0.001 days. In the first phase of the integration, the integrator is repeatedly called, using h = dt, until the planetary configuration at your chosen value of delta_t is reached. The default for delta_t is 1 day. Therefore, it takes a thousand cycles just to finish this first step. You'll need a bit of patience! While you are waiting, you can figure out how to add a progress indicator for this first phase. PHASE2 refers to the primary time interval, specified in the dialog box. By default, this contains 365 steps/cycles, with the progress indicator counting down by tens. The planetary configuration is appended to file SolarSys.out every 30 days.

Our custom dialog box, created by cmd_New, is built using the resource editor Constructor. Constructor allows you to create GUI elements for PowerPlant applications without any coding whatsoever. Not only that, these resources may be subsequently edited, e.g. translated into Chinese, without requiring any access to the source code.

Open the file PPStarterResources.ppob and delete the window titled <replace me>. Now, with Windows & Views selected, create a New Layout Resource of Kind LDialogBox. Name it Setup and give it a Resource ID of 1000. (The latter constant, and all the rest that we will need, are defined in the file CSetup.h.) Double click the Setup dialog that appears and modify its Property Inspector window as follows: Change the Class ID to our custom ID, SetU (case-sensitive). Position the top at 60 and the left at 20. Change the height to 180 and the Window proc to Moveable modal, with a title of Setup/Status. Last, change the Default (OK) and Cancel button IDs to 1001 and 1002, respectively.

Save As SolarSysApp.ppob and exit Constructor. At this point, the project should Run without error.

Finish the project by reopening the .ppob file and inserting the remaining dialog items:

Drag-and-drop the two LStdButtons (pbut) from the Controls submenu in the Catalog window. Enter Pane IDs of 1001 and 1002 along with the OK and Cancel titles. Likewise, change their Value Messages to msg_OK and msg_Cancel.

The rest of the items are in the Panes submenu, beginning with the ten checkboxes (cbox). Add the appropriate titles and be sure to disable the cbox for the Sun. (See Figure 1.) Make all of them Initially Checked, and give them the Pane IDs listed in CSetup.h.

Similarly, insert the three LEditField (edit) items with their corresponding IDs. Enter the initial text and, if you wish, the maximum (number of) characters. For the two integer editfields, set the key filter to Integer.

Finally, add the four LCaption (capt, formerly, statText) items. Since these are never referenced, their default IDs are adequate. If you want to match the styles in Figure 1, create a new Text Traits resource for a centered and underlined SystemFont and use it for the captions.

That's it! If everything has gone according to plan, the application should behave just like any normal Macintosh app. Double-click the icon and you should see the new dialog box. Cancel, and the app should quit, just as CSetup::ListenToMessage() says it should.

If you click OK or hit Return instead, then, about 80 seconds later (depending upon the speed of your computer), there will be a double beep and the progress indicator will start to count down to zero. Half a minute later, the application will beep and quit.

Ideally, the output file should contain the exact coordinates for the Sun and planets on 7/4/97 but, as you might expect, there will be some error. The correct answers are given in Table 1.

Table 1. Solar System Coordinates (4 July 1997)

Most of planetary motion occurs in the X,Y plane so most of the error occurs there as well. Also, planets that move faster will have more error than those that move more slowly. Were this a serious attempt to get an accurate answer, we would have made the time step adaptive and specific to each body-coordinate pair, getting smaller as error estimates increased, and vice versa. We also would have used a much more accurate technique, with adaptive order and, perhaps, even adaptive degree in our integration.

You can experiment with all of the options given in the Setup/Status dialog. It's fun to see what happens when planets are ignored (by unchecking them) or when step sizes change.

...or, perhaps, you had an idea of your own that you would like to try?

SimpleArrays.h
// A very simple pair of template classes for vectors and matrices

#pragma once

template<class Type>
class CVector {
public:
   // Constructors
   CVector();
   CVector(long firstIndex, long lastIndex);
   CVector(const CVector<Type> &original);
   // Destructor
   ~CVector();
   // Delete the old vector and make a new one.
   void   resize(long firstIndex, long lastIndex);
   // Subscript implementation with array-bounds check
   Type   &operator[] (long index);
   const Type   &operator[] (long index) const;
   // Basic accessors
   long   getFirstIndex() {return lower;};
   long   getLastIndex() {return upper;};
   unsigned long   size() {return Nelements;};
private:
   Type   *head;
   unsigned long   Nelements;
   long   upper, lower;
};

template <class Type>
class CMatrix {
public:
   // Constructors
   CMatrix();
   CMatrix(long firstRowIndex, long lastRowIndex,
           long firstColIndex, long lastColIndex);
   CMatrix(const CMatrix<Type> &original);
   // Destructor
   ~CMatrix();
   // Delete the old matrix and make a new one.
   void   resize(long firstRowIndex, long lastRowIndex,
                 long firstColIndex, long lastColIndex);
   // Subscript implementation with array-bounds check
   CVector<Type>   &operator[] (long index);
   const CVector<Type>   &operator[] (long index) const;
   // Basic accessors
   long   getFirstRowIndex() {return lowerVec;};
   long   getLastRowIndex() {return upperVec;};
   long   getFirstColIndex() {return lowerColNdx;};
   long   getLastColIndex() {return upperColNdx;};
   unsigned long   size() {return Nelements;};
private:
   // Use a vector of vectors to store 2-D data
   CVector< CVector<Type> * >   vectorPtr;
   unsigned long   Nelements;
   long   lowerVec, upperVec;
   long   lowerColNdx, upperColNdx;
};

const long   DEFAULT_LOW = 1;
const long   DEFAULT_HIGH = 10;

SimpleArrays.cp
// Implementation of CVector and CMatrix

#include <stdlib.h>
#include "SimpleArrays.h"

#ifndef Assert_            // when not using PowerPlant
#include <assert.h>
#define   Assert_ assert
#endif

template<class Type>
CVector<Type>::CVector()
{
   lower = DEFAULT_LOW; upper = DEFAULT_HIGH;
   Nelements = upper - lower + 1;
   head = new Type[Nelements];
   Assert_(head != nil);
}

template<class Type>
CVector<Type>::CVector(long firstIndex, long lastIndex)
{
   Assert_(lastIndex >= firstIndex);
   lower = firstIndex; upper = lastIndex;
   Nelements = upper - lower + 1;
   head = new Type[Nelements];
   Assert_(head != nil);
}

template<class Type>
CVector<Type>::CVector(const CVector<Type> &original)
{
   lower = original.lower; upper = original.upper;
   Nelements = original.Nelements;
   head = new Type[Nelements];
   Assert_(head != nil);
   ::BlockMoveData(&original[lower], head, 
         Nelements*sizeof(Type));
}

template<class Type>
CVector<Type>::~CVector()
{
   if (head) delete []head;
   head = nil;
}

//   Note: If lastIndex < firstIndex, all dynamic memory will be freed.
template<class Type>
void CVector<Type>::resize(long firstIndex, long lastIndex)
{
   if (head) delete []head;
   head = nil;
   lower = firstIndex;
   upper = lastIndex;
   if (upper >= lower) {
      Nelements = upper - lower + 1;
      head = new Type[Nelements];
      Assert_(head != nil);
   }
   else Nelements = 0;
}

template<class Type>
Type &CVector<Type>::operator[] (long index)
{
   Assert_((index >= lower) && (index <= upper));
   return head[index - lower];
}

template<class Type>
const Type &CVector<Type>::operator[] (long index) const
{
   Assert_((index >= lower) && (index <= upper));
   return head[index - lower];
}

template <class Type>
CMatrix<Type>::CMatrix()
{
   lowerVec = lowerColNdx = DEFAULT_LOW; 
   upperVec = upperColNdx = DEFAULT_HIGH;
   // Allocate one row
   for (long r = lowerVec; r <= upperVec; r++) {
      vectorPtr[r] = 
            new CVector<Type> (lowerColNdx, upperColNdx);
      Assert_(vectorPtr[r] != nil);
   }
   Nelements = (upperVec - lowerVec + 1)*
      (upperColNdx - lowerColNdx + 1);
}

template <class Type>
CMatrix<Type>::CMatrix(long firstRowIndex, long lastRowIndex,
                                 long firstColIndex, long lastColIndex)
{
   vectorPtr.resize(firstRowIndex, lastRowIndex);
   lowerVec = firstRowIndex;
   upperVec = lastRowIndex;
   lowerColNdx = firstColIndex;
   upperColNdx = lastColIndex;
   for(long r = lowerVec; r <= upperVec; r++) {
      vectorPtr[r] = 
         new CVector<Type> (lowerColNdx, upperColNdx);
      Assert_(vectorPtr[r] != nil);
   }
   Nelements = (upperVec - lowerVec + 1)*
                     (upperColNdx - lowerColNdx + 1);
}

template <class Type>
CMatrix<Type>::CMatrix(const CMatrix<Type> &original)
{
   lowerVec = original.lowerVec;
   upperVec = original.upperVec;
   lowerColNdx = original.lowerColNdx;
   upperColNdx = original.upperColNdx;
   vectorPtr.resize(lowerVec, upperVec);
   for(long r = lowerVec; r <= upperVec; r++) {
      vectorPtr[r] = new CVector<Type> (original[r]);
      Assert_(vectorPtr[r] != nil);
   }
   Nelements = (upperVec - lowerVec + 1)*
      (upperColNdx - lowerColNdx + 1);
}

template <class Type>
CMatrix<Type>::~CMatrix()
{
   for (long r = lowerVec; r <= upperVec; r++) {
      delete vectorPtr[r];
      vectorPtr[r] = nil;
   }
}

template<class Type>
void CMatrix<Type>::resize(long firstRowIndex,
                                       long lastRowIndex,
                                       long firstColIndex,
                                       long lastColIndex)
{
   // Delete old vectors (rows)
   for (long r = lowerVec; r <= upperVec; r++) {
      delete vectorPtr[r];
      vectorPtr[r] = nil;
   }
   // Resize Nrows
   vectorPtr.resize(firstRowIndex, lastRowIndex);
   lowerVec = firstRowIndex;
   upperVec = lastRowIndex;
   lowerColNdx = firstColIndex;
   upperColNdx = lastColIndex;
   // Allocate new vectors (rows)
   for (long r = lowerVec; r <= upperVec; r++) {
      vectorPtr[r] = 
         new CVector<Type> (firstColIndex, lastColIndex);
      Assert_(vectorPtr[r] != nil);
   }
   Nelements = (upperVec - lowerVec + 1)*
                     (upperColNdx - lowerColNdx + 1);
}

template<class Type>
CVector<Type> &CMatrix<Type>::operator[] (long row)
{
   Assert_((row >= lowerVec) && (row <= upperVec));
   return *vectorPtr[row];
}

template<class Type>
const 
CVector<Type> &CMatrix<Type>::operator[] (long row) const
{
   Assert_((row >= lowerVec) && (row <= upperVec));
   return *vectorPtr[row];
}

CSolver.h

#pragma once

#include "SimpleArrays.h"

// Solver Output modes
enum   {NORMAL, ERROR, IGNORING};

// CSolver is an "idler" LPeriodical which broadcasts status messages.
class   CSolver : public LPeriodical, public LBroadcaster {
public:
// These members are generic to all solvers:
   CSolver(LListener *setup);
   virtual   ~CSolver() {};
   void      SpendTime(const EventRecord &inEvt);
      void   doInput();
      void   doOutput(int mode = NORMAL);
private:
   void   Initialize();   // allocate arrays, etc.
   void   oneCycle();      // called by SpendTime()
   Uint32   *flags;
// These members are problem-specific:
   void   getSun(CMatrix<double> &State);
   void   getAcceleration(CMatrix<double> &State,
                                    CMatrix<double> &accel);
   Boolean   Integrate(double h);
   Boolean   Converged(CMatrix<double> State1,
                              CMatrix<double> State2);
   CVector< CMatrix<double> > state, accels;
   CMatrix<double> initial_state;
   CVector<double> GM;   // G*mass
   double   time, delta_t;
   short   last, current, next;
   enum   {Sun, Mercury, Venus, E_M, Mars,
               Jupiter, Saturn, Uranus, Neptune, Pluto};
   enum   {X, Y, Z};
};
const double TOLERANCE = 1e-12,   // precision
                  dt = 1e-3;               // a tiny time interval (days)
const   int   LINESIZE = 80;             // for doInput()

CSolver.cp

#include   <math.h>
#include   "CSolver.h"
#include   "CSetup.h"
#include   "Flags.h"
// Do NOT include template code in the project.
#include   "SimpleArrays.cp"

// the constructor
CSolver::CSolver(LListener* setup)
{
   AddListener(setup);
   // copy info from Setup
   CSetup   *S = (CSetup*) setup;
   flags = S->flags;
   delta_t = S->delta_t;
   Initialize();
}
// called when CPU is idle
void
CSolver::SpendTime(const EventRecord &inEvt) {
#pragma unused (inEvt)
   if (*flags & READY) {
      *flags &= ~READY;   // turn flag off
      oneCycle();
   }
}
// Get ready to begin.
void
CSolver::Initialize()
{
   GM.resize(Sun, Pluto);
   initial_state.resize(Sun, Pluto, X, Z);
   state.resize(-1, 1);
   accels.resize(-1, 1);
   for (int s = -1;s <= 1;s++) {
      state[s].resize(Sun, Pluto, X, Z);
      accels[s].resize(Mercury, Pluto, X, Z);
   }
   last = -1; current = 0; next = 1;
   doInput();
   getAcceleration(state[last], accels[last]);
   getAcceleration(state[current], accels[current]);
   time = 0.0;
   *flags |= READY;
   StartIdling();
}

// This is the basic control function.
void
CSolver::oneCycle()
{
   double   h = (*flags & PHASE2) ? delta_t : dt;
   
   if ((h == dt) && ((time + h) > delta_t))
      h = delta_t - time;
   if (Integrate(h)) {   // converged
      time += h;
      if (!(*flags & PHASE2))
         if (fabs(time - delta_t) < TOLERANCE) {
            // restore initial state
            for (long p = Sun;p <= Pluto;p++)
               for (long q = X;q <= Z;q++)
                  state[last][p][q] = initial_state[p][q];
            getAcceleration(initial_state, accels[last]);
            BroadcastMessage(msg_DONE1);
            return;
         }
      BroadcastMessage(msg_CONT);
   }
   else BroadcastMessage(msg_DVRG);
}

// Use Numerov's method to solve the differential equations, viz.,
//   q[k+1] = 2*q[k] - q[k-1] + (h^2/12)*(F[k+1] + 10*F[k] + F[k-1]) + R
//      where q'' = F(t, q [but NOT q']) and R = -h^6*q''''''/240.
Boolean
CSolver::Integrate(double h)
{
   CMatrix<double>   trial_state[2] = 
         {state[current], state[next]};   // copies
   double   fac = h*h/12;
   const   long   MaxIter = 10;
   long   iter = 0;
   short   ndx = 0;   // flip-flop

   for (long p = Mercury;p <= Pluto;p++)
      for (long q = X;q <= Z;q++)
         accels[next][p][q] = accels[current][p][q];
   // Iterate to convergence.
   do {
      for (long p = Mercury;p <= Pluto;p++)
         for (long q = X;q <= Z;q++)
            trial_state[1-ndx][p][q] = 
                  2*state[current][p][q] - state[last][p][q] +
                  fac*(accels[next][p][q] + 
                  10*accels[current][p][q] + accels[last][p][q]);
      getSun(trial_state[1-ndx]);
      getAcceleration(trial_state[1-ndx], accels[next]);
      // Divergence check
      if (++iter > MaxIter) return false;
      // Get ready to try again.
      ndx = 1 - ndx;
   } while (!Converged(trial_state[0], trial_state[1]));
   // Save result.
   for (long p = Sun;p <= Pluto;p++)
      for (long q = X;q <= Z;q++)
         state[next][p][q] = trial_state[ndx][p][q];
   last = current; current = next++;   // circular queue
   if (next > 1) next = -1;
   return true;   // converged
}

// Compute RHS of differential equation.
void
CSolver::getAcceleration(CMatrix<double> &State, 
                                    CMatrix<double> &accel)
{
   CVector<double>   factor(Sun, Pluto);
   double   dx, dy, dz, radius;
   long      p, pp;
   Uint32   ignore = MERCURY;

   for (p = Mercury;p <= Pluto;p++, ignore <<= 1) {
      //   factor[pp] == GM[pp]/radius[p][pp]^3
      for (pp = Sun;pp <= Pluto;pp++) {
         if ((*flags & ignore) || (pp == p))
            factor[pp] = 0.0;
         else {
            dx = State[pp][X] - State[p][X];
            dy = State[pp][Y] - State[p][Y];
            dz = State[pp][Z] - State[p][Z];
            radius = sqrt(dx*dx + dy*dy + dz*dz); 
            factor[pp] = GM[pp]/(radius*radius*radius);
         }
      }
      for (long q = X;q <= Z;q++) {
         accel[p][q] = 0.0;
         for (long pp = Sun;pp <= Pluto;pp++)
            accel[p][q] += 
                  factor[pp]*(State[pp][q] - State[p][q]);
      }
   }
}

// Keep the center-of-mass at the origin.
void
CSolver::getSun(CMatrix<double> &State)
{
   Uint32   ignore;

   for (long q = X;q <= Z;q++) {
      State[Sun][q] = 0.0;
      ignore = MERCURY;
      for (long p = Mercury;p <= Pluto;p++, ignore <<= 1)
         if (!(*flags & ignore))
            State[Sun][q] -= GM[p]*State[p][q];
      State[Sun][q] /= GM[Sun];
   }
}

// Are we done yet?
Boolean
CSolver::Converged(CMatrix<double> State1, 
                           CMatrix<double> State2)
{
   Uint32   ignore = MERCURY;
   Boolean   converged = true;
   
   for (long p = Mercury;converged && (p <= Pluto);
            p++, ignore <<= 1) {
      if (*flags & ignore) continue;
      for (long q = X;q <= Z;q++)
         if (fabs(State1[p][q] - State2[p][q]) > TOLERANCE)
            converged = false;
   }
   return converged;
}

CSolver_I_O.cp
// Input/output functions, part of CSolver

#include   <fstream>
#include   "CSetup.h"
#include   "CSolver.h"

// Read the input file.
void
CSolver::doInput()
{
   char   label[LINESIZE];
   ifstream   input("init.dat");

   if (!input.is_open()) {
      ::SysBeep(1);
      ::ExitToShell();
   }
   input.getline(label, LINESIZE);
   for (int p = Sun;p <= Pluto;p++)
      for (int q = X;q <= Z;q++) {
         input >> initial_state[p][q];
         state[last][p][q] = initial_state[p][q];
      }
   input.getline(label, LINESIZE);   // get the <CR>
   input.getline(label, LINESIZE);   // get the next line
   for (int p = Sun;p <= Pluto;p++)
      for (int q = X;q <= Z;q++)
         input >> state[current][p][q];
   input.getline(label, LINESIZE);
   input.getline(label, LINESIZE);
   for (int p = Sun;p <= Pluto;p++)
      input >> GM[p];
   input.close();
}


// Append the current configuration, or a message.
void
CSolver::doOutput(int mode)
{
   ofstream
      output("SolarSys.out", ios_base::out | ios::app);
   output.setf(ios_base::scientific | ios_base::showpos);
   static char   body[10][8] =
               {"Sun", "Mercury", "Venus", "E_M", "Mars",
                "Jupiter","Saturn", "Uranus", "Neptune", "Pluto"};

   switch (mode) {
      case NORMAL:
         output << time << endl;
         for (long p = Sun;p <= Pluto;p++) {
            for (long q = X;q <= Z;q++)
               output << state[current][p][q] << " ";
            output << endl;
         }
         break;
      case ERROR:
         output << time << endl;
         output << "Diverging!!" << endl;
         break;
      case IGNORING:
         Uint32   ignore = SUN;
         output << "Ignoring ";
         if (*flags < ignore) output << "none";
         else
            for (long p = Sun;p <= Pluto;p++, ignore <<= 1)
               if (*flags & ignore)
                  output << body[p] << '\t';
         output << endl << endl;
   }
   output.close();
}

CSetup.h

#pragma once

#include   "SolarSysApp.h"
#include   "CSolver.h"

class CSetup : public LDialogBox {
public:
   enum { class_ID = 'SetU' };
   static    CSetup*   CreateSetupStream( LStream *inStream );
               CSetup( LStream *inStream );
               ~CSetup();
   // CSolver contains all necessary
   //   knowledge for solving this particular problem
   CSolver   *Solver;
   double   delta_t;
   Uint32   *flags;
private:
   virtual   void   FinishCreateSelf();
   virtual   void
               ListenToMessage( MessageT inMessage, 
                                       void *ioParam);
   Int32   Nsteps, Nreport, Step;
// checkboxes
enum
{SU = 1010, ME, VE, EM, MA, JU, SA, UR, NE, PL,/* checkboxes */
   DT, NS, RP};   // editfields
};
// Resource constants
const   ResIDT   rRidL_Setup   = 1000;
const   PaneIDT   Setup_OK = 1001,
                     Setup_Cancel = 1002;
// Messages
const   MessageT   msg_CONT = 2000,
                        msg_DONE1 = 2001,
                        msg_DVRG = 2002;

CSetup.cp
// CSetup.cp -- a custom dialog (modified LDialogBox code)

#include <fp.h>
#include "CSetup.h"

CSetup*
CSetup::CreateSetupStream(LStream *inStream)
{
   return new CSetup(inStream);
}

CSetup::CSetup(LStream *inStream):LDialogBox(inStream)
{
   // called by CreateSetupStream()
}

CSetup::~CSetup()
{
}

void
CSetup::FinishCreateSelf()
{
   // get a pointer to the "global" Flags
   SolarSysApp   *theApp = (SolarSysApp*) mSuperCommander;
   flags = &theApp->Flags;
   
   UReanimator::LinkListenerToControls(this, this, 
                                                      rRidL_Setup);

   // Call inherited LDialogBox's FinishCreateSelf()
   //    sets up the default and cancel buttons.
   LDialogBox::FinishCreateSelf();
}

void
CSetup::ListenToMessage(MessageT inMessage, void *ioParam)
{
   const      Int32   statusRpt = 10;   // to update display
   static   Int32   disp;
   switch (inMessage) {
      case msg_OK:
         LButton   *OK_Button = 
               (LButton*) FindPaneByID(Setup_OK);
         OK_Button->Disable();
         // What planets are ignored?
         Uint32   temp_flag = MERCURY;
         for (PaneIDT k = ME;k <= PL;k++) {
            LControl   *ckbox = (LControl*) FindPaneByID(k);
            if (ckbox->GetValue() == 0) *flags |= temp_flag;
            temp_flag <<= 1;
         }
         // get delta_t, Nsteps, and reporting interval
         LEditField   *edt = (LEditField*) FindPaneByID(DT);
         LStr255      aString;
         edt->GetDescriptor(aString);
         delta_t = (double_t) aString;
         if (delta_t < dt) {   // override
            delta_t = dt;
            aString.Assign(dt, FLOATDECIMAL, 2);
            edt->SetDescriptor(aString);
         }
         edt = (LEditField*) FindPaneByID(NS);
         Nsteps = edt->GetValue();
         edt = (LEditField*) FindPaneByID(RP);
         Nreport = edt->GetValue();
         // create solver and launch
         Solver = new CSolver(this);
         Solver->doOutput(IGNORING);
         *flags |= TRYING;
         break;
      case msg_Cancel:
         mSuperCommander->ObeyCommand(cmd_Quit);
         break;
      case msg_CONT:
         if (*flags & PHASE2) {
            if (((++Step % Nreport) == 0) || (Step == Nsteps))
               Solver->doOutput();
            if (Step == disp) {
               LStr255      aString(Nsteps - disp);
               LEditField   *edt = (LEditField*) FindPaneByID(NS);
               edt->SetDescriptor(aString);
               disp += statusRpt;
               if (disp > Nsteps) disp = Nsteps;
            }
            if (Step == Nsteps) {   // All done!
               ::SysBeep(1);
               mSuperCommander->ObeyCommand(cmd_Quit);
            }
         }
         *flags |= READY;
         break;
      case msg_DONE1:
         *flags |= (PHASE2 | READY);
         Step = 1; disp = statusRpt;
         Solver->doOutput();
         ::SysBeep(1); ::SysBeep(1);
         break;
      case msg_DVRG:
         Solver->doOutput(ERROR);
         mSuperCommander->ObeyCommand(cmd_Quit);
         break;
      default:
         // Pass any others to inherited ListenToMessage.
         LDialogBox::ListenToMessage(inMessage, ioParam);
   }
}

SolarSysApp.h
// SolarSysApp.h -- the main module

#pragma once

#include   <LApplication.h>
#include   "Flags.h"

class   SolarSysApp : public LApplication {

public:
   SolarSysApp();   // constructor registers all PPobs
   virtual   ~SolarSysApp();
   
   // this overriding function performs application functions
   virtual Boolean   ObeyCommand(CommandT inCommand,
                                             void* ioParam);
   
   // this overriding function returns the status of menu items
   virtual void   FindCommandStatus(CommandT inCommand,
                     Boolean &outEnabled, Boolean &outUsesMark,
                     Char16 &outMark, Str255 outName);
   Uint32   Flags;
protected:
   virtual void      StartUp();   // overriding startup function
};

Flags.h
//========================================================
//   Flags.h -- each flag is one bit
//
//   To set a flag:         Flags |= thisFlag;
//   To reset a flag:      Flags &= ~thisFlag;
//   To toggle a flag:      Flags ^= thisFlag;
//   To test a flag:         enabled = Flags & thisFlag;
//   To initialize Flags:   Flags = 0;
//                        or
//                        Flags = thisFlag1 [ | thisFlag2 ...]
//========================================================

#pragma once

// App flags
const Uint32   WORKING   = 0x00000001;   // Setup dialog exists
const Uint32   TRYING   = 0x00000002;   // finished Setup
const Uint32   DONE      = 0x00000004;   // can Save now

// Solver flags
const Uint32   READY      = 0x00000010;   // for new cycle
const Uint32   PHASE2   = 0x00000020;   // in main time interval

// Ignore flags
const Uint32   SUN         = 0x00100000;
const Uint32   MERCURY   = 0x00200000;
const Uint32   VENUS      = 0x00400000;
const Uint32   EMBARY   = 0x00800000;
const Uint32   MARS      = 0x01000000;
const Uint32   JUPITER   = 0x02000000;
const Uint32   SATURN   = 0x04000000;
const Uint32   URANUS   = 0x08000000;
const Uint32   NEPTUNE   = 0x10000000;
const Uint32   PLUTO      = 0x20000000;

init.dat
initial state (JD2442963.5)
-7.826853112272E-04 -3.154354825996E-03 -1.347883896562E-03
+2.363722282085E-01 +1.925795112940E-01 +7.858575760206E-02
-2.814426096269E-01 +5.934177402851E-01 +2.847619979908E-01
+2.188361276660E-01 -9.139369488643E-01 -3.962761751021E-01
-1.650071503098E+00 +1.178289968881E-01 +9.882482675934E-02
+3.577768183934E+00 +3.214315993866E+00 +1.290688586589E+00
-5.339976482385E+00 +6.719990402690E+00 +3.004433799763E+00
-1.492977308373E+01 -1.013274418856E+01 -4.226360918526E+00
-8.802641438020E+00 -2.690623936169E+01 -1.079386542938E+01
-2.870849532338E+01 -8.815842320049E+00 +5.898219074445E+00
JD2442963.501 (t0 + dt)
-7.826790930626E-04 -3.154358772120E-03 -1.347885783055E-03
+2.363479684713E-01 +1.925982738720E-01 +7.859829766764E-02
-2.814612891282E-01 +5.934099950463E-01 +2.847596969258E-01
+2.188526517012E-01 -9.139336022900E-01 -3.962747240815E-01
-1.650072244682E+00 +1.178174001393E-01 +9.881952799675E-02
+3.577762846591E+00 +3.214321253281E+00 +1.290690971220E+00
-5.339981299311E+00 +6.719987284652E+00 +3.004432719300E+00
-1.492977078308E+01 -1.013274724759E+01 -4.226362290905E+00
-8.802638454393E+00 -2.690624016261E+01 -1.079386583145E+01
-2.870849425395E+01 -8.815845339411E+00 +5.898217810075E+00
mass factors
2.959122082856E-04
4.912547451451E-11
7.243456209633E-10
8.997011658557E-10
9.549528942224E-11
2.825342103446E-07
8.459468504831E-08
1.288816238138E-08
1.532112481284E-08
2.276247751864E-12

Michael McLaughlin, mpmcl@mitre.org, a former chemistry professor and Peace Corps volunteer, currently does R&D for future Air Traffic Control systems. He has been programming computers since 1965 but has long since forsaken Fortran, PLI, and Lisp in favor of C++ and assembly.