From: jpc@tauon.ph.unimelb.edu.au (John Costella)
Newsgroups: sci.virtual-worlds
Subject: TECH: Galcube source code
Date: 24 Dec 92 07:54:20 GMT
Article-I.D.: shelley.1hp6r7INNeel
Organization: University of Washington
NNTP-Posting-Host: stein.u.washington.edu
Originator: hlab@stein.u.washington.edu



I have received a number of queries about the source code for Galcube:
whether it would be made available, when, etc.

I wanted to clean up the code a little before posting it up for public
digestion, but with the amount of spare time I have to do this
cleaning up (=1e-15 s :) I have decided that it would be more useful
now-and-dirty than later-and-cleaner. And, since it ain't broke, maybe
I'd better not fix it. :)

So here it is (following). No responsibility taken by me if it blows
up your computer, etc, etc, blah, blah. It is written for the
Microsoft C compiler 6.0 or 7.0. This is the original cheap-and-nasty
program that, when compiled, yielded the executable posted here some
time back. Routines have been pulled out of all manner of musty places
of my hard disk, so conventions in some of them vary a little. :)

John Costella

----------------------------------------------------------------------------
John P. Costella              School of Physics, The University of Melbourne
jpc@tauon.ph.unimelb.edu.au         Tel: +61 3 543-7795, Fax: +61 3 347-4783
----------------------------------------------------------------------------

// GALCUBE.C: Rotate a cube using software simulation of Galilean
//            antialiasing.
//
//            Copyright (C) 1992 John P. Costella. May be copied,
//            used or modified freely for research or development
//            purposes, but notices of modifications should
//            accompany any executable or source code.


// Included files.

#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include <malloc.h>
#include <time.h>
#include <graph.h>



// Definitions.

#define FALSE 0
#define TRUE !FALSE

#define AN_RATE_REQ 25.0
#define RATE_EPSILON 0.5

#define NUM_PAL 16

#define NUM_TRIAL_FRAMES 50

#define START_EDGE_GUESS 50
#define MIN_EDGE_SIZE 20
#define MAX_EDGE_SIZE 80

#define NUM_TEXT_COLUMNS 80

#define DESIRED_REVS_PER_SECOND 0.25

#define SCREEN_WIDTH 640
#define SCREEN_HEIGHT 480

#define SCREEN_ROW_USED 23

#define MAX_XY_VEL (+16383)
#define MIN_XY_VEL (-16384)

#define MAX_XY_ACCEL (+2047)
#define MIN_XY_ACCEL (-2048)

#define MAX_Z_VEL (+127)
#define MIN_Z_VEL (-128)

#define MAX_Z_ACCEL (+31)
#define MIN_Z_ACCEL (-32)

#define MIN_PHYSICAL_Z 0
#define MAX_PHYSICAL_Z 255

#define MAX_FRAC_Z (+1)
#define MIN_FRAC_Z (-2)

#define MAX_GREY_FULL 28
#define MIN_GREY_FULL 0

#define MIN_GREY_DERIV 0
#define MAX_GREY_DERIV 6

#define GALILEAN_MODE_STRING "Galilean-antialiased display"
#define NORMAL_MODE_STRING "Standard (sample-and-hold) display"

#define STDUTILS_KEY_ALTERNATE   0x00
#define STDUTILS_KEY_ALTERNATE_2 0xe0

#define KEY_ESC 0x1b

#define PI_ON_2 1.570796327

#define JPC_CODE  "JPC"

#define WIREFRAME 'w'
#define AMBIENT   'a'
#define DIFFUSE   'd'

#define OPTION_WIREFRAME 0
#define OPTION_AMBIENT   1
#define OPTION_DIFFUSE   2

#define AMBIENT_PAL       8
#define AMBIENT_PAL_EDGE 11

#define BITS_PER_BYTE 8

#define POISSON_EPSILON 96
#define POISSON_MAP_OFFSET 150

#define BR_MULT 10

#define KEY_UP    0x4800
#define KEY_DOWN  0x5000
#define KEY_LEFT  0x4b00
#define KEY_RIGHT 0x4d00

#define KEY_CTRL_UP    0x8d00
#define KEY_CTRL_DOWN  0x9100
#define KEY_CTRL_LEFT  0x7300
#define KEY_CTRL_RIGHT 0x7400

#define KEY_ENTER 0x0d

#define KEY_CTRL_Q 0x11

#define DIFFUSE_LUM_MIN 0.1
#define DIFFUSE_LUM_MAX 1.0

#define DIFFUSE_THETA 35.0
#define DIFFUSE_PHI   20.0



// Custom types.

typedef int boolean;
typedef unsigned char byte;
typedef long RecI;
typedef unsigned int extkey;



// Structure definitions.

typedef struct XYMOTIONTAG
{
  boolean bDebris : 1;

  int iPosFrac : 4;
  int iVel : 15;
  int iAccel : 12;
}
xymotion;


typedef struct ZMOTIONTAG
{
  int iPosFrac : 2;
  int iVel : 8;
  int iAccel : 6;
}
zmotion;


typedef struct GREYSTRUCTTAG
{
  byte byPal : 4;
  int iFrac : 1;
  int iDeriv : 3;
}
greystruct;


typedef struct GALPIXELTAG
{
  byte byZ;
  greystruct gs;
  zmotion zm;
  xymotion xymX;
  xymotion xymY;
}
galpixel;


typedef struct GALPOINTTAG
{
  int iX;
  int iY;
  galpixel gal;
}
galpoint;


typedef struct PHYSPTTAG
{
  double dX;
  double dY;
  double dZ;
}
physpt;


typedef struct SCRNPTTAG
{
  int iX;
  int iY;
  byte byZ;
  byte iFracZ;
}
scrnpt;



// Function prototypes.

int main( int argc, char *argv[] );
int GetAnimationRate( double dRateRequired, long *plImageSize );
void PrintString( char sz[] );
void SetupGreyPalette( void );
void SetupDummyPalette( void );
int ClosestShade( double dShade );
void SetShade( double dShade );
galpixel **RecGalMatrix( RecI riRowFirst, RecI riRowLast, RecI riColumnFirst,
  RecI riColumnLast );
void ClearG2Image( galpixel **aagal, int iPal );
void ComputeG2Image( int iFrame, int iTotalFrames, galpixel **aagal );
void GalContinueG2Image( galpixel **aagal1, galpixel **aagal2, byte
  **aabyPoissonFlags );
void DisplayG2Image( galpixel **aagal, int iTop, int iLeft );
void PlotG2Pixel( galpixel **aagalPixmap, galpoint gpt );
void PlotG2Line( galpixel **aagalPixmap, galpoint gptPt0, galpoint gptPt1 );
int DeltaXYFrac( xymotion xymOld, xymotion *pxyNew );
void DeltaZFrac( byte byOldZ, zmotion zmOld, byte *pbyNewZ, zmotion *pzmNew );
void DeltaGS( greystruct gsOld, greystruct *pgsNew );
byte *RecBYVector( RecI riIndexFirst, RecI riIndexLast );
void RecFreeBYVector( byte *abyVector, RecI riIndexFirst, RecI riIndexLast );
byte **RecBYMatrix( RecI riRowFirst, RecI riRowLast, RecI riColumnFirst, RecI
  riColumnLast );
void RecFreeBYMatrix( byte **aabyMatrix, RecI riRowFirst, RecI riRowLast,
  RecI riColumnFirst, RecI riColumnLast );
physpt SetPhysPt( double dX, double dY, double dZ );
physpt RotatePoint( physpt pptOrig, double dAngle );
scrnpt PhysToScreen( physpt ppt, double dScale, int iXCentre, int iYCentre,
  int iZCentre );
void DrawLine( int iPt1, int iPt2, scrnpt *aspt, xymotion *axymX, xymotion
  *axymY, zmotion *azm, greystruct gs, galpixel **aagal );
void DrawQuad( int iIndexA, int iIndexB, int iIndexC, int iIndexD, scrnpt
  *aspt, xymotion *axymX, xymotion *axymY, zmotion *azm, greystruct gs,
  galpixel **aagal );
void PhysMotionToScreen( physpt pptVel, physpt pptAccel, xymotion *pxymX,
  xymotion *pxymY, zmotion *pzm, double dScale );
void GraftPosAndMotion( galpoint *pgpt, scrnpt spt, xymotion xymX, xymotion
  xymY, zmotion zm );
void DrawStaticBackground( galpixel **aagal, int iBackColour, byte byZ );
void SetPoissonFlag( boolean bValue, byte **aabyFlags, int iRow, int iColumn
  );
boolean GetPoissonFlag( byte **aabyFlags, int iRow, int iColumn );
void ClearPoissonMapDisplay( void );
boolean AllowBrowse( int iFrame, int iNumFramesGal, char szMsg[], galpixel
  **aagal );
void BrowsePos( galpixel **aagal, int iRow, int iColumn, int *piMainPal,
  int *piPoissonPal, char szMsg[] );
void RepairBrowse( galpixel **aagal, int iRow, int iColumn, int iMainPal, int
  iPoissonPal );
int QuadPtComp( const void *pvElem1, const void *pvElem2 );
double Interp( double dProp, double dTop, double dBottom );
galpixel AverageTwoGalpixels( galpixel gal1, galpixel gal2 );
galpixel AverageFourGalpixels( galpixel gal1, galpixel gal2, galpixel gal3,
  galpixel gal4 );
void SetDebrisFlags( galpixel *pgal, boolean bDebris, boolean bOld );
boolean IsDebris( galpixel gal );
boolean IsOld( galpixel gal );
char *RecCHVector( RecI riIndexFirst, RecI riIndexLast );
extkey getkey( void );
int stdkey( extkey key );
void vres_graphics( void );
void reset_screen( void );



// Macros.

#define round( x ) ( floor( (double) (x) + 0.5 ) )
#define booltoggle( x ) ( x = !x )
#define bound( lo, req, hi ) ( max( (lo), min( (hi), (req) ) ) )
#define RGB( r, g, b ) ( 0x3F3F3FL & ( (long) (b) << 16 | (long) (g) << 8 | \
  (long) (r) ) )
#define sqr( x ) ( (x) * (x) )



// Debastardisation.

#define _setcolour( x ) _setcolor( (x) )



// Global variables.

boolean gbJpcMode;
int gaiGreyRaw[ 15 ] = { 0, 5, 6, 7, 9, 11, 13, 16, 20, 24, 29, 35, 43, 52,
  63 }, giBackPal = 0, giNumScreenRows, giNumScreenColumns, giOverdriveRate,
  giShadingOption, giScreenLeft, giScreenTop, gaiXForPts[ 4 ], gaiYForPts[ 4
  ];
double gadGreyLum[ 15 ] = { 0.0, 3.79e-3, 5.67e-3, 7.96e-3, 0.0138, 0.0215,
  0.0311, 0.0490, 0.0801, 0.120, 0.181, 0.274, 0.432, 0.656, 1.0 },
  gadGreyCut[ 14 ] = { 3.10e-3, 4.64e-3, 6.72e-3, 0.0105, 0.0172, 0.0259,
  0.0390, 0.0626, 0.0980, 0.147, 0.223, 0.344, 0.532, 0.810 };



// main.

