TweetFollow Us on Twitter

Feb 99 Prog Challenge

Volume Number: 15 (1999)
Issue Number: 2
Column Tag: Programmers Challenge

Feb 99 Programmers Challenge

by Bob Boonstra, Westford, MA

Chinese Checkers

For this month's Challenge, I've reached once again into my stash of board games, this time coming up with one I haven't played in a very long time: Chinese Checkers. The prototype for the code you should write is:

#if defined(__cplusplus)
extern "C" {
#endif

typedef struct Position {
   long row;
   long col;
} Position;

void InitChineseCheckers(
   long numPlayers,   /* 2..6  */
   long gameSize,   /* base of home triangle, 3..63, you have size*(size+1)/2 pieces */
   long playerPosition[6],   /* 0..5, entries 0..numPlayers-1 are valid */
   long yourIndex   /* 0..numPlayers-1, your position is playerPosition[yourIndex] */
);

void YourMove(
   Position *from,   /* originating position */
   Position *to      /* destination position */
);

void OpponentMove(
   long opponent,    /* index in playerPosition[] of the player making move */
   Position from,    /* originating position */
   Position to       /* destination position */
   
);

void TermChineseCheckers(void);

#if defined(__cplusplus)
}
#endif

Chinese checkers is played on a six-pointed star shaped board by 2 to 6 players. The board can be viewed as six equilateral triangles attached to a hexagon. Each Position on the board is connected to six adjacent positions (ignoring boundary conditions), one each to the left, right, upper left, upper right, lower left, and lower right. At the beginning of a game, each player's pieces occupy one of the six triangular areas, referred to as that player's home position. The triangular area opposite the player's home position is called that player's goal position. The triangular areas that are not part of the home or goal positions are called neutral zones for that player. The objective of the game is to move your pieces from your home to your goal before any other player does the same.

Standard Chinese Checker boards use a board where each player has 10 pieces, initially placed in a home triangle with 4 positions at its base. We'll refer to the base of the home triangle as the gameSize. I've also seen games with 6 pieces per player (a gameSize of 3). In our game, the gameSize can be anywhere from 3 to 64, resulting in anywhere from 6 to 2016 pieces per player.

To communicate game play between my test code and your solution, we need a numbering system for the board. The positions on the board are described using row and column values. The board will have 4*gameSize rows, numbered from 0 to 4*gameSize-1. The star shape dictates that row r has the following number of positions:

Rows# of columns
0togameSize-1r+1
gameSizeto2*gameSize-14*gameSize-r+1
2*gameSizeto3*gameSize-1r+1
3*gameSizeto4*gameSize-1r+1-4*gameSize

The columns in row r with c columns are numbered from -(c-1)/2 to +(c-1)/2, inclusive, for even numbered rows, and from -(c/2)+1 to (c/2), inclusive, for odd numbered rows. For a position (r,c) in an even numbered row, the adjacent positions are (r,c-1), (r,c+1), (r+1,c), (r+1,c+1), (r-1,c), and (r-1,c+1). For the same position in an odd-numbered row, the adjacent positions are (r,c-1), (r,c+1), (r+1,c), (r+1,c-1), (r-1,c), and (r-1,c-1). (Provided, of course, that those positions exist on the board).

 0:                           0
 1:                         0   1 
 2:                      -1   0   1 
 3:                    -1   0   1   2
 4:  -6  -5  -4  -3  -2  -1   0   1   2   3   4   5   6 
 5:    -5  -4  -3  -2  -1   0   1   2   3   4   5   6
 6:      -5  -4  -3  -2  -1   0   1   2   3   4   5  
 7:        -4  -3  -2  -1   0   1   2   3   4   5
 8:          -4  -3  -2  -1   0   1   2   3   4  
 9:        -4  -3  -2  -1   0   1   2   3   4   5
10:      -5  -4  -3  -2  -1   0   1   2   3   4   5  
11:    -5  -4  -3  -2  -1   0   1   2   3   4   5   6
12:  -6  -5  -4  -3  -2  -1   0   1   2   3   4   5   6 
13:                    -1   0   1   2
14:                      -1   0   1 
15:                         0   1 
16:                           0

A move in Chinese Checkers can be either a non-jump move or a sequence of contiguous jumps. Relocating a piece from one position to an unoccupied adjacent position is a non-jump move. Jumping over an adjacent occupied position into an unoccupied position, collinear with the first two, is considered a jump. A single jump move can consist of jumping a single piece over an occupied position, followed by a jump from the new location over an adjacent piece to still another location, repeated as many times as the board situation allows. No pieces are removed from the board at any time. A move may not begin or end in a neutral zone triangles, although a multi-jump move may include an intermediate position in the neutral zone. Some Chinese Checker rules allow a player to pass; this is not allowed in our game.

Your InitChineseCheckers routine will be called at the beginning of each game with the game parameters, the number of players in the game (numPlayers), the size of the board (gameSize, described above), the board position for each player (playerPosition[]), and your index as a player in the game (yourIndex). The playerPosition is numbered clockwise starting from the (0,0) position.

After every player has been initialized, the player at playerPosition[0] will be asked to play by calling his YourMove routine. That player should return the Position of the piece he is moving in *from, and the position he is moving to in *to. For a non-jump move, the two positions must be adjacent. For a jump move, there must be a valid sequence of jumps between *from and *to. If you cannot make a legal move, set *from to the location of one of your pieces and *to to that same location. Illegal moves will result in disqualification.

After a player moves, the other players will be notified of the move by calling their OpponentMove routines, providing the index of the moving player in opponent, and the positions moved from and to in those parameters, respectively.

The game is over when one player has moved all of his pieces from the home area to the goal area, provided all players have had the same number of turns to move. At the end of the game, TermChineseCheckers will be called, where you can deallocate any memory allocated during the game.

The evaluation will consist of solutions playing against one another in combinations of between 2 and 6 players per game. Each solution will play the same number of games against each other solution, with the same number of opportunities to play first, with fair seating arrangements that will be determined after I see how many solutions are submitted. The winner will be the solution that wins the most games in the least amount of time. Specifically, the winner will be the solution with the greatest score, where the score is two times the number of games won, plus the number of games drawn, minus a time penalty of 1 game per second of execution time. Execution time will include the time used by all of your routines: InitChineseCheckers, YourMove, OpponentMove, and TermChineseCheckers.

This will be a native PowerPC Challenge, using the latest CodeWarrior environment. Solutions may be coded in C, C++, or Pascal.

Three Months Ago Winner

Congratulations to Tom Saxton for submitting the winning solution to the December, 1998, Rendezvous and Docking Challenge. Contestants were required to maneuver a spacecraft through a collection of orbiting objects to a successful docking with a target. The objective was to minimize both the fuel expended, measured as the total change in the velocity vector throughout the flight, and the amount of execution time used to perform the rendezvous. This was an extraordinarily difficult problem to solve. One person wrote to say that the final docking maneuver would have been a very difficult Challenge all by itself. Tom was the only one persistent enough to complete and submit a correct solution, or any solution at all, for that matter.

The test code consisted of four scenarios in a simulated solar system like our own. The first scenario, which I dubbed "Shuttle-Zarya" after the first space station assembly flight, involved maneuvering a ship in earth orbit to a target in a different earth orbit. The second required a ship in earth orbit to land on the moon. At least I called it a "landing" - the requirements for differential velocity at docking were very generous, so one might better refer to it as a controlled crash. The most demanding scenario, both in terms of the total delta-V required to accomplish the rendezvous and in terms of the amount of wall clock time needed to complete the simulation, was the third, where the ship starting in earth orbit was required to land on the Jovian moon Ganymede. The fourth and final scenario required the ship to leave earth orbit and land on the tiny Martian moon Phobos. In all cases the ship had an infinite supply of fuel available, with an engine that delivered between 10 and 30 meters/sec/sec of acceleration, or roughly 1 to 3 earth gravities..

For transfer between two orbits around the same object, the orbital technique requiring the least amount of velocity change is an elliptical transfer orbit tangent to the original and the final orbits. These orbits are called Hohmann transfer orbits, named after the person first recognizing its usefulness. Of course, these orbits allow a rendezvous only if the target object is at the point of intersection when the ship arrives at that point. This is why launch opportunities for real life Mars probes, for example, only occur at certain times, when the planets are properly aligned. In our problem, the ship and the target were not aligned in any particular way, so the Hohmann transfer did not always apply directly. One could have waited around until the alignment was correct, but Tom chose a different approach.

Tom's approach approximates the N-body problem presented by the Challenge as a collection of 2-body problems, recognizing the fact that the gravitational effect of one body is usually dominant. His solution includes four steps. Step 1 is to break the ship out of the orbit of any object that the target was not orbiting. Step 2 is to enter orbit around the object that the target is orbiting. Step 3 is to intersect the target orbit in the vicinity of the target. The final step is to close on the target and match velocities. The rendezvous state is stored in an ROK ("Rendezvous Orbit Kind") variable, and the state transitions occur in the _Plan routine.

The method used for the first step (escape orbit) is less efficient than it might be, in that the resulting orbit can overshoot the intended transfer orbit, although it gets the job done. This is evident in the Phobos test case, where the escape orbit actually takes the ship well outside the orbit of Mars. Tom comments, in private email, "The problem there is the inefficiency with which I escape the Earth's orbit. I didn't have time to come up with a general way to determine an optimal Hohmann transfer from the starting Earth orbit (or more generally, a starting orbit around a moon of an arbitrary planet). So, I just blast out of orbit when the ship is moving parallel to the parent's velocity (in the same direction if the ultimate target is in a larger orbit, or opposite if the target is in a smaller orbit)... This does yield a stupid transfer to Mars because the ship's escape orbit turns out to cross Mars' orbit. It would be less stupid transferring to one of the outer planets."

Once the ship and the target are in orbit about the same body, Tom computes a transfer orbit in the _FFindHohmannTransfer2 routine. This routine iterates over possible orbital starting positions, and when the iteration detects a sign change in the angular position error, it uses the secant method to converge on the starting point for the transfer. The routine actually computes two transfer orbits, one in the orbital plane of the ship and another in the orbital plane of the target. A plane change maneuver is executed where those orbits intersect, at the 90° point.

At the end of the transfer orbit, when the ship and the target are "close", Tom linearly decreases the relative velocity with the target, where "close" is defined as twice the distance needed to decelerate using maximum thrust. In the Ganymede test case, Tom actually "docked" with Ganymede a little harder than the problem permitted, which I attribute to a combination of the gravitational influence of Jupiter and numerical integration limitations in the test code related to the way the time step interval was calculated. Anyway, when you are the only one to solve an extraordinarily difficult Challenge, some liberties in scoring are justified.

Tom's code refers in places to "Fundamentals of Astrodynamics", by Bate, Mueller, and White. This 1971 Dover Publications text is an excellent reference, one that I have used myself on occasion and can recommend highly for anyone interested in astrodynamics. As I finish this article, I'm actually watching the shuttle Endeavor grapple Zarya and dock it to the Unity module of the International Space Station. Most exciting!

Tom was kind enough to submit his test code in addition to his solution, and allow it to be published on the Challenge ftp site. I've also included a version of my test code, which includes a graphic representation of the progress toward a solution that can be interesting to watch.

