/* #define E_FILE */
/*        Dynamic simulation of points on a sphere.

   Revision:  Option of retaining plotting information
	      for later animation.  --- Accomplished via
	      change to include file Gr-Disp.c

	      Color information changed to saving an index into
	      the palette.

   Language:  Borland C++, restricted to standard C subset

   Author:    Timothy Rolfe
	      Dakota State University
	      Madison  SD  57042-1799
*/
#include <stdio.h>
#include <stdlib.h>
#include "malloc.h"   /* Borland environment --- will ref. alloc.h */
#include <ctype.h>
#include <time.h>     /* Use time() for random number seed */
#include <math.h>
#include <graphics.h>
#include <conio.h>

/* For random numbers:  Borland C++'s implementation of random ---
			adjusted for where rand() is 32 bits         */
#define MAX_RAND 0x7FFFU
#ifdef random
#undef random
#endif
#define random(num)(int)(((long)(rand()&MAX_RAND)*(num))/(MAX_RAND+1))
#define randomize()     srand((unsigned)time(NULL))

#define f_rand() ( (double) random(MAX_RAND) / (MAX_RAND+1) )

/* NOTE:  scanf calls need to have format adjusted for REAL */
/*        Also the printf calls in Display                  */

typedef double REAL;

typedef struct Point
   {
      REAL Q[3];    /* Current Cartesian position     */
      REAL P[3];    /* Current momentum (unit mass)   */
      REAL dP[3];   /* Derivative of momentum         */
      REAL L[3];    /* Last position for integrator   */
      REAL M[3];    /* Last momentum for integrator   */
      REAL X[3];    /* Last point for plotting        */
      char Color;   /* Circle color in graphic display*/
   } Cell;

typedef Cell * CellPtr;

CellPtr Point = NULL;    /* Size will be done by realloc() */

REAL  Tot_Kin, Tot_Pot;  /* Total energies */

short Size = 0, MaxSize = 0;

void Initialize (int, int);  /* Passed the range for Point[]     */
void Read_Config (void);     /* Read configuration from file     */
void Display (REAL, int);    /* Displays up to global Size       */
void Gr_Disp (REAL, int);    /* Graphical display --- circles    */

#include "ANSIterm.c"        /* clrscr, gotoxy, etc.             */

void Rotate_X(REAL,int,int); /* Rotate YZ plane about the X axis */
void Rotate_Y(REAL,int,int); /* Rotate XZ plane about the Y axis */
void Rotate_Z(REAL,int,int); /* Rotate XY plane about the Z axis */

void Back_Prop(REAL,REAL);   /* Generate T-h state by back-diff. */
void Zero_Moment(REAL, REAL);/* Zero out CONSISTENTLY            */
void Derivatives(REAL);      /* Pass potential energy constant   */
void Time_Step(REAL);        /* Take one step of size h          */

void User_Rotate (REAL);     /* Under user control, rotate system*/
			     /* Calls to Display() need T        */

void Find_Eoln (void);       /* Discard char.s up to \n          */

