// ----------------------
// OpenGL Black Hole Simulator.
//
// Written by and Copyright Chris Halsall (chalsall@chalsall.com).
// First published on the O'Reilly Network on Linux.com
// (oreilly.linux.com).  September 2000.  All rights reserved.
//
// This code is licensed under the GNU GPL Version 2.0.
// See (URL: http://www.gnu.org/copyleft/gpl.html ) for details.
//
// Coded to the groovy tunes of Fluke: Risotto.
//
// Dedicated to Stephen W. Hawking, one of the greatest explorers 
// of our time.

#define PROGRAM_TITLE "O'Reilly Net: Black Hole -- C.Halsall"


#include <stdio.h>   // Always a good idea.
#include <stdlib.h>  // For RAND_MAX.
#include <time.h>    // For our FPS stats.
#include <math.h>    // For M_PI
#include <GL/gl.h>   // OpenGL itself.
#include <GL/glu.h>  // GLU support library.
#include <GL/glut.h> // GLUT support library.


// A quick macro to populate a 1x3 vector array. 
#define ourVectInit(a,x,y,z) { (a)[0]=x; (a)[1]=y; (a)[2]=z; }

// Structure to hold all the data for each particle.
typedef struct {
   int Running;
   float Pos[3];  // Position.
   float Vel[3];  // Velocity.
   float Grav[3]; // Acceleration.
   float Color[3];
} Particle;


// ------
// Some global variables.

// Window and Texture IDs, Window sizes.
int Texture_ID;
int Window_ID;
int Window_Width=640;
int Window_Height=480;

// We'll request a sphere from the GLU library at run-time.
struct GLUquadric *Black_Hole;

// "Name" for first (and only) OpenGL display list.
#define STAR_FIELD 1


// Pointer for allocated array of Particles.
Particle *Parts;

// Particle status variables.
int Parts_Running=0;
int Parts_Allocated=0;
int Parts_LastUnused=1;

// Number of parts initially in the system.  Make sure there's at least
// one over a hundred (101, 801), so our Root particle doesn't get deleted.
int Parts_Num = 801;

float Parts_Brightness = 0.15;


// Drawing flags.
int Draw_Axis = 0;
int Draw_Vectors = 0;
int Draw_Stars = 1;
int Heads_Up = 0;
int Texture_On = 1;


// Particle Gun variables.
float Gun_PX = 0;
float Gun_PY = 0;
float Gun_PZ = 4;

float Gun_VX =-0.135;
float Gun_VY = 0.0125;
float Gun_VZ = 0.0;
float Gun_R = 0.005;
float Gun_OffR = 0.050;

// Backwards firing and off-target probabilities.
float Gun_Back = 0.9, Gun_Off = 0.9;

// '.' key toggels between keypad adjusting gun position and eject vect.
int Gun_Control = 1;


// Orbit and motion settings.
int On_Orbit=1;
int Move_Enable=1;
int Move_Step=0;

// Observer initial placement.
float Obs_Angle=114.0;
float Obs_Height=.2;
float Obs_Dist=4.6;

// Calculated observer location.
float Obs[3];

// How quickly do the orbits decay? 
// Lower number (limit 1) is faster decay.

int Decay_Factor = 6000;

// The force of gravity exterted by the black hole.
float Grav = 0.075;

// Somewhat arbitrary values for Event and Escape horizons.

#define EVENT_HORIZON_GRAV .5

#define ESCAPE_HORIZON 30 
#define ESCAPE_HORIZON_SQR (ESCAPE_HORIZON * ESCAPE_HORIZON)


// ------
// Frames per second (FPS) statistic variables and routine.

#define FRAME_RATE_SAMPLES 50
int FrameCount=0;
float FrameRate=0;

static void ourDoFPS(
   void
)
{
   static clock_t last=0;
   clock_t now;
   float delta;

   if (++FrameCount >= FRAME_RATE_SAMPLES) {
      now  = clock();
      delta= (now - last) / (float) CLOCKS_PER_SEC;
      last = now;

      FrameRate = FRAME_RATE_SAMPLES / delta;
      FrameCount = 0;
   }
}