The table below lists the total delta-V expended for each of the four scenarios used to evaluate Tom's solution. It also lists the execution time in milliseconds, the score for each scenario, the number of numerical integration steps and the simulated flight time required to complete the rendezvous. For the Ganymede scenario, Tom's ship actually collided with Jupiter before being able to rendezvous with the target. Had that not happened, the delta-V would actually have been larger.

Test CaseTotal delta-VExecution Time (msec)ScoreIntegration StepsFlight Time
Shuttle-Zarya1.037E+036.793E+011.105E+032574.5 hours
Lunar landing1.142E+043.977E+021.182E+0414975.0 days
Ganymede (*)3.873E+051.485E+055.357E+055781835.5 years
Phobos1.462E+052.281E+053.754E+058851922.8 years

Top Contestants

Listed here are the Top Contestants for the Programmer's Challenge, including everyone who has accumulated 20 or more points during the past two years. The numbers below include points awarded over the 24 most recent contests, including points earned by this month's entrants.

  1. Munter, Ernst 204
  2. Saxton, Tom 79
  3. Boring, Randy 56
  4. Mallett, Jeff 50
  5. Rieken, Willeke 47
  6. Cooper, Greg 34
  7. Maurer, Sebastian 40
  8. Heithcock, JG 37
  9. Murphy, ACC 34
  10. Nicolle, Ludovic 27
  11. Lewis, Peter 31
  12. Hart, Alan 21
  13. Antoniewicz, Andy 20
  14. Brown, Pat 20
  15. Day, Mark 20
  16. Higgins, Charles 20
  17. Hostetter, Mat 20
  18. Studer, Thomas 20

There are three ways to earn points: (1) scoring in the top 5 of any Challenge, (2) being the first person to find a bug in a published winning solution or, (3) being the first person to suggest a Challenge that I use. The points you can win are:

1st place20 points
2nd place10 points
3rd place7 points
4th place4 points
5th place2 points
finding bug2 points
suggesting Challenge2 points

Here is Tom's winning Rendezvous solution:

Rendezvous.cpp

Copyright © 1998 Tom Saxton

#define USE_DOUBLES
#include "Rendezvous.h"

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

enum
{
   fFalse=0,
   fTrue=1
};

// disable asserts
#define Assert(f)

// basic constants
#define sEpsilon   (1.0e-6)
#define pi         (3.141592653589793233)
//
// Stupid Vector Tricks
//

static const Vector3D s_vecZero = { 0, 0, 0 };

_ZeroVec
// set a vector to zero
static void _ZeroVec(Vector3D *pvec)
{
   pvec->x = pvec->y = pvec->z = 0;
}

_NegateVec
// negate a vector
static void _NegateVec(Vector3D *pvec)
{
   pvec->x = -pvec->x;
   pvec->y = -pvec->y;
   pvec->z = -pvec->z;
}

_DotVec
// compute dot product of two vectors
static double _DotVec(
   const Vector3D &vec1,
   const Vector3D &vec2
   )
{
   return vec1.x*vec2.x + vec1.y*vec2.y + vec1.z*vec2.z;
}

_LengthVec
// compute length of a vector
static double _LengthVec(const Vector3D &vec)
{
   return sqrt(_DotVec(vec, vec));
}

_ScaleVec
// multiply a vector by a scale into another vector
static void _ScaleVec(
   const double g,
   const Vector3D &vec,
   Vector3D *pvecResult
   )
{
   pvecResult->x = g*vec.x;
   pvecResult->y = g*vec.y;
   pvecResult->z = g*vec.z;
}

_NormalizeVec
// create a normal vector (unit length) in the same direction as vec
static void _NormalizeVec(const Vector3D &vec, 
                                    Vector3D *pvecResult)
{
   double len = _LengthVec(vec);
   Assert(len > sEpsilon);
   _ScaleVec(1.0/len, vec, pvecResult);
}

_SubVec
// subtract two vectors
static void _SubVec(
   const Vector3D &vec1,
   const Vector3D &vec2,
   Vector3D *pvecResult
   )
{
   pvecResult->x = vec1.x - vec2.x;
   pvecResult->y = vec1.y - vec2.y;
   pvecResult->z = vec1.z - vec2.z;
}

_AddVec
// add two vectors
static void _AddVec(
   const Vector3D &vec1,
   const Vector3D &vec2,
   Vector3D *pvecResult
   )
{
   pvecResult->x = vec1.x + vec2.x;
   pvecResult->y = vec1.y + vec2.y;
   pvecResult->z = vec1.z + vec2.z;
}

_AddScaleVec
// add a vector to a scaled vector
static void _AddScaleVec(
   const Vector3D &vec1,
   const double g,
   const Vector3D &vec2,
   Vector3D *pvecResult
   )
{
   pvecResult->x = vec1.x + g*vec2.x;
   pvecResult->y = vec1.y + g*vec2.y;
   pvecResult->z = vec1.z + g*vec2.z;
}

_CrossVec
// compute the cross product of two vectors
static void _CrossVec(const Vector3D &vec1, const Vector3D &vec2, 
                              Vector3D *pvecResult)
{
   pvecResult->x = vec1.y*vec2.z - vec2.y*vec1.z;
   pvecResult->y = vec1.z*vec2.x - vec2.z*vec1.x;
   pvecResult->z = vec1.x*vec2.y - vec2.x*vec1.y;
}

_LinearComboVec
// compute a linear combination of two vectors
static void _LinearComboVec(
   const double g1,
   const Vector3D &vec1,
   const double g2,
   const Vector3D &vec2,
   Vector3D *pvecResult
   )
{
   pvecResult->x = g1*vec1.x + g2*vec2.x;
   pvecResult->y = g1*vec1.y + g2*vec2.y;
   pvecResult->z = g1*vec1.z + g2*vec2.z;
}

_DistBetweenPos
// compute the distance between two vectors
static double _DistBetweenPos(const Vector3D &pos1, 
                                 const Vector3D &pos2)
{
   Vector3D vecSep;
   _SubVec(pos1, pos2, &vecSep);
   return _LengthVec(vecSep);
}

#ifdef DEBUG
_FNormalVec
// does this vector have unit length?
static int _FNormalVec(const Vector3D &vec)
{
   double len = _LengthVec(vec);
   return fabs(len - 1.0) < sEpsilon;
}

_FEqualVec
// are these two vectors equal to within the specified relative tolerance?
static int _FEqualVec(const Vector3D &vec1, const Vector3D &vec2, 
                                       double epsilon)
{
   Vector3D dvec;
   
   _SubVec(vec1, vec2, &dvec);
   return _LengthVec(dvec)/fmax(_LengthVec(vec1)+
                     _LengthVec(vec2), 1.0) < epsilon;
}
#endif

_RadNormalize
// normalize a radian angle into the range [-pi, +pi]
static double _RadNormalize(double rad)
{
   while (rad < -pi)
      rad += 2*pi;
   while (rad > pi)
      rad -= 2*pi;
   Assert(-pi <= rad && rad <= pi);
   return rad;
}

//
// Numerical Integration routines
//

#define cbodyMax   255

#define secHour      (3600.0)
#define secDay      (24*secHour)
#define secYear      (365.26*secDay)

#define speedLight   (2.997925e8) // speed of light in m/s

typedef struct OE OE;

_FDone
// is the rendezvous successfully completed?
static Boolean _FDone(const Body abody[], int ibodyShip, 
                                       int ibodyTarget)
{
   Vector3D vecSeparation, vecApproach;
   
   _SubVec(abody[ibodyTarget].position, abody[ibodyShip].position, 
                  &vecSeparation);
   double dist = _LengthVec(vecSeparation);
   if (dist > abody[ibodyTarget].radius + 
                        10.0*abody[ibodyShip].radius)
      return FALSE;
      
   _SubVec(abody[ibodyTarget].velocity, abody[ibodyShip].velocity, 
                  &vecApproach);
   double speedRel = _LengthVec(vecApproach);
   
   if (speedRel > 1.0 && speedRel > 
                  0.0001*_LengthVec(abody[ibodyShip].velocity))
      return FALSE;
   
   return TRUE;
}

_CalcGravAccel
// calculate the gravitation acceleration on one body from another
static void _CalcGravAccel(
   const double G,
   int cbody,
   const Body pabody[],
   int ibodyTarget,
   int ibodySource,
   Vector3D *pvecAccel
)
{
   double len;
   
   _SubVec(pabody[ibodySource].position, 
                  pabody[ibodyTarget].position, pvecAccel);
   len = _LengthVec(*pvecAccel);
   Assert(len > sEpsilon);
   _ScaleVec(G*pabody[ibodySource].mass/(len*len*len), *pvecAccel, 
                        pvecAccel);
}

_GetSystemAccel
// calculate the gravitation acceleration for every body in the system
static void _GetSystemAccel(
   const double G,
   int cbody,
   const Body pabody[],
   Vector3D pavecAccel[])
{
   int ibodyTarget, ibodySource;
      
   for (ibodyTarget = 0; ibodyTarget < cbody; ++ibodyTarget)
   {
      Vector3D *pvecAccel = &pavecAccel[ibodyTarget];
      
      _ZeroVec(pvecAccel);
      
      for (ibodySource = 0; ibodySource < cbody; ++ibodySource)
      {
         Vector3D accel;
         
         if (ibodySource == ibodyTarget)
            continue;
         
         _CalcGravAccel(G, cbody, pabody, ibodyTarget, ibodySource, 
                                       &accel);
         _AddVec(*pvecAccel, accel, pvecAccel);
      }
   }
}