main (int argc, char *argv[])
{
   short Idx, Jdx, Kdx, Display_Freq;
   char  Ans[10];
   FILE *fptr;
   REAL  T = 0, T_final;
   REAL  h, Scale, Old_Kin;

#ifdef E_FILE
   fptr = fopen ("RUNTRACE.PRN", "a");
   if (fptr == NULL)
   {
      perror ("Failed to open the trace file");
      exit(1);
   }
#endif
   clrscr();

   if (argc < 2)
   {
      fputs ("Number of points:  ", stdout);
      scanf ("%hd", &Size);
      Find_Eoln();
   }
   else
   {
      sscanf (argv[1], "%hd", &Size);
      printf ("Number of points:  %d\n", Size);
   }

   Point = (CellPtr) realloc (Point, Size * (sizeof *Point));
   if (!Point)
   {
      puts ("Insufficient memory to allocate the array.");
      exit (61);
   }
   Initialize (MaxSize, Size-1);
   MaxSize = Size;

   fputs ("Initialize positions from file?  (Y/N) ", stdout);
   gets  (Ans);
   if (toupper(Ans[0]) == 'Y')
      Read_Config();

   clrscr();
   do
   {
      printf("Current T:  %2.2f\n", T);
      fputs ("T-final:  ", stdout);
      scanf ("%lf", &T_final);
      Find_Eoln();
      fputs ("Time step:  ", stdout);
      scanf ("%lf", &h);
      Find_Eoln();
      fputs ("Force scale factor:  ", stdout);
      scanf ("%lf", &Scale);
      Find_Eoln();
      fputs ("Display every ___ time steps:  ", stdout);
      scanf ("%hd", &Display_Freq);
      Find_Eoln();

      clrscr();
/*    Insure that the T-h image is consistent with the T image */
      Back_Prop(h, Scale);
/*    Force calculation of energy terms */
      Derivatives(Scale);
      Gr_Disp(T, -1);
#ifdef E_FILE
      fprintf (fptr, "%3f,%8g,%8g,%8g\n", T, Tot_Kin, Tot_Pot,
					  Tot_Kin + Tot_Pot);
#endif

      Kdx = 0;               /* Counts steps since last Display */

      Old_Kin = 0;
      while (T < T_final)
      {
	 if (Tot_Kin < Old_Kin)
	 {
	    Zero_Moment(h, Scale);
	 }
	 Old_Kin = Tot_Kin;
	 Derivatives(Scale);
	 Time_Step (h);
	 T += h;
	 if (++Kdx == Display_Freq)
	 {
/*          Display(T, -1);  */
	    Gr_Disp(T, 0);
	    Kdx = 0;
#ifdef E_FILE
	    fprintf (fptr, "%3f,%8g,%8g,%8g\n", T, Tot_Kin, Tot_Pot,
						Tot_Kin + Tot_Pot);
#endif
	 }
      }
/*    Display(T, -1);*/ /* Force the last display regardless */
      Gr_Disp(T, 0);
      do
      {
	 gotoxy (1,1);
	 fputs ("Options:\n", stdout);
	 fputs ("  0:  exit\n", stdout);
	 fputs ("  1:  Adjust constants\n", stdout);
	 fputs ("  2:  Adjust size as well\n", stdout);
	 fputs ("  3:  Restart, random configuration\n", stdout);
	 fputs ("  4:  System rotation dialog\n", stdout);
	 fputs ("\nYour choice:  ", stdout);
	 scanf ("%hd", &Idx);
	 Find_Eoln();
	 if (Idx == 4)
	    User_Rotate(T);
      }  while (Idx < 0 || Idx > 3);
      if (Idx > 1)
      {
	 printf ("Current size:  %d\n", Size);
	 fputs  ("Desired size:  ", stdout);
	 scanf  ("%hd", &Size);
	 Find_Eoln();
	 if (Size > MaxSize)
	 {
	    Point = (CellPtr) realloc (Point, Size * (sizeof *Point));
	    if (!Point)
	    {
	       puts ("Insufficient memory to allocate the array.");
	       exit (61);
	    }
	    Initialize (MaxSize, Size-1);
	    MaxSize = Size;
	 }
	 if (Idx == 3)
	 {
	    Initialize (0, Size-1);
	    T = 0;
	    fputs ("Initialize positions from file?  (Y/N) ", stdout);
	    gets  (Ans);
	    if (toupper(Ans[0]) == 'Y')
	       Read_Config();
	 }
      }
   } while (Idx > 0);

#ifdef E_FILE
   fclose(fptr);
#endif
   gotoxy (1,20);
   fputs ("                       ", stdout);
   gotoxy (1,21);
   fputs ("                       ", stdout);
   gotoxy (1,23);
   fputs ("                       ", stdout);
   gotoxy (1,22);
   fputs ("Save result to a file?  (Y/N) ", stdout);
   do
      Ans[0]  = getchar ();
   while (isspace(Ans[0]));
   gets (Ans+1);
   if (toupper(Ans[0]) == 'Y')
   {
      gotoxy (1,22);
      fputs ("                                        ", stdout);
      gotoxy (1,22);
      fputs ("File name:  ", stdout);
      gets  (Ans);
      fptr = fopen (Ans, "w");
      if (!fptr)
      {
         perror ("File open failed");
      }
      else
      {
         for (Idx = 0; Idx < Size; Idx++)
         {
            for (Jdx = 0; Jdx < 3; Jdx++)
      	    fprintf (fptr, "%15.8g", Point[Idx].Q[Jdx]);
      	    fputc ('\n', fptr);
         }
         fclose (fptr);
         fputs ("File save succeeded\n", stdout);
      }
   }
   closegraph();
   return 0;
}

void Find_Eoln (void)
{
   int CharIn;

   do
      CharIn = getchar();
   while (CharIn != '\n' && CharIn != EOF);
}