int main( int argc, char *argv[] )
{
  boolean bGalMode, bFinished, bQuit;
  byte **aabyImages, **aabyPoissonFlags;
  char *szMsg, *szPropTime, *szImageType;
  int iFrame, iSquareEdge, iNumFramesGal;
  long lSingleImageSize;
  clock_t clkEnd, clkFrame;
  extkey key;
  galpixel **aagal1, **aagal2;


  // Check the overdrive rate.

  if( argc < 3 || argc > 4 )
  {
    printf( "\n\nERROR: Command line call is incorrect\n\n" );

    exit( 1 );
  }

  if( argc == 4 )
  {
    if( !strcmp( argv[ 1 ], JPC_CODE ) )
    {
      gbJpcMode = TRUE;

      szPropTime = argv[ 2 ];

      szImageType = argv[ 3 ];
    }
    else
    {
      printf( "\n\nERROR: Too many parameters\n\n" );

      exit( 1 );
    }
  }
  else
  {
    gbJpcMode = FALSE;

    szPropTime = argv[ 1 ];

    szImageType = argv[ 2 ];
  }

  giOverdriveRate = bound( 1, atoi( szPropTime ), (int) floor( AN_RATE_REQ )
    );

  if( tolower( *szImageType ) == WIREFRAME )
    giShadingOption = OPTION_WIREFRAME;
  else if( tolower( *szImageType ) == AMBIENT )
    giShadingOption = OPTION_AMBIENT;
  else if( tolower( *szImageType ) == DIFFUSE )
    giShadingOption = OPTION_DIFFUSE;
  else
  {
    printf( "\n\nERROR: Shading mode unrecognised\n\n" );

    exit( 1 );
  }


  // Test the animation rate until we find 25 frames per second.

  iSquareEdge = GetAnimationRate( AN_RATE_REQ, &lSingleImageSize );

  giNumScreenRows = giNumScreenColumns = iSquareEdge;


  // Check to see if we are at maximum allowed edge size. If so, then
  // compute the number of C clock ticks in a single frame.

  clkFrame = (clock_t) ( (double) CLOCKS_PER_SEC / AN_RATE_REQ );


  // Compute the top and left edges of the screen display of the image.

  giScreenTop = ( SCREEN_HEIGHT - iSquareEdge ) / 2;
  giScreenLeft = ( SCREEN_WIDTH - iSquareEdge ) / 2;


  // Set up the video mode: full mode for JPC operation; dummy mode (no
  // visibility) for others.

  vres_graphics();

  if( gbJpcMode )
    SetupGreyPalette();
  else
    SetupDummyPalette();


  // Draw a bounding rectangle.

  SetShade( 0.1 );

  _rectangle( _GBORDER, (short) ( giScreenLeft - 1 ), (short) ( giScreenTop -
    1 ), (short) ( giScreenLeft + giNumScreenColumns ), (short) ( giScreenTop
    + giNumScreenRows ) );


  // Figure out how many frames are required to display the desired
  // revolution rate at the Galilean rate. Round off to a multiple of the
  // overdrive rate.

  iNumFramesGal = (int) round( AN_RATE_REQ / ( DESIRED_REVS_PER_SECOND * 4.0
    * (double) giOverdriveRate ) ) * giOverdriveRate;


  // Allocate a matrix of images.

  if( ( aabyImages = RecBYMatrix( (RecI) 0, (RecI) ( iNumFramesGal - 1 ),
    (RecI) 0, (RecI) ( lSingleImageSize - 1L ) ) ) == NULL )
  {
    reset_screen();

    printf( "\n\nERROR: Insufficient memory\n\n" );

    exit( 1 );
  }


  // Allocate memory for the two Galilean image planes, etc.

  if( ( aagal1 = RecGalMatrix( (RecI) 0, (RecI) ( giNumScreenRows - 1 ),
      (RecI) 0, (RecI) ( giNumScreenColumns - 1 ) ) ) == NULL ||
    ( aagal2 = RecGalMatrix( (RecI) 0, (RecI) ( giNumScreenRows - 1 ), (RecI)
      0, (RecI) ( giNumScreenColumns - 1 ) ) ) == NULL ||
    ( aabyPoissonFlags = RecBYMatrix( (RecI) 0, (RecI) ( giNumScreenRows - 1
      ), (RecI) 0, (RecI) ( ( giNumScreenColumns - 1 ) / BITS_PER_BYTE + 1 )
      ) ) == NULL ||
    ( szMsg = RecCHVector( (RecI) 0, (RecI) ( 4 * NUM_TEXT_COLUMNS ) ) ) ==
      NULL )
  {
    reset_screen();

    printf( "\n\nERROR: Insufficient memory\n\n" );

    exit( 1 );
  }


  // Loop through each of the Galilean frames.

  for( iFrame = 0, bQuit = FALSE; iFrame < iNumFramesGal && !bQuit; iFrame++
    )
  {
    // Print a message.

    sprintf( szMsg, "Computing frame %d of %d; please stand by...", iFrame +
      1, iNumFramesGal );

    PrintString( szMsg );


    // Optional clearing of the screen area.

    /*
    _setcolour( 0 );

    _rectangle( _GFILLINTERIOR, (short) giScreenLeft, (short) giScreenTop,
      (short) ( giScreenLeft + iSquareEdge - 1 ), (short) ( giScreenTop +
      iSquareEdge - 1 ) );
    */


    // If the current frame is a fully-computed one, compute it. Otherwise,
    // it is Galilean antialiased based on the previous G^(2) image.

    if( iFrame % giOverdriveRate == 0 )
    {
      if( iFrame % 2 == 0 )
      {
        ComputeG2Image( iFrame, iNumFramesGal, aagal1 );

        if( gbJpcMode )
          ClearPoissonMapDisplay();

        DisplayG2Image( aagal1, giScreenTop, giScreenLeft );
      }
      else
      {
        ComputeG2Image( iFrame, iNumFramesGal, aagal2 );

        if( gbJpcMode )
          ClearPoissonMapDisplay();

        DisplayG2Image( aagal2, giScreenTop, giScreenLeft );
      }
    }
    else
    {
      if( iFrame % 2 == 0 )
      {
        GalContinueG2Image( aagal2, aagal1, aabyPoissonFlags );

        DisplayG2Image( aagal1, giScreenTop, giScreenLeft );
      }
      else
      {
        GalContinueG2Image( aagal1, aagal2, aabyPoissonFlags );

        DisplayG2Image( aagal2, giScreenTop, giScreenLeft );
      }
    }


    // Suck up this image into the buffer.

    _getimage( (short) giScreenLeft, (short) giScreenTop, (short) (
      giScreenLeft + giNumScreenColumns - 1 ), (short) ( giScreenTop +
      giNumScreenRows - 1 ), aabyImages[ iFrame ] );


    // Allow a quit using ctrl-Q.

    while( kbhit() )
    {
      if( getkey() == KEY_CTRL_Q )
        bQuit = TRUE;
    }


    // If in JPC mode, give a chance for browsing. Skip if ctrl-Q'ed.

    if( gbJpcMode && !bQuit )
    {
      if( iFrame % 2 == 0 )
        bQuit = AllowBrowse( iFrame, iNumFramesGal, szMsg, aagal1 );
      else
        bQuit = AllowBrowse( iFrame, iNumFramesGal, szMsg, aagal2 );
    }
  }


  // Don't show the animation if we wanted to quit.

  if( !bQuit )
  {
    // Ask for a key to continue.

    PrintString( "Press any key for the simulation..." );

    getkey();


    // Clear the screen; set proper palette if necessary.

    _clearscreen( _GCLEARSCREEN );

    if( !gbJpcMode )
      SetupGreyPalette();


    // Display the animated frames until a key is hit.

    PrintString( NORMAL_MODE_STRING );

    for( iFrame = 0, bGalMode = FALSE, bFinished = FALSE; !bFinished;
      iFrame = ( iFrame == iNumFramesGal - 1 ) ? 0 : iFrame + 1 )
    {
      // Compute the end time.

      clkEnd = clock() + clkFrame;


      // Output the relevant image.

      if( bGalMode )
        _putimage( giScreenLeft, giScreenTop, aabyImages[ iFrame ], _GPSET );
      else
        _putimage( giScreenLeft, giScreenTop, aabyImages[ ( iFrame /
          giOverdriveRate ) * giOverdriveRate ], _GPSET );


      // Check for a keypress; esc exits, all others toggle mode.

      if( kbhit() )
      {
        key = getkey();

        if( key == KEY_ESC )
          bFinished = TRUE;
        else
        {
          booltoggle( bGalMode );

          if( bGalMode )
            PrintString( GALILEAN_MODE_STRING );
          else
            PrintString( NORMAL_MODE_STRING );
        }
      }

      // Wait until the end time is past.

      while( clock() < clkEnd )
        ;
    }
  }


  // Exit OK.

  RecFreeBYMatrix( aabyImages, (RecI) 0, (RecI) ( iNumFramesGal - 1 ), (RecI)
    0, (RecI) ( lSingleImageSize - 1L ) );

  RecFreeBYMatrix( aabyPoissonFlags, (RecI) 0, (RecI) ( giNumScreenRows - 1
    ), (RecI) 0, (RecI) ( ( giNumScreenColumns - 1 ) / BITS_PER_BYTE + 1 ) );

  reset_screen();

  if( bQuit )
    printf( "\nProgram terminated early by user...\n" );

  return 0;
}



// BrowsePos: Browse information at a given position in the image.

void BrowsePos( galpixel **aagal, int iRow, int iColumn, int *piMainPal,
  int *piPoissonPal, char szMsg[] )
{
  int iXInt, iYInt, iZInt, iXFrac, iYFrac, iZFrac, iGreyInt, iGreyFrac,
    iXVelInt, iYVelInt, iZVelInt, iXVelFrac, iYVelFrac, iZVelFrac,
    iXAccelInt, iYAccelInt, iZAccelInt, iXAccelFrac, iYAccelFrac,
    iZAccelFrac, iGreyDerivInt, iGreyDerivFrac;
  galpixel gal;


  // Grab current screen pixels.

  *piMainPal = (int) _getpixel( giScreenLeft + iColumn, giScreenTop + iRow );

  *piPoissonPal = (int) _getpixel( giScreenLeft - POISSON_MAP_OFFSET +
    iColumn, giScreenTop + iRow );


  // Change to the better of black or white for both.

  if( *piMainPal < 7 )
    _setcolour( 14 );
  else
    _setcolour( 0 );

  _setpixel( giScreenLeft + iColumn, giScreenTop + iRow );

  if( *piPoissonPal < 7 )
    _setcolour( 14 );
  else
    _setcolour( 0 );

  _setpixel( giScreenLeft - POISSON_MAP_OFFSET + iColumn, giScreenTop + iRow
    );


  // Retrieve the current galpixel.

  gal = aagal[ iRow ][ iColumn ];


  // Print details.

  iGreyInt = (int) gal.gs.byPal;
  iGreyFrac = gal.gs.iFrac;

  iGreyDerivInt = gal.gs.iDeriv >> 1;
  iGreyDerivFrac = gal.gs.iDeriv - ( iGreyDerivInt << 1 );

  iXInt = iColumn;
  iYInt = iRow;
  iZInt = (int) gal.byZ;

  iXFrac = gal.xymX.iPosFrac;
  iYFrac = gal.xymY.iPosFrac;
  iZFrac = gal.zm.iPosFrac;

  iXVelInt = ( gal.xymX.iVel + 64 ) >> 7;
  iYVelInt = ( gal.xymY.iVel + 64 ) >> 7;
  iZVelInt = ( gal.zm.iVel + 4 ) >> 3;

  iXVelFrac = gal.xymX.iVel - ( iXVelInt << 7 );
  iYVelFrac = gal.xymY.iVel - ( iYVelInt << 7 );
  iZVelFrac = gal.zm.iVel - ( iZVelInt << 3 );

  iXAccelInt = ( gal.xymX.iAccel + 64 ) >> 7;
  iYAccelInt = ( gal.xymY.iAccel + 64 ) >> 7;
  iZAccelInt = ( gal.zm.iAccel + 4 ) >> 3;

  iXAccelFrac = gal.xymX.iAccel - ( iXAccelInt << 7 );
  iYAccelFrac = gal.xymY.iAccel - ( iYAccelInt << 7 );
  iZAccelFrac = gal.zm.iAccel - ( iZAccelInt << 3 );

  sprintf( szMsg, "%sPos=(%d%+d/16,%d%+d/16,%d%+d/4); "
    "Vel=(%d%+d/128,%d%+d/128,%d%+d/8);\n"
    "Accel=(%d%+d/128,%d%+d/128,%d%+d/8); "
    "Grey=%d%+d/2; GreyDeriv=%d%+d/2",
    IsDebris( gal ) ? ( IsOld( gal ) ? "OLD DEBRIS: " : "NEW DEBRIS: " ) :
      "",
    iXInt, iXFrac, iYInt, iYFrac, iZInt, iZFrac,
    iXVelInt, iXVelFrac, iYVelInt, iYVelFrac, iZVelInt, iZVelFrac,
    iXAccelInt, iXAccelFrac, iYAccelInt, iYAccelFrac, iZAccelInt,
      iZAccelFrac,
    iGreyInt, iGreyFrac, iGreyDerivInt, iGreyDerivFrac );

  PrintString( szMsg );

}



// RepairBrowse: Repair the screen at the browse location.

void RepairBrowse( galpixel **aagal, int iRow, int iColumn, int iMainPal, int
  iPoissonPal )
{
  // Redraw the relevant pixels.

  _setcolour( (short) iMainPal );

  _setpixel( giScreenLeft + iColumn, giScreenTop + iRow );

  _setcolour( (short) iPoissonPal );

  _setpixel( giScreenLeft - POISSON_MAP_OFFSET + iColumn, giScreenTop + iRow
    );
}



// AllowBrowse: Allow browsing of information if in JPC mode. Return TRUE
//              if ctrl-q was hit; FALSE otherwise.

boolean AllowBrowse( int iFrame, int iNumFramesGal, char szMsg[], galpixel
  **aagal )
{
  boolean bFinished;
  int iRow, iColumn, iMainPal, iPoissonPal, iMult;
  extkey key;


  // Place a message on the screen.

  sprintf( szMsg, "Frame %d of %d, with previous Poisson flag map on left\n"
    "Press B to browse image, or Enter to continue calculations...", iFrame +
     1, iNumFramesGal );

  PrintString( szMsg );


  // Get a key and check for browse option.

  for( bFinished = FALSE; !bFinished; )
  {
    key = getkey();

    if( tolower( stdkey( key ) ) == 'b' )
      bFinished = TRUE;
    else if( key == KEY_ENTER )
      return FALSE;
    else if( key == KEY_CTRL_Q )
      return TRUE;
  }


  // Do until we are finished.

  for( iRow = iColumn = 0, iMult = 1; ; iMult = 1 )
  {
    // Show the current pos info and get pixel info.

    BrowsePos( aagal, iRow, iColumn, &iMainPal, &iPoissonPal, szMsg );


    // Get a key;

    key = getkey();


    // Repair the screen.

    RepairBrowse( aagal, iRow, iColumn, iMainPal, iPoissonPal );


    // Move or quit as appropriate.

    switch( key )
    {
      case KEY_CTRL_UP:
        iMult = BR_MULT;
      case KEY_UP:
        iRow = bound( 0, iRow - iMult, giNumScreenRows - 1 );
        break;

      case KEY_CTRL_DOWN:
        iMult = BR_MULT;
      case KEY_DOWN:
        iRow = bound( 0, iRow + iMult, giNumScreenRows - 1 );
        break;

      case KEY_CTRL_LEFT:
        iMult = BR_MULT;
      case KEY_LEFT:
        iColumn = bound( 0, iColumn - iMult, giNumScreenColumns - 1 );
        break;

      case KEY_CTRL_RIGHT:
        iMult = BR_MULT;
      case KEY_RIGHT:
        iColumn = bound( 0, iColumn + iMult, giNumScreenColumns - 1 );
        break;

      case KEY_ENTER:
        return FALSE;
        break;

      case KEY_CTRL_Q:
        return TRUE;
        break;

      default:
        break;
    }
  }
}



// ClearPoissonMapDisplay: Clear the Poisson map display area on the screen,
//                         when in JPC mode.

void ClearPoissonMapDisplay( void )
{
  _setcolour( 0 );

  _rectangle( _GFILLINTERIOR, giScreenLeft - POISSON_MAP_OFFSET, giScreenTop,
    giScreenLeft - POISSON_MAP_OFFSET + giNumScreenColumns - 1, giScreenTop
    + giNumScreenRows - 1 );
}