_UpdateStateCore
// advance the system state
static void _UpdateStateCore(
   const SimTime dt,
   const Vector3D &vecShipAccel,
   const int cbody,
   Body pabody[],
   Vector3D pavecGrav[],
   const int fDoCollisions
)
{
   int ibody;
   
   for (ibody = 0; ibody < cbody; ++ibody)
   {
      // compute new position
      _AddScaleVec(pabody[ibody].position, dt, 
                  pabody[ibody].velocity, &pabody[ibody].position);
      _AddScaleVec(pabody[ibody].position, dt*dt, pavecGrav[ibody], 
                  &pabody[ibody].position);

      // compute new velocity
      _AddScaleVec(pabody[ibody].velocity, dt, pavecGrav[ibody], 
                  &pabody[ibody].velocity);
      
      // add in ship's acceleration
      if (pabody[ibody].bodyType == kShip)
      {
         _AddScaleVec(pabody[ibody].position, dt*dt/2.0, 
                  vecShipAccel, &pabody[ibody].position);
         _AddScaleVec(pabody[ibody].velocity, dt, vecShipAccel, 
                  &pabody[ibody].velocity);
      }
   }
   
   // check for collisions
   if (fDoCollisions)
   {
      int ibodyOuter, ibodyInner;
   
      for (ibodyOuter = 0; ibodyOuter < cbody; ++ibodyOuter)
      {
         if (pabody[ibodyOuter].bodyType == kDestroyed)
            continue;
            
         for (ibodyInner = 0; ibodyInner < ibodyOuter; ++ibodyInner)
         {
            Vector3D dpos;
            double distSquared, radius;
            
            if (pabody[ibodyInner].bodyType == kDestroyed)
               continue;
            
            _SubVec(pabody[ibodyOuter].position, 
                  pabody[ibodyInner].position, &dpos);
            distSquared = _DotVec(dpos, dpos);
            radius = pabody[ibodyOuter].radius + 
                  pabody[ibodyInner].radius;
            
            if (distSquared < radius*radius)
            {
               Vector3D I;
               int ibodyKeep, ibodyDelete;

#ifdef DEBUG
               double dist, speedRel, speedAbs;
               Vector3D vRel;
               _SubVec(pabody[ibodyOuter].velocity, 
                  pabody[ibodyInner].velocity, &vRel);   
               dist = _LengthVec(dpos);
               speedRel = _LengthVec(vRel);
               speedAbs = _LengthVec(pabody[ibodyOuter].velocity);
               Assert(dist + speedRel + speedAbs != 0.0);
#endif
               Assert(fFalse);

               ibodyKeep = pabody[ibodyOuter].radius > 
            pabody[ibodyInner].radius ? ibodyOuter : ibodyInner;
               ibodyDelete = ibodyOuter + ibodyInner - ibodyKeep;
#if 0
               DrawBodyAtPos(pabody[ibodyKeep].position, 
                  pabody[ibodyKeep].radius, g_clGreen);
               DrawBodyAtPos(pabody[ibodyDelete].position, 
                  pabody[ibodyDelete].radius, g_clRed);
#endif
#ifdef VERBOSE
            printf("objects %d and %d intersect\n", ibodyOuter, 
                  ibodyInner);
               printf("object %d absorbs object %d\n", ibodyKeep, 
                  ibodyDelete);
#endif
               
               // calculate the net inertia of the two bodies
               _ScaleVec(pabody[ibodyKeep].mass, 
                  pabody[ibodyKeep].velocity, &I);
               _AddScaleVec(I, pabody[ibodyDelete].mass, 
                  pabody[ibodyDelete].velocity, &I);

               // the new body gets the combined mass and conserves momentum
            pabody[ibodyKeep].mass += pabody[ibodyDelete].mass;
               _ScaleVec(1.0/pabody[ibodyKeep].mass, I, 
                  &pabody[ibodyKeep].velocity);
               
               // REVIEW does the radius of the new body change?
               
               // destroy the old body
               pabody[ibodyDelete].bodyType = kDestroyed;
            }
         }
      }
   }
}

//
// Orbital Mechanics
//
// The basic orbital mechanics used here are taken from
//     "Fundamentals of Astrodynamics"
//     by Roger R. Bate et al,(c) 1971
//

typedef enum
{
   okElliptical,
   okParabolic,
   okHyperbolic,
   
   okMax
} OK;

struct OE // Orbital Elements
{
   // orbit kind
   OK ok;
   
   // known position and velocity at start
   Vector3D r0;
   Vector3D v0;
   
   // orbital parameters
   double mu;                     // mass*G
   double energy;            // orbital energy
   double sh;                     // scalar angular momentum
   double tPeriod;         // period (valid only for elliptical orbits
   double a, b, c, e;   // conic section parameters
   
   // unit vectors for perifocal coordinate system
   Vector3D normalP;      // direction from focus to periapsis
   Vector3D normalQ;      // direction of velocity vector at periapsis
   Vector3D normalW;      // direction of angular momentum vector
   
   // hint for Solving for "x"
   double tPrev;
   double xPrev;
};

_SolveOrbit
// given "body" structs for an object and the parent object it orbits, 
//   compute orbital elements
static void _SolveOrbit(const Body &body, const Body &bodyParent, 
                                 int fNeedEllipse, double G, OE *poe)
{
   Vector3D h;
   
   _SubVec(body.position, bodyParent.position, &poe->r0);
   _SubVec(body.velocity, bodyParent.velocity, &poe->v0);
   
   double dist = _LengthVec(poe->r0);
   double speed = _LengthVec(poe->v0);
   
   // compute angular momentum and energy
   poe->mu = G*bodyParent.mass;
   poe->energy = speed*speed/2.0 - poe->mu/dist;
   _CrossVec(poe->r0, poe->v0, &h);
   poe->sh = _LengthVec(h);
   _ScaleVec(1.0/poe->sh, h, &poe->normalW);
   
   // periapsis
   _LinearComboVec(
      speed*speed - poe->mu/dist,
      poe->r0,
      -_DotVec(poe->r0, poe->v0),
      poe->v0,
      &poe->normalP
      );
   _NormalizeVec(poe->normalP, &poe->normalP);
   
   // the third axis, perp to periapsis and angular momentum vectors
   _CrossVec(poe->normalW, poe->normalP, &poe->normalQ);
   
   // conic geometry
   poe->e = 1.0 + 2.0*poe->energy*poe->sh*poe->sh/
                                             (poe->mu*poe->mu);
   Assert(poe->e > -sEpsilon);
   poe->e = sqrt(fmax(0.0, poe->e));
   poe->a = -poe->mu/(2.0*poe->energy);
   poe->c = poe->e*poe->a;
   
   // determine kind of orbit
   if (poe->energy < -sEpsilon)
   {
      poe->ok = okElliptical;
      poe->b = sqrt(poe->a*poe->a - poe->c*poe->c);
      poe->tPeriod = 2*pi*poe->a*sqrt(poe->a/poe->mu);
   }
   else if (poe->energy > sEpsilon)
   {
      Assert(!fNeedEllipse);
      poe->ok = okHyperbolic;
      poe->b = sqrt(poe->c*poe->c - poe->a*poe->a);
   }
   else
   {
      Assert(fFalse);
      poe->ok = okParabolic;
   }
   
   // init "x" cache
   poe->tPrev = poe->xPrev = 0.0;
}

_OrbitFromEllipse
// from initial position and ellipse geometry, compute orbital elements
static void _OrbitFromEllipse(
   double a,
   double e,
   double mu,
   const Vector3D &normalP,
   const Vector3D &normalQ,
   const Vector3D &normalW,
   const Vector3D &r0,
   OE *poe
)
{
   poe->ok = okElliptical;
   poe->r0 = r0;
   poe->normalP = normalP;
   poe->normalQ = normalQ;
   poe->normalW = normalW;
   poe->mu = mu;
   poe->a = a;
   poe->e = e;
   poe->c = a*e;
   poe->b = sqrt(poe->a*poe->a - poe->c*poe->c);
   poe->tPeriod = 2*pi*poe->a*sqrt(poe->a/poe->mu);
   
   // solve for energy and angular momentum
   poe->energy = -mu/(2.0*poe->a);
   poe->sh = sqrt((e*e - 1.0)*mu*mu/(2.0*poe->energy));

   // solve for initial velocity
   double dist = _LengthVec(r0);
   double energyKinetic = poe->energy + mu/dist;
   Assert(energyKinetic > sEpsilon);
   double speed = sqrt(energyKinetic*2.0);
   double sinTheta = poe->sh/(speed*dist);
   Assert(fabs(sinTheta) < 1.0+sEpsilon);
   double cosTheta = sqrt(fmax(1.0 - sinTheta*sinTheta, 0.0));
   if (_DotVec(poe->r0, normalQ) < 0)
      cosTheta = -cosTheta;
   Vector3D rPerp;
   _CrossVec(normalW, r0, &rPerp);
   _LinearComboVec(
         cosTheta,
         r0,
         sinTheta,
         rPerp,
         &poe->v0);
   _ScaleVec(speed/dist, poe->v0, &poe->v0);
   Assert(_LengthVec(poe->v0) < speedLight);
   
   // init "x" cache
   poe->tPrev = poe->xPrev = 0.0;
}

_ProjectOrbit
// Project an orbit to time "t" relative to starting position
//
// see F of A, p 199
static void _ProjectOrbit(OE *poe, double t, Vector3D *pr, 
                                    Vector3D *pv)
{
   // solve for "x"
   double r0 = _LengthVec(poe->r0);
   double sqrtMu = sqrt(poe->mu);
   double r_dot_v = _DotVec(poe->r0, poe->v0);
   
   Assert(t >= 0);
   double x, z, C, S;
   
   if (poe->ok == okElliptical)
   {
      if (fabs(t) > poe->tPeriod)
      {
         int cOrbit = t/poe->tPeriod;
         t -= cOrbit * poe->tPeriod;
      }
      Assert(fabs(t) <= poe->tPeriod);
   }
   
   if (fabs(t) < fabs(t - poe->tPrev) )
   {
      poe->tPrev = 0.0;
      poe->xPrev = 0.0;
   }
   x = poe->xPrev + (t - poe->tPrev)*sqrt(poe->mu)/poe->a;
   for(int cIter = 0;;++cIter)
   {
      double tn;
      
      Assert(cIter < 100);
      
      z = x*x/poe->a;
      
      if (z > sEpsilon)
      {
         double sqrtZ = sqrt(z);
         C = (1.0 - cos(sqrtZ))/z;
         S = (sqrtZ - sin(sqrtZ))/(z*sqrt(z));
      }
      else if (z < -sEpsilon)
      {
         double sqrtNegZ = sqrt(-z);
         
         C = (1.0 - cosh(sqrtNegZ))/z;
         S = (sinh(sqrtNegZ) - sqrtNegZ)/(-z*sqrtNegZ);
      }   
      else
      {
         Assert(fabs(z) <= sEpsilon);
         C = 0.5 - z*(1.0/24.0 + z*(1.0/720.0 - z/40320.0));
         S = 1.0/6.0 - z*(1.0/120.0 + z*(1.0/5040.0 - z/362880.0));
      }
      
      tn = x*(x*(r_dot_v*C/sqrtMu + x*(1.0 - r0/poe->a)*S) + r0) 
                                          / sqrtMu;
      
      if (t < 1.0 || fabs(t-tn)/t < sEpsilon)
         break;
      
      double slope = (x*(x*C + r_dot_v*(1.0 - z*S)/sqrtMu) + 
                                          r0*(1.0 - z*C))/sqrtMu;
      
      Assert(fabs(slope) > sEpsilon);
      x += (t - tn)/slope;
   }
   
   poe->tPrev = t;
   poe->xPrev = x;
   
   // compute f and g
   double f = 1.0 - x*x*C/r0;
   double g = t - x*x*x*S/sqrtMu;
   
   // compute r and its length
   _LinearComboVec(f, poe->r0, g, poe->v0, pr);
   double r = _LengthVec(*pr);
#ifdef DEBUG
   {
      double rTheory = r0 + x*(x*(C*(1.0-r0/poe->a) - 
                     r_dot_v*x*S/(poe->a*sqrtMu)) + r_dot_v/sqrtMu);
      Assert(fabs((r - rTheory)/r) < (poe->e < 0.9 ? 
                                 1.0e-4 : 1.0e-3));
   }
#endif
   
   // compute time derivatives of f and g
   double gDot = 1.0 - x*x*C/r;
   double fDot = sqrtMu*x*(z*S - 1.0)/(r0*r);
   
   // compute v
   _LinearComboVec(fDot, poe->r0, gDot, poe->v0, pv);
   
   // verify our computation...
   Assert(1.0 - (f*gDot - fDot*g) < (poe->e < 0.9 
                              ? 1.0e-5 : 1.0e-3));
}