void Initialize (int Lo, int Hi)
/*
   Initialize a range of points:  random position on the sphere
   and initially unmoving --- dX = 0.  Last position for the
   integrator is set to the origin as a flag that back-propagation
   will be required.
*/
{
   int   Idx;
   REAL X, Y, Z, Scale;
   static int First = 1, Mono_Chrome;

   if (First)
   {
      First = 0;
      fputs ("Color (Y/N)? ", stdout);
      Idx = toupper(getchar());
      Mono_Chrome = (Idx != 'Y');
      while ( (Idx = getchar()) != '\n')
	 if (Idx == EOF) break;
   }

   randomize();

   for (Idx = Lo; Idx <= Hi; Idx++)
   {
      if (Mono_Chrome)
	 Point[Idx].Color = -Idx-1; /* Gr_Disp() will generate white */
      else
	 Point[Idx].Color = Idx+1;  /* color --- resolved by Gr_Disp */

      do  /* Choose random positions in the "unit" box, but reject   */
	  /* any positions not within the unit circle.               */
      {  X = 2 * f_rand() - 1;
	 Y = 2 * f_rand() - 1;
	 Z = 2 * f_rand() - 1;
	 Scale = X*X + Y*Y + Z*Z;  /* Right now, length squared */
      }  while (Scale > 1);

      Scale = 1.0 / sqrt (Scale);  /* Scale factor is 1/length  */

      Point[Idx].Q[0] = Point[Idx].X[0] = X * Scale;
      Point[Idx].L[0] = 0;
      Point[Idx].P[0] = Point[Idx].M[0] = Point[Idx].dP[0] = 0;
      Point[Idx].Q[1] = Point[Idx].X[1] = Y * Scale;
      Point[Idx].L[1] = 0;
      Point[Idx].P[1] = Point[Idx].M[1] = Point[Idx].dP[1] =  0;
      Point[Idx].Q[2] = Point[Idx].X[2] = Z * Scale;
      Point[Idx].L[2] = 0;
      Point[Idx].P[2] = Point[Idx].M[2] = Point[Idx].dP[2] =  0;
   }
}

void Read_Config (void)
{
   short Idx, Jdx;
   char  Ans[40];
   FILE *fptr;

   fputs ("File name:  ", stdout);
   gets  (Ans);
   fptr = fopen (Ans, "r");
   if (fptr == 0)
   {
      fputs ("Open failed; using random positions.\n", stdout);
      fputs ("Press <ENTER> to resume.  ", stdout);
      Find_Eoln();
   }
   else
   {
      for (Idx = 0; Idx < Size; Idx++)
      {
	 for (Jdx = 0; Jdx < 3; Jdx++)
	 {
	    if (feof(fptr)) break;
	    fscanf (fptr, "%lg", &Point[Idx].Q[Jdx]);
	 }
	 if (Jdx < 3)
	    break;
      }
      if (Idx < Size)
	 printf ("Only %d points read in.\n", Idx);
      Size = Idx;
   }
}


void Back_Prop(REAL h, REAL Scale)
{
   short Idx, Jdx;

   Derivatives(Scale); /* Does not depend on previous position    */
/* Newly generated points will have T-h image unset; generate one */
   for (Idx = 0; Idx < Size; Idx++)
      if (Point[Idx].L[0] == 0
      &&  Point[Idx].L[1] == 0
      &&  Point[Idx].L[2] == 0)
	 for (Jdx = 0; Jdx < 3; Jdx++)
	    Point[Idx].M[Jdx] = Point[Idx].Q[Jdx] - h * Point[Idx].P[Jdx];
}

void Zero_Moment(REAL h, REAL Scale)
{
   short Jdx, Kdx;

   for (Jdx = 0; Jdx < Size; Jdx++)
   {
      for (Kdx = 0; Kdx < 3; Kdx++)
      {
	 Point[Jdx].P[Kdx] = 0;
	 Point[Jdx].M[Kdx] = 0;
	 Point[Jdx].L[Kdx] = 0;
      }
   }
/* Generate a phony history consistent with present setup */
   Back_Prop(h, Scale);
}

int Compare(CellPtr Left, CellPtr Right)
/*
   Compare two cells --- based on the X component of the current
   Cartesian position (.Q[0]).
*/
{
   if (Left->Q[0] < Right->Q[0])
      return -1;
   else if (Left->Q[0] == Right->Q[0])
      return 0;
   else
      return 1;
}

void InSort (Cell X[], short Size)
/*
   Standard insertion sort.  Note the use of an external Compare
   function to hide what it means to compare two Cells.
*/
{
   short Lim, Hole, Test;
   Cell  Safe;

   for (Lim = 1; Lim < Size; Lim++)
   {
      Hole = Lim;
      Safe = X[Hole];
      Test = Hole - 1;
      while (Test >= 0 && Compare(&X[Test], &Safe) > 0)
      {
	 X[Hole] = X[Test];
	 Hole = Test;
	 Test = Hole - 1;
      }
      X[Hole] = Safe;
   }
}

