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

   Change 4:  Retain plotting information for later animation
              --- revision to Display().

   Language:  Standard C, avoiding vendor-specific features

   Author:    Timothy Rolfe
              Dakota State University
              Madison  SD  57042-1799
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.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>

/* 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        */
   } 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       */

#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,REAL);/* Make a T-h state by back-diff. */
void Derivatives(REAL,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, Fric_K;

#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 ("Friction constant:  ", stdout);
      scanf ("%lf", &Fric_K);
      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, Fric_K);
/*    Force calculation of energy terms */
      Derivatives(Scale, Fric_K);
      Display(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 */

      while (T < T_final)
      {
         Derivatives(Scale, Fric_K);
         Time_Step (h);
         T += h;
         if (++Kdx == Display_Freq)
         {
            Display(T, -1);
            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 */
      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') return 0;
   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);
   }
   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;

   randomize();

   for (Idx = Lo; Idx <= Hi; Idx++)
   {
      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, REAL Fric_K)
{
   short Idx, Jdx;

   Derivatives(Scale, Fric_K);
/* 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];
}

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;
   }
}

void Display(REAL T, int Sort_Flag)
/*
   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, Jdx;
   int Y, Z;
   static int First = 1;
   static FILE *fp = NULL;

/* Set up for the save file */
   if (First)
   {
      char Line_In[80];

      First = 0;
      fputs ("To disable output file, press just <Enter>\n", stdout);
      do
      {
         fputs ("\nOutput file name:  ", stdout);
         gets  (Line_In);
         if (Line_In[0] == '\0')
            break;
         fp = fopen (Line_In, "w");
      }  while (fp == NULL);
   }

   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];
      if (fp != NULL)
      {
         for (Jdx = 0; Jdx < 3; Jdx++)
            fprintf (fp, "%5.3f ", Point[Idx].Q[Jdx]);
         fprintf (fp, "%d\n", Idx+1);
      }
   }
   if (fp != NULL)
      fputc('\n', fp);    /* Blank line to delimit images */
   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, REAL Fric_K)
{
   short  Idx, Jdx, Kdx;
   REAL   Rhat[3], R, Force, Speed;
   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 */
      Speed = Rcalc (Origin, Point[Idx].P, Rhat);
      Tot_Kin += square (Speed) / 2;
      for (Kdx = 0; Kdx < 3; Kdx++)
/*       Frictional force is OPPOSED to the current motion */
         Point[Idx].dP[Kdx] -= Fric_K * Speed * Rhat[Kdx];

/*    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);
      Find_Eoln();
      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);
      Find_Eoln();
      return;
   }

   if ( Cross(V1, V2, Norm) == 0 )
   {
      fputs ("The points are co-linear\n", stdout);
      fputs ("Press <ENTER> to resume:  ", stdout);
      Find_Eoln();
      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(); Display(T, 0);
   for (Idx = 0; Idx < Lim; Idx++)
   {
      Angle = Angle + Phi;
      Rotate(Phi, 0, Size-1);
      Display(T, 0);
      if (Flag) continue;
      gotoxy (1,1);
      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 */
   Display(T, 0);
}

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

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

   do
   {
      clrscr();
      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(); Display(T, 0);
}