FSolveEllipse
// given distance to periapsis or apoapsis and one other point,
// determine ellipse geometry
static int _FSolveEllipse(
   double x1, // ship position at peripasis or apoapsis, y1==0
   double x2, // target position
   double y2,
   double *pa, // ellipse parameters
   double *pb,
   double *pc,
   int *pfPeriapsis
)
{
   int fPeriapsis;
   double r1, r2;
   double cos1, cos2;
   
   Assert(x1 > 0.0);
   r1 = x1;
   r2 = sqrt(x2*x2 + y2*y2);
   fPeriapsis = r2 > r1;
   cos1 = fPeriapsis ? 1.0 : -1.0;
   cos2 = cos1*x2/r2;
   
   // Find e and p such that r = p/(1 + e*cos(nu)) for both points
   double eLow = 0.0;
   double eHigh = 0.99;
   
   double rLow = r1*(1.0 + eLow*cos1)/(1.0 + eLow*cos2);
   double rHigh = r1*(1.0 + eHigh*cos1)/(1.0 + eHigh*cos2);
   
   double errorLow = fPeriapsis ? rLow - r2 : r2 - rLow;
   double errorHigh = fPeriapsis ? rHigh - r2 : r2 - rHigh;
   
   if (errorHigh < 0.0)
      return fFalse;
   
   Assert(errorLow <= 0.0);
   Assert(errorHigh >= 0.0);
   
   for(int cIter = 0; eHigh - eLow > sEpsilon; ++cIter)
   {
      double eTry, rTry, errorTry;
      
      Assert(cIter < 20);
      Assert(eHigh > eLow);
      Assert(errorLow < 0.0);
      Assert(errorHigh > 0.0);
      
//      eTry = eLow - errorLow*(eHigh - eLow)/(errorHigh - errorLow);
      eTry = 0.5*(eLow + eHigh);
      rTry = r1*(1.0 + eTry*cos1)/(1.0 + eTry*cos2);
      errorTry = fPeriapsis ? rTry - r2 : r2 - rTry;
      
      if (errorTry > 0.0)
      {
         eHigh = eTry;
         errorHigh = errorTry;
      }
      else
      {
         eLow = eTry;
         errorLow = errorTry;
      }
   }
   
   double e = (eLow + eHigh)*0.5;
   Assert(0.01 + sEpsilon < e && e < 0.99 - sEpsilon);
   
   double p = r1*(1.0 + e*cos1);
   *pa = p/(1.0 - e*e);
   *pc = e * *pa;
   *pb = sqrt(*pa * *pa - *pc * *pc);
   *pfPeriapsis = fPeriapsis;
   
   double c = fPeriapsis ? -*pc : *pc;
   Assert(fabs((x1-c)*(x1-c)/(*pa * *pa) - 1.0) < sEpsilon);
   Assert(fabs(((x2-c)*(x2-c)/(*pa * *pa) + 
                                 y2*y2/(*pb * *pb))- 1.0) < 1.0e-4);
   
   return fTrue;
}

_IbodyFindParent
// determine the closest body the specified body is orbiting (the "parent" body)
static int _IbodyFindParent(int cbody, const Body abody[], 
         int ibodyChild)
{
   const Body *pbodyChild = &abody[ibodyChild];
   
   // REVIEW: make the caller do this and pass in G
   Vector3D avecAccel[cbodyMax];
   _GetSystemAccel(1.0, cbody, abody, avecAccel);
   
   int ibodyParent = -1;
   double invDistParent = 0.0;
   for (int ibody = 0; ibody < cbody; ++ibody)
   {      
      if (ibody == ibodyChild)
         continue;
      
      const Body *pbody = &abody[ibody];
      if (pbody->mass < pbodyChild->mass)
         continue;
         
      Vector3D vecSep;
      _SubVec(pbody->position, pbodyChild->position, &vecSep);
      double dist = _LengthVec(vecSep);
      
      Vector3D dv;
      _SubVec(pbodyChild->velocity, pbody->velocity, &dv);
      double speedSq = _DotVec(dv, dv);
      
      double energy = speedSq*0.5 - pbody->mass/dist;
      
      if (energy < 0.0 && (1.0/dist) > invDistParent)
      {
         Vector3D accel2Body, accelNet, daccel;
         
         _ScaleVec(pbody->mass/(dist*dist*dist), vecSep, 
                  &accel2Body);
   _SubVec(avecAccel[ibodyChild], avecAccel[ibody], &accelNet);
         _SubVec(accel2Body, accelNet, &daccel);
         
         if (_LengthVec(daccel) < _LengthVec(accel2Body)*0.5)
         {
            ibodyParent = ibody;
            invDistParent = 1.0/dist;
         }
      }
   }
   
   return ibodyParent;
}

_IbodyFindCommonParent
// find the common parent body of two bodies
static int _IbodyFindCommonParent(int cbody, const Body pabody[], 
                                    int ibodyChild1, int ibodyChild2)
{
   int mp_ibody_level[cbodyMax];
   int level;
   int ibody, ibodyParent;
   
   for (ibodyParent = 0; ibodyParent < cbody; ++ibodyParent)
      mp_ibody_level[ibodyParent] = 0;
      
   for (ibody = ibodyChild1, level = 0; ibody != -1; 
                  ibody = ibodyParent)
   {
      ibodyParent = _IbodyFindParent(cbody, pabody, ibody);
      if (ibodyParent == -1)
         break;
      mp_ibody_level[ibodyParent] = ++level;
   }
   
   for (ibody = ibodyChild2; ibody != -1; ibody = ibodyParent)
   {
      ibodyParent = _IbodyFindParent(cbody, pabody, ibody);
      if (ibodyParent == -1)
         break;
      if (mp_ibody_level[ibodyParent] > 0)
         return ibodyParent;
   }
   
   Assert(fFalse);
   return -1;
}

_EFromNu
// compute the eccentric anomaly from the true anomaly
static double _EFromNu(const OE &oe, double nu)
{
   double sr;
   
   double cosNu = cos(nu);
   double sinNu = sin(nu);
   sr = oe.a*(1.0 - oe.e*oe.e)/(1.0 + oe.e*cosNu);
   
   double E = acos((oe.e + cosNu)/(1.0 + oe.e*cosNu));
   if (sinNu < 0.0)
      E = -E;
#if 0 // unstable method for finding E from nu
   double EBad = atan2(sr*sinNu, sr*cosNu + oe.c);
   Assert(fabs(E1 - E) < sEpsilon);
#endif
   while (E < 0.0)
      E += 2*pi;
   Assert(0.0 <= E && E <= 2*pi);
   
   return E;
}

_NuFromE
// compute the true anomaly from the eccentric anomaly
static double _NuFromE(const OE &oe, double E)
{
   double x, y, scale, nu;
   
   x = cos(E);
   y = sin(E);
   
   scale = 1.0/sqrt(x*x/(oe.a*oe.a) + y*y/(oe.b*oe.b));
   
   nu = atan2(scale*y, scale*x - oe.c);
   
   Assert(fabs(_EFromNu(oe, nu) - E) < sEpsilon);
   
   return nu;
}

_NuFromR
// compute r vector from the true anomaly
static double _NuFromR(const OE &oe, const Vector3D &r)
{
   double x = _DotVec(r, oe.normalP);
   double y = _DotVec(r, oe.normalQ);
   double z = _DotVec(r, oe.normalW);
   
#ifdef DEBUG
   {
      double sr = _LengthVec(r);
      Assert(fabs(z)/sr < sEpsilon);
      double sUnityCheck = (x+oe.c)*(x+oe.c)/(oe.a*oe.a) + 
                                                y*y/(oe.b*oe.b);
      Assert(fabs(sUnityCheck - 1.0) < 1.0e-5);
   }
#endif

   return atan2(y, x);
}

_TimeOfFlight
// compute time of flight along an orbit between to true anomaly measurements
static double _TimeOfFlight(const OE &oe, double nuStart, 
                                                      double nuEnd)
{
   int fNegate;
   
   // compute the eccentric anomaly of the start and target positions
   double EStart = _EFromNu(oe, nuStart);
   double EEnd = _EFromNu(oe, nuEnd);
   
   // normalize
   if ((fNegate = EEnd < EStart) != 0)
   {
      double ET = EStart;
      EStart = EEnd;
      EEnd = ET;
   }
   Assert(0.0 <= EStart && EStart <= EEnd && EEnd <= 2*pi);
   
   // compute time required to go from start to end
   double tFlight = oe.a*sqrt(oe.a/oe.mu)*
      ((EEnd - oe.e*sin(EEnd)) - (EStart - oe.e*sin(EStart)));
   
   // compensate for normalization
   if (fNegate)
      tFlight = oe.tPeriod - tFlight;

#ifdef DEBUG
   if (1)
   {
      OE oeCopy = oe;
      Vector3D r, v;
      double nuCheck;
      
      nuCheck = _NuFromR(oeCopy, oeCopy.r0);
      if (fabs(_RadNormalize(nuCheck - nuStart)) < sEpsilon)
      {
         _ProjectOrbit(&oeCopy, tFlight, &r, &v);
         nuCheck = _NuFromR(oeCopy, r);
         Assert(fabs(_RadNormalize(nuCheck - nuEnd)) < 1.0e-5);
      }
   }
#endif
   
   return tFlight;
}