// ComputeG2Image: Compute the G^(2) image for fully-computed frames.

void ComputeG2Image( int iFrame, int iTotalFrames, galpixel **aagal )
{
  static boolean bInitialised = FALSE;
  static int iHalfDiag, iXCentre, iYCentre;
  static double dOmega, dOmegaSqr, dScale, dDiffLum0, dDiffTheta, dDiffPhi,
    dDiffCosTheta, dDiffSinTheta;
  static physpt apptPos[ 7 ], apptVel[ 7 ], apptAccel[ 7 ];

  int iPt;
  double dAngle, dDiffSum, dDiffLumTop, dDiffLumRight, dDiffLumFront,
    dDiffSumSin, dDiffSumCos, dDiffDerivTop, dDiffDerivRight,
    dDiffDerivFront;
  physpt apptRotPos[ 7 ], apptRotVel[ 7 ], apptRotAccel[ 7 ];
  scrnpt aspt[ 7 ];
  greystruct gs, gsTop, gsRight, gsFront;
  xymotion axymX[ 7 ], axymY[ 7 ];
  zmotion azm[ 7 ];


  // Compute some things first time through.

  if( !bInitialised )
  {
    // Compute the centre of the square.

    iXCentre = giNumScreenColumns >> 1;
    iYCentre = giNumScreenRows >> 1;


    // Compute the edge size and scale.

    iHalfDiag = (int) floor( 0.5 * (double) min( giNumScreenRows,
      giNumScreenColumns ) ) - 1;

    dScale = (double) iHalfDiag * 0.62;


    // Compute the angular velocity of one of the corners.

    dOmega = PI_ON_2 / (double) iTotalFrames;

    dOmegaSqr = sqr( dOmega );


    // Compute points that we are starting from, in physical space.

    apptPos[ 0 ] = SetPhysPt( -1.0, +1.0, -1.0 );
    apptPos[ 1 ] = SetPhysPt( +1.0, +1.0, -1.0 );
    apptPos[ 2 ] = SetPhysPt( +1.0, +1.0, +1.0 );
    apptPos[ 3 ] = SetPhysPt( -1.0, +1.0, +1.0 );
    apptPos[ 4 ] = SetPhysPt( +1.0, -1.0, -1.0 );
    apptPos[ 5 ] = SetPhysPt( +1.0, -1.0, +1.0 );
    apptPos[ 6 ] = SetPhysPt( -1.0, -1.0, +1.0 );

    apptVel[ 0 ] = SetPhysPt( +dOmega, 0.0, -dOmega );
    apptVel[ 1 ] = SetPhysPt( +dOmega, 0.0, +dOmega );
    apptVel[ 2 ] = SetPhysPt( -dOmega, 0.0, +dOmega );
    apptVel[ 3 ] = SetPhysPt( -dOmega, 0.0, -dOmega );
    apptVel[ 4 ] = apptVel[ 1 ];
    apptVel[ 5 ] = apptVel[ 2 ];
    apptVel[ 6 ] = apptVel[ 3 ];

    apptAccel[ 0 ] = SetPhysPt( +dOmegaSqr, 0.0, +dOmegaSqr );
    apptAccel[ 1 ] = SetPhysPt( -dOmegaSqr, 0.0, +dOmegaSqr );
    apptAccel[ 2 ] = SetPhysPt( -dOmegaSqr, 0.0, -dOmegaSqr );
    apptAccel[ 3 ] = SetPhysPt( +dOmegaSqr, 0.0, -dOmegaSqr );
    apptAccel[ 4 ] = apptAccel[ 1 ];
    apptAccel[ 5 ] = apptAccel[ 2 ];
    apptAccel[ 6 ] = apptAccel[ 3 ];


    // In diffuse mode, compute relevant parameters.

    if( giShadingOption == OPTION_DIFFUSE )
    {
      dDiffTheta = DIFFUSE_THETA / 90.0 * PI_ON_2;
      dDiffPhi = DIFFUSE_PHI / 90.0 * PI_ON_2;

      dDiffCosTheta = cos( dDiffTheta );
      dDiffSinTheta = sin( dDiffTheta );

      dDiffLum0 = ( DIFFUSE_LUM_MAX - DIFFUSE_LUM_MIN ) / max( dDiffSinTheta,
        dDiffCosTheta );
    }


    // Initialised now.

    bInitialised = TRUE;
  }


  // Clear the image.

  ClearG2Image( aagal, 0 );


  // Compute the angle of rotation where angle is zero for iFrame = -1 and
  // angle is pi / 2 for iFrame = iTotalFrames - 1.

  dAngle = (double) ( iFrame + 1 ) / (double) iTotalFrames * PI_ON_2;


  // Rotate the original points to the rotated points, still in physical
  // space. Then convert them to display space.

  for( iPt = 0; iPt < 7; iPt++ )
  {
    apptRotPos[ iPt ] = RotatePoint( apptPos[ iPt ], dAngle );
    apptRotVel[ iPt ] = RotatePoint( apptVel[ iPt ], dAngle );
    apptRotAccel[ iPt ] = RotatePoint( apptAccel[ iPt ], dAngle );

    aspt[ iPt ] = PhysToScreen( apptRotPos[ iPt ], dScale, iXCentre,
      iYCentre, 128 );

    PhysMotionToScreen( apptRotVel[ iPt ], apptRotAccel[ iPt ], &axymX[ iPt
      ], &axymY[ iPt ], &azm[ iPt ], dScale );
  }


  // Draw a cube, either in wireframe, with ambient lighting, or with
  // diffuse reflections. Include extra edge highlights for ambient
  // lighting, to make the edges visible.

  if( giShadingOption == OPTION_AMBIENT )
  {
    gs.byPal = (byte) AMBIENT_PAL;
    gs.iFrac = 0;
    gs.iDeriv = 0;

    DrawStaticBackground( aagal, 0, MIN_PHYSICAL_Z );

    DrawQuad( 0, 1, 3, 2, aspt, axymX, axymY, azm, gs, aagal );
    DrawQuad( 2, 3, 5, 6, aspt, axymX, axymY, azm, gs, aagal );
    DrawQuad( 1, 2, 4, 5, aspt, axymX, axymY, azm, gs, aagal );
  }
  else if( giShadingOption == OPTION_DIFFUSE )
  {
    // Add the phi component of the light source to the cube's angle.

    dDiffSum = dAngle + dDiffPhi;

    dDiffSumSin = sin( dDiffSum );
    dDiffSumCos = cos( dDiffSum );


    // Compute the luminosity of each face.

    dDiffLumTop = DIFFUSE_LUM_MIN + dDiffLum0 * dDiffSinTheta;

    dDiffLumRight = DIFFUSE_LUM_MIN + dDiffLum0 * dDiffCosTheta * max( 0,
      dDiffSumSin );

    dDiffLumFront = DIFFUSE_LUM_MIN +  dDiffLum0 * dDiffCosTheta * max( 0,
      dDiffSumCos );

    dDiffDerivTop = 0.0;

    if( dDiffSumSin > 0.0 )
      dDiffDerivRight = dDiffLum0 * dOmega * dDiffCosTheta * dDiffSumCos;
    else
      dDiffDerivRight = 0.0;

    if( dDiffSumCos > 0.0 )
      dDiffDerivFront = - dDiffLum0 * dOmega * dDiffCosTheta * dDiffSumSin;
    else
      dDiffDerivFront = 0.0;


    // Convert these to grey scale.

    gsTop.byPal = (byte) bound( 0, round( dDiffLumTop * 14.0 ), 14 );
    gsTop.iFrac = bound( -1, (int) round( ( dDiffLumTop * 14.0 - (double)
      gsTop.byPal ) * 2.0 ), 0 );
    gsTop.iDeriv = bound( -4, (int) round( dDiffDerivTop * 28.0 ), +3 );

    gsRight.byPal = (byte) bound( 0, round( dDiffLumRight * 14.0 ), 14 );
    gsRight.iFrac = bound( -1, (int) round( ( dDiffLumRight * 14.0 - (double)
      gsRight.byPal ) * 2.0 ), 0 );
    gsRight.iDeriv = bound( -4, (int) round( dDiffDerivRight * 28.0 ), +3 );

    gsFront.byPal = (byte) bound( 0, round( dDiffLumFront * 14.0 ), 14 );
    gsFront.iFrac = bound( -1, (int) round( ( dDiffLumFront * 14.0 - (double)
      gsFront.byPal ) * 2.0 ), 0 );
    gsFront.iDeriv = bound( -4, (int) round( dDiffDerivFront * 28.0 ), +3 );

    DrawStaticBackground( aagal, 0, MIN_PHYSICAL_Z );

    DrawQuad( 0, 1, 3, 2, aspt, axymX, axymY, azm, gsTop, aagal );
    DrawQuad( 2, 3, 5, 6, aspt, axymX, axymY, azm, gsFront, aagal );
    DrawQuad( 1, 2, 4, 5, aspt, axymX, axymY, azm, gsRight, aagal );
  }

  if( giShadingOption == OPTION_WIREFRAME || giShadingOption ==
    OPTION_AMBIENT )
  {
    if( giShadingOption == OPTION_WIREFRAME )
      gs.byPal = (byte) 14;
    else
      gs.byPal = (byte) AMBIENT_PAL_EDGE;

    gs.iFrac = 0;
    gs.iDeriv = 0;

    DrawLine( 0, 1, aspt, axymX, axymY, azm, gs, aagal );
    DrawLine( 1, 2, aspt, axymX, axymY, azm, gs, aagal );
    DrawLine( 2, 3, aspt, axymX, axymY, azm, gs, aagal );
    DrawLine( 3, 0, aspt, axymX, axymY, azm, gs, aagal );
    DrawLine( 1, 4, aspt, axymX, axymY, azm, gs, aagal );
    DrawLine( 4, 5, aspt, axymX, axymY, azm, gs, aagal );
    DrawLine( 5, 2, aspt, axymX, axymY, azm, gs, aagal );
    DrawLine( 5, 6, aspt, axymX, axymY, azm, gs, aagal );
    DrawLine( 6, 3, aspt, axymX, axymY, azm, gs, aagal );
  }
}



// DrawStaticBackground: Draw a static background at a finite distance.

void DrawStaticBackground( galpixel **aagal, int iBackColour, byte byZ )
{
  int iRow, iColumn;
  galpixel galPlot;


  // Set up the plottable pixel.

  galPlot.xymX.iPosFrac = 0;
  galPlot.xymX.iVel = 0;
  galPlot.xymX.iAccel = 0;

  galPlot.xymY.iPosFrac = 0;
  galPlot.xymY.iVel = 0;
  galPlot.xymY.iAccel = 0;

  galPlot.zm.iPosFrac = 0;
  galPlot.zm.iVel = 0;
  galPlot.zm.iAccel = 0;

  galPlot.gs.byPal = (byte) iBackColour;
  galPlot.gs.iFrac = 0;
  galPlot.gs.iDeriv = 0;

  galPlot.byZ = byZ;


  // Place the pixels all over the display.

  for( iRow = 0; iRow < giNumScreenRows; iRow++ )
    for( iColumn = 0; iColumn < giNumScreenColumns; iColumn++ )
      aagal[ iRow ][ iColumn ] = galPlot;
}



// QuadPtComp: Comparison routine for qsort, for point y-values.

int QuadPtComp( const void *pvElem1, const void *pvElem2 )
{
  int iX1, iX2, iY1, iY2;


  // Retrieve the two pairs of x and y values.

  iX1 = gaiXForPts[ *(int *)pvElem1 ];
  iX2 = gaiXForPts[ *(int *)pvElem2 ];

  iY1 = gaiYForPts[ *(int *)pvElem1 ];
  iY2 = gaiYForPts[ *(int *)pvElem2 ];


  // Compare them.

  if( iY1 > iY2 )
    return +1;
  else if( iY1 == iY2 )
  {
    if( iX1 > iX2 )
      return +1;
    else if( iX1 == iX2 )
      return 0;
    else
      return -1;
  }
  else
    return -1;
}



// DrawQuad: Draw a parallelogram quadrilateral.
//
//           Completely rewritten 31 October 1992.