#include "Gr-Disp.c"

void Display(REAL T, int Sort_Flag)
/*
   Character-graphics display module

   Assume a screen of 23 lines by 80 columns, using all lines and
   the rightmost 59 columns.

   Display will be of the y-z plane, with X information used to determine
   bright (positive x) or dim (negative x) in the display.

   To insure that two points differing only in X-coordinate are displayed
   correctly, the array is first SORTED based on the .X value.  A simple
   insertion sort is adequate:  size is probably less than 20, and the
   are will also mostly stay ordered.

   Previous position is blanked, then new position is plotted ---
   and new position is plotted a second time, should any point have been
   masked by an erase for a later point.
*/
{
   int Idx;
   int Y, Z;

   if (Sort_Flag)
      InSort (Point, Size);

/* First, the blank/plot pass:  */
   for (Idx = 0; Idx < Size; Idx++)
   {
      Y = 50 + (int) (29 * Point[Idx].X[1]);
/*    Note:  row designations are reversed compared to Y coordinates */
      Z = 12 - (int) (11 * Point[Idx].X[2]);
      lowvideo();          /* Blank is low video */
      gotoxy (Y, Z);
      putchar(' ');
      Y = 50 + (int) (29 * Point[Idx].Q[1]);
      Z = 12 - (int) (11 * Point[Idx].Q[2]);
      gotoxy (Y, Z);
      if (Point[Idx].Q[0] < 0)
	 lowvideo();
      else
	 highvideo();
      putchar(Idx < 10 ? '0' + Idx : 'a' + Idx - 10);
   }
/* Then replot the new positions */
   for (Idx = 0; Idx < Size; Idx++)
   {
      Y = 50 + (int) (29 * Point[Idx].Q[1]);
      Z = 12 - (int) (11 * Point[Idx].Q[2]);
      gotoxy (Y, Z);
      if (Point[Idx].Q[0] < 0)
	 lowvideo();
      else
	 highvideo();
      putchar(Idx < 10 ? '0' + Idx : 'a' + Idx - 10);
/*    Update "previous" position for the next plot */
      Point[Idx].X[0] = Point[Idx].Q[0];
      Point[Idx].X[1] = Point[Idx].Q[1];
      Point[Idx].X[2] = Point[Idx].Q[2];
   }
   normvideo();
   gotoxy (1,20);
   printf ("     Time:  %3.3lf   \n", T);
   printf ("  Kinetic:  %7lg   \n", Tot_Kin);
   printf ("Potential:  %7lg   \n", Tot_Pot);
   printf ("    Total:  %7lg   \n", Tot_Kin + Tot_Pot);
}

/* Note: The rotations are of the body with respect to the axes.  */
/*       Goldstein's formulas for Euler rotations are rotations   */
/*       of the axes with respect to the body (exactly reversed). */

void Rotate_X (REAL Angle, int Low, int High)
/*
   Rotate points between Low and High by Angle (in radians)
   about the X axis.
*/
{
   double  Sine, Cosine, Y, Z;
   short  Idx;

   Cosine = cos(Angle);
   Sine = sin(Angle);

   for (Idx = Low; Idx <= High; Idx++)
   {
      Y = Cosine * Point[Idx].Q[1] - Sine * Point[Idx].Q[2];
      Z = Sine * Point[Idx].Q[1] + Cosine * Point[Idx].Q[2];
      Point[Idx].Q[1] = Y;
      Point[Idx].Q[2] = Z;

      Y = Cosine * Point[Idx].L[1] - Sine * Point[Idx].L[2];
      Z = Sine * Point[Idx].L[1] + Cosine * Point[Idx].L[2];
      Point[Idx].L[1] = Y;
      Point[Idx].L[2] = Z;

      Y = Cosine * Point[Idx].P[1] - Sine * Point[Idx].P[2];
      Z = Sine * Point[Idx].P[1] + Cosine * Point[Idx].P[2];
      Point[Idx].P[1] = Y;
      Point[Idx].P[2] = Z;

      Y = Cosine * Point[Idx].M[1] - Sine * Point[Idx].M[2];
      Z = Sine * Point[Idx].M[1] + Cosine * Point[Idx].M[2];
      Point[Idx].M[1] = Y;
      Point[Idx].M[2] = Z;

      Y = Cosine * Point[Idx].dP[1] - Sine * Point[Idx].dP[2];
      Z = Sine * Point[Idx].dP[1] + Cosine * Point[Idx].dP[2];
      Point[Idx].dP[1] = Y;
      Point[Idx].dP[2] = Z;
   }
}