_SolveHohmannTransfer
//
// Given:
//   orbital parameters for the ship and target
//   true anomaly of ship at start of simulation
//   true anomaly from which ship is to begin elliptical transfer orbit
//
// Determine:
//   time ship must coast before reaching departure true anomaly
//   time needed to complete transfer orbit (half period)
//   how far away the target is at end of transfer orbit,
//     expressed as delta in true anomaly in target's orbit
//
static void _SolveHohmannTransfer(
   const Body &bodyParent,
   OE *poeShip,
   OE *poeTarget,
   double nuShipStart,
   double nuShip,
   OE *poeTransfer,
   double *pdtCoast,
   double *pdtTransfer,
   double *pdnuError
)
{
   // get starting point for elliptical transfer
   double cosShip = cos(nuShip);
   double sinShip = sin(nuShip);
   double srShip = poeShip->a*(1.0 - poeShip->e*poeShip->e)/
                                    (1.0 + poeShip->e*cosShip);
   Vector3D normalShip;
   _LinearComboVec(
         cosShip,
         poeShip->normalP,
         sinShip,
         poeShip->normalQ,
         &normalShip
         );
   Assert(_FNormalVec(normalShip));
   
   double nuTarget, srTarget;
   double sinTarget, cosTarget;
   // get ending point for elliptical transfer
   {
      Vector3D normalStartOnPlane;
      
      // project starting point vector onto target orbital plane
      double dzOffPlane = _DotVec(normalShip, poeTarget->normalW);
      double dzCorrection = _DotVec(poeShip->normalW, 
                                                         poeTarget->normalW);
      _AddScaleVec(normalShip, -dzOffPlane/dzCorrection, 
                                 poeShip->normalW, &normalStartOnPlane);
      Assert(fabs(_DotVec(normalStartOnPlane, poeTarget->normalW)) < 
                                                         sEpsilon);
      _NormalizeVec(normalStartOnPlane, &normalStartOnPlane);
      
      // point on target plane opposite to projected starting point
      cosTarget = -_DotVec(normalStartOnPlane, poeTarget->normalP);
      sinTarget = -_DotVec(normalStartOnPlane, poeTarget->normalQ);
      
      // solve for distance and true anomaly of point on target orbit opposite to 
      // projected starting point
      nuTarget = atan2(sinTarget, cosTarget);
      srTarget = poeTarget->a*(1.0 - poeTarget->e*poeTarget->e)/
                                       (1.0 + poeTarget->e*cosTarget);
   }
   // compute elements of transfer orbit
   poeTransfer->ok = okElliptical;
   poeTransfer->mu = poeShip->mu;
   poeTransfer->energy = -poeTransfer->mu/(srShip + srTarget);
   poeTransfer->normalP = normalShip;
   poeTransfer->normalW = poeShip->normalW;
   poeTransfer->a = (srShip + srTarget)*0.5;
   poeTransfer->c = poeTransfer->a - srShip;
   if (poeTransfer->c < 0.0)
   {
      // ship starts at apoapsis, rendezvous is at periapsis
      poeTransfer->c = -poeTransfer->c;
      _NegateVec(&poeTransfer->normalP);
   }
   _CrossVec(poeTransfer->normalW, poeTransfer->normalP, 
                        &poeTransfer->normalQ);
   Assert(_FNormalVec(poeTransfer->normalQ));
   poeTransfer->e = poeTransfer->c/poeTransfer->a;
   poeTransfer->b = sqrt(poeTransfer->a*poeTransfer->a - 
                        poeTransfer->c*poeTransfer->c);
   poeTransfer->tPeriod = 2.0*pi*
         sqrt(poeTransfer->a*poeTransfer->a*poeTransfer->a/
                        poeTransfer->mu);
   
   // solve for starting point of transfer orbit
   // (ship is at apoapsis or periapsis, so r0 and v0 are perpendicular)
   double energyKinetic = poeTransfer->energy + 
                        poeTransfer->mu/srShip;
   Assert(energyKinetic > sEpsilon);
   double sv0 = sqrt(energyKinetic*2.0);
   _ScaleVec(srShip, normalShip, &poeTransfer->r0);
_CrossVec(poeTransfer->normalW, normalShip, &poeTransfer->v0);
   _ScaleVec(sv0, poeTransfer->v0, &poeTransfer->v0);
   poeTransfer->sh = srShip*sv0;
   
   poeTransfer->tPrev = poeTransfer->xPrev = 0.0;
   
   // compute flight times and position of target at arrival
   *pdtCoast = _TimeOfFlight(*poeShip, nuShipStart, nuShip);
   *pdtTransfer = poeTransfer->tPeriod * 0.5;
   Vector3D rTarget, vTarget;
   _ProjectOrbit(poeTarget, *pdtCoast + *pdtTransfer, &rTarget, 
               &vTarget);
   double nuTargetArrive = atan2(_DotVec(rTarget, 
   poeTarget->normalQ), _DotVec(rTarget, poeTarget->normalP));
   
   // how much did we miss by?
   *pdnuError = _RadNormalize(nuTargetArrive - nuTarget);
}

_FFindHohmannTransfer2
// Given the ship, the target, and their common parent body, determine a valid
// Hohmann (elliptical) transfer orbit between the two, including possible change
// of plane maneuver.
//
// Returns:
//    time to coast before entering transfer orbit
//    first leg of transfer orbit, in ship's orbital plane,
//        and time to spend on first leg
//    second leg of transfer orbit, in target's orbital plane,
//        and time to complete second leg
//    total delta-v need to complete all maneuvers
static int _FFindHohmannTransfer2(
      const SimConstants &simConstants,
      const Body &bodyShip,
      const Body &bodyTarget,
      const Body &bodyParent,
      double *pdtCoast,
      double *pdtChangePlanes,
      double *pdtEndTransfer,
      OE *poeStartTransfer,
      OE *poeEndTransfer,
      double *psdv
)
{
   OE oeShip, oeTarget;
   int fHaveTransfer = fFalse;
   
   // solve for ship and target orbits
   _SolveOrbit(bodyShip, bodyParent, fTrue   /*fNeedEllipse*/, 
            simConstants.G, &oeShip);
   _SolveOrbit(bodyTarget, bodyParent, fTrue /*fNeedEllipse*/, 
            simConstants.G, &oeTarget);
   
   double nuShipStart = atan2(_DotVec(oeShip.r0, oeShip.normalQ), 
            _DotVec(oeShip.r0, oeShip.normalP));
   
   // determine length of time to examine, before orbital angles start to repeat,
   // assuming basically circular orbits
   double nuFirstShip, nuLimShip, dnuStepShip, dnuLimShip;
   double wShip = 2*pi/oeShip.tPeriod;
   double wTarget = 2*pi/oeTarget.tPeriod;
   double wMin = fmin(wShip, wTarget);
   double dw = fabs(wShip - wTarget);
   if (dw > 0.01*fmax(wShip, wTarget))
   {
      dnuLimShip = 2*pi*wShip/dw;
      dnuStepShip = dnuLimShip/20.0;
   }
   else
   {
      dnuLimShip = 4*pi*wShip/wMin;
      dnuStepShip = dnuLimShip/40.0;
   }
   nuFirstShip = _NuFromR(oeShip, oeShip.r0);
   nuLimShip = nuFirstShip + dnuLimShip*(1.0+sEpsilon);
   
   // loop over possible orbital starting positions, looking for sign change in 
   // target position error
   double dnuErrorPrev = 0.0;
   int fFirst = fTrue;
   for (double nuShip = nuFirstShip; 
               nuShip <= nuLimShip && !fHaveTransfer; 
               nuShip += dnuStepShip)
   {
      OE oeTransfer;
      double dnuError;
      double dtCoast, dtTransfer;
      
      _SolveHohmannTransfer(
         bodyParent,
         &oeShip,
         &oeTarget,
         nuShipStart,
         nuShip,
         &oeTransfer,
         &dtCoast,
         &dtTransfer,
         &dnuError
         );

      // when this fires, we have an interval which contains a valid transfer orbit
      if (!fFirst
&& (dnuErrorPrev < 0.0 ? dnuError >= 0.0 : dnuError <= 0.0)
         && fabs(dnuError - dnuErrorPrev) < pi
         )
      {
         double nu1 = nuShip - dnuStepShip;
         double nu2 = nuShip;
         double dnuError1 = dnuErrorPrev;
         double dnuError2 = dnuError;
      
         // use secant method to find correct transfer orbit
         for (int cIter = 0; ; ++cIter)
         {
            Assert(cIter < 20);
            
            double dnuErrorSecant;
            double nuShipSecant = nu1 - dnuError1 * 
                        (nu2 - nu1)/(dnuError2 - dnuError1);
            
            if (nuShipSecant < nuShip - dnuStepShip || 
                           nuShip < nuShipSecant)
            {
               Assert(fFalse);
               break;
            }
            
            _SolveHohmannTransfer(
               bodyParent,
               &oeShip,
               &oeTarget,
               nuShipStart,
               nuShipSecant,
               &oeTransfer,
               &dtCoast,
               &dtTransfer,
               &dnuErrorSecant
               );
            
            nu1 = nu2;
            dnuError1 = dnuError2;
            nu2 = nuShipSecant;
            dnuError2 = dnuErrorSecant;
            
            if (fabs(dnuError2 - dnuError1) <= sEpsilon)
            {
               fHaveTransfer = fTrue;
               break;
            }
         }
         
         if (!fHaveTransfer)
            continue;
         
         // we have the transfer orbit, rotated into the plane of the ship's current orbit
         *poeStartTransfer = oeTransfer;
         *pdtCoast = dtCoast;
         
         // compute correct transfer orbit in plane of destination
         {
            Vector3D rTargetEnd, vTargetEnd;
            Vector3D normalP, normalQ, normalW;
            double a, b, c, e;
            int fPeriapsis;
            
            // compute rendezvous position
      _ProjectOrbit(&oeTarget, dtCoast+dtTransfer, &rTargetEnd, 
                        &vTargetEnd);
            
            // find the ellipse that fits the two points 
            // (with start and rendezvous at peri/apoapsis)
            _FSolveEllipse(
                  _LengthVec(poeStartTransfer->r0),
                  -_LengthVec(rTargetEnd),
                  0.0,
                  &a,
                  &b,
                  &c,
                  &fPeriapsis
                  );
            e = c/a;
   Assert(fabs((a - poeStartTransfer->a)/poeStartTransfer->a) 
                                                < 2.0e-5);
   Assert(fabs((b - poeStartTransfer->b)/poeStartTransfer->b) 
                                                < 2.0e-5);
   Assert(fabs((c - poeStartTransfer->c)/poeStartTransfer->c) 
                                                < 2.0e-5);
            Assert(fabs(e - poeStartTransfer->e) < 1.0e-6);
            
            _NormalizeVec(rTargetEnd, &normalP);
            if (fPeriapsis)
               _NegateVec(&normalP);
      _CrossVec(poeStartTransfer->normalW, normalP, &normalQ);
            _NormalizeVec(normalQ, &normalQ);
            _CrossVec(normalP, normalQ, &normalW);
            Assert(_FNormalVec(normalW));
            
         // compute time to plane change, and time from plane change to rendezvous
         int fStartAtPeriapsis = _DotVec(poeStartTransfer->r0, 
                                    poeStartTransfer->normalP) > 0.0;
            double nuStart, nuPlaneChange;
            nuStart = fStartAtPeriapsis ? 0.0 : pi;
            nuPlaneChange = nuStart + pi*0.5;
            *pdtChangePlanes = _TimeOfFlight(*poeStartTransfer, 
                                             nuStart, nuPlaneChange);
            *pdtEndTransfer = dtTransfer - *pdtChangePlanes;
            
            // compute rendezvous position
            Vector3D rPlaneChange;
            _ScaleVec(
                  poeStartTransfer->b * sqrt(1.0 - 
                        poeStartTransfer->c*poeStartTransfer->c/
                        (poeStartTransfer->a*poeStartTransfer->a)),
                  poeStartTransfer->normalQ,
                  &rPlaneChange
                  );
            if (!fStartAtPeriapsis)
               _NegateVec(&rPlaneChange);
   double nuCheck = _NuFromR(*poeStartTransfer, rPlaneChange);
         Assert(fabs(_RadNormalize(nuCheck - nuPlaneChange)) < 
                                             sEpsilon);

            // compute orbital elements for second leg of transfer
         _OrbitFromEllipse(a, e, poeStartTransfer->mu, normalP, 
                  normalQ, normalW, rPlaneChange, poeEndTransfer);
            Assert(fabs((poeEndTransfer->a - 
            poeStartTransfer->a)/poeStartTransfer->a) < 2.0e-5);
            Assert(fabs((poeEndTransfer->b - 
            poeStartTransfer->b)/poeStartTransfer->b) < 2.0e-5);
         }
                  
         // what's the total delta-v for the transfer?
         Vector3D r, vShipStart, vShipEnd, vTargetEnd, dv;
         
         *psdv = 0.0;
         
         _ProjectOrbit(&oeShip, dtCoast, &r, &vShipStart);
         _SubVec(vShipStart, oeTransfer.v0, &dv);
         *psdv += _LengthVec(dv);
         
         _ProjectOrbit(&oeTransfer, dtTransfer, &r, &vShipEnd);
         _ProjectOrbit(&oeTarget, dtTransfer+dtCoast, &r, 
                        &vTargetEnd);
         _SubVec(vShipEnd, vTargetEnd, &dv);
         *psdv += _LengthVec(dv);
      }
      
      fFirst = fFalse;
      dnuErrorPrev = dnuError;
   }
   
   return fHaveTransfer;
}