void DrawQuad( int iIndexA, int iIndexB, int iIndexC, int iIndexD, scrnpt
  *aspt, xymotion *axymX, xymotion *axymY, zmotion *azm, greystruct gs,
  galpixel **aagal )
{
  int iY, iPtStartTop, iPtStartBottom, iPtEndTop, iPtEndBottom, iPt,
    iIndexStartTop, iIndexStartBottom, iIndexEndTop, iIndexEndBottom,
    aiPtsSorted[ 4 ], aiIndices[ 4 ], aiAllocIndices[ 4 ], iTemp;
  double dStartProp, dEndProp, dStartZ, dEndZ;
  galpoint gptStart, gptEnd;
  scrnpt sptStartTop, sptStartBottom, sptEndTop, sptEndBottom;
  xymotion xymXStartTop, xymXStartBottom, xymXEndTop, xymXEndBottom,
    xymYStartTop, xymYStartBottom, xymYEndTop, xymYEndBottom;
  zmotion zmStartTop, zmStartBottom, zmEndTop, zmEndBottom;


  // Sort the points into y-order.

  gaiXForPts[ 0 ] = aspt[ iIndexA ].iX;
  gaiXForPts[ 1 ] = aspt[ iIndexB ].iX;
  gaiXForPts[ 2 ] = aspt[ iIndexC ].iX;
  gaiXForPts[ 3 ] = aspt[ iIndexD ].iX;

  gaiYForPts[ 0 ] = aspt[ iIndexA ].iY;
  gaiYForPts[ 1 ] = aspt[ iIndexB ].iY;
  gaiYForPts[ 2 ] = aspt[ iIndexC ].iY;
  gaiYForPts[ 3 ] = aspt[ iIndexD ].iY;

  aiAllocIndices[ 0 ] = iIndexA;
  aiAllocIndices[ 1 ] = iIndexB;
  aiAllocIndices[ 2 ] = iIndexC;
  aiAllocIndices[ 3 ] = iIndexD;

  for( iPt = 0; iPt < 4; iPt++ )
    aiPtsSorted[ iPt ] = iPt;

  qsort( (void *) aiPtsSorted, (size_t) 4, sizeof( int ), QuadPtComp );

  for( iPt = 0; iPt < 4; iPt++ )
    aiIndices[ iPt ] = aiAllocIndices[ aiPtsSorted[ iPt ] ];


  // Run through all y-values spanned by the quad.

  for( iY = gaiYForPts[ aiPtsSorted[ 0 ] ]; iY <= gaiYForPts[ aiPtsSorted[ 3
    ] ]; iY++ )
  {
    // Depending on what the current scan line is, change the starting
    // and ending edges, by specifying their top and bottom points.

    if( iY <= gaiYForPts[ aiPtsSorted[ 1 ] ] )
    {
      // In top region. Start on 0-1 edge and finish on 0-2 edge.

      iPtStartTop = 0;
      iPtStartBottom = 1;

      iPtEndTop = 0;
      iPtEndBottom = 2;
    }
    else if( iY <= gaiYForPts[ aiPtsSorted[ 2 ] ] )
    {
      // In middle region. Start on 1-3 edge and finish on 0-2 edge.

      iPtStartTop = 1;
      iPtStartBottom = 3;

      iPtEndTop = 0;
      iPtEndBottom = 2;
    }
    else
    {
      // In bottom region. Start on 1-3 edge and finish on 2-3 edge.

      iPtStartTop = 1;
      iPtStartBottom = 3;

      iPtEndTop = 2;
      iPtEndBottom = 3;
    }


    // Retrieve the original index value for the starting and ending edges'
    // top and bottom points.

    iIndexStartTop = aiIndices[ iPtStartTop ];
    iIndexStartBottom = aiIndices[ iPtStartBottom ];

    iIndexEndTop = aiIndices[ iPtEndTop ];
    iIndexEndBottom = aiIndices[ iPtEndBottom ];


    // Retrieve quantities for the top and bottom points of each of the
    // starting and ending edges.

    sptStartTop = aspt[ iIndexStartTop ];
    sptStartBottom = aspt[ iIndexStartBottom ];
    sptEndTop = aspt[ iIndexEndTop ];
    sptEndBottom = aspt[ iIndexEndBottom ];

    xymXStartTop = axymX[ iIndexStartTop ];
    xymXStartBottom = axymX[ iIndexStartBottom ];
    xymXEndTop = axymX[ iIndexEndTop ];
    xymXEndBottom = axymX[ iIndexEndBottom ];

    xymYStartTop = axymY[ iIndexStartTop ];
    xymYStartBottom = axymY[ iIndexStartBottom ];
    xymYEndTop = axymY[ iIndexEndTop ];
    xymYEndBottom = axymY[ iIndexEndBottom ];

    zmStartTop = azm[ iIndexStartTop ];
    zmStartBottom = azm[ iIndexStartBottom ];
    zmEndTop = azm[ iIndexEndTop ];
    zmEndBottom = azm[ iIndexEndBottom ];


    // Compute how far down the starting and ending edges we are.

    if( ( iTemp = aspt[ iIndexStartBottom ].iY - aspt[ iIndexStartTop ].iY )
        != 0 )
      dStartProp = (double) ( iY - aspt[ iIndexStartTop ].iY ) / (double)
        iTemp;
    else
      dStartProp = 0.0;

    if( ( iTemp = aspt[ iIndexEndBottom ].iY - aspt[ iIndexEndTop ].iY ) !=
        0 )
      dEndProp = (double) ( iY - aspt[ iIndexEndTop ].iY ) / (double) iTemp;
    else
      dEndProp = 0.0;


    // Interpolate starting and ending galpixel values between the
    // end-points; place these values into gptStart and gptEnd.
    //
    // Note: assumes constant greyscake across quad, in present form.

    gptStart.gal.xymX.iPosFrac = gptStart.gal.xymY.iPosFrac = 0;
    gptEnd.gal.xymX.iPosFrac = gptEnd.gal.xymY.iPosFrac = 0;

    gptStart.iY = gptEnd.iY = iY;

    gptStart.gal.gs = gptEnd.gal.gs = gs;

    gptStart.iX = (int) round( Interp( dStartProp, (double) sptStartTop.iX,
      (double) sptStartBottom.iX ) );
    gptEnd.iX = (int) round( Interp( dEndProp, (double) sptEndTop.iX,
      (double) sptEndBottom.iX ) );

    dStartZ = Interp( dStartProp, (double) sptStartTop.byZ + 0.25 * (double)
      zmStartTop.iPosFrac, (double) sptStartBottom.byZ + 0.25 * (double)
      zmStartBottom.iPosFrac );
    dEndZ = Interp( dEndProp, (double) sptEndTop.byZ + 0.25 * (double)
      zmEndTop.iPosFrac, (double) sptEndBottom.byZ + 0.25 * (double)
      zmEndBottom.iPosFrac );

    gptStart.gal.byZ = (byte) round( dStartZ );
    gptEnd.gal.byZ = (byte) round( dEndZ );

    gptStart.gal.zm.iPosFrac = bound( -2, (int) round( ( dStartZ - (double)
      gptStart.gal.byZ ) * 4.0 ), +1 );
    gptEnd.gal.zm.iPosFrac = bound( -2, (int) round( ( dEndZ - (double)
      gptEnd.gal.byZ ) * 4.0 ), +1 );

    gptStart.gal.xymX.iVel = (int) round( Interp( dStartProp, (double)
      xymXStartTop.iVel, (double) xymXStartBottom.iVel ) );
    gptStart.gal.xymX.iAccel = (int) round( Interp( dStartProp, (double)
      xymXStartTop.iAccel, (double) xymXStartBottom.iAccel ) );

    gptEnd.gal.xymX.iVel = (int) round( Interp( dEndProp, (double)
      xymXEndTop.iVel, (double) xymXEndBottom.iVel ) );
    gptEnd.gal.xymX.iAccel = (int) round( Interp( dEndProp, (double)
      xymXEndTop.iAccel, (double) xymXEndBottom.iAccel ) );

    gptStart.gal.xymY.iVel = (int) round( Interp( dStartProp, (double)
      xymYStartTop.iVel, (double) xymYStartBottom.iVel ) );
    gptStart.gal.xymY.iAccel = (int) round( Interp( dStartProp, (double)
      xymYStartTop.iAccel, (double) xymYStartBottom.iAccel ) );

    gptEnd.gal.xymY.iVel = (int) round( Interp( dEndProp, (double)
      xymYEndTop.iVel, (double) xymYEndBottom.iVel ) );
    gptEnd.gal.xymY.iAccel = (int) round( Interp( dEndProp, (double)
      xymYEndTop.iAccel, (double) xymYEndBottom.iAccel ) );

    gptStart.gal.zm.iVel = (int) round( Interp( dStartProp, (double)
      zmStartTop.iVel, (double) zmStartBottom.iVel ) );
    gptStart.gal.zm.iAccel = (int) round( Interp( dStartProp, (double)
      zmStartTop.iAccel, (double) zmStartBottom.iAccel ) );

    gptEnd.gal.zm.iVel = (int) round( Interp( dEndProp, (double)
      zmEndTop.iVel, (double) zmEndBottom.iVel ) );
    gptEnd.gal.zm.iAccel = (int) round( Interp( dEndProp, (double)
      zmEndTop.iAccel, (double) zmEndBottom.iAccel ) );


    // Draw the horizontal line.

    PlotG2Line( aagal, gptStart, gptEnd );
  }
}



// Interp: Interpolate a floating-point value.

double Interp( double dProp, double dTop, double dBottom )
{
  return ( 1.0 - dProp ) * dTop + dProp * dBottom;
}



// DrawLine: Draw a line from ComputeImage.

void DrawLine( int iPt1, int iPt2, scrnpt *aspt, xymotion *axymX, xymotion
  *axymY, zmotion *azm, greystruct gs, galpixel **aagal )
{
  galpoint gpt1, gpt2;


  // Graft the position and motion information for the end-points.

  GraftPosAndMotion( &gpt1, aspt[ iPt1 ], axymX[ iPt1 ], axymY[ iPt1 ], azm[
    iPt1 ] );
  GraftPosAndMotion( &gpt2, aspt[ iPt2 ], axymX[ iPt2 ], axymY[ iPt2 ], azm[
    iPt2 ] );


  // Finish setting up of the galpoints.

  gpt1.gal.xymX.iPosFrac = gpt2.gal.xymX.iPosFrac = 0;
  gpt1.gal.xymY.iPosFrac = gpt2.gal.xymY.iPosFrac = 0;

  gpt1.gal.gs = gpt2.gal.gs = gs;


  // Plot the line.

  PlotG2Line( aagal, gpt1, gpt2 );
}



// SetDebrisFlags: Set the debris flags of a galpixel.

void SetDebrisFlags( galpixel *pgal, boolean bDebris, boolean bOld )
{
  if( bDebris )
    pgal->xymX.bDebris = (byte) 1;
  else
    pgal->xymX.bDebris = (byte) 0;

  if( bOld )
    pgal->xymY.bDebris = (byte) 1;
  else
    pgal->xymY.bDebris = (byte) 0;
}



// IsDebris: Determines whether a galpixel is debris or not.

boolean IsDebris( galpixel gal )
{
  if( (int) gal.xymX.bDebris == 0 )
    return FALSE;
  else
    return TRUE;
}



// IsOld: Determines whether a galpixel is old debris or not.

boolean IsOld( galpixel gal )
{
  if( (int) gal.xymY.bDebris == 0 )
    return FALSE;
  else
    return TRUE;
}



// GraftPosAndMotion: Graft the position and motion values onto an existing
//                    Galilean point.

void GraftPosAndMotion( galpoint *pgpt, scrnpt spt, xymotion xymX, xymotion
  xymY, zmotion zm )
{
  pgpt->iX = spt.iX;
  pgpt->iY = spt.iY;
  pgpt->gal.byZ = spt.byZ;

  pgpt->gal.xymX = xymX;
  pgpt->gal.xymY = xymY;
  pgpt->gal.zm = zm;
  pgpt->gal.zm.iPosFrac = spt.iFracZ;
}



// PhysMotionToScreen: Convert physical motion to fractional screen
//                     coordinates. Have hard-wired number of bits.

void PhysMotionToScreen( physpt pptVel, physpt pptAccel, xymotion *pxymX,
  xymotion *pxymY, zmotion *pzm, double dScale )
{
  int iXVel, iYVel, iZVel, iXAccel, iYAccel, iZAccel;
  double dScrnXVel, dScrnYVel, dScrnZVel, dScrnXAccel, dScrnYAccel,
    dScrnZAccel;


  // Compute screen coordinate values of vel and accel.

  dScrnXVel = dScale * pptVel.dX;
  dScrnYVel = -dScale * 0.8944 * ( pptVel.dY - 0.5 * pptVel.dZ );
  dScrnZVel = dScale * pptVel.dZ;

  dScrnXAccel = dScale * pptAccel.dX;
  dScrnYAccel = -dScale * 0.8944 * ( pptAccel.dY - 0.5 * pptAccel.dZ );
  dScrnZAccel = dScale * pptAccel.dZ;


  // Convert to integral values.

  iXVel = bound( MIN_XY_VEL, (int) round( dScrnXVel * 128.0 ), MAX_XY_VEL );
  iYVel = bound( MIN_XY_VEL, (int) round( dScrnYVel * 128.0 ), MAX_XY_VEL );
  iZVel = bound( MIN_Z_VEL, (int) round( dScrnZVel * 8.0 ), MAX_Z_VEL );

  iXAccel = bound( MIN_XY_ACCEL, (int) round( dScrnXAccel * 128.0 ),
    MAX_XY_ACCEL );
  iYAccel = bound( MIN_XY_ACCEL, (int) round( dScrnYAccel * 128.0 ),
    MAX_XY_ACCEL );
  iZAccel = bound( MIN_Z_ACCEL, (int) round( dScrnZAccel * 8.0 ),
    MAX_Z_ACCEL );


  // Store in the appropriate structures.

  pxymX->iVel = iXVel;
  pxymY->iVel = iYVel;
  pzm->iVel = iZVel;

  pxymX->iAccel = iXAccel;
  pxymY->iAccel = iYAccel;
  pzm->iAccel = iZAccel;
}



// PhysToScreen: Convert a physical space point to screen coords.

scrnpt PhysToScreen( physpt ppt, double dScale, int iXCentre, int iYCentre,
  int iZCentre )
{
  double dZ;
  scrnpt sptRes;


  // Compute z value.

  dZ = bound( (double) MIN_PHYSICAL_Z, (double) iZCentre + ppt.dZ * dScale,
    (double) MAX_PHYSICAL_Z );

  sptRes.byZ = (byte) round( dZ );
  sptRes.iFracZ = bound( -2, (int) round( ( dZ - (double) sptRes.byZ ) * 4.0
    ), +1 );

  sptRes.iX = iXCentre + (int) round( ppt.dX * dScale );

  sptRes.iY = iYCentre - (int) round( 0.8944 * dScale * ( ppt.dY - 0.5 *
    ppt.dZ ) );


  // Return value.

  return sptRes;
}



// SetPhysPt: Set a physical space point.

physpt SetPhysPt( double dX, double dY, double dZ )
{
  physpt pptRes;


  pptRes.dX = dX;
  pptRes.dY = dY;
  pptRes.dZ = dZ;

  return pptRes;
}


// RotatePoint: Rotate a point around the z-axis.