void Rotate_Y (REAL Angle, int Low, int High)
/*
   Rotate points between Low and High by Angle (in radians)
   about the Y axis.
*/
{
   double  Sine, Cosine, X, Z;
   short  Idx;

   Cosine = cos(Angle);
   Sine = sin(Angle);

   for (Idx = Low; Idx <= High; Idx++)
   {
      Z = Cosine * Point[Idx].Q[2] - Sine * Point[Idx].Q[0];
      X = Sine * Point[Idx].Q[2] + Cosine * Point[Idx].Q[0];
      Point[Idx].Q[2] = Z;
      Point[Idx].Q[0] = X;

      X = Sine * Point[Idx].L[2] + Cosine * Point[Idx].L[0];
      Z = Cosine * Point[Idx].L[2] - Sine * Point[Idx].L[0];
      Point[Idx].L[2] = Z;
      Point[Idx].L[0] = X;

      X = Sine * Point[Idx].P[2] + Cosine * Point[Idx].P[0];
      Z = Cosine * Point[Idx].P[2] - Sine * Point[Idx].P[0];
      Point[Idx].P[2] = Z;
      Point[Idx].P[0] = X;

      X = Sine * Point[Idx].M[2] + Cosine * Point[Idx].M[0];
      Z = Cosine * Point[Idx].M[2] - Sine * Point[Idx].M[0];
      Point[Idx].M[2] = Z;
      Point[Idx].M[0] = X;

      X = Sine * Point[Idx].dP[2] + Cosine * Point[Idx].dP[0];
      Z = Cosine * Point[Idx].dP[2] - Sine * Point[Idx].dP[0];
      Point[Idx].dP[2] = Z;
      Point[Idx].dP[0] = X;
   }
}

void Rotate_Z (REAL Angle, int Low, int High)
/*
   Rotate points between Low and High by Angle (in radians)
   about the Z axis.
*/
{
   double Sine, Cosine;
   double X, Y;
   short  Idx;

   Cosine = cos(Angle);
   Sine = sin(Angle);

   for (Idx = Low; Idx <= High; Idx++)
   {
      X = Cosine * Point[Idx].Q[0] - Sine * Point[Idx].Q[1];
      Y = Sine * Point[Idx].Q[0] + Cosine * Point[Idx].Q[1];
      Point[Idx].Q[0] = X;
      Point[Idx].Q[1] = Y;

      X = Cosine * Point[Idx].L[0] - Sine * Point[Idx].L[1];
      Y = Sine * Point[Idx].L[0] + Cosine * Point[Idx].L[1];
      Point[Idx].L[0] = X;
      Point[Idx].L[1] = Y;

      X = Cosine * Point[Idx].P[0] - Sine * Point[Idx].P[1];
      Y = Sine * Point[Idx].P[0] + Cosine * Point[Idx].P[1];
      Point[Idx].P[0] = X;
      Point[Idx].P[1] = Y;

      X = Cosine * Point[Idx].M[0] - Sine * Point[Idx].M[1];
      Y = Sine * Point[Idx].M[0] + Cosine * Point[Idx].M[1];
      Point[Idx].M[0] = X;
      Point[Idx].M[1] = Y;

      X = Cosine * Point[Idx].dP[0] - Sine * Point[Idx].dP[1];
      Y = Sine * Point[Idx].dP[0] + Cosine * Point[Idx].dP[1];
      Point[Idx].dP[0] = X;
      Point[Idx].dP[1] = Y;
   }
}

REAL square(REAL X)
{
   return X * X;
}

REAL Rcalc (REAL X1[], REAL X2[], REAL Rhat[])
{
   short Idx;
   REAL  R_squared, R, Work[3];    /* Allow for aliasing of Rhat */

   for (Idx = R_squared = 0; Idx < 3; Idx++)
   {
      Work[Idx] = X2[Idx] - X1[Idx];
      R_squared += square(Work[Idx]);
   }
   R = sqrt(R_squared);
   for (Idx = 0; Idx < 3; Idx++)
      Rhat[Idx] = R != 0 ? Work[Idx] / R : 0;
   return R;
}

REAL Cross (REAL A[], REAL B[], REAL C[])
{
   REAL  N[3], Sum;          /* Allow for aliasing */
   short Idx, Jdx, Kdx;

   for (Idx = 0; Idx < 3; Idx++)
   {
      Jdx = (Idx + 1) % 3;
      Kdx = (Jdx + 1) % 3;
      N[Idx] = A[Jdx]*B[Kdx] - A[Kdx]*B[Jdx];
   }

   for (Idx = Sum = 0; Idx < 3; Idx++)
     Sum = square( C[Idx] = N[Idx] );
   return sqrt(Sum);
}