_OptimizeDock
// for a simple docking maneuver between nearly equal orbits, determine timing to 
// minimize delta-v
static void _OptimizeDock(OE *poeShip, OE *poeTarget, 
                           const Body &bodyParent, double *ptWait)
{
   double tPeriod = fmax(poeShip->tPeriod, poeTarget->tPeriod);
   double dt = fmin(poeShip->tPeriod, poeTarget->tPeriod)/20.0;
   
   double tBest = 0.0;
   double sdvMin = 10.0*speedLight;
   
   double tLim = tPeriod;
   for (double t = 0.0; t <= tLim; t += dt)
   {
      Vector3D rShip, vShip;
      Vector3D rTarget, vTarget;
      Vector3D dv;
      double sdv;
      
      _ProjectOrbit(poeShip, t, &rShip, &vShip);
      _ProjectOrbit(poeTarget, t, &rTarget, &vTarget);
      
      _SubVec(vTarget, vShip, &dv);
      sdv = _LengthVec(dv);

      if (sdv < sdvMin)
      {
         tLim = t + tPeriod;
         sdvMin = sdv;
         tBest = t;
      }
   }
   
   *ptWait = tBest;
}

// Rendezvous Operation Kind
//
// opcodes for simple orbital maneuvers
typedef enum
{
   rokPlan,                     // go back into plan mode to compute next action
   rokFindEscape,            // pick a good time to blast out of orbit
   rokEscapeOrbit,            // escape orbit from current parent object
   rokAssumeStandardOrbitMrSulu,   // duh
   rokCoast,                  // let Newton have his way
   rokBeginTransfer,         // jump into a transfer orbit
   rokMonitorTransfer,      // monitor transfer orbit, adjust as needed
   rokDock,                     // perform final docking maneuver
   
   rokMax
} ROK;

typedef struct RO RO;
struct RO // Rendezvous Operation
{
   ROK      rok;
   double   dt;
   
   // valid on a per-rok basis
   OE      oe;
   int      ibodyTarget;
   int      ibodyParent;
   int      ibodyGrandParent;
   int      fWantInner;
};

static int s_cbody;
static int s_ibodyShip;
static int s_ibodyTargetFinal;
static SimTime s_tCur;
static SimTime s_tStartOp;
static SimTime s_tNextOp;

static RO s_aro[cbodyMax+1];
static int s_iroCur, s_cro;
static double s_sdv;
static double s_speedMax;

static SimConstants s_simconst;

static void _Plan(const Body abody[]);

InitSim
// callback to begin simulation
pascal void InitSim(
  const SimConstants simConstants, /* simulation constants */
  const SInt8 cbody,  /* number of bodies in the simulation */
  const Body abody[]  /* parameters of simulation bodies */
)
{
   s_simconst = simConstants;
   s_cbody = cbody;
   
   Assert(s_cbody <= cbodyMax);
   
#ifdef DEBUG
   s_ibodyShip = s_ibodyTargetFinal = -1;
#endif
   for (int ibody = 0; ibody < s_cbody; ++ibody)
   {
      switch (abody[ibody].bodyType)
      {
      case kShip:
         Assert(s_ibodyShip == -1);
         s_ibodyShip = ibody;
         break;
      case kTarget:
         Assert(s_ibodyTargetFinal == -1);
         s_ibodyTargetFinal = ibody;
         break;
      }
   }
   Assert(s_ibodyShip != -1);
   Assert(s_ibodyTargetFinal != -1);

   _Plan(abody);
}

_Plan
// plan next maneuver
static void _Plan(const Body abody[])
{
   int ibodyTargetCur, ibodyShipParent, ibodyTargetParent, ibodyCommonParent;
   
   s_speedMax = 1000.0; // maximum closing speed
   
   s_cro = 0; // no ops in the list yet
   
   ibodyTargetCur = s_ibodyTargetFinal; // until proven otherwise
   ibodyShipParent = _IbodyFindParent(s_cbody, abody, 
                              s_ibodyShip);
   ibodyTargetParent = _IbodyFindParent(s_cbody, abody, 
                              ibodyTargetCur);
   ibodyCommonParent = _IbodyFindCommonParent(s_cbody, abody, 
                              s_ibodyShip, ibodyTargetCur);
   
   if (ibodyShipParent == -1 || ibodyTargetParent == -1 || 
                              ibodyCommonParent == -1)
   {
      // situation defies analysis, just drive there
      Assert(fFalse);
      goto LDock;
   }
   else if (ibodyShipParent != ibodyCommonParent)
   {
      int ibodyShipGrandparent = _IbodyFindParent(s_cbody, abody, 
                              ibodyShipParent);
      Assert(ibodyShipGrandparent != -1);
      
      double srShip =                               
            _DistBetweenPos(abody[ibodyTargetCur].position, 
                           abody[ibodyCommonParent].position);
      double srTarget = 
            _DistBetweenPos(abody[ibodyTargetCur].position, 
                           abody[ibodyCommonParent].position);
      
      // escape current orbit to orbit parent of current parent
      s_aro[s_cro].rok = rokFindEscape;
      s_aro[s_cro].ibodyParent = ibodyShipParent;
      s_aro[s_cro].ibodyGrandParent = ibodyShipGrandparent;
      s_aro[s_cro].fWantInner = srShip > srTarget;
      s_aro[s_cro++].dt = -1.0;
      
      s_aro[s_cro].rok = rokEscapeOrbit;
      s_aro[s_cro].ibodyParent = ibodyShipParent;
      s_aro[s_cro].ibodyGrandParent = ibodyShipGrandparent;
      s_aro[s_cro++].dt = -1.0;
      
      s_aro[s_cro].rok = rokPlan;
      s_aro[s_cro++].dt = -1.0;
   }
   else
   {
      Assert(ibodyShipParent == ibodyCommonParent);
      
      // if the target is orbiting something orbiting the common parent,
      // then, for now, make our target the body orbiting the common parent
      if (ibodyTargetParent != ibodyCommonParent)
      {
         do
         {
            ibodyTargetCur = ibodyTargetParent;
            ibodyTargetParent = _IbodyFindParent(s_cbody, abody, 
                     ibodyTargetCur);
            Assert(ibodyTargetParent != -1);
         } while (ibodyTargetParent != ibodyCommonParent);
      }

      OE oeShip;
      OE oeTarget;
      double sdvBest, dtCoastBest, dtChangePlanes, dtDockBest;
      OE oeStartTransferBest, oeEndTransferBest;
      int fHaveTransfer;
      
      _SolveOrbit(abody[s_ibodyShip], abody[ibodyShipParent], 
            fFalse/*fNeedEllipse*/, s_simconst.G, &oeShip);
   _SolveOrbit(abody[ibodyTargetCur], abody[ibodyTargetParent], 
            fFalse/*fNeedEllipse*/, s_simconst.G, &oeTarget);

#ifdef DEBUG_CHART
      DrawOrbit(abody[ibodyShipParent], &oeShip);
      DrawOrbit(abody[ibodyTargetParent], &oeTarget);
#endif

      if (oeShip.ok != okElliptical)
      {
         // ship has the right parent, but has too much relative velocity.
         // get it into an elliptical orbit, the call Plan again
         s_aro[s_cro].rok = rokAssumeStandardOrbitMrSulu;
         s_aro[s_cro].ibodyParent = ibodyShipParent;
         s_aro[s_cro].oe = oeTarget;
         s_aro[s_cro++].dt = -1.0;
         
         s_aro[s_cro].rok = rokPlan;
         s_aro[s_cro++].dt = -1.0;
         
         goto LRet;
      }
      
      if (oeTarget.ok != okElliptical)
      {
         Vector3D dv;
         
         // target is in non-elliptical orbit. As a terrible solution,
         // just drive toward the target
         ibodyTargetCur = s_ibodyTargetFinal;
         _SubVec(abody[s_ibodyShip].velocity, 
                        abody[ibodyTargetCur].velocity, &dv);
         s_speedMax = fmax(s_speedMax, _LengthVec(dv));
         goto LDock;
      }
      
      fHaveTransfer = fFalse;
      if (ibodyShipParent == ibodyTargetParent)
      {         
         fHaveTransfer = _FFindHohmannTransfer2(
               s_simconst,
               abody[s_ibodyShip],
               abody[ibodyTargetCur],
               abody[ibodyTargetParent],
               &dtCoastBest,
               &dtChangePlanes,
               &dtDockBest,
               &oeStartTransferBest,
               &oeEndTransferBest,
               &sdvBest
               );
      }
#if 0
      if (!fHaveTransfer)
      {
         double tCoast, dtCoast;
         dtCoast = oeShip.tPeriod/20.0;
for (tCoast = 0.0; tCoast < oeShip.tPeriod; tCoast += dtCoast)
         {
            Vector3D rShip, vShip, posShip, rTarget, vTarget;
            double tDock, sdv;
            OE oeTransfer;
            Body bodyTarget;
            
            bodyTarget = abody[ibodyTargetCur];
            
            _ProjectOrbit(&oeShip, tCoast, &rShip, &vShip);
      _ProjectOrbit(&oeTarget, tCoast, &rTarget, &vTarget);
            
            // REVIEW - project correct positions for parents?
   _AddVec(abody[ibodyShipParent].position, rShip, &posShip);
      _AddVec(abody[ibodyShipParent].velocity, vShip, &vShip);
            _AddVec(abody[ibodyTargetParent].position, rTarget, 
                  &bodyTarget.position);
            _AddVec(abody[ibodyTargetParent].velocity, vTarget, 
                  &bodyTarget.velocity);

   #ifdef DEBUG_CHART
            DotAtPos(posShip, 2, *g_pclShip);
            DotAtPos(bodyTarget.position, 2, *g_pclTarget);
   #endif
            
            if (_FFindHohmannTransfer(
               s_simconst,
               bodyTarget,
               abody[ibodyTargetParent],
               posShip,
               vShip,
               &oeTransfer,
               &tDock,
               &sdv))
            {
               if (!fHaveTransfer || sdv < sdvBest)
               {
                  fHaveTransfer = fTrue;
                  sdvBest = sdv;
                  dtCoastBest = tCoast;
                  dtDockBest = tDock;
                  oeTransferBest = oeTransfer;
               }
            }
         }
      }
#endif
      
      if (fHaveTransfer)
      {
         s_aro[s_cro].rok = rokCoast;
         s_aro[s_cro++].dt = dtCoastBest;
         
         s_aro[s_cro].rok = rokBeginTransfer;
         s_aro[s_cro].oe = oeStartTransferBest;
         s_aro[s_cro].ibodyTarget = ibodyTargetCur;
         s_aro[s_cro].ibodyParent = ibodyShipParent;
         s_aro[s_cro++].dt = dtChangePlanes;
         
         s_aro[s_cro].rok = rokBeginTransfer;
         s_aro[s_cro].oe = oeEndTransferBest;
         s_aro[s_cro].ibodyTarget = ibodyTargetCur;
         s_aro[s_cro].ibodyParent = ibodyShipParent;
         s_aro[s_cro++].dt = dtDockBest;
         
         if (ibodyTargetCur != s_ibodyTargetFinal)
         {
            s_aro[s_cro].rok = rokPlan;
            s_aro[s_cro++].dt = -1.0;
            goto LRet;
         }
      }
      else if (ibodyTargetParent == ibodyShipParent && 
                        oeTarget.ok == okElliptical)
      {
   _OptimizeDock(&oeShip, &oeTarget, abody[ibodyTargetParent], 
                        &dtCoastBest);
         
         s_aro[s_cro].rok = rokCoast;
         s_aro[s_cro++].dt = dtCoastBest;
      }
      else
      {
         // bad news, nothing smart to do. Try just driving there...
         Assert(fFalse);
      }

LDock:
      // finally, dock
      s_aro[s_cro].rok = rokDock;
      s_aro[s_cro++].dt = -1.0;
   }

LRet:
   Assert(s_cro > 0 && 
         (s_aro[s_cro-1].rok == rokDock || 
          s_aro[s_cro-1].rok == rokPlan));
   s_tNextOp = 0.0;
   s_iroCur = -1;

   s_tCur = 0.0;
}