physpt RotatePoint( physpt pptOrig, double dAngle )
{
  double dCos, dSin;
  physpt pptRes;


  // Compute sin and cos.

  dSin = sin( dAngle );
  dCos = cos( dAngle );


  // Compute the rotation.

  pptRes.dX = pptOrig.dX * dCos - pptOrig.dZ * dSin;

  pptRes.dY = pptOrig.dY;

  pptRes.dZ = pptOrig.dZ * dCos + pptOrig.dX * dSin;


  // Return point.

  return pptRes;
}



// PlotG2Line: Plot a line in (flat) G^(2) space, assuming that the line
//             is rigid. Only physical screen coordinates currently
//             supported. Linearly interpolate everything.

void PlotG2Line( galpixel **aagalPixmap, galpoint gptPt0, galpoint gptPt1 )
{
  boolean bHor;
  byte byZ;
  int iX, iY, iXDist, iYDist, iFracZ, iGrey, iGreyDeriv, iFracGrey;
  double dXInc, dYInc, dZInc, dX, dY, dZ, dXVelInc, dYVelInc, dZDist,
    dZVelInc, dXVel, dYVel, dZVel, dXAccelInc, dYAccelInc, dZAccelInc,
    dXAccel, dYAccel, dZAccel, dOneOnXDist, dOneOnYDist, dGrey, dGreyDeriv,
    dGreyInc, dGreyDerivInc;
  galpoint gptPlot, gptStart, gptEnd;


  // Set up useful quantities.

  iXDist = gptPt1.iX - gptPt0.iX;
  iYDist = gptPt1.iY - gptPt0.iY;
  dZDist = (double) gptPt1.gal.byZ - (double) gptPt0.gal.byZ + 0.25 *
    (double) gptPt1.gal.zm.iPosFrac - 0.25 * (double) gptPt0.gal.zm.iPosFrac;


  // If line is more horizontal, run along x values; else y values.

  if( abs( iXDist ) >= abs( iYDist ) )
    bHor = TRUE;
  else
    bHor = FALSE;


  // Set up quantities and plot.

  if( bHor )
  {
    // Start at the more left-most x value. Swap if necessary.

    if( iXDist > 0 )
    {
      gptStart = gptPt0;
      gptEnd = gptPt1;
    }
    else
    {
      gptStart = gptPt1;
      gptEnd = gptPt0;

      iXDist = -iXDist;
      iYDist = -iYDist;
      dZDist = -dZDist;
    }


    // Set up starting floating-point quantities. Y is measured in pixels.
    // Z is measured in integral z-buffer values, and (for historical
    // reasons) includes quarters of units, stored in the zmotion struct.
    // All velocities and accelerations are integral; the fractional bits
    // are not separated off for special treatment here. The grey level is
    // measured in integral screen units for gs.byPal (0 through 14), as well
    // as halves of screen units (stored in gs.byFrac). The grey-level time
    // derivative is integral; its one fractional bit is not separated off
    // for special treatment here.

    dY = (double) gptStart.iY;
    dZ = (double) gptStart.gal.byZ + 0.25 * (double)
      gptStart.gal.zm.iPosFrac;

    dXVel = (double) gptStart.gal.xymX.iVel;
    dYVel = (double) gptStart.gal.xymY.iVel;
    dZVel = (double) gptStart.gal.zm.iVel;

    dXAccel = (double) gptStart.gal.xymX.iAccel;
    dYAccel = (double) gptStart.gal.xymY.iAccel;
    dZAccel = (double) gptStart.gal.zm.iAccel;

    dGrey = (double) gptStart.gal.gs.byPal + 0.5 * (double)
      gptStart.gal.gs.iFrac;

    dGreyDeriv = (double) gptStart.gal.gs.iDeriv;


    // Set up incremental quantities. Not necessary if iXDist is zero (and
    // would indeed bugger up the division). Same units as noted above.

    if( iXDist != 0 )
    {
      dOneOnXDist = 1.0 / (double) iXDist;

      dYInc = (double) iYDist * dOneOnXDist;
      dZInc = dZDist * dOneOnXDist;

      dXVelInc = (double) ( gptEnd.gal.xymX.iVel - gptStart.gal.xymX.iVel )
        * dOneOnXDist;
      dYVelInc = (double) ( gptEnd.gal.xymY.iVel - gptStart.gal.xymY.iVel )
        * dOneOnXDist;
      dZVelInc = (double) ( gptEnd.gal.zm.iVel - gptStart.gal.zm.iVel )
        * dOneOnXDist;

      dXAccelInc = (double) ( gptEnd.gal.xymX.iAccel -
        gptStart.gal.xymX.iAccel ) * dOneOnXDist;
      dYAccelInc = (double) ( gptEnd.gal.xymY.iAccel -
        gptStart.gal.xymY.iAccel ) * dOneOnXDist;
      dZAccelInc = (double) ( gptEnd.gal.zm.iAccel - gptStart.gal.zm.iAccel )
        * dOneOnXDist;

      dGreyInc = ( (double) gptEnd.gal.gs.byPal + 0.5 * (double)
        gptEnd.gal.gs.iFrac - (double) gptStart.gal.gs.byPal - 0.5 *
        (double) gptStart.gal.gs.iFrac ) * dOneOnXDist;

      dGreyDerivInc = ( (double) gptEnd.gal.gs.iDeriv - (double)
        gptStart.gal.gs.iDeriv ) * dOneOnXDist;
    }


    // Run along x values along the line.

    for( iX = gptStart.iX; iX <= gptEnd.iX; iX++ )
    {
      // Compute rounded values for the current quantised x-value. Z buffer
      // and grey scale information is also computed down to fractions.

      iY = (int) round( dY );

      byZ = (byte) round( dZ );
      iFracZ = bound( -2, (int) round( ( dZ - (double) byZ ) * 4.0 ), +1 );

      iGrey = (int) round( dGrey );
      iFracGrey = bound( -1, (int) round( ( dGrey - (double) iGrey ) * 2.0 ),
        0 );

      iGreyDeriv = (int) round( dGreyDeriv );


      // Set up the plottable galpixel. Fractional x-y positions are set
      // to zero to avoid extra galpixel attrition; this could be set
      // to the true fractional position, to partially spatio-temporally
      // antialias the line at the pixel level, but this is probably only
      // of use in wireframe; shaded surfaces are less useful. Also would
      // carry a greater risk of attrition. Clear debris flags.

      gptPlot.iX = iX;
      gptPlot.iY = iY;
      gptPlot.gal.byZ = byZ;

      gptPlot.gal.xymX.iPosFrac = 0;
      gptPlot.gal.xymY.iPosFrac = 0;
      gptPlot.gal.zm.iPosFrac = iFracZ;

      gptPlot.gal.xymX.iVel = (int) round( dXVel );
      gptPlot.gal.xymY.iVel = (int) round( dYVel );
      gptPlot.gal.zm.iVel = (int) round( dZVel );

      gptPlot.gal.xymX.iAccel = (int) round( dXAccel );
      gptPlot.gal.xymY.iAccel = (int) round( dYAccel );
      gptPlot.gal.zm.iAccel = (int) round( dZAccel );

      gptPlot.gal.gs.byPal = (byte) iGrey;
      gptPlot.gal.gs.iFrac = iFracGrey;
      gptPlot.gal.gs.iDeriv = iGreyDeriv;

      SetDebrisFlags( &gptPlot.gal, FALSE, FALSE );

      PlotG2Pixel( aagalPixmap, gptPlot );

      dY += dYInc;
      dZ += dZInc;

      dXVel += dXVelInc;
      dYVel += dYVelInc;
      dZVel += dZVelInc;

      dXAccel += dXAccelInc;
      dYAccel += dYAccelInc;
      dZAccel += dZAccelInc;

      dGrey += dGreyInc;
      dGreyDeriv += dGreyDerivInc;
    }
  }
  else
  {
    // Line is more vertical. Start at the more top-most (smaller) y value.
    // Swap if necessary.

    if( iYDist > 0 )
    {
      gptStart = gptPt0;
      gptEnd = gptPt1;
    }
    else
    {
      gptStart = gptPt1;
      gptEnd = gptPt0;

      iXDist = -iXDist;
      iYDist = -iYDist;
      dZDist = -dZDist;
    }


    // Set up starting floating-point quantities. See above.

    dX = (double) gptStart.iX;
    dZ = (double) gptStart.gal.byZ + 0.25 * (double)
      gptStart.gal.zm.iPosFrac;

    dXVel = (double) gptStart.gal.xymX.iVel;
    dYVel = (double) gptStart.gal.xymY.iVel;
    dZVel = (double) gptStart.gal.zm.iVel;

    dXAccel = (double) gptStart.gal.xymX.iAccel;
    dYAccel = (double) gptStart.gal.xymY.iAccel;
    dZAccel = (double) gptStart.gal.zm.iAccel;

    dGrey = (double) gptStart.gal.gs.byPal + 0.5 * (double)
      gptStart.gal.gs.iFrac;

    dGreyDeriv = (double) gptStart.gal.gs.iDeriv;


    // Set up incremental quantities. Not necessary if iYDist is zero (and
    // would indeed bugger up the division). Same units as noted above.

    if( iYDist != 0 )
    {
      dOneOnYDist = 1.0 / (double) iYDist;

      dXInc = (double) iXDist * dOneOnYDist;
      dZInc = dZDist * dOneOnYDist;

      dXVelInc = (double) ( gptEnd.gal.xymX.iVel - gptStart.gal.xymX.iVel )
        * dOneOnYDist;
      dYVelInc = (double) ( gptEnd.gal.xymY.iVel - gptStart.gal.xymY.iVel )
        * dOneOnYDist;
      dZVelInc = (double) ( gptEnd.gal.zm.iVel - gptStart.gal.zm.iVel )
        * dOneOnYDist;

      dXAccelInc = (double) ( gptEnd.gal.xymX.iAccel -
        gptStart.gal.xymX.iAccel ) * dOneOnYDist;
      dYAccelInc = (double) ( gptEnd.gal.xymY.iAccel -
        gptStart.gal.xymY.iAccel ) * dOneOnYDist;
      dZAccelInc = (double) ( gptEnd.gal.zm.iAccel - gptStart.gal.zm.iAccel )
        * dOneOnYDist;

      dGreyInc = ( (double) gptEnd.gal.gs.byPal + 0.5 * (double)
        gptEnd.gal.gs.iFrac - (double) gptStart.gal.gs.byPal - 0.5 *
        (double) gptStart.gal.gs.iFrac ) * dOneOnYDist;

      dGreyDerivInc = ( (double) gptEnd.gal.gs.iDeriv - (double)
        gptStart.gal.gs.iDeriv ) * dOneOnYDist;
    }


    // Run along y values along the line.

    for( iY = gptStart.iY; iY <= gptEnd.iY; iY++ )
    {
      // Compute rounded values for the current quantised y-value. Z buffer
      // and grey scale information is also computed down to fractions.

      iX = (int) round( dX );

      byZ = (byte) round( dZ );
      iFracZ = bound( -2, (int) round( ( dZ - (double) byZ ) * 4.0 ), +1 );

      iGrey = (int) round( dGrey );
      iFracGrey = bound( -1, (int) round( ( dGrey - (double) iGrey ) * 2.0 ),
        0 );

      iGreyDeriv = (int) round( dGreyDeriv );


      // Set up the plottable galpixel. See above.

      gptPlot.iX = iX;
      gptPlot.iY = iY;
      gptPlot.gal.byZ = byZ;

      gptPlot.gal.xymX.iPosFrac = 0;
      gptPlot.gal.xymY.iPosFrac = 0;
      gptPlot.gal.zm.iPosFrac = iFracZ;

      gptPlot.gal.xymX.iVel = (int) round( dXVel );
      gptPlot.gal.xymY.iVel = (int) round( dYVel );
      gptPlot.gal.zm.iVel = (int) round( dZVel );

      gptPlot.gal.xymX.iAccel = (int) round( dXAccel );
      gptPlot.gal.xymY.iAccel = (int) round( dYAccel );
      gptPlot.gal.zm.iAccel = (int) round( dZAccel );

      gptPlot.gal.gs.byPal = (byte) iGrey;
      gptPlot.gal.gs.iFrac = iFracGrey;
      gptPlot.gal.gs.iDeriv = iGreyDeriv;

      SetDebrisFlags( &gptPlot.gal, FALSE, FALSE );

      PlotG2Pixel( aagalPixmap, gptPlot );

      dX += dXInc;
      dZ += dZInc;

      dXVel += dXVelInc;
      dYVel += dYVelInc;
      dZVel += dZVelInc;

      dXAccel += dXAccelInc;
      dYAccel += dYAccelInc;
      dZAccel += dZAccelInc;

      dGrey += dGreyInc;
      dGreyDeriv += dGreyDerivInc;
    }
  }
}



// PlotG2Pixel: Plot a single G^(2) galpixel.

void PlotG2Pixel( galpixel **aagalPixmap, galpoint gpt )
{
  // Check for out of range.

  if( gpt.iX < 0 || gpt.iX >= giNumScreenColumns || gpt.iY < 0 || gpt.iY >=
      giNumScreenRows )
    return;


  // If all OK, plot the pixel.

  aagalPixmap[ gpt.iY ][ gpt.iX ] = gpt.gal;
}



// ClearG2Image: Clear a G^(2) image structure by filling it with stationary
//               galpixels of the colour specified (often black).

void ClearG2Image( galpixel **aagal, int iPal )
{
  int iRow, iColumn;
  galpixel gal;


  // Set up the pixel that we will store throughout the whole structure.
  // Set the z-buffer distance to be minimum. Make sure debris and old
  // flags are set.

  gal.byZ = (byte) MIN_PHYSICAL_Z;

  gal.gs.byPal = (byte) iPal;
  gal.gs.iFrac = 0;
  gal.gs.iDeriv = 0;

  gal.zm.iPosFrac = MIN_FRAC_Z;
  gal.zm.iVel = 0;
  gal.zm.iAccel = 0;

  gal.xymX.iPosFrac = 0;
  gal.xymX.iVel = 0;
  gal.xymX.iAccel = 0;

  gal.xymY.iPosFrac = 0;
  gal.xymY.iVel = 0;
  gal.xymY.iAccel = 0;

  SetDebrisFlags( &gal, TRUE, TRUE );


  // Store this throughout the G^(2) structure.

  for( iRow = 0; iRow < giNumScreenRows; iRow++ )
    for( iColumn = 0; iColumn < giNumScreenColumns; iColumn++ )
      aagal[ iRow ][ iColumn ] = gal;
}