void Derivatives(REAL Scale)
{
   short  Idx, Jdx, Kdx;
   REAL   Rhat[3], R, Force;
   static REAL Origin[3] = {0, 0, 0};

/* First:  zero out the derivative components */
   for (Idx = 0; Idx < Size; Idx++)
      Point[Idx].dP[0] = Point[Idx].dP[1] = Point[Idx].dP[2] = 0;

   Tot_Pot = Tot_Kin = 0;

   for (Idx = 0; Idx < Size; Idx++)
   {
/*    Quantities based on Idx alone:  friction and kinetic energy */
      Tot_Kin += square (Rcalc (Origin, Point[Idx].P, Rhat)) / 2;

/*    Pair-wise quantities */
      for (Jdx = Idx + 1; Jdx < Size; Jdx++)
      {
	 R = Rcalc(Point[Jdx].Q, Point[Idx].Q, Rhat);
/*       Potential energy accumulation  */
	 Tot_Pot += Scale / R;
/*       Force calculation:  */
	 Force = Scale / square(R);
/*       As REPULSIVE, the force on Idx is along the vector and */
/*       the force on Jdx is opposite to the vector.            */
	 for (Kdx = 0; Kdx < 3; Kdx++)
	 {
	    Point[Idx].dP[Kdx] += Force * Rhat[Kdx];
	    Point[Jdx].dP[Kdx] -= Force * Rhat[Kdx];
	 }
      }
   }
}

void Time_Step (REAL h)
/*
   For each point, rotate it into the XZ plane (angle phi), then
   rotate it to the Z axis (angle theta).

   Discard the Z component to the force, then apply a central-differences
   time-step.  Fine-tune the resulting position to lie on the unit sphere.

   Finally, reverse the theta and phi rotations.
*/
{
   REAL  Phi, Theta, R, X, Y, Z;
   short Idx;

   for (Idx = 0; Idx < Size; Idx++)
   {
/*    Rotation onto the Z axis */
      Phi = atan2 (Point[Idx].Q[1], Point[Idx].Q[0]);
      Rotate_Z (-Phi, Idx, Idx);
      Theta = atan2(Point[Idx].Q[0], Point[Idx].Q[2]);
      Rotate_Y (-Theta, Idx, Idx);

/*    First, update position based on momentum (P = dQ)  */
      X = Point[Idx].L[0] + 2*h * Point[Idx].P[0];
      Y = Point[Idx].L[1] + 2*h * Point[Idx].P[1];
      Z = Point[Idx].Q[2];   /* I.e., NO change --- and should still be 1 */
      Point[Idx].P[2] = 0;
      Point[Idx].L[0] = Point[Idx].Q[0];
      Point[Idx].L[1] = Point[Idx].Q[1];
      Point[Idx].L[2] = Point[Idx].Q[2];
/*    Force the new position onto the unit sphere.  */
      R = sqrt (X*X + Y*Y + Z*Z);
      Point[Idx].Q[0] = X/R;
      Point[Idx].Q[1] = Y/R;
      Point[Idx].Q[2] = Z/R;

/*    Then update momentum based on force --- dP, given unit mass */
      X = Point[Idx].M[0] + 2*h * Point[Idx].dP[0];
      Y = Point[Idx].M[1] + 2*h * Point[Idx].dP[1];
      Z = Point[Idx].P[2];   /* I.e., NO change */
      Point[Idx].dP[2] = 0;
      Point[Idx].M[0] = Point[Idx].P[0];
      Point[Idx].M[1] = Point[Idx].P[1];
      Point[Idx].M[2] = Point[Idx].P[2];
      Point[Idx].P[0] = X;
      Point[Idx].P[1] = Y;
      Point[Idx].P[2] = Z;

/*    Undo the earlier rotation */
      Rotate_Y (Theta, Idx, Idx);
      Rotate_Z (Phi, Idx, Idx);
      R = 0;   /* Code for "Run to Cursor" in Borland */
   }
}

REAL Pi = 0;   /* Several functions require Pi, so make it global */