_SetAccel
// set ship's acceleration for the current time step such that ship matches
// velocity with target object, or other position/velocity. If ship is more
// than "distDecel" away from target position, try to set relative velocity
// to speedMax in direction of target. As distance drops below distDecel,
// reduce closing speed linearly.
static void _SetAccel(
   int cbody,
   const Body abody[],
   const Body abodyNew[],
   const Vector3D avecGrav[],
   int ibodyShip,
   double G,
   const Vector3D &posWant,
   const Vector3D &vWant,
   double radiusTarget,
   double dt,
   double accelMax,
   double speedMax,
   double distDecel,
   int fDocking,
   Vector3D *pvecAccel
   )
{
   Vector3D vecSepCur;
   Vector3D vecSepNext;
   Vector3D vecRelVelocity;

   if (fDocking)
   {
      // projected separation and velocity
      _SubVec(abody[s_ibodyTargetFinal].position, 
               abody[ibodyShip].position, &vecSepCur);
      _SubVec(abodyNew[s_ibodyTargetFinal].position, 
               abodyNew[ibodyShip].position, &vecSepNext);
      _SubVec(abodyNew[ibodyShip].velocity, 
      abodyNew[s_ibodyTargetFinal].velocity, &vecRelVelocity);
   }
   else
   {
      // projected separation and velocity
      _SubVec(posWant, abody[ibodyShip].position, &vecSepCur);
_SubVec(posWant, abodyNew[ibodyShip].position, &vecSepNext);
_SubVec(abodyNew[ibodyShip].velocity, vWant, &vecRelVelocity);
   }
   
   // what would we like our relative velocity to be?
   Vector3D vecRelVelocityWant;
   double distSep = _LengthVec(vecSepNext);
   Assert(distSep > sEpsilon);
   
   double dist = distSep - radiusTarget - 
                        5*abodyNew[ibodyShip].radius;
   double speedWant = speedMax;
   double speedMin = 0.5*fmin(1.0, _LengthVec(vWant)*0.0001);
   if (_DotVec(vecSepNext, vecSepCur) < 0.0)
   {
      speedWant = 0.0;
   }
   else if (dist < distDecel)
   {
      speedWant *= dist/distDecel;
      if (speedWant < speedMin)
         speedWant = speedMin;
   }
   _ScaleVec(speedWant/distSep, vecSepNext, &vecRelVelocityWant);
   
   // what delta-v do we need to get the right velocity?
   Vector3D dv;
   _SubVec(vecRelVelocityWant, vecRelVelocity, &dv);
   double sdv = _LengthVec(dv);
   
   // apply acceleration, subject to limit
   double accel = sdv/dt;
   if (accel > accelMax)
      accel = accelMax;
   if (accel < 0.000001 * accelMax)
      *pvecAccel = s_vecZero;
   else
      _ScaleVec(accel/sdv, dv, pvecAccel);

}

AdvanceTime
// callback to advance time
pascal Boolean AdvanceTime(  /* return TRUE if docking complete or giving up */
  const SimTime dt,    /* propagate forward this amount of time */
  Vector3D *pvecShipAccel,
    /* return ship acceleration vector applied this step */
  Body abody[]             
    /* return body parameters at the end of this time step */
)
{
   static Vector3D avecGrav[cbodyMax];
   static Body abodyNew[cbodyMax];

   Vector3D posWant, vWant;
   double radiusTarget = 1.0;
   
   // project system forward
   memcpy(abodyNew, abody, s_cbody*sizeof(abody[0]));
   _GetSystemAccel(s_simconst.G, s_cbody, abody, avecGrav);
   _UpdateStateCore(dt, s_vecZero, s_cbody, abodyNew, avecGrav, 
            fFalse/*fDoCollisions*/);
   
   // default action is no thrust   
   _ZeroVec(pvecShipAccel);

   if (s_tNextOp >= 0.0 && s_tCur + dt*0.5 >= s_tNextOp)
   {
LAdvanceOp:
      ++s_iroCur;
      Assert(s_iroCur < s_cro);
      s_tStartOp = s_tCur;
      if (s_aro[s_iroCur].dt > 0.0)
         s_tNextOp = s_tCur + s_aro[s_iroCur].dt;
      else
         s_tNextOp = -1.0;
   }
   switch (s_aro[s_iroCur].rok)
   {
   default:
      Assert(fFalse);
   case rokCoast:
      goto LUpdate;
      
   case rokPlan:
      _Plan(abody);
      goto LAdvanceOp;
   
   case rokAssumeStandardOrbitMrSulu:
      {
         OE oeCur;
         int ibodyParent = s_aro[s_iroCur].ibodyParent;
         
         _SolveOrbit(abody[s_ibodyShip], abody[ibodyParent], 
               fFalse/*fNeedEllipse*/, s_simconst.G, &oeCur);
         if (oeCur.ok == okElliptical && oeCur.e < 0.5)
            goto LAdvanceOp;
         
         Vector3D vecSep, vRel, vRelWant;
         _SubVec(abody[s_ibodyShip].position, 
                  abody[ibodyParent].position, &vecSep);
         _SubVec(abody[s_ibodyShip].velocity, 
                  abody[ibodyParent].velocity, &vRel);
         double dist = _LengthVec(vecSep);
         
         Vector3D dv;
#if 0
         if (oeCur.energy > s_aro[s_iroCur].oe.energy*0.5)
         {
            // slow down!
            dv = vRel;
            _NegateVec(&dv);
         }
         else
#endif
         {
            // increase angular momentum!
      _CrossVec(s_aro[s_iroCur].oe.normalW, vecSep, &vRelWant);
            _NormalizeVec(vRelWant, &vRelWant);
            _ScaleVec(sqrt(oeCur.mu/dist), vRelWant, &vRelWant);
            
            _SubVec(vRelWant, vRel, &dv);
         }
         
         double sdv = _LengthVec(dv);
         double accel = sdv/dt;
         if (accel > s_simconst.maxAcceleration)
            accel = s_simconst.maxAcceleration;
         _ScaleVec(accel/sdv, dv, pvecShipAccel);
      }
      goto LUpdate;
   
   case rokFindEscape:
      {
         int ibodyParent = s_aro[s_iroCur].ibodyParent;
      int ibodyGrandParent = s_aro[s_iroCur].ibodyGrandParent;
         
         Vector3D vecGrandToParent, vRel, vRelNew, vRelParent;
         _SubVec(abody[ibodyParent].position, 
         abody[ibodyGrandParent].position, &vecGrandToParent);
         _SubVec(abody[s_ibodyShip].velocity, 
                  abody[ibodyParent].velocity, &vRel);
         _SubVec(abodyNew[s_ibodyShip].velocity, 
                  abodyNew[ibodyParent].velocity, &vRelNew);
         _SubVec(abody[ibodyParent].velocity, 
                  abody[ibodyGrandParent].velocity, &vRelParent);
         
         if (_DotVec(vRel, vecGrandToParent) * 
                  _DotVec(vRelNew, vecGrandToParent) <= 0.0)
         {
            // begin maneuver
            double dot = _DotVec(vRel, vRelParent);
      if (s_aro[s_iroCur].fWantInner ? dot < 0.0 : dot > 0.0)
               goto LAdvanceOp;
         }
      }
      goto LUpdate;
      
   case rokEscapeOrbit:
      {
         int ibodyParent = s_aro[s_iroCur].ibodyParent;
      int ibodyGrandParent = s_aro[s_iroCur].ibodyGrandParent;
         
         if (_IbodyFindParent(s_cbody, abody, s_ibodyShip) != 
                                       ibodyParent)
         {
            // we have escaped!
            goto LAdvanceOp;
         }
         
         // compute gravitational acceleration on ship from parent
         Vector3D accelGrav1, accelGrav2;
      _CalcGravAccel(s_simconst.G, s_cbody, abody, s_ibodyShip, 
                                             ibodyParent, &accelGrav1);
   _CalcGravAccel(s_simconst.G, s_cbody, abodyNew, s_ibodyShip, 
                                             ibodyParent, &accelGrav2);
         
         // set ship's accel to negate average gravity
         _AddVec(accelGrav1, accelGrav2, pvecShipAccel);
         _ScaleVec(-0.5, *pvecShipAccel, pvecShipAccel);
         
         // limit to maximum acceleration
         double accel = _LengthVec(*pvecShipAccel);
         if (accel > s_simconst.maxAcceleration)
   _ScaleVec(s_simconst.maxAcceleration/accel, *pvecShipAccel, 
                                             pvecShipAccel);
      }
      goto LUpdate;
      
   case rokMonitorTransfer:
   case rokBeginTransfer:
      {
         Vector3D rWant, vecSep, vRel;
         int ibodyTarget = s_aro[s_iroCur].ibodyTarget;
      _ProjectOrbit(&s_aro[s_iroCur].oe, s_tCur+dt-s_tStartOp, 
                                 &rWant, &vWant);
         _AddVec(abodyNew[s_aro[s_iroCur].ibodyParent].position, 
                        rWant, &posWant);
         _AddVec(abodyNew[s_aro[s_iroCur].ibodyParent].velocity, 
                        vWant, &vWant);
         
         _SubVec(abodyNew[s_ibodyShip].velocity, vWant, &vRel);
   _SubVec(posWant, abodyNew[s_ibodyShip].position, &vecSep);
         
         double errorTolerance = s_aro[s_iroCur].rok == 
                        rokMonitorTransfer ? 1.0e-3 : 2.0e-5;
   if (_LengthVec(vecSep)/_LengthVec(posWant) < errorTolerance
      && _LengthVec(vRel)/_LengthVec(vWant) < errorTolerance)
         {
            s_aro[s_iroCur].rok = rokMonitorTransfer;
            goto LUpdate;
         }
         
         // set max closing speed to current relative speed
         _SubVec(abodyNew[s_ibodyShip].velocity, 
                        abodyNew[ibodyTarget].velocity, &vRel);
         s_speedMax = _LengthVec(vRel);
         
         // check to see if we need to early exit from the transfer orbit
         if (s_aro[s_iroCur+1].rok == rokDock)
         {
            Assert(s_iroCur + 1 < s_cro);
            int ibodyParent = _IbodyFindParent(s_cbody, abody, 
                        s_ibodyShip);
            if (ibodyParent != s_aro[s_iroCur].ibodyParent)
            {
               _SubVec(abodyNew[ibodyTarget].position, 
                        abodyNew[s_ibodyShip].position, &vecSep);
      double distCollide = 10.0*abodyNew[s_ibodyShip].radius + 
                        abodyNew[ibodyTarget].radius;
               double dist = _LengthVec(vecSep);
               double speedClose = _DotVec(vRel, vecSep) / dist;
               double dtClose = fmax(0.0, dist - distCollide) / 
                                                            speedClose;
         double dtStop = speedClose/s_simconst.maxAcceleration;
   if (dtClose < 0.0 || dtClose > fmax(10 * dt, 2 * dtStop))
                  goto LUpdate;
               // we've changed parent, and are moving toward the target,
               // go into docking mode
               goto LAdvanceOp;
            }
         }
      }
      goto LMatch;

   case rokDock:
LDock:
      posWant = abodyNew[s_ibodyTargetFinal].position;
      vWant = abodyNew[s_ibodyTargetFinal].velocity;
      radiusTarget = abodyNew[s_ibodyTargetFinal].radius;
      
LMatch:
      double distDecel = 2.0*fmax(
         0.5*(s_speedMax*s_speedMax)/s_simconst.maxAcceleration,
            2.0*dt*s_speedMax);
            
      // set ship's accel to give right velocity at end of time frame
      _SetAccel(
            s_cbody,
            abody,
            abodyNew,
            avecGrav,
            s_ibodyShip,
            s_simconst.G,
            posWant,
            vWant,
            radiusTarget,
            dt,
            s_simconst.maxAcceleration,
            s_speedMax,
            distDecel,
            s_aro[s_iroCur].rok == rokDock,
            pvecShipAccel);

LUpdate:
      // add in ship's accel, and recompute future state
      _UpdateStateCore(dt, *pvecShipAccel, s_cbody,
                                                abody, avecGrav, fTrue);
      break;
   }
   
   s_tCur += dt;
   return _FDone(abody, s_ibodyShip, s_ibodyTargetFinal) || 
                                 abody[s_ibodyShip].bodyType != kShip;
}
 