// DisplayG2Image: Display a G^(2) image on the screen.

void DisplayG2Image( galpixel **aagal, int iTop, int iLeft )
{
  int iOldPal = -1, iPal, iRow, iColumn;


  // Loop through the pixels and plot the colour information.

  for( iRow = 0; iRow < giNumScreenRows; iRow++ )
    for( iColumn = 0; iColumn < giNumScreenColumns; iColumn++ )
    {
      // Don't bother setting colour unless it's different from the last.

      if( ( iPal = (byte) aagal[ iRow ][ iColumn ].gs.byPal ) != iOldPal )
        _setcolour( iOldPal = iPal );


      // Plot the screen pixel.

      _setpixel( (short) ( iLeft + iColumn ), (short) ( iTop + iRow ) );
    }
}



// GalContinueG2Image: Continue the image from aagal1 to aagal2 based on the
//                     galpixel information.

void GalContinueG2Image( galpixel **aagal1, galpixel **aagal2, byte
  **aabyPoissonFlags )
{
  boolean bInSurface, bLeftIn, bAboveIn, bRightIn, bBelowIn, bOrigDebris,
    bDestDebris, bDestOld, bReplace, bFlag;
  byte byNewZ;
  int iRow, iColumn, iOffsetX, iOffsetY, iNewRow, iNewColumn, iLapl,
    iNumAverageable, iAv;
  xymotion xymNewX, xymNewY;
  zmotion zmNew;
  greystruct gsNew;
  xymotion xymX, xymY;
  galpixel galSource, galDest, galRevStart, galNew, galAbove, galBelow,
    galLeft, galRight, agalAv[ 2 ];


  // Compute the Poisson flags, unless we are in non-JPC wireframe mode.

  if( gbJpcMode || giShadingOption != OPTION_WIREFRAME )
  {
    for( iRow = 0; iRow < giNumScreenRows; iRow++ )
      for( iColumn = 0; iColumn < giNumScreenColumns; iColumn++ )
      {
        // Check for edges.

        if( iRow == 0 || iRow == giNumScreenRows - 1 || iColumn == 0 ||
            iColumn == giNumScreenColumns - 1 )
          SetPoissonFlag( FALSE, aabyPoissonFlags, iRow, iColumn );
        else
        {
          // Compute the Laplacian.

          iLapl = ( (int) aagal1[ iRow ][ iColumn + 1 ].byZ << 2 )
            + ( (int) aagal1[ iRow ][ iColumn - 1 ].byZ << 2 )
            + ( (int) aagal1[ iRow - 1 ][ iColumn ].byZ << 2 )
            + ( (int) aagal1[ iRow + 1 ][ iColumn ].byZ << 2 )
            - ( (int) aagal1[ iRow ][ iColumn ].byZ << 4 )
            + (int) aagal1[ iRow ][ iColumn + 1 ].zm.iPosFrac
            + (int) aagal1[ iRow ][ iColumn - 1 ].zm.iPosFrac
            + (int) aagal1[ iRow + 1 ][ iColumn ].zm.iPosFrac
            + (int) aagal1[ iRow - 1 ][ iColumn ].zm.iPosFrac
            - ( (int) aagal1[ iRow - 1 ][ iColumn ].zm.iPosFrac << 2 );


          // Check for being close to zero.

          if( iLapl < POISSON_EPSILON && iLapl > -POISSON_EPSILON )
            SetPoissonFlag( FALSE, aabyPoissonFlags, iRow, iColumn );
          else
            SetPoissonFlag( TRUE, aabyPoissonFlags, iRow, iColumn );
        }
      }
  }


  // Display Poisson flag information to left of actual display, if in
  // JPC mode.

  if( gbJpcMode )
    for( iRow = 0; iRow < giNumScreenRows; iRow++ )
      for( iColumn = 0; iColumn < giNumScreenColumns; iColumn++ )
      {
        if( GetPoissonFlag( aabyPoissonFlags, iRow, iColumn ) )
          _setcolour( 14 );
        else
          _setcolour( 8 );

        _setpixel( giScreenLeft + iColumn - POISSON_MAP_OFFSET, giScreenTop +
           iRow );
      }



  // Start off by simply copying the information from aagal1 to aagal2 as
  // "old debris". For wireframe images, however, clear the galpixmap to
  // the background colour, rather than using debris.

  if( giShadingOption == OPTION_WIREFRAME )
    ClearG2Image( aagal2, giBackPal );
  else
    for( iRow = 0; iRow < giNumScreenRows; iRow++ )
      for( iColumn = 0; iColumn < giNumScreenColumns; iColumn++ )
      {
        aagal2[ iRow ][ iColumn ] = aagal1[ iRow ][ iColumn ];

        SetDebrisFlags( &aagal2[ iRow ][ iColumn ], TRUE, TRUE );
      }


  // Go through the galpixels and propagate them. Any debris (whether marked
  // as old or not) in the source is propagated as new debris; it beats
  // old-debris in the destination, but loses to non-debris; against other
  // new-debris in the destination, it fights off via the z-buffer.
  // Non-debris in the source is propagated as such; it beats any type of
  // debris in the destination.
  //
  // (For *hardware* implementations, the above copying procedure would be
  // done by setting the debris and old flags of the destination en masse,
  // and then propagating the old debris below, in parallel. In this case,
  // old debris beats old debris, but loses out to everything else. This is
  // not needed in this pre-computed simulation.)
  //
  // This debris algorithm supercedes that listed in the paper UM-P-92/105;
  // it was developed on 2 November 1992, after submission of the paper.

  for( iRow = 0; iRow < giNumScreenRows; iRow++ )
    for( iColumn = 0; iColumn < giNumScreenColumns; iColumn++ )
    {
      // Get the current source galpixel.

      galSource = aagal1[ iRow ][ iColumn ];


      // Check for debris nature of the old galpixel.

      bOrigDebris = IsDebris( galSource );


      // Update the x and y motion strcutures and get the pixel offsets.

      iOffsetX = DeltaXYFrac( galSource.xymX, &xymNewX );
      iOffsetY = DeltaXYFrac( galSource.xymY, &xymNewY );

      DeltaZFrac( galSource.byZ, galSource.zm, &byNewZ, &zmNew );

      DeltaGS( galSource.gs, &gsNew );


      // Check if the pixel has moved off the screen.

      iNewRow = iRow + iOffsetY;
      iNewColumn = iColumn + iOffsetX;

      if( iNewRow < giNumScreenRows && iNewRow >= 0 && iNewColumn <
        giNumScreenColumns && iNewColumn >= 0 )
      {
        galDest = aagal2[ iNewRow ][ iNewColumn ];

        bDestDebris = IsDebris( galDest );
        bDestOld = IsOld( galDest );


        // Perform the above debris algorithm checks.

        if( bOrigDebris )
        {
          // Original galpixel is itself some form of debris. Beats old-
          // debris in the destination, loses to non-debris, and fights it
          // off with other new-debris.

          if( bDestDebris )
          {
            if( bDestOld )
              bReplace = TRUE;
            else
            {
              if( byNewZ > galDest.byZ || ( byNewZ == galDest.byZ &&
                  zmNew.iPosFrac > galDest.zm.iPosFrac ) )
                bReplace = TRUE;
              else
                bReplace = FALSE;
            }
          }
          else
            bReplace = FALSE;
        }
        else
        {
          // Original galpixel was non-debris. Beats any debris, but must
          // have a z-buffer duel with any non-debris.

          if( bDestDebris )
            bReplace = TRUE;
          else
          {
            if( byNewZ > galDest.byZ || ( byNewZ == galDest.byZ &&
                zmNew.iPosFrac > galDest.zm.iPosFrac ) )
              bReplace = TRUE;
            else
              bReplace = FALSE;
          }
        }


        // If the above algorithm so specifies, replace the destination
        // galpixel with the source one.

        if( bReplace )
        {
          galDest.byZ = byNewZ;
          galDest.gs = gsNew;
          galDest.zm = zmNew;
          galDest.xymX = xymNewX;
          galDest.xymY = xymNewY;


          // Source debris of any age is propagated as new debris; otherwise
          // it is non-debris.

          if( bOrigDebris )
            SetDebrisFlags( &galDest, TRUE, FALSE );
          else
            SetDebrisFlags( &galDest, FALSE, FALSE );


          // Put the propagated galpixel in place.

          aagal2[ iNewRow ][ iNewColumn ] = galDest;
        }
      }
    }


  // Run through the pixels of the new frame. If we find a piece of debris,
  // reverse-propagate the galpixels surrounding it. Check the Poisson flags
  // of those reverse-propagated galpixels to determine whether the
  // expansion algorithm should be used. If so, average the surrounding
  // galpixels. Only do all this if not in wireframe mode.

  if( giShadingOption != OPTION_WIREFRAME )
  {
    for( iRow = 1; iRow < giNumScreenRows - 1; iRow++ )
      for( iColumn = 1; iColumn < giNumScreenColumns - 1; iColumn++ )
      {
        // Check for debris status.

        if( IsDebris( aagal2[ iRow ][ iColumn ] ) )
        {
          // If any of the galpixels surrounding are in a surface, then
          // assume we are in a surface. Start off assuming not. Note that
          // debris in surrounding pixels may still be used; new debris must
          // be checked for the *same* source position, whereas old debris,
          // like non-debris, must be reverse-propagated.

          bInSurface = FALSE;
          iNumAverageable = 0;


          // Start by getting the galpixels in question.

          galAbove = aagal2[ iRow - 1 ][ iColumn ];
          galBelow = aagal2[ iRow + 1 ][ iColumn ];
          galLeft = aagal2[ iRow ][ iColumn - 1 ];
          galRight = aagal2[ iRow ][ iColumn + 1 ];

          galRevStart = galAbove;

          if( !IsDebris( galRevStart ) || ( IsDebris( galRevStart ) && IsOld(
            galRevStart ) ) )
          {
            xymX = galRevStart.xymX;
            xymY = galRevStart.xymY;

            xymX.iVel = -xymX.iVel;
            xymY.iVel = -xymY.iVel;

            iOffsetX = DeltaXYFrac( xymX, &xymNewX );
            iOffsetY = DeltaXYFrac( xymY, &xymNewY );

            bFlag = GetPoissonFlag( aabyPoissonFlags, iRow - 1 + iOffsetY,
              iColumn + iOffsetX );
          }
          else
            bFlag = GetPoissonFlag( aabyPoissonFlags, iRow - 1, iColumn );

          if( bFlag )
            bAboveIn = FALSE;
          else
          {
            bAboveIn = TRUE;

            bInSurface = TRUE;

            iNumAverageable++;
          }


          galRevStart = galBelow;

          if( !IsDebris( galRevStart ) || ( IsDebris( galRevStart ) && IsOld(
            galRevStart ) ) )
          {
            xymX = galRevStart.xymX;
            xymY = galRevStart.xymY;

            xymX.iVel = -xymX.iVel;
            xymY.iVel = -xymY.iVel;

            iOffsetX = DeltaXYFrac( xymX, &xymNewX );
            iOffsetY = DeltaXYFrac( xymY, &xymNewY );

            bFlag = GetPoissonFlag( aabyPoissonFlags, iRow + 1 + iOffsetY,
              iColumn + iOffsetX );
          }
          else
            bFlag = GetPoissonFlag( aabyPoissonFlags, iRow + 1, iColumn );

          if( bFlag )
            bBelowIn = FALSE;
          else
          {
            bBelowIn = TRUE;

            bInSurface = TRUE;

            iNumAverageable++;
          }


          galRevStart = galLeft;

          if( !IsDebris( galRevStart ) || ( IsDebris( galRevStart ) && IsOld(
            galRevStart ) ) )
          {
            xymX = galRevStart.xymX;
            xymY = galRevStart.xymY;

            xymX.iVel = -xymX.iVel;
            xymY.iVel = -xymY.iVel;

            iOffsetX = DeltaXYFrac( xymX, &xymNewX );
            iOffsetY = DeltaXYFrac( xymY, &xymNewY );

            bFlag = GetPoissonFlag( aabyPoissonFlags, iRow + iOffsetY, iColumn
              - 1 + iOffsetX );
          }
          else
            bFlag = GetPoissonFlag( aabyPoissonFlags, iRow, iColumn - 1 );

          if( bFlag )
            bLeftIn = FALSE;
          else
          {
            bLeftIn = TRUE;

            bInSurface = TRUE;

            iNumAverageable++;
          }


          galRevStart = galRight;

          if( !IsDebris( galRevStart ) || ( IsDebris( galRevStart ) && IsOld(
            galRevStart ) ) )
          {
            xymX = galRevStart.xymX;
            xymY = galRevStart.xymY;

            xymX.iVel = -xymX.iVel;
            xymY.iVel = -xymY.iVel;

            iOffsetX = DeltaXYFrac( xymX, &xymNewX );
            iOffsetY = DeltaXYFrac( xymY, &xymNewY );

            bFlag = GetPoissonFlag( aabyPoissonFlags, iRow + iOffsetY, iColumn
              + 1 + iOffsetX );
          }
          else
            bFlag = GetPoissonFlag( aabyPoissonFlags, iRow, iColumn + 1 );

          if( bFlag )
            bRightIn = FALSE;
          else
          {
            bRightIn = TRUE;

            bInSurface = TRUE;

            iNumAverageable++;
          }


          // Check to see if we are in the surface.

          if( bInSurface )
          {
            // Now average galpixel attributes, depending on how many
            // surrounding galpixels were deemed to be in the surface. If all
            // four were, then we can average (in hardware) by adding them
            // and shifting right two bits, for all relevant attributes.
            // Three, however, cannot be done without a maths coprocessor; in
            // this case, choose the two best for averaging. Two can be done
            // as per four, but shifting right one bit. Finally, if only one
            // surrounding galpixel is in the surface, then we must simply
            // copy it. All of this can be done very simply on a monolithic
            // circuit, and a little more messily with discrete chips; it is
            // all very fast when in hardware.

            switch( iNumAverageable )
            {
              case 1:
                // Find which one is in, and simply copy it.

                if( bAboveIn )
                  galNew = galAbove;
                else if( bBelowIn )
                  galNew = galBelow;
                else if( bLeftIn )
                  galNew = galLeft;
                else
                  galNew = galRight;

                break;

              case 2:
                // Find out which two are in the surface, and average them.

                iAv = 0;

                if( bAboveIn )
                  agalAv[ iAv++ ] = galAbove;

                if( bBelowIn )
                  agalAv[ iAv++ ] = galBelow;

                if( bLeftIn )
                  agalAv[ iAv++ ] = galLeft;

                if( bRightIn )
                  agalAv[ iAv++ ] = galRight;

                if( iAv != 2 )
                {
                  // We've buggered it somewhere, and overwritten memory to
                  // boot. Quit before we trash the machine.

                  reset_screen();

                  printf( "\n\nERROR: Fatal error in Poisson average 2\n\n"
                    );

                  exit( 1 );
                }


                // Average the two that we've found.

                galNew = AverageTwoGalpixels( agalAv[ 0 ], agalAv[ 1 ] );

                break;

              case 3:
                // We have to choose two of the three points to average, for
                // hardware reasons. Choose the two that are diametrically
                // opposite each other with the chosen pixel as centre.
                // Identified uniquely by which one is *not* in the surface.

                if( !bAboveIn || !bBelowIn )
                  galNew = AverageTwoGalpixels( galLeft, galRight );
                else
                  galNew = AverageTwoGalpixels( galAbove, galBelow );

                break;

              case 4:
                // Average all four, since they are all in the surface.

                galNew = AverageFourGalpixels( galAbove, galBelow, galLeft,
                  galRight );

                break;

              default:
                // Should not be here.

                break;
            }


            // Zero the fractional position in x and y directions.

            galNew.xymX.iPosFrac = galNew.xymY.iPosFrac = 0;


            // Store the new (copied) galpixel.

            aagal2[ iRow ][ iColumn ] = galNew;
          }
        }
    }
  }
}