void To_Z (REAL Q[], REAL T)
{
   REAL Phi, Theta, Phi2, R;
   int  Idx2;
   char Opt2;

/* Final point of dialog --- which item in the XZ plane? */
   gotoxy (1,18);
   fputs  ("Which point to XZ plane? ", stdout);
   scanf  ("%d", &Idx2);
   Find_Eoln();
   gotoxy (1,18);
   fputs  ("                                  ", stdout);

/* If the point is EXACTLY on the Z axis, atan2() will have      */
/* indigestion --- but then Phi might as well be zero anyway.    */
   Phi = (Q[0] == 0 && Q[1] == 0) ? Phi = 0 : atan2 (Q[1], Q[0]);
   R = sqrt (Q[0]*Q[0] + Q[1]*Q[1] + Q[2]*Q[2]);
   Theta = acos (Q[2] / R);

/* Note: The rotations are of the body with respect to the axes.  */
/*       Goldstein's formulas for Euler rotations are rotations   */
/*       of the axes with respect to the body (exactly reversed). */
   Rotate_Z (-Phi, 0, Size-1);
   Rotate_Y (-Theta, 0, Size-1);
/* If the point is EXACTLY on the Z axis, atan2() will have      */
/* indigestion --- but then Phi2 might as well be zero anyway.   */
   Phi2 = (Point[Idx2].Q[0] == 0 && Point[Idx2].Q[1] == 0) ?
	  Phi2 = 0 :
	  atan2 (Point[Idx2].Q[1], Point[Idx2].Q[0]);
   Rotate_Z (-Phi2, 0, Size-1);

   Display(T, 0);
   gotoxy (1,18); printf ("  Phi = %7lg\n", -180*Phi/Pi);
		  printf ("Theta = %7lg\n", -180*Theta/Pi);
		  printf ("Phi-2 = %7lg\n", -180*Phi2/Pi);
   gotoxy (1,16);
   fputs ("Retain?  (Y/N) ", stdout);
   do
      scanf ("%c", &Opt2);
   while (isspace(Opt2));
   if (toupper(Opt2) != 'Y')
   {
      Rotate_Z (Phi2, 0, Size-1);
      Rotate_Y (Theta, 0, Size-1);
      Rotate_Z (Phi, 0, Size-1);
      Display(T, 0);
   }
}

void To_YZ (REAL Q[], REAL T)
{
   REAL Phi;
   char Opt2;

/* If the point is EXACTLY on the Z axis, atan2() will have      */
/* indigestion --- but then Phi might as well be zero anyway.    */
   Phi = (Q[0] == 0 && Q[1] == 0) ? Phi = 0 : atan2 (Q[1], Q[0]);
/* -Phi will take us to the XZ plan, then Pi/2 takes us to YZ    */
   Rotate_Z (Pi/2 - Phi, 0, Size-1);
   Display(T, 0);
   gotoxy (1,18); printf ("  Phi = %7lg\n", -180*Phi/Pi);
   gotoxy (1,16);
   fputs ("Retain?  (Y/N) ", stdout);
   do
      scanf ("%c", &Opt2);
   while (isspace(Opt2));
   if (toupper(Opt2) != 'Y')
   {
      Rotate_Z (Phi - Pi/2, 0, Size-1);
      Display(T, 0);
   }
}

void Surface_Normal (REAL T)
{
   REAL  V1[3], V2[3], Norm[3];
   short P1, P2, P3;
   char  CharIn;

   gotoxy (1,20); fputs("                        ", stdout);
   gotoxy (1,21); fputs("                        ", stdout);
   gotoxy (1,22); fputs("                        ", stdout);
   gotoxy (1,23); fputs("                        ", stdout);

   gotoxy (1,20);
   fputs ("Index of 1st:  ", stdout);
   scanf ("%d", &P1);
   Find_Eoln();
   gotoxy (1,21);
   fputs ("Index of 2nd:  ", stdout);
   scanf ("%d", &P2);
   Find_Eoln();
   gotoxy (1,22);
   fputs ("Index of 3rd:  ", stdout);
   scanf ("%d", &P3);
   Find_Eoln();

   if ( Rcalc(Point[P1].Q, Point[P2].Q, V1) == 0 )
   {
      fputs ("1st and 2nd points are identical\n", stdout);
      fputs ("Press <ENTER> to resume:  ", stdout);
      do
	 CharIn = getchar();
      while (CharIn != '\n');
      return;
   }

   if ( Rcalc(Point[P2].Q, Point[P3].Q, V2) == 0 )
   {
      fputs ("2nd and 3rd points are identical\n", stdout);
      fputs ("Press <ENTER> to resume:  ", stdout);
      do
	 CharIn = getchar();
      while (CharIn != '\n');
      return;
   }

   if ( Cross(V1, V2, Norm) == 0 )
   {
      fputs ("The points are co-linear\n", stdout);
      fputs ("Press <ENTER> to resume:  ", stdout);
      do
	 CharIn = getchar();
      while (CharIn != '\n');
      return;
   }

   To_Z (Norm, T);
}