// ------
// String rendering routine; leverages on GLUT routine.

static void ourPrintString(
   void *font,
   char *str
)
{
   int i,l=strlen(str);

   for(i=0;i<l;i++)
      glutBitmapCharacter(font,*str++);
}


// ------
// Reallocates our array of particles, adding to or removing as
// requested.

Particle *ourAllocParticles(
   int Num
)
{
   int i;
   Particle *P;

   P = realloc(Parts, sizeof(Particle) * Num);

   if (!P)
      return 0;

   if (Parts_Allocated < Num)
      memset( &P[Parts_Allocated],
         0,sizeof(Particle) * (Num-Parts_Allocated));

   if (Parts_LastUnused > Num)
      Parts_LastUnused = Num;

   Parts_Running = 0;
   for (i = 0; i < Num; i++)
      if (P[i].Running) 
         Parts_Running++;

   Parts_Allocated = Num;
   Parts = P;

   return P;
}


// ------
// Function to return a random floating number between 0 and the passed
// parameter.

float ourRand(
   float Max
)
{
   return( (Max * rand()) / RAND_MAX );
}


// ------
// Builds a Display List containing a random star field.
//
// Note: this could also be done by calculating the star points in this 
// routine, which would be faster than having OpenGL perform two 
// rotations (matrix multiplications) for each star.  However, this 
// technique is simpler and faster for the programmer, and demonstrates 
// how successive transformations can be a powerful tool.

void ourBuildStarfield(
   int Stars
)
{
   int Cnt;

   glNewList(STAR_FIELD, GL_COMPILE);

   glMatrixMode(GL_MODELVIEW);
   glPushMatrix();

   for ( Cnt = 0; Cnt < Stars; Cnt++) {

      // Vary the color for each star.
      glColor4f(
         0.8 + ourRand(0.2),
         0.8 + ourRand(0.2),
         0.8 + ourRand(0.2),
         .95); 

      // Vary the size.  Ensure integer sizes to avoid alias shimmering.
      glPointSize(ourRand(2) > 1 ? 1.0 : 2.0);

      // Spin your Universe, round and round....
      glRotatef(ourRand(100),0.0f,1.0f,0.0f);
      glRotatef(ourRand(100),1.0f,0.0f,0.0f);

      glBegin(GL_POINTS);
         glVertex3f(15.0, 0.0f, 0.0f);
      glEnd();
   }

   glPopMatrix();
   glEndList();
}


// ------
// Function builds a simple alpha channel texture of a dot,
// and then creates mipmaps.  This could instead load textures from
// graphics files from disk, or render textures based on external
// input.

void ourBuildTextures(
   void
)
{
   GLenum gluerr;
   GLubyte tex[128][128];
   int x,y,t;
   int hole_size = 3300; // ~ == 57.45 ^ 2.

   // Generate a texture index, then bind it for future operations.
   glGenTextures(1,&Texture_ID);
   glBindTexture(GL_TEXTURE_2D,Texture_ID);

   // Iterate across the texture array.
  
   for(y=0;y<128;y++) {
      for(x=0;x<128;x++) {

        // Make a round dot in the texture's alpha-channel.

         // Calculate distance to center (squared).
         t = (x-64)*(x-64) + (y-64)*(y-64) ;

         if ( t < hole_size) // Don't take square root; compare squared.
            tex[x][y]= 240 - (240 * t) / hole_size + ourRand(15);
         else
            tex[x][y]=0;   // Outside of the dot, it's transparent.

      }
   }

   // The GLU library helps us build MipMaps for our texture.

   if ((gluerr=gluBuild2DMipmaps(GL_TEXTURE_2D, 1, 128, 128, GL_ALPHA,
                 GL_UNSIGNED_BYTE, (void *)tex))) {

      fprintf(stderr,"GLULib%s\n",gluErrorString(gluerr));
      exit(-1);
   }

   // Some pretty standard settings for wrapping and filtering.
   glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_WRAP_S,GL_REPEAT);
   glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_WRAP_T,GL_REPEAT);

   glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_MIN_FILTER,GL_LINEAR_MIPMAP_LINEAR);
   glTexParameterf(GL_TEXTURE_2D,GL_TEXTURE_MAG_FILTER,GL_LINEAR);

   // We start with GL_MODULATE mode.
   glTexEnvf(GL_TEXTURE_ENV,GL_TEXTURE_ENV_MODE,GL_MODULATE);
}