// AverageTwoGalpixels: Average the relevant entries in two galpixel
//                      structures. If both galpixels are debris, the result
//                      is new debris.

galpixel AverageTwoGalpixels( galpixel gal1, galpixel gal2 )
{
  int iZ, iGrey;
  galpixel galRes;


  // Average the relevant quantities.

  iZ = ( ( (int) gal1.byZ << 2 ) + (int) gal1.zm.iPosFrac
    + ( (int) gal2.byZ << 2 ) + (int) gal2.zm.iPosFrac + 1 ) >> 1;

  galRes.byZ = (byte) ( ( iZ + 2 ) >> 2 );
  galRes.zm.iPosFrac = iZ - ( (int) galRes.byZ << 2 );

  iGrey = ( ( (int) gal1.gs.byPal << 1 ) + gal1.gs.iFrac
    + ( (int) gal2.gs.byPal << 1 ) + gal2.gs.iFrac + 1 ) >> 1;

  galRes.gs.byPal = (int) ( ( iGrey + 1 ) >> 1 );
  galRes.gs.iFrac = iGrey - ( (int) galRes.gs.byPal << 1 );

  galRes.gs.iDeriv = ( gal1.gs.iDeriv + gal2.gs.iDeriv + 1 ) >> 1;

  galRes.zm.iVel = ( gal1.zm.iVel + gal2.zm.iVel + 1 ) >> 1;
  galRes.zm.iAccel = ( gal1.zm.iAccel + gal2.zm.iAccel + 1 ) >> 1;

  galRes.xymX.iVel = ( gal1.xymX.iVel + gal2.xymX.iVel + 1 ) >> 1;
  galRes.xymX.iAccel = ( gal1.xymX.iAccel + gal2.xymX.iAccel + 1 ) >> 1;

  galRes.xymY.iVel = ( gal1.xymY.iVel + gal2.xymY.iVel + 1 ) >> 1;
  galRes.xymY.iAccel = ( gal1.xymY.iAccel + gal2.xymY.iAccel + 1 ) >> 1;

  if( IsDebris( gal1 ) && IsDebris( gal2 ) )
    SetDebrisFlags( &galRes, TRUE, FALSE );
  else
    SetDebrisFlags( &galRes, FALSE, FALSE );


  // Return the averaged galpixel.

  return galRes;
}



// AverageFourGalpixels: Average the relevant entries in four galpixel
//                       structures.

galpixel AverageFourGalpixels( galpixel gal1, galpixel gal2, galpixel gal3,
  galpixel gal4 )
{
  int iZ, iGrey;
  galpixel galRes;


  // Average the relevant quantities.

  iZ = ( ( (int) gal1.byZ << 2 ) + (int) gal1.zm.iPosFrac
    + ( (int) gal2.byZ << 2 ) + (int) gal2.zm.iPosFrac
    + ( (int) gal3.byZ << 2 ) + (int) gal3.zm.iPosFrac
    + ( (int) gal4.byZ << 2 ) + (int) gal4.zm.iPosFrac + 2 ) >> 2;

  galRes.byZ = (byte) ( ( iZ + 2 ) >> 2 );
  galRes.zm.iPosFrac = iZ - ( (int) galRes.byZ << 2 );

  iGrey = ( ( (int) gal1.gs.byPal << 1 ) + gal1.gs.iFrac
    + ( (int) gal2.gs.byPal << 1 ) + gal2.gs.iFrac
    + ( (int) gal3.gs.byPal << 1 ) + gal3.gs.iFrac
    + ( (int) gal4.gs.byPal << 1 ) + gal4.gs.iFrac + 2 ) >> 2;

  galRes.gs.byPal = (int) ( ( iGrey + 1 ) >> 1 );
  galRes.gs.iFrac = iGrey - ( (int) galRes.gs.byPal << 1 );

  galRes.gs.iDeriv = ( gal1.gs.iDeriv + gal2.gs.iDeriv + gal3.gs.iDeriv
    + gal4.gs.iDeriv + 2 ) >> 2;

  galRes.zm.iVel = ( gal1.zm.iVel + gal2.zm.iVel + gal3.zm.iVel +
    gal4.zm.iVel + 2 ) >> 2;
  galRes.zm.iAccel = ( gal1.zm.iAccel + gal2.zm.iAccel + gal3.zm.iAccel +
    gal4.zm.iAccel + 2 ) >> 2;

  galRes.xymX.iVel = ( gal1.xymX.iVel + gal2.xymX.iVel + gal3.xymX.iVel +
    gal4.xymX.iVel + 2 ) >> 2;
  galRes.xymX.iAccel = ( gal1.xymX.iAccel + gal2.xymX.iAccel +
    gal3.xymX.iAccel + gal4.xymX.iAccel + 2 ) >> 2;

  galRes.xymY.iVel = ( gal1.xymY.iVel + gal2.xymY.iVel + gal3.xymY.iVel +
    gal4.xymY.iVel + 2 ) >> 2;
  galRes.xymY.iAccel = ( gal1.xymY.iAccel + gal2.xymY.iAccel +
    gal3.xymY.iAccel + gal4.xymY.iAccel + 2 ) >> 2;

  if( IsDebris( gal1 ) && IsDebris( gal2 ) && IsDebris( gal3 ) && IsDebris(
      gal4 ) )
    SetDebrisFlags( &galRes, TRUE, FALSE );
  else
    SetDebrisFlags( &galRes, FALSE, FALSE );


  // Return the averaged galpixel.

  return galRes;
}





// SetPoissonFlag: Set a Poisson flag entry.

void SetPoissonFlag( boolean bValue, byte **aabyFlags, int iRow, int iColumn
  )
{
  static byte byMask[ 8 ] = { 0x01, 0x02, 0x04, 0x08, 0x10, 0x20, 0x40, 0x80
    };

  byte byRead;
  int iByteNum, iBytePos;


  // Compute byte entry.

  iByteNum = iColumn / BITS_PER_BYTE;
  iBytePos = iColumn % BITS_PER_BYTE;


  // Set the bit in question.

  byRead = aabyFlags[ iRow ][ iByteNum ];

  if( bValue )
    byRead |= byMask[ iBytePos ];
  else
    byRead &= ~byMask[ iBytePos ];

  aabyFlags[ iRow ][ iByteNum ] = byRead;
}



// GetPoissonFlag: Get a Poisson flag entry.

boolean GetPoissonFlag( byte **aabyFlags, int iRow, int iColumn )
{
  static byte byMask[ 8 ] = { 0x01, 0x02, 0x04, 0x08, 0x10, 0x20, 0x40, 0x80
    };

  byte byRead;
  int iByteNum, iBytePos;


  // Compute byte entry.

  iByteNum = iColumn / BITS_PER_BYTE;
  iBytePos = iColumn % BITS_PER_BYTE;


  // Get the bit in question.

  byRead = aabyFlags[ iRow ][ iByteNum ];

  if( byRead & byMask[ iBytePos ] )
    return TRUE;
  else
    return FALSE;
}



// DeltaGS: Get the change in grey-scale structure and update it.

void DeltaGS( greystruct gsOld, greystruct *pgsNew )
{
  int iNewGrey;


  // Add the fractional part of the old grey level to the rate of change.

  iNewGrey = (int) gsOld.iFrac + ( (int) gsOld.byPal << 1 ) + (int)
    gsOld.iDeriv;

  if( iNewGrey > MAX_GREY_FULL )
    iNewGrey = MAX_GREY_FULL;
  else if( iNewGrey < MIN_GREY_FULL )
    iNewGrey = MIN_GREY_FULL;


  // Copy the new info in.

  pgsNew->byPal = (byte) ( ( iNewGrey + 1 ) >> 1 );
  pgsNew->iFrac = iNewGrey - ( (int) iNewGrey << 1 );
  pgsNew->iDeriv = gsOld.iDeriv;
}



// DeltaZFrac: Get the change in z-buffer from a zmotion structure. The new
//             zmotion structure is returned in pzNew, and the new z-buffer
//             position in pbyNewZ.

void DeltaZFrac( byte byOldZ, zmotion zmOld, byte *pbyNewZ, zmotion *pzmNew )
{
  int iDeltaPos, iIntegralZ, iNewVel, iNewZ, iPosFrac;


  // Add the fractional part of the position to vel and accel / 2, in the
  // full precision of the vel and accel, and then shift it to remove
  // unwanted bits. (Add 1 to round correctly.)

  iDeltaPos = ( ( ( zmOld.iPosFrac << 1 ) + zmOld.iVel + ( ( zmOld.iAccel + 1
    ) >> 1 ) ) + 1 ) >> 1;


  // Slice off the integral number of pixels. First add one half (or 2 in
  // the integral form) so that we round off correctly.

  iIntegralZ = ( iDeltaPos + 2 ) >> 2;


  // Find the remaining fractional position (positive or negative) by
  // subtracting iIntegralPos pixels (or iIntegralPos * 16 in integral form).

  iPosFrac = iDeltaPos - ( iIntegralZ << 2 );


  // Update the z-buffer depth, and check for overflows.

  iNewZ = (int) byOldZ + iIntegralZ;

  if( iNewZ > MAX_PHYSICAL_Z )
  {
    iNewZ = MAX_PHYSICAL_Z;

    iPosFrac = MAX_FRAC_Z;
  }
  else if( iNewZ < MIN_PHYSICAL_Z )
  {
    iNewZ = MIN_PHYSICAL_Z;

    iPosFrac = MIN_FRAC_Z;
  }

  pzmNew->iPosFrac = iPosFrac;
  *pbyNewZ = iNewZ;


  // Update the velocity value in the new structure. Need to check for
  // overflows.

  iNewVel= zmOld.iVel + zmOld.iAccel;

  if( iNewVel > MAX_Z_VEL )
    iNewVel = MAX_Z_VEL;
  else if( iNewVel < MIN_Z_VEL )
    iNewVel = MIN_Z_VEL;

  pzmNew->iVel = iNewVel;


  // Copy the acceleration over.

  pzmNew->iAccel = zmOld.iAccel;
}



// DeltaXYFrac: Get the change in position from an xymotion structure.
//              The new xymotion structure is returned in pxyNew. The
//              number of complete pixels to move is returned by the function
//              itself.

int DeltaXYFrac( xymotion xymOld, xymotion *pxymNew )
{
  int iDeltaPos, iIntegralPos, iNewVel;


  // Add the fractional part of the position to vel and accel / 2, in the
  // full precision of the vel and accel, and then shift it to remove
  // unwanted bits. (Add 4 to round correctly.)

  iDeltaPos = ( ( ( xymOld.iPosFrac << 3 ) + xymOld.iVel + ( ( xymOld.iAccel
    + 1 ) >> 1 ) ) + 4 ) >> 3;


  // Slice off the integral number of pixels. First add one half (or 8 in
  // the integral form) so that we round off correctly.

  iIntegralPos = ( iDeltaPos + 8 ) >> 4;


  // Find the remaining fractional position (positive or negative) by
  // subtracting iIntegralPos pixels (or iIntegralPos * 16 in integral form).

  pxymNew->iPosFrac = iDeltaPos - ( iIntegralPos << 4 );


  // Update the velocity value in the new structure. Need to check for
  // overflows.

  iNewVel= xymOld.iVel + xymOld.iAccel;

  if( iNewVel > MAX_XY_VEL )
    iNewVel = MAX_XY_VEL;
  else if( iNewVel < MIN_XY_VEL )
    iNewVel = MIN_XY_VEL;

  pxymNew->iVel = iNewVel;


  // Copy the acceleration over.

  pxymNew->iAccel = xymOld.iAccel;


  // Return the integral number of pixels to shift.

  return iIntegralPos;
}