AAPL
$118.93
Apple Inc.
-0.07
MSFT
$48.12
Microsoft Corpora
+0.37
GOOG
$541.03
Google Inc.
+0.66

MacTech Search:
Community Search:

Software Updates via MacUpdate

Skype 7.2.0.412 - Voice-over-internet ph...
Skype allows you to talk to friends, family and co-workers across the Internet without the inconvenience of long distance telephone charges. Using peer-to-peer data transmission technology, Skype... Read more
HoudahSpot 3.9.6 - Advanced file search...
HoudahSpot is a powerful file search tool built upon MacOS X Spotlight. Spotlight unleashed Create detailed queries to locate the exact file you need Narrow down searches. Zero in on files Save... Read more
RapidWeaver 6.0.3 - Create template-base...
RapidWeaver is a next-generation Web design application to help you easily create professional-looking Web sites in minutes. No knowledge of complex code is required, RapidWeaver will take care of... Read more
iPhoto Library Manager 4.1.10 - Manage m...
iPhoto Library Manager lets you organize your photos into multiple iPhoto libraries. Separate your high school and college photos from your latest summer vacation pictures. Or keep some photo... Read more
iExplorer 3.5.1.9 - View and transfer al...
iExplorer is an iPhone browser for Mac lets you view the files on your iOS device. By using a drag and drop interface, you can quickly copy files and folders between your Mac and your iPhone or... Read more
MacUpdate Desktop 6.0.3 - Discover and i...
MacUpdate Desktop 6 brings seamless 1-click installs and version updates to your Mac. With a free MacUpdate account and MacUpdate Desktop 6, Mac users can now install almost any Mac app on macupdate.... Read more
SteerMouse 4.2.2 - Powerful third-party...
SteerMouse is an advanced driver for USB and Bluetooth mice. It also supports Apple Mighty Mouse very well. SteerMouse can assign various functions to buttons that Apple's software does not allow,... Read more
iMazing 1.1 - Complete iOS device manage...
iMazing (was DiskAid) is the ultimate iOS device manager with capabilities far beyond what iTunes offers. With iMazing and your iOS device (iPhone, iPad, or iPod), you can: Copy music to and from... Read more
PopChar X 7.0 - Floating window shows av...
PopChar X helps you get the most out of your font collection. With its crystal-clear interface, PopChar X provides a frustration-free way to access any font's special characters. Expanded... Read more
OneNote 15.4 - Free digital notebook fro...
OneNote is your very own digital notebook. With OneNote, you can capture that flash of genius, that moment of inspiration, or that list of errands that's too important to forget. Whether you're at... Read more

Latest Forum Discussions

See All

7 tips to get ahead of the competition i...
7 tips to get ahead of the competition in Dynasty of Dungeons Posted by Simon Reed on November 28th, 2014 [ permalink ] Playcrab has launched their action-packed new dungeon crawler, Dynasty of Dungeons, today. | Read more »
Monster Strike Review
Monster Strike Review By Campbell Bird on November 28th, 2014 Our Rating: :: BILLIARD STRATEGYUniversal App - Designed for iPhone and iPad Collect monsters and battle by flinging them across the battlefield in this strangely... | Read more »
Proun+ Review
Proun+ Review By Jennifer Allen on November 28th, 2014 Our Rating: :: TWITCHY RACINGUniversal App - Designed for iPhone and iPad Twitchy racing aplenty in Proun+, an enjoyably tricky title.   | Read more »
Lucha Amigos (Games)
Lucha Amigos 1.0 Device: iOS Universal Category: Games Price: $1.99, Version: 1.0 (iTunes) Description: Forget Ninja Turtles, and meet Wrestlers Turtles! Crazier, Spicier and…Bouncier! Sling carapaces of 7 Luchadores to knock all... | Read more »
Record of Agarest War Zero (Games)
Record of Agarest War Zero 1.0 Device: iOS Universal Category: Games Price: $7.99, Version: 1.0 (iTunes) Description: HyperDevbox Holiday Turkey Black Friday Special Pricing! To celebrate the opening of the holiday season HyperDevbox... | Read more »
Raby (Games)
Raby 1.0.3 Device: iOS Universal Category: Games Price: $2.99, Version: 1.0.3 (iTunes) Description: ***WARNING - Raby runs on: iPhone 5, iPhone 5C, iPhone 5S, iPhone 6, iPhone 6 Plus, iPad Mini Retina, iPad Mini 3, iPad 4, iPad Air,... | Read more »
Oddworld: Stranger's Wrath (Games)
Oddworld: Stranger's Wrath 1.0 Device: iOS Universal Category: Games Price: $5.99, Version: 1.0 (iTunes) Description: ** PLEASE NOTE: Oddworld Stranger's Wrath requires at least an iPhone 4S, iPad 2, iPad Mini or iPod Touch 5th gen... | Read more »
Bounce On Back (Games)
Bounce On Back 1.0 Device: iOS Universal Category: Games Price: $2.99, Version: 1.0 (iTunes) Description: | Read more »
Dwelp (Games)
Dwelp 1.0 Device: iOS Universal Category: Games Price: $.99, Version: 1.0 (iTunes) Description: === 50% off for a limited time, to celebrate release === Dwelp is an elegant little puzzler with a brand new game mechanic. To complete a... | Read more »
Make Way for Fat Chicken, from the Maker...
Make Way for Fat Chicken, from the Makers of Scrap Squad Posted by Jessica Fisher on November 26th, 2014 [ permalink ] Relevant Games has announced they will be releasing their reverse tower defense game, | Read more »

Price Scanner via MacPrices.net

Black Friday: Up to $60 off Mac minis, NY tax...
 B&H Photo has new 2014 Mac minis on sale for up to $60 off MSRP as part of their Black Friday sale. Shipping is free, and B&H charges NY sales tax only: - 1.4GHz Mac mini: $449.99 $50 off... Read more
Black Friday: 27-inch 5K iMac for $2299, save...
 B&H Photo continues to offer Black Friday sale prices on the 27″ 3.5GHz 5K iMac, in stock today and on sale for $2299 including free shipping plus NY sales tax only. Their price is $200 off MSRP... Read more
Karalux Announces 24K Gold-Plated iPhone 6
Karalux, a Vietnam-based jewellery firm, has launched a unique 24 karat gold-plated iPhone 6 version with gold-cast monolithic dragon on its back panel. The real 24 karat gold plated enclosure doesn’... Read more
Black Friday: 13-inch 2.6GHz Retina MacBook P...
 B&H Photo has lowered their price for the 13″ 2.6GHz/128GB Retina MacBook Pro to $1159 for Black Friday. That’s $140 off MSRP, and it’s the lowest price for this model (except for Apple’s $1099... Read more
View all the Black Friday sales on our Mac Pr...
We’ve updated our Mac Price Trackers with the latest information on prices, bundles, and availability on systems from Apple’s authorized internet/catalog resellers. View Black Friday sale prices at a... Read more
Black Friday: 11-inch MacBook Air for $779, s...
 Best Buy has lowered their price for the 2014 11″ 1.4GHz/128GB MacBook Air to $779.99 for Black Friday. That’s $120 off MSRP. Choose free shipping or free local store pickup (if available). Sale... Read more
Apple Store Black Friday sale for 2014: $100...
BLACK FRIDAY The Apple Store has posted their Black Friday deals for 2014. Receive a $100 PRODUCT(RED) branded iTunes gift card with the purchase of select Macs, $50 with iPads, and $25 with iPods,... Read more
Black Friday: 15% off iTunes Gift Cards
Staples is offering 15% off $50 and $100 iTunes Gift Cards on their online store as part of their Black Friday sale. Click here for more information. Shipping is free. Best Buy is offering $100... Read more
BEVL Releases Dock Tailored for iPhone 6 and...
Seattle based BEVL has released their first product: an iPhone dock that is divergent in build quality, rock-solid function and visual simplicity to complement the iPhone. BEVL is now accepting... Read more
Black Friday: $150 off 13-inch Retina MacBook...
 Best Buy has 13-inch 2.6GHz Retina MacBook Pros on sale for $150 off MSRP on their online store as part of their Black Friday sale. Choose free shipping or free local store pickup (if available).... Read more

Jobs Board

*Apple* Solutions Consultant (ASC) - Apple (...
**Job Summary** The ASC is an Apple employee who serves as an Apple brand ambassador and influencer in a Reseller's store. The ASC's role is to grow Apple Read more
Senior Event Manager, *Apple* Retail Market...
…This senior level position is responsible for leading and imagining the Apple Retail Team's global event strategy. Delivering an overarching brand story; in-store, 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
*Apple* Solutions Consultant (ASC) - Apple (...
**Job Summary** The ASC is an Apple employee who serves as an Apple brand ambassador and influencer in a Reseller's store. The ASC's role is to grow Apple Read more
*Apple* Solutions Consultant (ASC) - Apple (...
**Job Summary** The ASC is an Apple employee who serves as an Apple brand ambassador and influencer in a Reseller's store. The ASC's role is to grow Apple Read more
All contents are Copyright 1984-2011 by Xplain Corporation. All rights reserved. Theme designed by Icreon.