// ------
// Callback routine executed whenever our window is resized.  Lets us
// request the newly appropriate perspective projection matrix for
// our needs.  Try removing the gluPerspective() call to see what happens.

void cbResizeScene(
   int Width,
   int Height
)
{
   // Let's not core dump, no matter what.
   if (Height == 0)
      Height = 1;

   glViewport(0, 0, Width, Height);

   glMatrixMode(GL_PROJECTION);
   glLoadIdentity();
   gluPerspective(45.0f,(GLfloat)Width/(GLfloat)Height,0.05f,100.0f);

   glMatrixMode(GL_MODELVIEW);

   Window_Width  = Width;
   Window_Height = Height;
}


// ------
// Fires the Particle Gun, or, sets up the passed Particle to be
// placed at the Particle Gun location, and fired in a direction
// specified in Gun_Va, with 'a' being "X", "Y", or "Z".
//
// The particles are normally fired with a randomness of Gun_R.
// Occationally (by default 10%) a larger randomness of Gun_OffR is
// added.  Also, 10% of the time, the particles are fired backwards.
// This is used to introduce a bit of non-uniformity.  Gun_R and
// Gun_OffR can be controled with the '3' and '6', and '/' and '*'
// keys, respectively.  If set to zero, effects are removed.

static void ourFireParticleGun(
   Particle *p
)
{
   float r;
   int c;
   int Dir = 1;

   if (!p->Running) {
      p->Running=1;
      Parts_Running++;
   }

   if (p == Parts) {  // Root part.
       ourVectInit(p->Color,1,1,0);                // Bright Yellow
       ourVectInit(p->Pos,5.0,0,0);                // Location
       ourVectInit(p->Vel,0,Gun_VY,-Gun_VX*0.95);  // Velocity
       return;
   }

   // Regular particle.

   ourVectInit(p->Pos,Gun_PX,Gun_PY,Gun_PZ);
   r = Gun_R;

// This creates a few a very unpredicatable trajectories.  It actually
// works out to be much less than a full 10 percent, as many are eatten
// or escape within a very short period of time.  Only a few enter a 
// stable orbit. 

   if (ourRand(1) > Gun_Off) r += Gun_OffR; 

   if (ourRand(1) > Gun_Back) Dir = -1;

   ourVectInit(p->Vel,
      (Gun_VX + r-ourRand(2*r)) * Dir,
      (Gun_VY + r-ourRand(2*r)) * Dir,
      (Gun_VZ + r-ourRand(2*r)) * Dir);

   c = (int)(ourRand(5) + 1.5); // Range of 1 to 6.

   // The last set of numbers bias the colors to blue.  Red is nice too.

   ourVectInit(p->Color,
      (c & 0x01 ? 0.9 : 1.0) * 0.7,
      (c & 0x02 ? 0.9 : 1.0) * 0.7,
      (c & 0x04 ? 0.9 : 1.0) * 1.0
      );
}


// ------
// Calculates the next position based on the current velocity vector,
// then calculates the new velocity vector based on the particle's 
// proximity to the black hole.
//
// We do the motion calculation before updating the velocity vector
// (and calculating the acceleration-because-of-gravity vector) so 
// that our Vector Display option will be correct.  If we didn't do
// this, the gravity vector would not point towards (0,0,0) when
// we drew it, outside of this function.