// GetAnimationRate: Find the animation rate of the current machine in the
//                   VGA 640 x 480 16-colour mode. Return the edge size of
//                   the largest square that achieves a RAM-to-screen rate
//                   (using _putimage) of iRateRequired updates per second.
//
//                   NOTE: Only valid on a non-multitasked machine (i.e run
//                         from plain DOS or Exclusive mode of Windows).

int GetAnimationRate( double dRateRequired, long *plImageSize )
{
  boolean bFinished = FALSE;
  byte *abyImage;
  int iFrame, iEdgeSize, iLastEdgeSize;
  long lImageSize;
  double dSeconds, dFramesPerSec;
  time_t timeOrig, timeFinal;


  // Set the vres graphics screen.

  vres_graphics();


  // Print a string.

  PrintString( "Testing machine animation speed..." );


  // Start off by having a square edge size given by the default.

  iEdgeSize = START_EDGE_GUESS;
  iLastEdgeSize = 0;


  // Loop through until we are satisfied.

  while( !bFinished )
  {
    // Find out how much memory is needed for a square image of edge size
    // iEdgeSize.

    lImageSize = _imagesize( (short) 0, (short) 0, (short) ( iEdgeSize - 1 ),
      (short) ( iEdgeSize - 1 ) );


    // Allocate a vector of the required number of bytes.

    if( ( abyImage = RecBYVector( (RecI) 0, (RecI) ( lImageSize - 1L ) ) )
      == NULL )
    {
      reset_screen();

      printf( "\n\nERROR: Insufficient memory\n\n" );

      exit( 1 );
    }


    // Get an image from the screen.

    _getimage( (short) 0, (short) 0, (short) ( iEdgeSize - 1 ), (short) (
       iEdgeSize - 1 ), abyImage );


    // Get the starting time.

    timeOrig = clock();


    // Display the image a number of times in succession.

    for( iFrame = 0; iFrame < NUM_TRIAL_FRAMES; iFrame++ )
    {
      _setcolour( 0 );

      _rectangle( _GFILLINTERIOR, (short) 0, (short) 0, (short) ( iEdgeSize -
        1 ), (short) ( iEdgeSize - 1 ) );

      _putimage( 0, 0, abyImage, _GPSET );
    }


    // Get the final time.

    timeFinal = clock();

    dSeconds = (double) ( timeFinal - timeOrig ) / CLOCKS_PER_SEC;


    // Compute the number of frames displayed per second.

    dFramesPerSec = (double) NUM_TRIAL_FRAMES / dSeconds;


    // Check if this rate is close enough to the desired rate.

    if( fabs( dFramesPerSec - dRateRequired ) <= RATE_EPSILON )
      bFinished = TRUE;
    else
    {
      // If we have hit a bound, or are stuck, then we are finished.

      if( iEdgeSize == MIN_EDGE_SIZE || iEdgeSize == MAX_EDGE_SIZE ||
          iEdgeSize == iLastEdgeSize )
        bFinished = TRUE;
      else
      {
        // Store the last edge size.

        iLastEdgeSize = iEdgeSize;


        // Compute the optimal size, assuming that the time to display
        // the image is proportional to its area.

        iEdgeSize = bound( MIN_EDGE_SIZE, (int) ( (double) iEdgeSize * sqrt(
          dFramesPerSec / dRateRequired ) ), MAX_EDGE_SIZE );
      }
    }

    // Free the memory used for the vector, to be reallocated (if necessary)
    // on the next iteration.

    RecFreeBYVector( abyImage, (RecI) 0, (RecI) ( lImageSize - 1L ) );
  }


  // Set text mode again.

  reset_screen();


  // Return the value of the edge size and image size.

  *plImageSize = lImageSize;

  return iEdgeSize;
}



// Print a string at the centre top of the graphics screen.

void PrintString( char sz[] )
{
  boolean bTwoRows;
  char szOne[] = " ";
  int iColumn, iStringLength, iChar, iSecondRowOffset;


  // Set the text colour.

  _settextcolor( (short) ( NUM_PAL - 1 ) );


  // Go to the left of the first row used.

  _settextposition( (short) SCREEN_ROW_USED, (short) 0 );


  // Erase one row's worth of characters.

  for( iColumn = 0; iColumn < NUM_TEXT_COLUMNS; iColumn++ )
    _outtext( " " );


  // Go to the left of the second row used.

  _settextposition( (short) ( SCREEN_ROW_USED + 1 ), (short) 0 );


  // Erase one row's worth of characters.

  for( iColumn = 0; iColumn < NUM_TEXT_COLUMNS; iColumn++ )
    _outtext( " " );


  // Look for a \n, which signifies the start of the second line.

  for( iChar = 0, bTwoRows = FALSE; iChar < (int) strlen( sz ) && !bTwoRows;
      iChar++ )
    if( sz[ iChar ] == '\n' )
    {
      bTwoRows = TRUE;

      iSecondRowOffset = iChar + 1;
    }


  // Check if one or two rows.

  if( bTwoRows )
  {
    // Find the (limited) length of the first line.

    iStringLength = bound( 0, iSecondRowOffset - 1, NUM_TEXT_COLUMNS );


    // Centre the text.

    _settextposition( (short) SCREEN_ROW_USED, (short) ( ( NUM_TEXT_COLUMNS -
      iStringLength ) / 2 + 1 ) );


    // Print the text.

    for( iChar = 0; iChar < iStringLength; iChar++ )
    {
      szOne[ 0 ] = sz[ iChar ];

      if( szOne[ 0 ] != '\r' && szOne[ 0 ] != '\t' )
        _outtext( szOne );
    }


    // Find the (limited) length of the second line.

    iStringLength = bound( 0, strlen( sz ) - iSecondRowOffset,
      NUM_TEXT_COLUMNS );


    // Centre the text.

    _settextposition( (short) SCREEN_ROW_USED + 1, (short) ( (
      NUM_TEXT_COLUMNS - iStringLength + 1 ) / 2 ) );


    // Print the text.

    for( iChar = 0; iChar < iStringLength; iChar++ )
    {
      szOne[ 0 ] = sz[ iChar + iSecondRowOffset ];

      if( szOne[ 0 ] != '\r' && szOne[ 0 ] != '\t' )
        _outtext( szOne );
    }
  }
  else
  {
    // Find the (limited) length of the string.

    iStringLength = bound( 0, strlen( sz ), NUM_TEXT_COLUMNS );


    // Centre the text.

    _settextposition( (short) SCREEN_ROW_USED, (short) ( ( NUM_TEXT_COLUMNS -
      iStringLength ) / 2 + 1 ) );


    // Print the text.

    for( iChar = 0; iChar < iStringLength; iChar++ )
    {
      szOne[ 0 ] = sz[ iChar ];

      if( szOne[ 0 ] != '\r' && szOne[ 0 ] != '\t' )
        _outtext( szOne );
    }
  }
}



// SetupGreyPalette: Setup a palette whereby palette entries 0 through
//                   NUM_PAL - 2 are grey-scale values, and NUM_PAL - 1
//                   is reserved for system use.

void SetupGreyPalette( void )
{
  int iPal, iGrey;


  // Loop through the palette values.

  for( iPal = 0; iPal < NUM_PAL - 1; iPal++ )
  {
    iGrey = gaiGreyRaw[ iPal ];

    _remappalette( (short) iPal, RGB( iGrey, iGrey, iGrey ) );
  }

  _remappalette( (short) NUM_PAL - 1, RGB( 63, 63, 63 ) );
}



// SetupDummyPalette: Setup a dummy palette, whereby palette entries 0
//                    through NUM_PAL - 2 are all black, and NUM_PAL - 1
//                    is reserved for system use (white).

void SetupDummyPalette( void )
{
  int iPal;


  // Loop through the palette values.

  for( iPal = 0; iPal < NUM_PAL - 1; iPal++ )
    _remappalette( (short) iPal, RGB( 0, 0, 0 ) );

  _remappalette( (short) NUM_PAL - 1, RGB( 63, 63, 63 ) );
}



// SetShade: Set the current colour to the palette entry closest to the
//           desired intensity. The intensity shade should be between 0.0
//           and 1.0 inclusive, although values outside this range are
//           mapped to black and white respectively anyway.

void SetShade( double dShade )
{
  _setcolour( (short) ClosestShade( dShade ) );
}



// ClosestShade: Find the shade of grey closest to the desired value. The
//               input intensity should be a value between 0.0 and 1.0. The
//               value returned is a palette index between 0 and 14.

int ClosestShade( double dShade )
{
  int iPal;


  // Go through the cuts in reverse order, because the latter palette entries
  // encompass a larger range of (linear) intensities.

  for( iPal = NUM_PAL - 3; iPal >= 0; iPal-- )
    if( dShade > gadGreyCut[ iPal ] )
      return iPal + 1;


  // Return black otherwise.

  return 0;
}



// RecGalMatrix: Allocates a matrix of Galilean pixels.

galpixel **RecGalMatrix( RecI riRowFirst, RecI riRowLast, RecI riColumnFirst,
  RecI riColumnLast )
{
  RecI ri;
  galpixel **aagal;


  // Allocate pointers to rows.

  aagal = (galpixel **) halloc( riRowLast - riRowFirst + (RecI) 1, sizeof(
    galpixel * ) );

  if( aagal == NULL )
    return NULL;

  aagal -= riRowFirst;


  // Allocate rows and set pointers to them.

  for( ri = riRowFirst; ri <= riRowLast; ri++ )
  {
    aagal[ ri ] = (galpixel *) halloc( riColumnLast - riColumnFirst + (RecI)
      1, sizeof( galpixel ) );

    if( aagal[ ri ] == NULL )
      return NULL;

    aagal[ ri ] -= riColumnFirst;
  }


  // Return pointer to array of pointers to rows.

  return aagal;
}


// RecBYVector: Allocates a byte vector with range [ riIndexFirst ..
//              riIndexLast ].

byte *RecBYVector( RecI riIndexFirst, RecI riIndexLast )
{
  byte *abyVector;


  abyVector = (byte *) halloc( riIndexLast - riIndexFirst + (RecI) 1, sizeof(
    byte ) );

  if( abyVector == NULL )
    return NULL;

  return abyVector - riIndexFirst;
}



// RecFreeBYVector: Frees a byte vector allocated by RecBYVector().

void RecFreeBYVector( byte *abyVector, RecI riIndexFirst, RecI riIndexLast )
{
  hfree( abyVector + riIndexFirst );
}



// RecBYMatrix: Allocates a byte matrix with range [ riRowFirst .. riRowLast
//              ][ riColumnFirst .. riColumnLast ]

byte **RecBYMatrix( RecI riRowFirst, RecI riRowLast, RecI riColumnFirst, RecI
  riColumnLast )
{
  RecI ri;
  byte **aabyMatrix;


  // Allocate pointers to rows.

  aabyMatrix = (byte **) halloc( riRowLast - riRowFirst + (RecI) 1, sizeof(
    byte * ) );

  if( aabyMatrix == NULL )
    return NULL;

  aabyMatrix -= riRowFirst;


  // Allocate rows and set pointers to them.

  for( ri = riRowFirst; ri <= riRowLast; ri++ )
  {
    aabyMatrix[ ri ] = (byte *) halloc( riColumnLast - riColumnFirst + (RecI)
      1, sizeof( byte ) );

    if( aabyMatrix[ ri ] == NULL )
      return NULL;

    aabyMatrix[ ri ] -= riColumnFirst;
  }


  // Return pointer to array of pointers to rows.

  return aabyMatrix;
}



// RecFreeBYMatrix: Frees a byte matrix allocated with RecBYMatrix().

void RecFreeBYMatrix( byte **aabyMatrix, RecI riRowFirst, RecI riRowLast,
  RecI riColumnFirst, RecI riColumnLast )
{
  RecI ri;


  for( ri = riRowLast; ri >= riRowFirst; ri-- )
    hfree( aabyMatrix[ ri ] + riColumnFirst );

  hfree( aabyMatrix + riRowFirst );
}



// getkey: Get a key from the keyboard, including extended codes. Return
//         the result in an unsigned int. If normal, the lowest eight bits
//         are the code. If extended, these eight bits are zero, and the
//         upper eight bits are the code.

extkey getkey( void )
{
  extkey key;


  key = (extkey) getch();

  if( key == STDUTILS_KEY_ALTERNATE || key == STDUTILS_KEY_ALTERNATE_2 )
  {
    key = (extkey) getch();

    key <<= 8;
  }


  return key;
}



// RecCHVector: Allocates a char vector with range [iIndexFirst..iIndexLast].

char *RecCHVector( RecI riIndexFirst, RecI riIndexLast )
{
  char *achVector;


  achVector = (char *) halloc( riIndexLast - riIndexFirst + (RecI) 1, sizeof(
    char ) );

  if( achVector == NULL )
    return NULL;

  return achVector - riIndexFirst;
}



// stdkey: Return an integer equal to the standard part of the extended key
//         entered.

int stdkey( extkey key )
{
  return (int) ( key & 0x00ff );
}



// vres_graphics: Set the very high resolution graphics mode.

void vres_graphics( void )
{
  _setvideomode( _VRES16COLOR );
}



// Set the usual text mode.

void reset_screen( void )
{
  _setvideomode( _TEXTC80 );
}




John

----------------------------------------------------------------------------
John P. Costella              School of Physics, The University of Melbourne
jpc@tauon.ph.unimelb.edu.au         Tel: +61 3 543-7795, Fax: +61 3 347-4783
----------------------------------------------------------------------------