void Rotate_One( void Rotate (REAL, int, int), REAL T )
{
   REAL Angle;
   char Opt2;

   fputs ("Size of rotation in degrees: ", stdout);
   scanf ("%lf", &Angle);
   Find_Eoln();
   Angle = Pi * Angle / 180;
   Rotate (Angle, 0, Size-1);
   Display (T, 0);
   gotoxy (1,16);
   fputs ("Retain?  (Y/N) ", stdout);
   do
      scanf ("%c", &Opt2);
   while (isspace(Opt2));
   if (toupper(Opt2) != 'Y')
   {
      Rotate (-Angle, 0, Size-1);
   }
}

void Rotate_Full ( void Rotate (REAL, int, int), REAL T )
{
   REAL  Phi, Angle;
   short Idx, Lim;
   char  Flag, CharIn;

   fputs ("How many steps for the full rotation?  ", stdout);
   scanf ("%hd", &Lim);
   Find_Eoln();
   fputs ("Wait for <ENTER> at each step?  (Y/N)  ", stdout);
   do
      CharIn = getchar();
   while (isspace(CharIn));
   Flag = toupper(CharIn) != 'Y';
   Phi =  2 * Pi / Lim;  Angle = 0;
   clrscr(); Gr_Disp(T, 1);
   for (Idx = 0; Idx < Lim; Idx++)
   {
      Angle = Angle + Phi;
      Rotate(Phi, 0, Size-1);
      Gr_Disp(T, 0);
      if (Flag) continue;
      do
      {
	 CharIn = getchar();
	 if (toupper(CharIn) == 'Q')   /* Early exit request */
	    Flag = 2;
      }  while (CharIn != '\n');
      if (Flag == 2)
	 break;
   }
   Rotate(-Angle, 0, Size-1);  /* Undo, especially on early exit */
   closegraph();
   Display(T, 0);
}

void User_Rotate(REAL T)
{
   short Option = 0, Idx;
   char  Opt2;

   if (Pi == 0)
      Pi = 4 * atan(1.0);

   do
   {
      clrscr();

      closegraph();  /* Exit graphical mode for dialog */

      Display(T, 0);
      gotoxy (1,1);
      fputs ("Options:\n", stdout);
      fputs ("  0:  exit\n", stdout);
      fputs ("  1:  move particle\n", stdout);
      fputs ("      to Z axis\n", stdout);
      fputs ("  2:  move particle\n", stdout);
      fputs ("      into the YZ plane\n", stdout);
      fputs ("  3:  set surface (3 points)\n", stdout);
      fputs ("      normal to Z axis\n", stdout);
      fputs ("  4:  single rotation\n", stdout);
      fputs ("  5:  Full rotation\n", stdout);
      fputs ("\nYour choice:  ", stdout);
      scanf ("%d", &Option);
      Find_Eoln();
      switch (Option)
      {
	 case 0:  break;

	 case 1:  fputs ("\nIndex:  ", stdout);
		  scanf ("%d", &Idx);
		  Find_Eoln();
		  To_Z (Point[Idx].Q, T);
		  break;
	 case 2:  fputs ("\nIndex:  ", stdout);
		  scanf ("%d", &Idx);
		  Find_Eoln();
		  To_YZ (Point[Idx].Q, T);
		  break;
	 case 3:  Surface_Normal(T);
		  break;
	 case 4:  fputs ("\nAbout which axis? (X, Y, Z) ", stdout);
		  do
		     scanf ("%c", &Opt2);
		  while (isspace(Opt2));
		  switch (toupper(Opt2))
		  {
		     case 'X': Rotate_One (Rotate_X, T); break;
		     case 'Y': Rotate_One (Rotate_Y, T); break;
		     case 'Z': Rotate_One (Rotate_Z, T); break;
		     default:  ;  /* I.e., ignore! */
		  }
		  break;
	 case 5:  fputs ("\nAbout which axis? (X, Y, Z) ", stdout);
		  do
		     scanf ("%c", &Opt2);
		  while (isspace(Opt2));
		  Opt2 = toupper(Opt2);
		  switch (toupper(Opt2))
		  {
		     case 'X': Rotate_Full(Rotate_X, T); break;
		     case 'Y': Rotate_Full(Rotate_Y, T); break;
		     case 'Z': Rotate_Full(Rotate_Z, T); break;
		     default:  ;  /* I.e., ignore! */
		  }
		  break;

	 default: printf ("Unrecognized option (%d)\n", Option);
      }
   } while (Option > 0);
   clrscr();
   Gr_Disp(T, -1);   /* Back to graphical display */
}