static void ourMoveParticle(
   Particle *p
)
{
   float dp2, dsx, dsy, dsz, G, d;

   //  Used to randomly kill and re-create particles.
   if (p != Parts)
      if (ourRand(1) >  0.9998) {
      ourFireParticleGun(p);
      return;
   }

   // We're actually going to move this particle...

   // We first move it based on the LAST iteration's 
   // calculation of the Velocity...
   p->Pos[0] += p->Vel[0];
   p->Pos[1] += p->Vel[1];
   p->Pos[2] += p->Vel[2];

   // ...and then proceed to calculate the force of gravity at the new 
   // position, and update our velocity vector.

   dsx = p->Pos[0] * p->Pos[0];
   dsy = p->Pos[1] * p->Pos[1];
   dsz = p->Pos[2] * p->Pos[2];

   // Calculate the square of the distance.
   dp2 = dsx + dsy + dsz;

   if (dp2) {
      // May wish to scale dp2 (change 1.0); effects gravity gradiant.
      G = Grav / (dp2 * 1.0);  
      d = sqrt(dp2);
   }

   // If the force of gravity is too strong, our algorithim breaks
   // down and conservation of energy isn't maintained.  We consider 
   // this the event horizon, and recycle the particle.
   if (G > EVENT_HORIZON_GRAV) {
      ourFireParticleGun(p);
      return;
   }

   if (dp2 > ESCAPE_HORIZON_SQR) { 
      // Particle escaped; lucky it.
      ourFireParticleGun(p);
      return;
   }

   // OK, this particle is staying in the system.  Calculate the 
   // vectors....

   // We store the components of the force of gravity for 
   // our Vectors display.  Note the negative magnitude; the vector
   // must point _from_ our particle _towards_ (0,0,0).
   p->Grav[0] = - G * p->Pos[0] / d;
   p->Grav[1] = - G * p->Pos[1] / d;
   p->Grav[2] = - G * p->Pos[2] / d;

   // Simply add the gravity vector to the current velocity vector.
   p->Vel[0] += p->Grav[0];
   p->Vel[1] += p->Grav[1];
   p->Vel[2] += p->Grav[2];

   if (p != Parts) {
      // This handles orbit decay; not correctly, but well enough.
      // (Decay should be a ratio to the vector length, applied to each 
      // vector component here, rather than each component being effected
      // based on its individual size.  The effect is to circlurize the
      // orbit, which we want anyway.)

      p->Vel[0] -= p->Vel[0] / Decay_Factor;
      p->Vel[1] -= p->Vel[1] / Decay_Factor;
      p->Vel[2] -= p->Vel[2] / Decay_Factor;
   }
}


// ------
// Angle to Radian conversion.

float ourA2R(
   float Angle
)
{
   return Angle * M_PI/180;
}


// ------
// Calculates the observer's XYZ position from their Distance from 
// the origin, the angle and the height.

static void ourCalcObs(void)
{
  Obs[0]=Obs_Dist * sin(ourA2R(Obs_Angle));
  Obs[1]=Obs_Height;
  Obs[2]=Obs_Dist * cos(ourA2R(Obs_Angle));
}


// ------
// Draws the X, Y, and Z axis lines.

void ourRenderAxis(
   void
)
{
   glBegin(GL_LINES);

   glColor4f(0.5,0.5,0.0,1.0);  // Mid-level yellow.

   // Three primary axis lines.
   glVertex3f(100,0,0);
   glVertex3f(-100,0,0);
   glVertex3f(0,100,0);
   glVertex3f(0,-100,0);
   glVertex3f(0,0,100);
   glVertex3f(0,0,-100);

   glColor4f(0.25,0.25,0.0,1.0);  // Low-level yellow.

   // Two pairs of secondary lines for X and Z axis.
   glVertex3f(100,1,0);
   glVertex3f(-100,1,0);
   glVertex3f(100,-1,0);
   glVertex3f(-100,-1,0);

   glVertex3f(0,1,100);
   glVertex3f(0,1,-100);
   glVertex3f(0,-1,100);
   glVertex3f(0,-1,-100);

   glColor4f(0.0,0.5,0.0,1.0); // Mid-level green.

   // Lable the X axis.
   glVertex3f(1.0,0.9,0);
   glVertex3f(1.1,0.8,0);
   glVertex3f(1.1,0.9,0);
   glVertex3f(1.0,0.8,0);

   // And the Z.
   glVertex3f(0,0.9,1.0);
   glVertex3f(0,0.9,1.1);
   glVertex3f(0,0.9,1.1);
   glVertex3f(0,0.8,1.0);
   glVertex3f(0,0.8,1.0);
   glVertex3f(0,0.8,1.1);

   glEnd();
}


