TweetFollow Us on Twitter

Scientific Prog in PP

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

An Array of Excuses

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

A Scientific Problem

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.

Building the Application

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.

Building the SetupWindow Resource

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.

Solar System, 4 July 1977

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)

BodyXYZDist. from Sun
Sun+1.864677e-03-3.988760e-03 -1.784555e-03n/a
Mercury-1.704775e-01+2.299152e-01 +1.410340e-013.237438e-01
Venus +6.740839e-01 -2.398561e-01 -1.504364e-01 7.277428e-01
Earth_Moon +2.172080e-01 -9.155892e-01 -3.970636e-01 1.016678e+00
Mars +1.355082e+00 +3.693063e-01 +1.327712e-01 1.410195e+00
Jupiter +1.182357e+00 +4.544504e+00 +1.919276e+00 5.076694e+00
Saturn -6.944600e+00 +5.417624e+00 +2.535628e+00 9.169827e+00
Uranus -1.404475e+01 -1.121683e+01 -4.713703e+00 1.858057e+01
Neptune -7.707629e+00 -2.717941e+01 -1.093293e+01 3.028910e+01
Pluto -2.829839e+01 -9.911559e+00 +5.432751e+00 3.047292e+01

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.

 

Community Search:
MacTech Search:

Software Updates via MacUpdate