// ------
// Draws the Gravity and Velocity Vectors for each active particle.

void ourRenderVectors(
   void
)
{
   int i;

   glBegin(GL_LINES);

   for (i=0; i<Parts_Num; i++) {

      if (!Parts[i].Running) continue;

      // Draw the velocity vector as green.
      glColor4f(0.0,1.0,0.0,Parts_Brightness + .5);

      glVertex3f(
         Parts[i].Pos[0],
         Parts[i].Pos[1],
         Parts[i].Pos[2]
         );

      glVertex3f(
         Parts[i].Pos[0] + Parts[i].Vel[0] ,
         Parts[i].Pos[1] + Parts[i].Vel[1] ,
         Parts[i].Pos[2] + Parts[i].Vel[2] 
         );

      // Draw the gravity vector as red.
      glColor4f(1.0,0.0,0.0,Parts_Brightness + .5);

      glVertex3f(
         Parts[i].Pos[0],
         Parts[i].Pos[1],
         Parts[i].Pos[2]
         );

      glVertex3f(
         Parts[i].Pos[0] + Parts[i].Grav[0] ,
         Parts[i].Pos[1] + Parts[i].Grav[1] ,
         Parts[i].Pos[2] + Parts[i].Grav[2] 
      );
    }

    glEnd();
}


// ------
// Draws the heads-up-display.
void ourRenderHeadsUp(
  void
)
{
   char buf[80];

   glLoadIdentity();
     // We need to change the projection matrix for the text rendering.
   glMatrixMode(GL_PROJECTION);

   // But we like our current view too; so we save it here.
   glPushMatrix();

   // Now we set up a new projection for the text.
   glLoadIdentity();
   glOrtho(0,Window_Width,0,Window_Height,-1.0,1.0);

   // No need for textured text.
   glDisable(GL_TEXTURE_2D);

   // We don't want depth-testing either.
   glDisable(GL_DEPTH_TEST);


   // Draw various variables for the user.

   glColor3f(1.0,1.0,0.0);

   sprintf(buf,"Parts: %d / %d  Bright:%.2f", 
      Parts_Running, Parts_Allocated, Parts_Brightness);
   glRasterPos2i(10,48);
   ourPrintString(GLUT_BITMAP_HELVETICA_12,buf);

   sprintf(buf,"Rnd - Normal:%.3f Extreme:%.3f",
      Gun_R, Gun_OffR);
   glRasterPos2i(10,6);
   ourPrintString(GLUT_BITMAP_HELVETICA_12,buf);

   if (Gun_Control)
      glColor3f(1.0,1.0,0.0);
   else
      glColor3f(0.5,1.0,0.0);

   sprintf(buf,"GunP: (%.3f,%.3f,%.3f)", Gun_PX, Gun_PY, Gun_PZ);
   glRasterPos2i(10,34);
   ourPrintString(GLUT_BITMAP_HELVETICA_12,buf);

   if (!Gun_Control)
      glColor3f(1.0,1.0,0.0);
   else
      glColor3f(0.5,1.0,0.0);

   sprintf(buf,"GunV: (%.3f,%.3f,%.3f)",
      Gun_VX, Gun_VY, Gun_VZ);
   glRasterPos2i(10,20);
   ourPrintString(GLUT_BITMAP_HELVETICA_12,buf);

   // Now we want to render the calulated FPS at the top.
  
   // To ease, simply translate up.  Note we're working in screen
   // pixels in this projection.
  
   glTranslatef(6.0f,Window_Height - 14,0.0f);

   glColor4f(0.9,0.2,0.2,.95);
   sprintf(buf,"FPS: %f F: %2d", FrameRate, FrameCount);
   glRasterPos2i(6,0);
   ourPrintString(GLUT_BITMAP_HELVETICA_12,buf);

   // Lets also show the current position of the Root Particle
   sprintf(buf,"RootP: ( %-.4f, %-.4f, %-.4f)", 
      Parts->Pos[0], Parts->Pos[1],Parts->Pos[2]);
   glRasterPos2i(6,-16);
   ourPrintString(GLUT_BITMAP_HELVETICA_12,buf);

   // And the Root Particle's velocity.
   sprintf(buf,"RootV: ( %-.4f, %-.4f, %-.4f)", 
      Parts->Vel[0], Parts->Vel[1],Parts->Vel[2]);
   glRasterPos2i(6,-32);
   ourPrintString(GLUT_BITMAP_HELVETICA_12,buf);

   // Done with this special projection matrix.  Throw it away.
   glPopMatrix();
}


// ------
// Routine which actually does the drawing

static void cbRenderScene(void)
{
   Particle *p;
   int i;

   // For the first few objects, we want full depth-buffer testing.
   glEnable(GL_DEPTH_TEST);
   glDepthMask(GL_TRUE);

   // Clear the screen.
   glClearColor(0.00,0.00,0.00,1.0);
   glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);

   // Ensure we're working with the model matrix.
   glMatrixMode(GL_MODELVIEW);

   // Reset to 0,0,0; no rotation, no scaling.
   glLoadIdentity();

   // Are we on-orbit, or wandering around?
   if (On_Orbit) {
      gluLookAt(Parts->Pos[0],Parts->Pos[1],Parts->Pos[2],
         0.0,0.0,0.0,
         0.0,1.0,0.0);
   } else {
      gluLookAt(Obs[0],Obs[1],Obs[2],
         0.0,0.0,0.0,
         0.0,1.0,0.0);
   }

   // No texturing.
   glDisable(GL_TEXTURE_2D);

   // Black holes are BLACK!
   glColor4f(0.0,0.0,0,1.0);
   gluSphere(Black_Hole,0.02,8,8);

   if (Draw_Stars)
      glCallList(STAR_FIELD);

   if (Draw_Vectors)
      ourRenderVectors();

   if (Draw_Axis)
      ourRenderAxis();

   // We don't want any of the particles to obscure any others, but
   // we DO want the black hole to block any particles behind it.
   // Note that GL_DEPTH_TEST is still enabled.
   glDepthMask(GL_FALSE);

   // Enable the dot texture.  "Oh no!  Not THE DOT!"
   if (Texture_On)
      glEnable(GL_TEXTURE_2D);

   // Iterate through the array of particles, drawing all that are
   // active.  For those that aren't active, 0.03% of the time, we 
   // introduce them into the system.

   for(i=0; i<Parts_Num; i++) {
      p = &Parts[i];

      if (!p->Running) {
         if (Move_Enable && ourRand(1) > 0.9997)
            ourFireParticleGun(p);
      } else {
         // Set the part's color.
         glColor4f(p->Color[0],p->Color[1],p->Color[2],
            Parts_Brightness);

         // Draw two intersecting quads, along XY and ZY axis.
         glBegin(GL_QUADS);

         glTexCoord2f(0.0,0.0);
         glVertex3f(p->Pos[0]-.00,p->Pos[1]-.10,p->Pos[2]-.10);
         glTexCoord2f(1.0,0.0);
         glVertex3f(p->Pos[0]-.00,p->Pos[1]+.10,p->Pos[2]-.10);
         glTexCoord2f(1.0,1.0);
         glVertex3f(p->Pos[0]-.00,p->Pos[1]+.10,p->Pos[2]+.10);
         glTexCoord2f(0.0,1.0);
         glVertex3f(p->Pos[0]-.00,p->Pos[1]-.10,p->Pos[2]+.10);

         glTexCoord2f(0.0,0.0);
         glVertex3f(p->Pos[0]-.10,p->Pos[1]-.10,p->Pos[2]-.00);
         glTexCoord2f(1.0,0.0);
         glVertex3f(p->Pos[0]-.10,p->Pos[1]+.10,p->Pos[2]-.00);
         glTexCoord2f(1.0,1.0);
         glVertex3f(p->Pos[0]+.10,p->Pos[1]+.10,p->Pos[2]+.00);
         glTexCoord2f(0.0,1.0);
         glVertex3f(p->Pos[0]+.10,p->Pos[1]-.10,p->Pos[2]+.00);

         glEnd();

         if (Move_Enable)
            ourMoveParticle(p);
      }
   }


   if (Heads_Up)
      ourRenderHeadsUp();

   // All done drawing.  Let's show it.
   glutSwapBuffers();

   // This handles our single-step function.
   if (Move_Step)
      Move_Step = Move_Enable = 0;

   // Collect the FPS statistics.
   ourDoFPS();
}