Firefox 37.0 - Fast, safe Web browser. (...
Firefox offers a fast, safe Web browsing experience. Browse quickly, securely, and effortlessly. With its industry-leading features, Firefox is the choice of Web development professionals and casual... Read more
Arq 4.11 - Online backup to Google Drive...
Arq is super-easy online backup for the Mac. Back up to your own Google Drive storage (15GB free storage), your own Amazon Glacier ($.01/GB per month storage) or S3, or any SFTP server. Arq backs up... Read more
MacFamilyTree 7.3.4 - Create and explore...
MacFamilyTree gives genealogy a facelift: it's modern, interactive, incredibly fast, and easy to use. We're convinced that generations of chroniclers would have loved to trade in their genealogy... Read more
Yummy FTP 1.10.2 - FTP/SFTP/FTPS client...
Yummy FTP is an FTP + SFTP + FTPS file transfer client which focuses on speed, reliability and productivity. Whether you need to transfer a few files or a few thousand, schedule automatic backups, or... Read more
VueScan 9.5.08 - Scanner software with a...
VueScan is a scanning program that works with most high-quality flatbed and film scanners to produce scans that have excellent color fidelity and color balance. VueScan is easy to use, and has... Read more
Iridient Developer 3.0.1 - Powerful imag...
Iridient Developer (was RAW Developer) is a powerful image conversion application designed specifically for OS X. Iridient Developer gives advanced photographers total control over every aspect of... Read more
Air Video Server HD 2.1.0 - Stream video...
Air Video Server HD streams videos instantly from your computer on your iPhone, iPad, iPod touch or Apple TV. No need to worry about converting or transferring files. We took everything that was... Read more
Duplicate Annihilator 5.7.5 - Find and d...
Duplicate Annihilator takes on the time-consuming task of comparing the images in your iPhoto library using effective algorithms to make sure that no duplicate escapes. Duplicate Annihilator... Read more
BusyContacts 1.0.2 - Fast, efficient con...
BusyContacts is a contact manager for OS X that makes creating, finding, and managing contacts faster and more efficient. It brings to contact management the same power, flexibility, and sharing... Read more
Capture One Pro 8.2.0.82 - RAW workflow...
Capture One Pro 8 is a professional RAW converter offering you ultimate image quality with accurate colors and incredible detail from more than 300 high-end cameras -- straight out of the box. It... Read more

2K Announces WWE 2K, Mobile's First...
It seems like this month has been pretty big for wrestling. First Wrestlemania, then 2K has announces that they're releasing  WWE 2K for iOS. It's a simulation-based WWE game where you'll get to play with several WWE superstars such as John Cena, ... | Read more »
How the Apple Watch Could Change the Fac...
The Apple Watch is still a ways out, but my previous musings on the wearable’s various features got me thinking: what might it be like a year after launch? Two years? Five years? What if it becomes a symbiotic part of the iOS framework to the point... | Read more »
Pie In The Sky: A Pizza Odyssey (Games)
Pie In The Sky: A Pizza Odyssey 1.0 Device: iOS Universal Category: Games Price: $2.99, Version: 1.0 (iTunes) Description: A game about delivering pizza. In space. | Read more »
Chosen Gives Hopeful Singers, Songwriter...
If YouTube videos and reality TV shows like The Voice have taught us one thing, it’s that there are a lot of people out there who are anxious to show the world their talents. And if they’ve taught us a second thing, it’s that there’s an almost... | Read more »
Android's Popular OfficeSuite Now A...
Once only available for Android devices, OfficeSuite has finally landed on the app store. The Mobile Systems app lets you view, edit, create, and share Word, Excel, and PowerPoint documents as well as convert them to/from PDFs. It's touted as being... | Read more »
Warhammer: Arcane Magic is Coming Soon,...
Turbo Tape Games has announced that they're joining forces with Games Workshop to bring the turn-based strategy board game, Warhammer: Arcane Magic, to life on the iOS. | Read more »
Fast & Furious: Legacy's Creati...
| Read more »
N-Fusion and 505's Ember is Totally...
| Read more »
These are All the Apple Watch Apps and G...
The Apple Watch is less than a month from hitting store shelves, and once you get your hands on it you're probably going to want some apps and games to install. Fear not! We've compiled a list of all the Apple Watch apps and games we've been able to... | Read more »
Appy to Have Known You - Lee Hamlet Look...
Being at 148Apps these past 2 years has been an awesome experience that has taught me a great deal, and working with such a great team has been a privilege. Thank you to Rob Rich, and to both Rob LeFebvre and Jeff Scott before him, for helping me... | Read more »

Price Scanner via MacPrices.net

Adobe Brings Powerful Layout-Design Capabilit...
Adobe today announced the availability of Adobe Comp CC, a free iPad app that enables rapid creation of layout concepts for mobile, Web and print projects. With Comp CC, designers can rough out and... Read more
Apple offering refurbished 27-inch 5K iMacs f...
The Apple Store is offering Apple Certified Refurbished 27″ 3.5GHz 5K iMacs for $2119 including free shipping. Their price is $380 off the price of new models, and it’s the lowest price available for... Read more
16GB iPad mini on sale for $199, save $50
Walmart has 16GB iPad minis (1st generation) available for $199.99 on their online store, including free shipping. Their price is $50 off MSRP. Online orders only. Read more
New 128GB MacBook Airs on sale for $50 off MS...
 B&H Photo has 128GB 11″ and 13″ 2015 MacBook Airs on sale today for $50 off MSRP including free shipping plus NY sales tax only: - 11″ 1.6GHz/128GB MacBook Air (Model #MJVM2LL/A): $849 $50 off... Read more
13-inch 2.6GHz Retina MacBook Pro (refurbishe...
The Apple Store has Apple Certified Refurbished 13″ 2.6GHz/128GB Retina MacBook Pros available for $979 including free shipping. Original MSRP for this model was $1299. Read more
Save up to $600 with Apple refurbished Mac Pr...
The Apple Store is offering Apple Certified Refurbished Mac Pros for up to $600 off the cost of new models. An Apple one-year warranty is included with each Mac Pro, and shipping is free. The... Read more
Samsung Galaxy S 6 and Galaxy S 6 edge U.S. P...
Samsung Electronics America, Inc. has announced the Galaxy S 6 and Galaxy S 6 edge will be available in the U.S. beginning April 10, with pre-orders being accepted now. “We have completely reimagined... Read more
13-inch 2.5GHz MacBook Pro (refurbished) avai...
The Apple Store has Apple Certified Refurbished 13″ 2.5GHz MacBook Pros available for $829, or $270 off the cost of new models. Apple’s one-year warranty is standard, and shipping is free: - 13″ 2.... Read more
Save up to $80 on iPad Air 2s, NY tax only, f...
 B&H Photo has iPad Air 2s on sale for $80 off MSRP including free shipping plus NY sales tax only: - 16GB iPad Air 2 WiFi: $469.99 $30 off - 64GB iPad Air 2 WiFi: $549.99 $50 off - 128GB iPad... Read more
iMacs on sale for up to $205 off MSRP
B&H Photo has 21″ and 27″ iMacs on sale for up to $205 off MSRP including free shipping plus NY sales tax only: - 21″ 1.4GHz iMac: $1019 $80 off - 21″ 2.7GHz iMac: $1189 $110 off - 21″ 2.9GHz... Read more

Jobs Board

DevOps Software Engineer - *Apple* Pay, iOS...
**Job Summary** Imagine what you could do here. At Apple , great ideas have a way of becoming great products, services, and customer experiences very quickly. Bring Read more
*Apple* Retail - Multiple Positions (US) - A...
Sales Specialist - Retail Customer Service and Sales Transform Apple Store visitors into loyal Apple customers. When customers enter the store, you're also the Read more
Sr. Technical Services Consultant, *Apple*...
**Job Summary** Apple Professional Services (APS) has an opening for a senior technical position that contributes to Apple 's efforts for strategic and transactional Read more
Lead *Apple* Solutions Consultant - Retail...
**Job Summary** Job Summary The Lead ASC is an Apple employee who serves as the Apple business manager and influencer in a hyper-business critical Reseller's store Read more
*Apple* Pay - Site Reliability Engineer - Ap...
**Job Summary** Imagine what you could do here. At Apple , great ideas have a way of becoming great products, services, and customer experiences very quickly. Bring Read more
All contents are Copyright 1984-2011 by Xplain Corporation. All rights reserved. Theme designed by Icreon.