// ------
// Callback function called when a normal key is pressed.

void cbKeyPressed(
   unsigned char key, 
   int x, int y
)
{
   int t;

   switch (key) {
   case 27:
   case 'q': case 'Q':
      exit(0);
      break;

   // Toggle drawing.
   case 'a': case 'A':
      Draw_Axis = !Draw_Axis;
      break;
   case 'v': case 'V':
      Draw_Vectors = !Draw_Vectors;
      break;
   case 's': case 'S':
     Draw_Stars = !Draw_Stars;
     break;

   // Adjust particle brightness.
   case 'b':
     Parts_Brightness+=0.01;
     break;
   case 'B':
     Parts_Brightness-=0.01;
     break;

   // Toggle being on-orbit.
   case 'o': case 'O':
      On_Orbit = ! On_Orbit;
      break;

   // Toggle Texture.
   case 't': case 'T':
      Texture_On = !Texture_On;
      break;

   // Impart an impulse on the Root particle.
   case 'x':
      Parts->Vel[0]+=.0001;
      break;
   case 'X':
      Parts->Vel[0]-=.0001;
      break;
   case 'c': case 'y':
      Parts->Vel[1]+=.0001;
      break;
   case 'C': case 'Y':
      Parts->Vel[1]-=.0001;
      break;
   case 'z':
      Parts->Vel[2]+=.0001;
      break;
   case 'Z':
      Parts->Vel[2]-=.0001;
      break;

   case 'r': case 'R':
      ourFireParticleGun(Parts);

   // Single-step through motion calculations.
   case 'm':
      Move_Step = 1;
      Move_Enable = 1;
      break;

   // Normal motion.
   case 'M':
      Move_Enable= !Move_Enable;
      break;

   // Heads up display.
   case 'h': case 'H':
      Heads_Up = !Heads_Up;
      break;

   // Inject one (or all) free particle(s).
   case 'i': case 'I':
      for (t=Parts_LastUnused; t<Parts_Num; t++) {
         if (!Parts[t].Running) {
            ourFireParticleGun(&Parts[t]);
            Parts_LastUnused = t;
            if (key == 'i') break;
         }
      }
      break;

   // Toggle gun control between ejection velocity and position.
   case '.':
      Gun_Control= !Gun_Control;
      break;

   // Control Particle Gun velocity vector.
   case '7':
      if (Gun_Control) Gun_VZ-=0.001; else Gun_PZ-=0.01;
      break;
   case '8':
      if (Gun_Control) Gun_VZ+=0.001; else Gun_PZ+=0.01;
      break;

   case '4':
      if (Gun_Control) Gun_VY-=0.001; else Gun_PY-=0.01;
      break;
   case '5':
      if (Gun_Control) Gun_VY+=0.001; else Gun_PY+=0.01;
      break;

   case '1':
      if (Gun_Control) Gun_VX-=0.001; else Gun_PX-=0.01;
      break;
   case '2':
      if (Gun_Control) Gun_VX+=0.001; else Gun_PX+=0.01;
      break;

   // Range of randomness which is added to each initial velocity vector.
   case '3':
      Gun_R-=0.001;
      break;
   case '6':
      Gun_R+=0.001;
      break;

   // Controls the large random shots which occur rarely. When set to 
   // 0, the system becomes highly controlled.
   case '/':
      Gun_OffR -= .001;
      break;
   case '*':
      Gun_OffR += .001;
      break;

   // Adds or removes particles to/from the system.
   case '-':
      if (Parts_Num > 100) {
         Parts_Num -= 100;
         ourAllocParticles(Parts_Num);
      }
      break;

   case '+':
      Parts_Num += 100;
      ourAllocParticles(Parts_Num);
      break;

   default:
      printf("No action for %d\n", key);

   }
}


// ------
// Callback Function called when a special key is pressed.

static void cbSpecialKeyPressed(int key, int x, int y)
{
   switch (key) {
   case GLUT_KEY_PAGE_UP:
      Obs_Dist -= 0.05f;
      break;

   case GLUT_KEY_PAGE_DOWN:
      Obs_Dist += 0.05f;
      break;

   case GLUT_KEY_LEFT:
      Obs_Angle-=2.0;
      break;

   case GLUT_KEY_RIGHT:
      Obs_Angle+=2.0;
      break;

   case GLUT_KEY_DOWN:
      Obs_Height-=0.05;
      break;

   case GLUT_KEY_UP:
      Obs_Height+=0.06;
      break;
   }

   // We don't know anything changed, but it never hurts.
   ourCalcObs();
}


// ------
// Does everything needed before losing control to the main
// OpenGL event loop

void ourInit(
   int Width,
   int Height
)
{
   ourBuildTextures();
   ourBuildStarfield(500);

   glEnable(GL_BLEND);
   glDisable(GL_ALPHA_TEST);

   // Enable flat shading -- no need for smooth.
   glShadeModel(GL_FLAT);

   // Blending mode used for fire, lit gas, etc.
   glBlendFunc(GL_SRC_ALPHA,GL_ONE);

   // Calculate the non-on-orbit observer's position.
   ourCalcObs();

   // Load up the correct perspective matrix; using a callback directly.
   cbResizeScene(Width, Height);

   if (!(Black_Hole = gluNewQuadric()))
      exit;

   // Allocate our first block of particles.
   ourAllocParticles(Parts_Num);

   // Fire off the first (Root) Particle.
   ourFireParticleGun(Parts);

}


// ------
// The main() function.  Inits OpenGL.  Calls our own init function,
// then passes control onto OpenGL.

int main(int argc,char **argv)
{
   glutInit(&argc,argv);

   // To see OpenGL drawing, take out the GLUT_DOUBLE request.
   glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH);
   glutInitWindowSize(Window_Width,Window_Height);

   // Open a window

   if (!(Window_ID=glutCreateWindow( PROGRAM_TITLE ))) {
      fprintf(stderr,"Error opening a window.\n");
      exit(-1);
   }
  
   // Register the callback function to do the drawing.
   glutDisplayFunc(&cbRenderScene);

   // If there's nothing to do, draw.
   glutIdleFunc(&cbRenderScene);

   // It's a good idea to know when our window's resized.
   glutReshapeFunc(&cbResizeScene);

   // And let's get some keyboard input.
   glutKeyboardFunc(&cbKeyPressed);
   glutSpecialFunc(&cbSpecialKeyPressed);

   // OK, OpenGL's ready to go.  Let's call our own init function.
   ourInit(Window_Width, Window_Height);

   // Print out a bit of help dialog.
   printf("\n" PROGRAM_TITLE "\n\n\
Use arrow keys to rotate around or move along Y axis.\n\
Page up/down will move observer towards/away from the Y axis.\n\n\
'O' toggles observer onto non-decaying orbit of yellow root particle.\n\
'H' toggles heads-up-display of various status variables.\n\n\
'x'/'X', 'c'/'C', 'z'/'Z' thrusts root particle along X, Y, Z axis\n\
in positive/negative directions, respectively; 'R' resets.\n\n\
Numerical Keypad controls Particle Gun parameters.  '.' key switches\n\
between effecting Gun Velocity Vector and Position.\n\n\
'+' and '-' add to or remove particles from the system.\n\n\
Use first letter of shown display mode settings to alter.\n\n\
Q or [Esc] to quit; OpenGL window must have focus for input.\n");

   // Pass off control to OpenGL.
   // Above functions are called as appropriate.
   glutMainLoop();

   // Free our allocated particle array.
   ourAllocParticles(0);

   return 1;
}