//               INTERGALACTIC GRAVITATIONAL MOTION SIMULATOR
//
//                    v3.2     Last Updated 16 April 2000
//                         Programmed by D. Perley
//          http://www.people.cornell.edu/pages/dap29/programs.html
//
//   You are free to perform any modifications to this code that you desire,
//   provided you retain my name above and credit me where appropriate.
//   The code may be redistributed free of restriction.
//
//   All components of the software - interface, 3D projection, and n-body
//   simulator are included below.  Only the graphics file EGAVGA.BGI and
//   the Borland C/C++ 3.1 header/library files are necessary outside this
//   file to recompile the program.

#include <stdio.h>
#include <graphics.h>
#include <stdlib.h>
#include <conio.h>
#include <math.h>
#include <ctype.h>
#include <alloc.h>
#include <string.h>
#include <process.h>

#define Gphys 6.673e-11          //universal gravitational constant
#define maxp 502                 //maximum particles supportable (lim 968)
#define sccx 320                 //x-coordinate of center of screen
#define sccy 240                 //y-coordinate of center of screen
#define zinc 1.5                 //zoom step factor
#define sinc 1.5                 //shift increment
#define rinc 5                   //rotation increment angle
#define dtor 0.0174532925        //conversion: degrees to radians

void simmodule(void);            //main module for simulator
void corecalc(void);             //calculates acceleration and moves particles
void input(void);                //finds user input and routes appropriately
void graphiccalc(int p1);        //calculates graphical gx and gy for p1
void drawparticle(int p1, int style, int cset);  //draws/erases specified p
void drawallparticles(int cset1, int cset2);//draws/erases/redraws all p
void pause(void);                //holds simulation during user pause
void imagezoom(int zdir);        //magnifies or minifies view space
void imageshift(int sdir);       //shifts image horizontally
void changedelay(int ddir);      //increases or decreases delay
void displayinfo(void);          //displays simulation statistics
void clearinfo(void);            //clears displayinfo() display
void printmode(void);            //special color/style for using PrintScrn
void imagerotate(int rdir);      //rotates image
void initialize(void);           //sets default conditions and init. grapics
void initgalaxies(void);         //creates galaxy disks
void transform(float rxy, float rxz, float ryz); //transforms over 3 rotations
void retroproject(void);         //find initial galaxy positions in rev. time
void snapshot(void);             //records visual graphic data to disk
void storepos(void);             //store particle positions in dynamic array
void storeinitial(void);         //store essential initial conditions to disk
void storesimulation(void);      //store simulation to disk
void loadsimulation(void);       //load initstats and simulation details to disk
void findpos(void);              //find current particle positions from databank
void initmemory(void);           //initializes dynamic memory
void terminate(void);            //terminates active simulation
int evendiv(float num, float den); //determines whether numbers divide evenly
void IFmodule(void);             //main module for interface
void IFprintheader(void);        //print interface header
void IFgetsellabel(char * selptr, int l); //retrieve selection label
int IFinput(void);               //find user input in interface
void IFchangesel(int dir, int l);//change selection in interface
void liminput(char * inputstring, int lim, int digok, int alpok, int buf);
void IFproperties(void);         //module for interface properties page
void IFgalaxies(void);           //module for interface galaxies page
int checkretro(void);            //check retroprojecton parameters for errors
int checkdisk(void);             //check galaxy particle parameters for errors
int checkprop(void);             //check property parameters for errors
void errmsg(char * errormessage);//get retroprojection error message
void printfunct(int f1, int f2, int f3, int f4, int f5); //print funct. footer
double fix(double rnum, int ndig);//round after ndig decimal place
void IFparticles(void);          //module for interface particles mage
void printerrmsg(void);          //prints error message in interface
void savesettings(void);         //saves interface settings to disk
void loadsettings(int loadall);  //loades interface settings from disk
void tracking(void);             //fixes camera on object
int dopplercolor(int p1);        //assigns p color based on z-axis velocity
int fadecolor(int p1);           //assigns p color based on z-axis distance
void adjcolors(int cdir);        //adjusts settings for color functions
void kcalc(void);                //calculates acceleration equation constants
void printtime(void);            //prints simulation name and current time


float t;                   //time
float timestep = 1;        //time increment
float maxt;                //auto-termination time
int dataentry = 0;         //storage index
int nump;                  //number of particles
char com;                  //user input variable
float centx = 0;           //x-coordinate of effective center of screen
float centy = 0;           //y-coordinate of effective center of screen
float centz = 0;           //z-coordinate of effective center of screen
float offx = 0;            //x-axis distance from center and tracked particle
float offy = 0;            //y-axis distance from center and tracked particle
float rota = 0;            //x-y rotation
float rotb = 0;            //x-z rotation
float rotc = 0;            //y-z rotation
float zoom = 1;            //zoom factor
int ptrack = 0;            //tracked particle ID
int exitflag = 0;          //flags exit command
int pauseflag = 0;         //flags pause
int displayflag = 0;       //flags display info
int store = 0;             //data storage 0=none 1=RAM 2=disk
int source = 0;            //source of particle data 0=simulation 1=file
int redraw = 2;            //pixel redraw metod 0=flash 1=sequence
int monitor = 1;           //flags view simulation
int style = 0;             //particle set appearance
int track = 0;             //tracking mode 0=none 1=g1 2=g2 3=center-of-mass
int doppler = 0;           //flags doppler mimicry
int fade = 0;              //depth coloration
long delay = 0;            //steps in timing loop delay
int tracers = 1;           //flags tracers
int recfreq = 30;          //store data every recfreq stu
int drawfreq = 2;          //redraw particles every drawfreq stu
int autosnap = 0;          //automatically snapshot every autosnap stu
int color;                 //color of drawn particle
int fileindex = 0;         //snapshot file index
int skip;                  //flags skip user input
float transx0, transy0, transz0;
float transx1, transy1, transz1;   //3-D transformation step variables
float transx2, transy2, transz2;
float igdis;               //perigalacticon distance
float igvel;               //perigalacticon velocity
float igecc;               //eccentricity
float iginc;               //kinematic inclination
float igarg;               //kinematic rotation
float igt0;                //initial time (0 = perigalacticon)
double smu = 1e+39;        //value of simulation mass unit
double sdu = 1e+19;        //value of simulation distance unit
double stu = 1e+13;        //value of simulation time unit
double G = 6.673e-3;       //value of simulation gravitational constant
int sduor = 0;
int stuor = 1;             //precedence in which units are rescaled
int smuor = 2;
int gor = 3;
int check = 0;             //flags acceptible parameters
int nextpage = 0;          //flags change pages
int err = 0;               //interface error code
int errl = 0;              //interface error location
int rerr = 0;              //retrotrajectory error code
float galr, perir;         //initial radii of encounter
int runsim = 0;            //flags initiate simulation
int loadselflag = 0;       //flags reprint all selections
int pnucleus1, pnucleus2;  //particle IDs of galactic nucleus particles
float contrast = 5;        //color contrast value
float restvel = 0;         //color base value
int autoterminate;         //flags automatically terminate simulation
int pageturn;              //flags cycle through interface

float x[maxp + 1];         //x-coordinate
float y[maxp + 1];         //y-coordinate
float z[maxp + 1];         //z-coordinate
float vx[maxp + 1];        //x-axis velocity
float vy[maxp + 1];        //y-axis velocity
float vz[maxp + 1];        //z-axis velocity
float ax[maxp + 1];        //x-axis acceleration
float ay[maxp + 1];        //y-axis acceleration
float az[maxp + 1];        //z-axis acceleration
float m[maxp + 1];         //mass
int c[maxp + 1];           //pixel base color
int c1[maxp + 1];          //pixel temporary color
int c2[maxp + 1];          //tracer color
int rx[maxp + 1];          //graphical x-coordinate
int ry[maxp + 1];          //graphical y-coordinate

float gx[3];               //galaxy initial x-coordinate
float gy[3];               //galaxy initial y-coordinate
float gz[3];               //galaxy initial z-coordinate
float gvx[3];              //galaxy initial x velocity
float gvy[3];              //galaxy initial y velocity
float gvz[3];              //galaxy initial z velocity
float gm[3];               //galaxy mass
int gc[3] = {0, 10, 13};   //galaxy disk color
int gring[3];              //number of rings in disk
int grmin[3];              //radius of innermost ring
int grmax[3];              //radius of outermost ring
int gtilt[3];              //rotation of galaxy about y axis
int gskew[3];              //rotation of galaxy about x axis
int gp[3];                 //particles per ring
int gpo[3];                //relative mass density at center
int gpmax[3];              //relative mass density at edge
int grmmax[3];             //radius of massive disk

float ki[3];
float kj[3];               //constants to abbreviate acceleratory equation
float kk[3];


int numentry = 65;         //total number of times stored in simulation file
float far * sx[maxp + 1];  //points to memory storage for x coordinates
float far * sy[maxp + 1];  //points to memory storage for y coordinates
float far * sz[maxp + 1];  //points to memory storage for z coordinates


FILE * fptr;                                     //miscellaneous file pointer
char simname[7] = "UNTITL";                      //simulation name
char capsname[7];                                //simulation name (caps)
char filepath[30] = "";  //file directory path
char filename[13];                               //file name
char filelocation[45];                           //file name and location


/****************************************************************************/
int main(void)
{
   //Main module for IGMS, calls major subroutine modules

   do
   {
      IFmodule();
      initmemory();
      initialize();
      simmodule();
      terminate();
   } while (1);
}
/****************************************************************************/
void simmodule(void)
{
   //Module for gravitational simulator, cycles time and calls sim. functions

   long n;                    //miscellaneous counting variable
   int p1;                    //particle ID

   if (track) tracking();
   drawallparticles(0, 1);
   displayflag = 1;
   displayinfo();

   pause();
   clearinfo();

   while (!exitflag)
   {
      if (store == 2) storepos();                   //record

      if (!source) corecalc();                      //calculate new particle
      if (source) findpos();                        //positions

      if (autosnap > 0) { if (evendiv(t, autosnap)) snapshot(); } //autosnap

      if (!source) t += timestep; else t += recfreq * timestep;   //timestep

      for (n = 1; n < (delay * (recfreq * source + 1)); n ++);    //delay                        //delay

      if (track) tracking();                                      //track
      if ((evendiv(t / timestep, drawfreq)) || (source))
      {
	 if (redraw == 0) {drawallparticles(2 + !tracers, 0); drawallparticles(0, 1);}
	 if (redraw == 1) drawallparticles(2 + !tracers, 1);                 //redraw
	 if (redraw == 2) {drawallparticles(2 + !tracers, 1); drawallparticles(0, 1);}
	 printtime();
      }


      if (kbhit()) input();      //check for user input

      if ((store || source) && (t >= maxt)) {exitflag = 1; autoterminate = 1;}

   }

   return;
}

/****************************************************************************/
void corecalc(void)
{
   //Calculates acceleration, accelerates, and moves all particles

   int p1, p2;                //particle IDs
   int galaxy;                //galaxy ID
   float dx, dy, dz;          //distance between particles along axes
   float D;                   //resultant distance between particles
   float A;                   //resultant acceleration

   for (p1 = 1; p1 <= nump; p1 ++)
   {
      ax[p1] = ay[p1] = az[p1] = 0;
      for (p2 = 1; p2 <= nump; p2 ++)            //determine acceleration
      {
	 if ((p1 != p2) && (m[p2] != 0))
	 {
	    dx = x[p2] - x[p1];
	    dy = y[p2] - y[p1];
	    dz = z[p2] - z[p1];
	    D = sqrt(pow(dx, 2) + pow(dy, 2) + pow(dz, 2));
	    galaxy = 0;
	    if (p2 == pnucleus1) galaxy = 1; if (p2 == pnucleus2) galaxy = 2;
	    if (D < grmmax[galaxy]) A = G * gm[galaxy] * (ki[galaxy] * D + kj[galaxy]) / kk[galaxy];
	    else A = G * m[p2] / pow(D, 2);
	    ax[p1] += dx * A / D;
	    ay[p1] += dy * A / D;
	    az[p1] += dz * A / D;
	 }
      }
   }

   for (p1 = 1; p1 <= nump; p1++)
   {
      vx[p1] += ax[p1] * timestep;                            //accelerate
      vy[p1] += ay[p1] * timestep;
      vz[p1] += az[p1] * timestep;

      x[p1] += vx[p1] * timestep;                             //move
      y[p1] += vy[p1] * timestep;
      z[p1] += vz[p1] * timestep;
   }

   galr = D;
   return;
}

/****************************************************************************/
void input(void)
{
   //Finds user input and routes to appropriate function
   com = getch();
   if (com == 27) {exitflag = 1; autoterminate = 0;} else exitflag = 0;
   if (com == '-') changedelay(0);
   if (com == '=') changedelay(1);
   if (com == 'i') imagezoom(0);
   if (com == 'o') imagezoom(1);
   if (com == 72) imageshift(1);
   if (com == 75) imageshift(2);
   if (com == 77) imageshift(3);
   if (com == 80) imageshift(4);
   if (com == 'r') imageshift(0);
   if (com == ',') imagerotate(0);
   if (com == '.') imagerotate(1);
   if (com == 'l') imagerotate(2);
   if (com == ';') imagerotate(3);
   if (com == 'p') imagerotate(4);
   if (com == '[') imagerotate(5);
   if (com == 'c') {clearviewport(); displayflag = 0; drawallparticles(0, 1);}
   if ((com == ' ') && (!pauseflag)) pause();
   if (com == '/') {if (displayflag) clearinfo(); displayflag = !displayflag; pause();}
   if (com == 't') tracers = !tracers;
   if (com == 'a') printmode();
   if (com == 's') snapshot();
   if (com == 'm') adjcolors(0);
   if (com == 'n') adjcolors(1);
   if (com == 'j') adjcolors(2);
   if (com == 'k') adjcolors(3);


   return;
}
/****************************************************************************/
void drawallparticles(int cset1, int cset2)
{
   //Draws all particles according to color/reprocess instructions

   int p1;     //particle ID
   int pstyle; //appearance of particles 0=point 1=cross 2,3=circle

   for (p1 = 1; p1 <= nump; p1 ++)
   {
      if (style == 0) pstyle = 0;
      if (style == 1) pstyle = 1;
      if (style == 2) pstyle = 2;
      if (style == 3) {if ((p1 == pnucleus1) || (p1 == pnucleus2)) pstyle = 1; else pstyle = 0;}
      if (style == 4) {if ((p1 == pnucleus1) || (p1 == pnucleus2)) pstyle = 3; else pstyle = 1;}
      if (cset1) drawparticle(p1, pstyle, cset1);
      if (cset2) graphiccalc(p1);
      if (cset2) drawparticle(p1, pstyle, cset2);
   }
   return;
}
/****************************************************************************/
void graphiccalc(int p1)
{
   //Recalculates graphical coordinates and colors of particle

   transx0 = x[p1] - centx;
   transy0 = y[p1] - centy;
   transz0 = z[p1] - centz;
   transform(rota, rotb, rotc);
   rx[p1] = (int)(zoom * (transx2) + sccx + offx);
   ry[p1] = (int)(zoom * (transy2) + sccy + offy);
   if (!doppler || !fade) c1[p1] = c[p1];
   if (doppler) c1[p1] = dopplercolor(p1);
   if (fade) c1[p1] = fadecolor(p1);
   return;
}

/****************************************************************************/
void drawparticle(int p1, int pstyle, int cset)
{
   //Draws specified particle with specified style and color option

   if (cset == 1) color = c1[p1];   //draw particle
   if (cset == 2) color = c2[p1];   //erase particle with tracer
   if (cset == 3) color = 0;        //erase
   if (cset == 4) color = 15;       //black (15) particle for printmode


   if (pstyle == 0) putpixel(rx[p1], ry[p1], color);  //0 = point
   if (pstyle == 1)
   {
      putpixel (rx[p1], ry[p1], color);               //1 = cross
      putpixel (rx[p1] + 1, ry[p1], color);
      putpixel (rx[p1] - 1, ry[p1], color);
      putpixel (rx[p1], ry[p1] + 1, color);
      putpixel (rx[p1], ry[p1] - 1, color);
   }
   if ((pstyle == 2) || (pstyle == 3))                //2 = circle
   {                                                  //3 = filled circle
      setcolor(color);
      setfillstyle(1, color);
      circle(rx[p1], ry[p1], 2);
      if (pstyle == 3) floodfill(rx[p1], ry[p1], color);
   }
   return;
}

/****************************************************************************/
void pause(void)
{
   //Pauses simulation and holds input function until unpaused

   pauseflag = 1;
   do
   {
      if (displayflag) displayinfo();
      input();
   } while ((com != ' ') && (exitflag == 0));

   if ((displayflag == 1) && (exitflag == 0)) clearinfo();
   if (tracers == 2) tracers = 1;
   pauseflag = 0;
   displayflag = 0;
   return;
}

/****************************************************************************/
void imagezoom(int zdir)
{
   //Changes zoom factor

   if ((!zdir) && (zoom < 1000)) {zoom *= zinc; offx *= zinc; offy *= zinc;}
   if ((zdir) && (zoom > 0.01)) {zoom /= zinc; offx /= zinc; offy /= zinc;}
   clearviewport();
   printtime();
   drawallparticles(0, 1);
   return;
}

/****************************************************************************/
void imageshift(int sdir)
{
   //Moves camera

   if (sdir == 0) {offy = 0; offx = 0;}
   if (sdir == 1) offy -= sinc; // zoom
   if (sdir == 2) offx -= sinc; // zoom
   if (sdir == 3) offx += sinc; // zoom
   if (sdir == 4) offy += sinc; // zoom

   if (tracers == 1) {tracers = 2; clearviewport(); printtime();}

   drawallparticles(3, 1); drawallparticles(0, 1);

   if ((!pauseflag) && (tracers == 2)) tracers = 1;

   return;
}
/****************************************************************************/
void changedelay(int ddir)
{
   //Adjusts delay value

   if (ddir)
   {
      if ((delay > 0) && (delay < 400000)) delay *= 2;
      if (delay == 0) delay = 100;
   }
   else
   {
      if (delay == 100) delay = 0;
      if (delay > 100) delay /= 2;
   }
   return;
}
/****************************************************************************/
void displayinfo(void)
{
   //Prints simulator information on screen

   int n;              //miscellaneous counting variable

   for (n = 0; n < 7; n ++) capsname[n] = toupper(simname[n]);
   gotoxy(1, 1); printf(" simname: %s", capsname);
   if (timestep >= 1)
   {
      gotoxy(1, 2); printf("    time: %.0f  ", t);
      gotoxy(1, 4); printf("timestep: %.0f  ", timestep);
   }
   else
   {
      gotoxy(1, 2); printf("    time: %.3f  ", t);
      gotoxy(1, 4); printf("timestep: %.3f  ", timestep);
   }
   gotoxy(1, 3); printf("          %.1f Myr  ", t * stu / 3.1557e+13);
   gotoxy(1, 5); printf("   delay: %ld  ", delay);
   gotoxy(1, 6); printf("    zoom: %.0f\% ", zoom * 100);
   gotoxy(1, 7); printf("x offset: %.2f ", offx);
   gotoxy(1, 8); printf("y offset: %.2f ", offy);
   gotoxy(1, 9); printf("a rotate: %.0f ", rota);
   gotoxy(1, 10); printf("b rotate: %.0f ", rotb);
   gotoxy(1, 11); printf("c rotate: %.0f ", rotc);
   gotoxy(1, 12); printf("contrast: %.2f ", contrast);
   gotoxy(1, 13); printf(" restvel: %.1f ", restvel);
   return;
}
/****************************************************************************/
void clearinfo(void)
{
   //Clears previously printed displayinfo information

   int n;              //miscellaneous counting variable

   gotoxy(1,1);
   for (n = 1; n <= 13; n++) printf("                      \n");
   printtime();
   return;
}
/****************************************************************************/
void printmode(void)
{
   //Reprints all particles as black on white background for easy printing

   drawallparticles(3, 0);
   setbkcolor(15);
   drawallparticles(0, 1);
   setpalette(15, 0);
   setpalette(8, 7);
   setcolor(15);
   drawallparticles(4, 0);
   getch();
   setbkcolor(0);
   setpalette(8, 56);
   setpalette(15, 63);
   drawallparticles(3, 0);
   drawallparticles(0, 1);

   return;
}

/****************************************************************************/
void imagerotate(int rdir)
{
   //changes 3-D viewing angle

   if (tracers == 1)
   {
      tracers = 2; clearviewport();
      printtime();
   }

   if (rdir == 0) rota += rinc;
   if (rdir == 1) rota -= rinc;
   if (rdir == 2) rotb += rinc;
   if (rdir == 3) rotb -= rinc;
   if (rdir == 4) rotc += rinc;
   if (rdir == 5) rotc -= rinc;

   if (rota < -179) rota = 180;
   if (rotb < -179) rotb = 180;
   if (rotc < -179) rotc = 180;

   if (rota > 184) rota = -175;
   if (rotb > 184) rotb = -175;
   if (rotc > 184) rotc = -175;

   drawallparticles(3, 0);
   drawallparticles(0, 1);

   return;
}

/****************************************************************************/
void initialize(void)
{
   //Module for preparing system and variables for simulation
   int gdriver; int gmode; int errorcode;

   dataentry = 0;
   if (!source)
   {
      if (!pageturn)
      {
	 if ((checkretro()) && (!rerr)) retroproject();
	 else t = 0;
	 initgalaxies();
      }
      if (store == 0) maxt = 50000;
      if (store > 0) maxt = numentry * recfreq * timestep + igt0;
   }
   else
   {
      loadsimulation();
      dataentry = 0;
      findpos();
      t = igt0;
      maxt = (numentry - 1) * recfreq * timestep + igt0;
   }

   x[0] = y[0] = z[0] = 0;
   centx = centy = centz = 0;

   if ((track == 1) || (track == 2)) ptrack = track; else ptrack = 0;
   if ((track == 3) || (source == 1)) {m[1] = gm[1]; m[2] = gm[2];}

   if (gm[1] > 0) pnucleus1 = 1; else pnucleus1 = 0;
   if (gm[2] > 0) pnucleus2 = 2; else pnucleus2 = 0;

   gdriver = DETECT;
   initgraph(&gdriver, &gmode, "EGAVGA.BGI");

   return;
}
/****************************************************************************/
void initgalaxies(void)
{
   //Creates galaxy disks based on disk parameters

   int galaxy;                       //galaxy ID
   float ringspace;                  //space between rings
   int ring;                         //ring ID
   float ringr;                      //radius of ring
   int np;                           //new particle ID
   int npmax;                        //ID of last particle of ring
   float ia;                         //orientation of ring particle to nucleus
   int ringp;                        //particles per ring
   float egm;                        //effective galactic mass

   np = 1;
   if (gm[1] > 0)
   {
      x[1] = gx[1]; y[1] = gy[1]; z[1] = gz[1]; m[1] = gm[1];
      vx[1] = gvx[1]; vy[1] = gvy[1]; vz[1] = gvz[1];
      c[1] = 15; c2[1] = 8;
      np ++;
   }
   if (gm[2] > 0)
   {
      x[2] = gx[2]; y[2] = gy[2]; z[2] = gz[2]; m[2] = gm[2];
      vx[2] = gvx[2]; vy[2] = gvy[2]; vz[2] = gvz[2];
      c[2] = 15; c2[2] = 8;
      np ++;
   }

   kcalc();

   for (galaxy = 1; galaxy <= 2; galaxy ++)
   {
      if ((gring[galaxy] > 0) && (gm[galaxy] > 0))
      {
	 if (gring[galaxy] > 1) ringspace = ((float)(grmax[galaxy] - grmin[galaxy]) / (float)(gring[galaxy] - 1)); else ringspace = 0;
	 ringr = grmin[galaxy];
	 ringp = gp[galaxy];

	 for (ring = 1; ring <= gring[galaxy]; ring ++)
	 {
	    ia = 0;
	    npmax = np + ringp - 1;
	    if (ringr < grmmax[galaxy]) egm = m[galaxy] * pow(ringr, 2) * (ki[galaxy] * ringr + kj[galaxy]) / kk[galaxy];
	    else egm = m[galaxy];
	    for (; np <= npmax; np ++)
	    {
	       m[np] = 0;
	       c[np] = gc[galaxy];
	       c2[np] = 0;

	       ia += 2 * 3.14159 / ringp;
	       transx0 = cos(ia) * ringr;
	       transy0 = sin(ia) * ringr;
	       transz0 = 0;
	       transform(0, gtilt[galaxy], gskew[galaxy]);
	       x[np] = transx2 + gx[galaxy];
	       y[np] = transy2 + gy[galaxy];
	       z[np] = transz2 + gz[galaxy];

	       transx0 = - sqrt(G * egm / ringr) * sin(ia);
	       transy0 =   sqrt(G * egm / ringr) * cos(ia);
	       transz0 = 0;
	       transform(0, gtilt[galaxy], gskew[galaxy]);
	       vx[np] = transx2 + gvx[galaxy];
	       vy[np] = transy2 + gvy[galaxy];
	       vz[np] = transz2 + gvz[galaxy];

	    }
	    ringr += ringspace;
	 }
      }
   }

   nump = np - 1;

   return;
}
/****************************************************************************/
void transform(float rxy, float rxz, float ryz)
{
   //Applies 3-D axis rotation to recalcuate x, y, z coordinates or vectors

   rxy = rxy * dtor;
   rxz = rxz * dtor;
   ryz = ryz * dtor;
   transx1 = transx0 * cos(rxy) + transy0 * sin(rxy);
   transy1 = transy0 * cos(rxy) - transx0 * sin(rxy);
   transx2 = transx1 * cos(rxz) + transz0 * sin(rxz);
   transz1 = transz0 * cos(rxz) - transx1 * sin(rxz);
   transy2 = transy1 * cos(ryz) + transz1 * sin(ryz);
   transz2 = transz1 * cos(ryz) - transy1 * sin(ryz);
   return;
}

/****************************************************************************/
void retroproject(void)
{
   //Simulates galactic nuclei in reverse time to find initial properties

   float egm1;           //effective galaxy 1 mass
   float egm2;           //effective galaxy 2 mass
   float igvel1;         //resultant velocity of galaxy 1
   float igvel2;         //resultant velocity of galaxy 2

   transx0 = igdis;
   transy0 = 0;
   transz0 = 0;
   transform(igarg - 90, iginc, 0);
   x[1] = -transx2 / 2;
   y[1] = -transy2 / 2;
   z[1] = -transz2 / 2;
   x[2] = transx2 / 2;
   y[2] = transy2 / 2;
   z[2] = transz2 / 2;

   kcalc();
   if (igdis < grmmax[1]) egm1 = gm[1] * pow(igdis, 2) * (ki[1] * igdis + kj[1]) / kk[1];
   else egm1 = gm[1];
   if (igdis < grmmax[2]) egm2 = gm[2] * pow(igdis, 2) * (ki[2] * igdis + kj[2]) / kk[2];
   else egm2 = gm[2];

   igvel1 = sqrt(G * egm2 * (1 + igecc) / ( (egm1 / egm2 + 1) * igdis));
   igvel2 = sqrt(G * egm1 * (1 + igecc) / ( (egm2 / egm1 + 1) * igdis));

   transx0 = 0;
   transy0 = igvel1;
   transz0 = 0;
   transform(igarg - 90, iginc, 0);
   vx[1] = -transx2;
   vy[1] = -transy2;
   vz[1] = -transz2;

   transx0 = 0;
   transy0 = igvel2;
   transz0 = 0;
   transform(igarg - 90, iginc, 0);
   vx[2] = transx2;
   vy[2] = transy2;
   vz[2] = transz2;

   m[1] = gm[1];
   m[2] = gm[2];

   timestep = -timestep;
   nump = 2;
   t = 0;

   galr = perir = igdis;

   while (t > igt0)
   {
      corecalc();
      t += timestep;
   }

   timestep = -timestep;

   gx[1] = x[1];   gy[1] = y[1];   gz[1] = z[1];
   gx[2] = x[2];   gy[2] = y[2];   gz[2] = z[2];
   gvx[1] = vx[1]; gvy[1] = vy[1]; gvz[1] = vz[1];
   gvx[2] = vx[2]; gvy[2] = vy[2]; gvz[2] = vz[2];

   return;
}

/****************************************************************************/
void snapshot(void)
{
   //Records graphical image of data to disk for later viewing with snapreader

   int n;                     //miscellaneous counting variable
   int p1;                    //particle ID

   strcpy(filelocation, filepath);   //hold directory path

   clearviewport();

   strcpy(filename, simname);
   if (fileindex < 10)
   {
      strcat(filename, "0");
      strcat(filename, ecvt(fileindex, 0, 0, 0));
   }
   else strcat(filename, ecvt(fileindex, 2, 0, 0));
   strcat(filename, ".SNP");

   gotoxy(40, 1);
   printf("Snapshot image created as %s", filename);

   strcat(filelocation, filename);

   fptr = fopen(filelocation, "w");
   if (!fptr) {printf("\a"); exit(1);}

   fprintf(fptr, "%d %f ", nump, t);
   for (p1 = 0; p1 <= nump; p1 ++)
   {
      fprintf(fptr, "%d %d %d ", rx[p1], ry[p1], c[p1]);
   }
   fclose(fptr);

   fileindex += 1;
   drawallparticles(0, 1);

   return;
}

/*******************************************************************/
void storepos(void)
{
   //Stores current particle positions in memory

   int p1;     //particle ID

   if ((evendiv(t / timestep, recfreq)) && (dataentry < numentry))
   {
      for (p1 = 1; p1 <= nump; p1 ++)
      {
	 sx[p1][dataentry] = x[p1];
	 sy[p1][dataentry] = y[p1];
	 sz[p1][dataentry] = z[p1];
      }
      gotoxy(1, 1); printf("%d", dataentry);
      dataentry ++;
   }

   return;
}
/****************************************************************************/
void storeinitial(void)
{
   //stores essential initial properties and interface settings to disk

   int p1;        //particle ID

   strcpy(filelocation, filepath);
   strcpy(filename, simname);
   strcat(filename, ".INT");
   strcat(filelocation, filename);
   fptr = fopen(filelocation, "w");
   if (!fptr) {printf("\a"); exit(1);}
   fprintf(fptr, "%d %d %f %d %f ", nump, recfreq, timestep, numentry, igt0);
   for (p1 = 1; p1 <= nump; p1 ++)
   {
      fprintf(fptr, "%d ", c[p1]);
   }
   fclose(fptr);

   strcpy(filelocation, filepath);
   strcpy(filename, simname);
   strcat(filename, "I.SET");
   strcat(filelocation, filename);
   savesettings();

   return;
}

/****************************************************************************/
void storesimulation(void)
{
   //Stores all historic particle positions to disk
   int p1;      //particle ID

   strcpy(filelocation, filepath);

   strcpy(filename, simname);
   strcat(filename, ".SIM");
   strcat(filelocation, filename);

   fptr = fopen(filelocation, "w");
   if (!fptr) {printf("\a"); exit(1);}

   for (dataentry = 0; dataentry < numentry; dataentry ++)
   {
      for (p1 = 1; p1 <= nump; p1 ++)
      {
	 fprintf(fptr, "%f %f %f ", sx[p1][dataentry], sy[p1][dataentry], sz[p1][dataentry]);
      }
   }

   fclose(fptr);

   return;
}

/****************************************************************************/
void loadsimulation(void)
{
   //Loads historic particle positions from disk into memory

   int p1;          //particle ID

   strcpy(filelocation, filepath);
   strcpy(filename, simname);
   strcat(filename, ".INT");
   strcat(filelocation, filename);

   fptr = fopen(filelocation, "r");
   if (!fptr) {printf("\a"); exit(1);}


   fscanf(fptr, "%d %d %f %d ", &nump, &recfreq, &timestep, &numentry);
   if (numentry == 0) numentry = 60;     //accomodate old files
   else fscanf(fptr, "%f ", &igt0);
   for (p1 = 1; p1 <= nump; p1 ++)
   {
     fscanf(fptr, "%d ", &c[p1]);
   }
   fclose(fptr);

   strcpy(filelocation, filepath);
   strcpy(filename, simname);
   strcat(filename, ".SIM");
   strcat(filelocation, filename);

   fptr = fopen(filelocation, "r");
   if (!fptr) {printf("\a"); exit(1);}    //should never happen

   for (dataentry = 0; dataentry < numentry; dataentry ++)
   {
      for (p1 = 1; p1 <= nump; p1 ++)
      {
	 fscanf(fptr, "%f %f %f ", &sx[p1][dataentry], &sy[p1][dataentry], &sz[p1][dataentry]);
      }
   }

   fclose(fptr);

   return;
}
/****************************************************************************/
void findpos(void)
{
   //Loads current position and calculates velocity from historic data

   int p1;             //particle ID

   for (p1 = 1; p1 <= nump; p1 ++)
   {
      x[p1] = sx[p1][dataentry];
      y[p1] = sy[p1][dataentry];
      z[p1] = sz[p1][dataentry];
      if (dataentry + 1 < numentry)
      {
	 vx[p1] = (sx[p1][dataentry + 1] - x[p1]) / recfreq / timestep;
	 vy[p1] = (sy[p1][dataentry + 1] - y[p1]) / recfreq / timestep;
	 vz[p1] = (sz[p1][dataentry + 1] - z[p1]) / recfreq / timestep;
      }
   }
   dataentry ++;

   return;
}

/****************************************************************************/
void initmemory(void)
{
   //Initializes dynamic memory to store historic coordinate data

   int index;                 //pointer array reference

   for (index = 0; index <= maxp; index ++)
   {
      sx[index] = (float far *) farmalloc(numentry * sizeof(float));
      if (!sx[index]) {printf("Out of Memory"); getch(); exit(1);}
   }
   for (index = 0; index <= maxp; index ++)
   {
      sy[index] = (float far *) farmalloc(numentry * sizeof(float));
      if (!sy[index]) {printf("Out of Memory"); getch(); exit(1);}
   }
   for (index = 0; index <= maxp; index ++)
   {
      sz[index] = (float far *) farmalloc(numentry * sizeof(float));
      if (!sz[index]) {printf("Out of Memory"); getch(); exit(1);}
   }

   return;
}
/****************************************************************************/
void terminate(void)
{
   //Frees memory and calls termination functions

   int index;         //pointer array reference
   int n;             //miscellaneous counting variable

   for (n = 0; n < 7; n ++) capsname[n] = toupper(simname[n]);

   if ((!autoterminate) && (!store)) {gotoxy(38, 1); printf("simulation terminated by user at %.0f", t);}
   if ((!autoterminate) && (store))  {gotoxy(41, 1); printf("simulation aborted by user at %.0f", t);}
   if ((autoterminate) && (store))   {gotoxy(34, 1); printf("simulation saved as %s.SIM at %.0f", capsname, t);}
   if ((autoterminate) && (!store))  {gotoxy(31, 1); printf("simulation terminated at end of data at %.0f", t);}

   if ((dataentry == numentry) && (store == 2))
   {
      storeinitial();
      storesimulation();
   }

   getch();
   closegraph();

   for (index = 0; index < maxp; index ++) farfree(sx[index]);
   for (index = 0; index < maxp; index ++) farfree(sy[index]);
   for (index = 0; index < maxp; index ++) farfree(sz[index]);

   t = 0;

   return;
}
/****************************************************************************/
int evendiv(float num, float den)
{
   //Determines whether one number is a multiple of another
   return ( (num / den) == floor(num / den));
}
/****************************************************************************/
void IFmodule(void)
{
   //Module for interface operations

   _setcursortype(_NOCURSOR);

   if (runsim) restorecrtmode();
   textmode(3);
   nextpage = 1;
   pageturn = 0;

   do
   {
      if (nextpage) IFproperties();
      if (nextpage) IFgalaxies();
      if (nextpage) IFparticles();
   } while ((!exitflag) && (!runsim));

   if (exitflag) exit(0);

   return;
}

/****************************************************************************/
void IFgalaxies(void)
{
   //Module for galaxy parameter page of interface

   int l;                         //current cursor location and option ID
   float holdfloat;               //holds generic floating-point data
   char holdstring[13] = "\0";    //holds generic string data

   struct galaxies
   {
      char label[16];
      char help[60];
      int * ivar;
      float * fvar;
      int ypos, xlpos, xspos;
   } ;
   struct galaxies galx[41];

   nextpage = 0;
   pageturn = 1;


   strcpy(galx[1].label, "        x:");
   strcpy(galx[1].help, "Initial x-coordinate of galaxy center");
   galx[1].fvar = &gx[1]; galx[1].ivar = 0;
   galx[1].ypos = 5; galx[1].xlpos = 2; galx[1].xspos = 12;

   strcpy(galx[2].label, "        y:");
   strcpy(galx[2].help, "Initial y-coordinate of galaxy center");
   galx[2].fvar = &gy[1]; galx[2].ivar = 0;
   galx[2].ypos = 6; galx[2].xlpos = 2; galx[2].xspos = 12;

   strcpy(galx[3].label, "        z:");
   strcpy(galx[3].help, "Initial z-coordinate of galaxy center");
   galx[3].fvar = &gz[1]; galx[3].ivar = 0;
   galx[3].ypos = 7; galx[3].xlpos = 2; galx[3].xspos = 12;

   strcpy(galx[4].label, "       vx:");
   strcpy(galx[4].help, "Initial x-axis velocity of galaxy");
   galx[4].fvar = &gvx[1]; galx[4].ivar = 0;
   galx[4].ypos = 8; galx[4].xlpos = 2; galx[4].xspos = 12;

   strcpy(galx[5].label, "       vy:");
   strcpy(galx[5].help, "Initial y-axis velocity of galaxy");
   galx[5].fvar = &gvy[1]; galx[5].ivar = 0;
   galx[5].ypos = 9; galx[5].xlpos = 2; galx[5].xspos = 12;

   strcpy(galx[6].label, "       vz:");
   strcpy(galx[6].help, "Initial z-axis velocity of galaxy");
   galx[6].fvar = &gvz[1]; galx[6].ivar = 0;
   galx[6].ypos = 10; galx[6].xlpos = 2; galx[6].xspos = 12;

   strcpy(galx[7].label, "     Mass:");
   strcpy(galx[7].help, "Total mass of galaxy");
   galx[7].fvar = &gm[1]; galx[7].ivar = 0;
   galx[7].ypos = 11; galx[7].xlpos = 2; galx[7].xspos = 12;

   strcpy(galx[8].label, "     Rmin:");
   strcpy(galx[8].help, "Radius of innermost ring");
   galx[8].fvar = 0; galx[8].ivar = &grmin[1];
   galx[8].ypos = 13; galx[8].xlpos = 2; galx[8].xspos = 12;

   strcpy(galx[9].label, "     Rmax:");
   strcpy(galx[9].help, "Radius of outermost ring");
   galx[9].fvar = 0; galx[9].ivar = &grmax[1];
   galx[9].ypos = 14; galx[9].xlpos = 2; galx[9].xspos = 12;

   strcpy(galx[10].label, "    Rings:");
   strcpy(galx[10].help, "Total number of rings");
   galx[10].fvar = 0; galx[10].ivar = &gring[1];
   galx[10].ypos = 15; galx[10].xlpos = 2; galx[10].xspos = 12;

   strcpy(galx[11].label, "Particles:");
   strcpy(galx[11].help, "Particles per ring");
   galx[11].fvar = 0; galx[11].ivar = &gp[1];
   galx[11].ypos = 16; galx[11].xlpos = 2; galx[11].xspos = 12;

   strcpy(galx[12].label, "     Tilt:");
   strcpy(galx[12].help, "x-z inclination of galaxy");
   galx[12].fvar = 0; galx[12].ivar = &gtilt[1];
   galx[12].ypos = 17; galx[12].xlpos = 2; galx[12].xspos = 12;

   strcpy(galx[13].label, "     Skew:");
   strcpy(galx[13].help, "y-z inclination of galaxy");
   galx[13].fvar = 0; galx[13].ivar = &gskew[1];
   galx[13].ypos = 18; galx[13].xlpos = 2; galx[13].xspos = 12;

   strcpy(galx[14].label, "    Color:");
   strcpy(galx[14].help, "Color of disk particles");
   galx[14].fvar = 0; galx[14].ivar = &gc[1];
   galx[14].ypos = 19; galx[14].xlpos = 2; galx[14].xspos = 12;

   strcpy(galx[15].label, "        x:");
   strcpy(galx[15].help, galx[1].help);
   galx[15].fvar = &gx[2]; galx[15].ivar = 0;
   galx[15].ypos = 5; galx[15].xlpos = 26; galx[15].xspos = 36;

   strcpy(galx[16].label, "        y:");
   strcpy(galx[16].help, galx[2].help);
   galx[16].fvar = &gy[2]; galx[16].ivar = 0;
   galx[16].ypos = 6; galx[16].xlpos = 26; galx[16].xspos = 36;

   strcpy(galx[17].label, "        z:");
   strcpy(galx[17].help, galx[3].help);
   galx[17].fvar = &gz[2]; galx[17].ivar = 0;
   galx[17].ypos = 7; galx[17].xlpos = 26; galx[17].xspos = 36;

   strcpy(galx[18].label, "       vx:");
   strcpy(galx[18].help, galx[4].help);
   galx[18].fvar = &gvx[2]; galx[18].ivar = 0;
   galx[18].ypos = 8; galx[18].xlpos = 26; galx[18].xspos = 36;

   strcpy(galx[19].label, "       vy:");
   strcpy(galx[19].help, galx[5].help);
   galx[19].fvar = &gvy[2]; galx[19].ivar = 0;
   galx[19].ypos = 9; galx[19].xlpos = 26; galx[19].xspos = 36;

   strcpy(galx[20].label, "       vz:");
   strcpy(galx[20].help, galx[6].help);
   galx[20].fvar = &gvz[2]; galx[20].ivar = 0;
   galx[20].ypos = 10; galx[20].xlpos = 26; galx[20].xspos = 36;

   strcpy(galx[21].label, "     Mass:");
   strcpy(galx[21].help, galx[7].help);
   galx[21].fvar = &gm[2]; galx[21].ivar = 0;
   galx[21].ypos = 11; galx[21].xlpos = 26; galx[21].xspos = 36;

   strcpy(galx[22].label, "     Rmin:");
   strcpy(galx[22].help, galx[8].help);
   galx[22].fvar = 0; galx[22].ivar = &grmin[2];
   galx[22].ypos = 13; galx[22].xlpos = 26; galx[22].xspos = 36;

   strcpy(galx[23].label, "     Rmax:");
   strcpy(galx[23].help, galx[9].help);
   galx[23].fvar = 0; galx[23].ivar = &grmax[2];
   galx[23].ypos = 14; galx[23].xlpos = 26; galx[23].xspos = 36;

   strcpy(galx[24].label, "    Rings:");
   strcpy(galx[24].help, galx[10].help);
   galx[24].fvar = 0; galx[24].ivar = &gring[2];
   galx[24].ypos = 15; galx[24].xlpos = 26; galx[24].xspos = 36;

   strcpy(galx[25].label, "Particles:");
   strcpy(galx[25].help, galx[11].help);
   galx[25].fvar = 0; galx[25].ivar = &gp[2];
   galx[25].ypos = 16; galx[25].xlpos = 26; galx[25].xspos = 36;

   strcpy(galx[26].label, "     Tilt:");
   strcpy(galx[26].help, galx[12].help);
   galx[26].fvar = 0; galx[26].ivar = &gtilt[2];
   galx[26].ypos = 17; galx[26].xlpos = 26; galx[26].xspos = 36;

   strcpy(galx[27].label, "     Skew:");
   strcpy(galx[27].help, galx[13].help);
   galx[27].fvar = 0; galx[27].ivar = &gskew[2];
   galx[27].ypos = 18; galx[27].xlpos = 26; galx[27].xspos = 36;

   strcpy(galx[28].label, "    Color:");
   strcpy(galx[28].help, galx[14].help);
   galx[28].fvar = 0; galx[28].ivar = &gc[2];
   galx[28].ypos = 19; galx[28].xlpos = 26; galx[28].xspos = 36;

   strcpy(galx[29].label, "      Distance:");
   strcpy(galx[29].help, "Distance between galactic centers at perigalacticon");
   galx[29].fvar = &igdis; galx[29].ivar = 0;
   galx[29].ypos = 5; galx[29].xlpos = 53; galx[29].xspos = 68;

   strcpy(galx[30].label, "  Eccentricity:");
   strcpy(galx[30].help, "Relative velocity of galactic centers at perigalacticon");
   galx[30].fvar = &igecc; galx[30].ivar = 0;
   galx[30].ypos = 6; galx[30].xlpos = 53; galx[30].xspos = 68;

   strcpy(galx[31].label, "   Inclination:");
   strcpy(galx[31].help, "Inclination of galaxy 2 orbital plane to galaxy 1 disk");
   galx[31].fvar = &iginc; galx[31].ivar = 0;
   galx[31].ypos = 7; galx[31].xlpos = 53; galx[31].xspos = 68;

   strcpy(galx[32].label, "      Argument:");
   strcpy(galx[32].help, "Rotation of galaxy 2 perigalacticon point in orbital plane");
   galx[32].fvar = &igarg; galx[32].ivar = 0;
   galx[32].ypos = 8; galx[32].xlpos = 53; galx[32].xspos = 68;

   strcpy(galx[33].label, "  Initial Time:");
   strcpy(galx[33].help, "Time of start of simulation relative to perigalacticon");
   galx[33].fvar = &igt0; galx[33].ivar = 0;
   galx[33].ypos = 9; galx[33].xlpos = 53; galx[33].xspos = 68;

   strcpy(galx[34].label, " Cent. Density:");
   strcpy(galx[34].help, "Relative area mass density at center of galactic disk");
   galx[34].fvar = 0; galx[34].ivar = &gpo[1];
   galx[34].ypos = 12; galx[34].xlpos = 53; galx[34].xspos = 68;

   strcpy(galx[35].label, "  Edge Density:");
   strcpy(galx[35].help, "Relative area mass density at edge of massive disk");
   galx[35].fvar = 0; galx[35].ivar = &gpmax[1];
   galx[35].ypos = 13; galx[35].xlpos = 53; galx[35].xspos = 68;

   strcpy(galx[36].label, "Massive Radius:");
   strcpy(galx[36].help, "Radius of massive disk");
   galx[36].fvar = 0; galx[36].ivar = &grmmax[1];
   galx[36].ypos = 14; galx[36].xlpos = 53; galx[36].xspos = 68;

   strcpy(galx[37].label, " Cent. Density:");
   strcpy(galx[37].help, galx[34].help);
   galx[37].fvar = 0; galx[37].ivar = &gpo[2];
   galx[37].ypos = 17; galx[37].xlpos = 53; galx[37].xspos = 68;

   strcpy(galx[38].label, "  Edge Density:");
   strcpy(galx[38].help, galx[35].help);
   galx[38].fvar = 0; galx[38].ivar = &gpmax[2];
   galx[38].ypos = 18; galx[38].xlpos = 53; galx[38].xspos = 68;

   strcpy(galx[39].label, "Massive Radius:");
   strcpy(galx[39].help, galx[36].help);
   galx[39].fvar = 0; galx[39].ivar = &grmmax[2];
   galx[39].ypos = 19; galx[39].xlpos = 53; galx[39].xspos = 68;




   textcolor(7); textbackground(0);

   clrscr();
   IFprintheader();

   gotoxy(72, 1); cprintf("Galaxies");

   gotoxy(2, 4); cprintf("GALAXY 1");
   gotoxy(26, 4); cprintf("GALAXY 2");
   gotoxy(54, 4); cprintf("RETROTRAJECTORY");
   gotoxy(54, 11); cprintf("GALAXY 1 MODEL");
   gotoxy(54, 16); cprintf("GALAXY 2 MODEL");

   textcolor(7);
   for (l = 1; l <= 39; l ++)
   {
      gotoxy(galx[l].xlpos, galx[l].ypos); cprintf("%s", galx[l].label);

      gotoxy(galx[l].xspos, galx[l].ypos);
      if (galx[l].fvar != 0)
      {
	 * galx[l].fvar = fix(* galx[l].fvar, 5);
	 gcvt(* galx[l].fvar, 7, holdstring);
	 cprintf(" %s ", holdstring);
      }
      if (galx[l].ivar != 0) cprintf(" %d ", * galx[l].ivar);
   }

   l = 1;

   do
   {
      textcolor(7);
      gotoxy(10, 22); cprintf("                                                            ");
      gotoxy(10, 22); cprintf("%s", galx[l].help);

      printfunct(0, 4, 1, checkretro(), 0);

      textcolor(15);
      gotoxy(galx[l].xlpos, galx[l].ypos); cprintf("%s", galx[l].label);

      textbackground(3);


      gotoxy(galx[l].xspos + 1, galx[l].ypos);

      if (galx[l].fvar != 0) gcvt(* galx[l].fvar, 7, holdstring);
      if (galx[l].ivar != 0) itoa(* galx[l].ivar, holdstring, 10);

      liminput(holdstring, 10, 1, 0, 1);

      holdfloat = strtod(holdstring, 0);
      if (((l >= 7) && (l <= 11)) || ((l >= 21) && (l <= 25)) || (l == 14) || (l == 28) || (l == 29) || (l == 30) || (l > 33)) holdfloat = fabs(holdfloat);
      if ((l == 33) && (holdfloat != 0)) holdfloat = -fabs(holdfloat);
      if (((l == 10) || (l == 24)) && (holdfloat > 30)) holdfloat = 30;
      if (((l == 11) || (l == 25)) && (holdfloat > 300)) holdfloat = 300;
      if ((l == 12) || (l == 13) || (l == 26) || (l == 27) || (l == 31) || (l == 32))
      {
	 if (holdfloat < -180) holdfloat += 360 * (int)(abs(holdfloat) / 360) + 360;
	 holdfloat += 179;
	 holdfloat = (int)(holdfloat) % 360;
	 holdfloat -= 179;
      }
      if ((l == 14) || (l == 28)) holdfloat = (int)(holdfloat) % 16;
      if ((l == 33) && (holdfloat < -1e+6)) holdfloat = 0;
      if (((l == 8) || (l == 9) || (l == 22) || (l == 23) || (l == 36) || (l == 39)) && (holdfloat > 1000)) holdfloat = 1000;



      if (galx[l].fvar != 0) * galx[l].fvar = holdfloat;
      if (galx[l].ivar != 0)
      {
	 if (fabs(holdfloat) < 32767) * galx[l].ivar = (int)(holdfloat);
	 else * galx[l].ivar = 0;
      }


      textcolor(7); textbackground(0);
      gotoxy(galx[l].xspos, galx[l].ypos); cprintf("             ");

      gotoxy(galx[l].xlpos, galx[l].ypos); cprintf("%s", galx[l].label);
      gotoxy(galx[l].xspos, galx[l].ypos);
      if (galx[l].fvar != 0) {gcvt(* galx[l].fvar, 7, holdstring); cprintf(" %s ", holdstring);}
      if (galx[l].ivar != 0) cprintf(" %d ", * galx[l].ivar);


      if (com == 72) l --;
      if (com == 80) l ++;
      if (l == 0) l = 39;
      if (l == 40) l = 1;
      if (com == 60) nextpage = -1;
      if (com == 61) exitflag = 1;
      if (com == 62)
      {
	 if ((checkretro()) && (!rerr))
	 {
	    err = 0;
	    retroproject();
	    if (!err)
	    {
	       for (l = 1; l <= 20; l ++)
	       {
		  if (l == 7) l = 15;
		  * galx[l].fvar = fix(* galx[l].fvar, 5);
		  gotoxy(galx[l].xspos, galx[l].ypos);
		  gcvt(* galx[l].fvar, 7, holdstring);
		  cprintf(" %s ", holdstring);
		  holdfloat = strtod(holdstring, 0);
		  * galx[l].fvar = holdfloat;
	       }
	    }
	    l = 1;
	 }
	 err = rerr;
	 printerrmsg();
	 if (err) l = errl;
      }

      if (nextpage)
      {
	 err = checkdisk();
	 printerrmsg();
	 if (err) { l = errl;  nextpage = 0; }
      }
   } while (!exitflag && !nextpage);
   com = 0;
}
/****************************************************************************/
void IFproperties(void)
{
   //Module for metaproperty parameter if simulator interface

   int n;                       //miscellaneous counting variable
   int l;                       //cursor position and option ID
   int dir;                     //selection change instructional code
   char sellabel[18];           //selection label

   struct properties
   {
      char label[18];
      char help[60];
      int ypos, xlpos, xspos;
   } ;
   struct properties prop[21];

   nextpage = 0;
   com = 0;
   check = 0;
   runsim = 0;



   strcpy(prop[1].label, " Simulation Name:");
   strcpy(prop[1].help, "Simulation title, max 6 characters");
   prop[1].ypos = 5; prop[1].xlpos = 5; prop[1].xspos = 22;

   strcpy(prop[2].label, "           Input:");
   strcpy(prop[2].help, "Source of coordinate data");
   prop[2].ypos = 6; prop[2].xlpos = 5; prop[2].xspos = 22;

   strcpy(prop[3].label, "          Output:");
   strcpy(prop[3].help, "Storage of coordinate data");
   prop[3].ypos = 7; prop[3].xlpos = 5; prop[3].xspos = 22;

   strcpy(prop[4].label, "Record Frequency:");
   strcpy(prop[4].help, "Coordinate data is recorded in memory every x timesteps");
   prop[4].ypos = 8; prop[4].xlpos = 5; prop[4].xspos = 22;

   strcpy(prop[5].label, "   Load Settings:");
   strcpy(prop[5].help, "Load previously recorded settings into interface");
   prop[5].ypos = 10; prop[5].xlpos = 5; prop[5].xspos = 22;

   strcpy(prop[6].label, "   Save Settings:");
   strcpy(prop[6].help, "Save current settings to disk");
   prop[6].ypos = 11; prop[6].xlpos = 5; prop[6].xspos = 22;

   strcpy(prop[7].label, " Coarse Timestep:");
   strcpy(prop[7].help, "Basic timestep unit (stu)");
   prop[7].ypos = 15; prop[7].xlpos = 5; prop[7].xspos = 22;

   strcpy(prop[8].label, "   Fine Timestep:");
   strcpy(prop[8].help, "Minimum timestep unit (stu)");
   prop[8].ypos = 16; prop[8].xlpos = 5; prop[8].xspos = 22;

   strcpy(prop[9].label, "           Delay:");
   strcpy(prop[9].help, "Duration of decelerating timing loop");
   prop[9].ypos = 17; prop[9].xlpos = 5; prop[9].xspos = 22;

   strcpy(prop[10].label, "      Monitoring:");
   strcpy(prop[10].help, "Display simulation on screen");
   prop[10].ypos = 5; prop[10].xlpos = 39; prop[10].xspos = 56;

   strcpy(prop[11].label, "   Redraw Method:");
   strcpy(prop[11].help, "Particle redraw: erase before or during particle draw");
   prop[11].ypos = 6; prop[11].xlpos = 39; prop[11].xspos = 56;

   strcpy(prop[12].label, "  Draw Frequency:");
   strcpy(prop[12].help, "Screen is updated every x timesteps");
   prop[12].ypos = 7; prop[12].xlpos = 39; prop[12].xspos = 56;

   strcpy(prop[13].label, "  Particle Style:");
   strcpy(prop[13].help, "Size and shape of particle images");
   prop[13].ypos = 8; prop[13].xlpos = 39; prop[13].xspos = 56;

   strcpy(prop[14].label, "        Tracking:");
   strcpy(prop[14].help, "Fix camera on specified object");
   prop[14].ypos = 9; prop[14].xlpos = 39; prop[14].xspos = 56;

   strcpy(prop[15].label, " Doppler Mimicry:");
   strcpy(prop[15].help, "Color particles by line-of-sight velocity");
   prop[15].ypos = 10; prop[15].xlpos = 39; prop[15].xspos = 56;

   strcpy(prop[16].label, "     Z-Axis Fade:");
   strcpy(prop[16].help, "Color particles by line-of-sight distance");
   prop[16].ypos = 11; prop[16].xlpos = 39; prop[16].xspos = 56;

   strcpy(prop[17].label, "            SDU =");
   strcpy(prop[17].help, "Simulation distance unit (meters)");
   prop[17].ypos = 15; prop[17].xlpos = 39; prop[17].xspos = 56;

   strcpy(prop[18].label, "            SMU =");
   strcpy(prop[18].help, "Simulation mass unit (kilograms)");
   prop[18].ypos = 16; prop[18].xlpos = 39; prop[18].xspos = 56;

   strcpy(prop[19].label, "            STU =");
   strcpy(prop[19].help, "Simulation time unit (seconds)");
   prop[19].ypos = 17; prop[19].xlpos = 39; prop[19].xspos = 56;

   strcpy(prop[20].label, "              G =");
   strcpy(prop[20].help, "Simulation gravitational constant (smu^3 * stu^2 / sdu)");
   prop[20].ypos = 18; prop[20].xlpos = 39; prop[20].xspos = 56;




   clrscr();
   textcolor(7);
   IFprintheader();

   gotoxy(70, 1); cprintf("Properties");
   gotoxy(2, 4); cprintf("FILE");
   gotoxy(2, 14); cprintf("TIME");
   gotoxy(36, 4); cprintf("GRAPHICS");
   gotoxy(36, 14); cprintf("SCALING");

   textcolor(7);
   for (l = 1; l <= 20; l ++)
   {
      gotoxy(prop[l].xlpos, prop[l].ypos); cprintf("%s", prop[l].label);
      IFgetsellabel(sellabel, l);
      gotoxy(prop[l].xspos, prop[l].ypos); cprintf(" %s ", sellabel);
   }

   l = 1;

   do
   {
      textcolor(7);
      gotoxy(10, 22); cprintf("                                                        ");
      gotoxy(10, 22); cprintf("%s", prop[l].help);

      printfunct((check + 1), !source * 3, 1, 0, 0);

      textcolor(15);
      gotoxy(prop[l].xlpos, prop[l].ypos); cprintf("%s", prop[l].label);


      if ((l != 1) && (l < 17)) textbackground(3);
      IFgetsellabel(sellabel, l);
      if (strlen(sellabel) > 0) {gotoxy(prop[l].xspos, prop[l].ypos); cprintf(" %s ", sellabel);}

      textcolor(7); textbackground(0);

      dir = IFinput();

      if (dir < 3)
      {
	 gotoxy(prop[l].xlpos, prop[l].ypos); cprintf("%s", prop[l].label);
	 IFgetsellabel(sellabel, l);
	 gotoxy(prop[l].xspos, prop[l].ypos); cprintf(" %s ", sellabel);
      }


      if (dir == 1) l++;
      if (dir == 2) l--;
      if (l > 20) l = 1;
      if (l < 1) l = 20;

      if (check == -1)
      {
	 err = 0;
	 err = checkprop();
	 printerrmsg();
	 if ((source == 1) && (!err))
	 {
	    strcpy(filelocation, filepath);
	    strcpy(filename, simname);
	    strcat(filename, "I.SET");
	    strcat(filelocation, filename);
	    loadsettings(0);
	 }
	 if (!err) check = 1;
	 if (err) {if (errl) l = errl; check = 0; }
      }

      gotoxy(prop[l].xspos + 1, prop[l].ypos);
      textcolor(15); textbackground(3);

      if (dir > 2) IFchangesel(dir, l);

      textcolor(7); textbackground(0);

      if (!runsim) {gotoxy(prop[l].xspos, prop[l].ypos); cprintf("                ");}


      if (l >= 17)                                //print other selections
      {
	 for (n = 17; n <= 20; n ++)
	 {
	    gotoxy(prop[n].xspos, prop[n].ypos); cprintf("               ");
	    IFgetsellabel(sellabel, n);
	    gotoxy(prop[n].xspos, prop[n].ypos); cprintf(" %s ", sellabel);
	 }
      }
      if ((l == 7) || (l == 8))
      {
	 if (l == 7) n = 8; else n = 7;
	 gotoxy(prop[n].xspos, prop[n].ypos); cprintf("               ");
	 IFgetsellabel(sellabel, n);
	 gotoxy(prop[n].xspos, prop[n].ypos); cprintf(" %s ", sellabel);
      }
      if (loadselflag)
      {
	 for (n = 1; n <= 20; n ++)
	 {
	    gotoxy(prop[n].xspos, prop[n].ypos); cprintf("               ");
	    IFgetsellabel(sellabel, n);
	    gotoxy(prop[n].xspos, prop[n].ypos); cprintf(" %s ", sellabel);
	 }
	 loadselflag = 0;
      }

   } while (!exitflag && !nextpage && !runsim);


   return;
}
/****************************************************************************/
void IFparticles(void)
{
   //Particle-viewing page of interface, displays properties for all particles

   int page = 0;         //particle page ID
   int p1;               //particle ID
   int endp;             //last particle on page ID
   int ypos;             //cursor y position
   int n;                //miscellaneous counting variable

   nextpage = 0;

   initgalaxies();

   textcolor(7); textbackground(0);
   clrscr();
   IFprintheader();
   gotoxy(71, 1); cprintf("Particles");
   gotoxy(1, 3); cprintf("#p");
   gotoxy(8, 3); cprintf("x");
   gotoxy(18, 3); cprintf("y");
   gotoxy(28, 3); cprintf("z");
   gotoxy(41, 3); cprintf("vx");
   gotoxy(51, 3); cprintf("vy");
   gotoxy(61, 3); cprintf("vz");
   gotoxy(72, 3); cprintf("m");
   gotoxy(77, 3); cprintf("c");

   printfunct(0, 5, 1, 0, 0);

   do
   {
     if (nump > (page + 1) * 16) endp = ((page + 1) * 16); else endp = nump;

     for (ypos = 4; ypos < 4 + 16; ypos ++)
     {
	 gotoxy(1, ypos); for (n = 1; n < 81; n ++) putch(' ');
     }
     ypos = 4;
     for (p1 = page * 16 + 1; p1 <= endp; p1 ++)
     {
	 gotoxy(1, ypos); cprintf("%d", p1);
	 gotoxy(8 - (x[p1] < 0), ypos); cprintf("%.3f", x[p1]);
	 gotoxy(18 - (y[p1] < 0), ypos); cprintf("%.3f", y[p1]);
	 gotoxy(28 - (z[p1] < 0), ypos); cprintf("%.3f", z[p1]);
	 gotoxy(41 - (vx[p1] < 0), ypos); cprintf("%.4f", vx[p1]);
	 gotoxy(51 - (vy[p1] < 0), ypos); cprintf("%.4f", vy[p1]);
	 gotoxy(61 - (vz[p1] < 0), ypos); cprintf("%.4f", vz[p1]);
	 gotoxy(72, ypos); cprintf("%.0f", m[p1]); //.4f for non-pt. mass
	 gotoxy(77, ypos); cprintf("%d", c[p1]);
	 ypos += 1;
      }
      com = getch();
      if (com == 61) exitflag = 1;
      if ((com == 80) && (endp < nump)) page ++;
      if ((com == 72) && (page > 0)) page --;
      if (com == 60) nextpage = 1;

   } while (!nextpage && !exitflag);

   return;
}
/****************************************************************************/
void IFprintheader(void)
{
   //Prints header on top of interface screen

   int n;           //miscellaneous counting variable

   gotoxy(1, 1); cprintf("INTERGALACTIC GRAVITATIONAL MOTION SIMULATOR 3.2");
   gotoxy(1, 2);   for (n = 1; n < 81; n ++) putch('-');
   gotoxy(1, 21);  for (n = 1; n < 81; n ++) putch('-');

   gotoxy(1, 22); cprintf("   Help:");
   gotoxy(1, 23); cprintf("  Error:");
   return;
}
/****************************************************************************/
void IFgetsellabel(char * selptr, int l)
{
   //Copies properties selection label into memory

   strcpy(selptr, "");
   if (l == 1) strcpy(selptr, simname);
   if ((l == 2) && (source == 0)) strcpy(selptr, "simulation");
   if ((l == 2) && (source == 1)) strcpy(selptr, "file");
   if ((l == 3) && (store == 0)) strcpy(selptr, "none");
   if ((l == 3) && (store == 1)) strcpy(selptr, "RAM");
   if ((l == 3) && (store == 2)) strcpy(selptr, "disk");
   if (l == 4) itoa(recfreq, selptr, 10);
   if (l == 7) gcvt(timestep, 5, selptr);
   if (l == 8) gcvt(timestep, 5, selptr);
   if (l == 9) ltoa(delay, selptr, 10);
   if ((l == 10) && (monitor == 0)) strcpy(selptr, "disabled");
   if ((l == 10) && (monitor == 1)) strcpy(selptr, "enabled");
   if ((l == 11) && (redraw == 0)) strcpy(selptr, "flash");
   if ((l == 11) && (redraw == 1)) strcpy(selptr, "sequence");
   if ((l == 11) && (redraw == 2)) strcpy(selptr, "sequence++");
   if (l == 12) itoa(drawfreq, selptr, 10);
   if ((l == 13) && (style == 0)) strcpy(selptr, "point");
   if ((l == 13) && (style == 1)) strcpy(selptr, "cross");
   if ((l == 13) && (style == 2)) strcpy(selptr, "circle");
   if ((l == 13) && (style == 3)) strcpy(selptr, "varsmall");
   if ((l == 13) && (style == 4)) strcpy(selptr, "varlarge");
   if ((l == 14) && (track == 0)) strcpy(selptr, "none");
   if ((l == 14) && (track == 1)) strcpy(selptr, "galaxy 1");
   if ((l == 14) && (track == 2)) strcpy(selptr, "galaxy 2");
   if ((l == 14) && (track == 3)) strcpy(selptr, "center of mass");
   if ((l == 15) && (doppler == 0)) strcpy(selptr, "disabled");
   if ((l == 15) && (doppler == 1)) strcpy(selptr, "enabled");
   if ((l == 16) && (fade == 0)) strcpy(selptr, "disabled");
   if ((l == 16) && (fade == 1)) strcpy(selptr, "chromatic");
   if (l == 17) gcvt(sdu, 4, selptr);
   if (l == 18) gcvt(smu, 4, selptr);
   if (l == 19) gcvt(stu, 4, selptr);
   if (l == 20) gcvt(G, 4, selptr);

   return;
}
/****************************************************************************/
int IFinput(void)
{
   //Gets user input and flags appropriately in interface

   int dir = 0;       //selection change instructional code

   do
   {
      if (!skip) com = getch();
      skip = 0;

      if (com == 80) dir = 1;
      if (com == 72) dir = 2;
      if (com == 75) dir = 3;
      if (com == 77) dir = 4;
      if (com == 32) dir = 4;
      if ((com == 59) && (check != 1)) check = -1;
      if ((com == 59) && (check == 1)) runsim = 1; else runsim = 0;
      if ((com == 60) && (!source)) nextpage = -1; else nextpage = 0;
      if (com == 61) exitflag = 1; else exitflag = 0;
      //if (com == 63) IFviewcat(); //viewcat = 1; else viewcat = 0;

   } while (!dir && !exitflag && !check && !nextpage);

   return(dir);
}

/****************************************************************************/
void IFchangesel(int dir, int l)
{
   char holdstring[13] = "\0";         //miscellaneous input/output string

   check = 0;
   if ((l == 1) && ((dir == 4) || (dir == 5))) liminput(simname, 6, 0, 1, 1);
   if (l == 2) source = abs(!source);
   if (l == 3)
   {
      if (dir == 3) store --;
      if (dir == 4) store ++;
      if (store == 3) store = 0;
      if (store == -1) store = 2;
   }
   if (l == 4)
   {
      if (dir == 3) recfreq -= (4 * (recfreq > 25) + 1) * (recfreq > 1);
      if (dir == 4) recfreq += (4 * (recfreq >= 25) + 1) * (recfreq < 100);
   }
   if (((l == 5) || (l == 6)) && (dir == 4))
   {
      cprintf("       .SET");
      gotoxy(wherex() - 10, wherey());
      liminput(holdstring, 6, 0, 1, 0);
      textcolor(7); textbackground(0);
      if ((strlen(holdstring) > 0) && (com == 13))
      {
	 strcpy(filelocation, filepath);
	 strcpy(filename, holdstring);
	 strcat(filename, ".SET");
	 strcat(filelocation, filename);
	 if (l == 5) loadsettings(1);
	 if (l == 6) savesettings();
	 printerrmsg();
      }

   }

   if ((l == 7) || (l == 8))
   {
      if ((dir == 3) && (timestep <= 1) && (timestep > 0.0625)) timestep /= 2;
      if ((dir == 3) && (timestep > 1)) timestep -= 1;
      if ((dir == 3) && (timestep > 10)) timestep -= 4;
      if ((dir == 4) && (timestep >= 10) && (timestep < 200)) timestep += 4;
      if ((dir == 4) && (timestep >= 1) && (timestep < 200)) timestep += 1;
      if ((dir == 4) && (timestep < 1)) timestep *= 2;
   }
   if (l == 9) changedelay(dir - 3);
   if (l == 10) monitor = abs(!monitor);
   if (l == 11)
   {
      if (dir == 3) redraw --;
      if (dir == 4) redraw ++;
      if (redraw == 3) redraw = 0;
      if (redraw == -1) redraw = 2;
   }
   if (l == 12)
   {
      if (dir == 3) drawfreq -= (4 * (drawfreq > 25) + 1) * (drawfreq > 1);
      if (dir == 4) drawfreq += (4 * (drawfreq >= 25) + 1) * (drawfreq < 100);
   }
   if (l == 13)
   {
      if (dir == 3) style --;
      if (dir == 4) style ++;
      if (style == 5) style = 0;
      if (style == -1) style = 4;
   }
   if (l == 14)
   {
      if (dir == 3) track --;
      if (dir == 4) track ++;
      if (track == 4) track = 0;
      if (track == -1) track = 3;
   }
   if (l == 15) doppler = abs(!doppler);
   if (l == 16) fade = abs(!fade);
   if (l >= 17)
   {
      IFgetsellabel(holdstring, l);
      liminput(holdstring, 10, 1, 0, 1);
      sduor ++;  smuor ++;  stuor ++;  gor ++;

      if (l == 17) {sdu = strtod(holdstring, 0); sduor = 0;}
      if (l == 18) {smu = strtod(holdstring, 0); smuor = 0;}
      if (l == 19) {stu = strtod(holdstring, 0); stuor = 0;}
      if (l == 20) {G = strtod(holdstring, 0); gor = 0;}

      if (sdu == 0) sdu = 1e+19;
      if (smu == 0) smu = 1e+39;
      if (stu == 0) stu = 1e+13;
      if (G == 0) G = 6.673e-3;


      if ((gor > sduor) && (gor > smuor) && (gor > stuor)) G = Gphys / pow(sdu, 3) * smu * pow(stu, 2);
      if ((stuor > sduor) && (stuor > smuor) && (stuor > gor)) stu = pow((G / Gphys * pow(sdu, 3) / smu), (0.5));
      if ((smuor > sduor) && (smuor > stuor) && (smuor > gor)) smu = G / Gphys * pow(sdu, 3) / pow(stu, 2);
      if ((sduor > smuor) && (sduor > stuor) && (sduor > gor)) sdu = pow((smu * pow(stu, 2) * Gphys / G), (0.33333333));

   }
   return;
}
/****************************************************************************/
void liminput(char * inputstring, int lim, int digok, int alpok, int buf)
{
   //Special input function to recieve limited string input from keyboard

   int pos;                //cursor position in string
   int inx;                //initial cursor x position on screen
   int iny;                //initial cursor y position on screen

   inx = wherex();
   iny = wherey();

   gotoxy(inx - 1, iny);
   for (pos = 0; pos <= lim; pos ++) putch(' ');
   if (buf) putch(' ');

   pos = strlen(inputstring);

   do
   {
      gotoxy(inx, iny); cprintf("%s", inputstring);

      if (pos < lim) {gotoxy(pos + inx, iny); putch('_');}

      do
      {
	 com = getch();
      } while (
      !(  ( islower(com) || isdigit(com) || (com == '_') ) && alpok )
      &&
      !(  ( isdigit(com) || (com == 'e') || (com == '.') || (com == '+') || (com == '-') ) && digok )
      &&
      !(  ( (com == 27) || (com == 13) || ((com >= 59) && (com <= 68)) || (com == 80) || (com == 72) || (com == 8) )  )
      );

      if (!(  ( (com == 27) || (com == 13) || ((com >= 59) && (com <= 68)) || (com == 80) || (com == 72) || (com == 8) )  ))
      {
	 if (islower(com) && (alpok == 1)) com -= 32;
	 if (pos == lim) pos --;
	 inputstring[pos] = com;
	 inputstring[pos + 1] = '\0';
	 pos ++;
	 if (isupper(com) && (alpok == 1)) com += 32;
      }

      if ((com == 8) && (pos > 0))
      {
	 inputstring[pos] = '\0';
	 gotoxy(pos + inx, iny); if (pos < lim) putch(' ');
	 pos --;
	 inputstring[pos] = '\0';
      }

      if ((com == 27) || ((com >= 59) && (com <= 68)) || (com == 80) || (com == 72)) skip = 1;

   } while (!( (com == 27) || (com == 13) || ((com >= 59) && (com <= 68)) || (com == 80) || (com == 72)));

   return;
}
/****************************************************************************/
int checkretro(void)
{
   //Checks validity of retroprojection parameters

   rerr = 0;
   if (((gpo[1] != 0) || (gpmax[1] != 0)) && (grmmax[1] == 0)) rerr = 19;
   if (((gpo[2] != 0) || (gpmax[2] != 0)) && (grmmax[2] == 0)) rerr = 20;
   if ((gpo[1] == 0) && (gpmax[1] == 0) && (grmmax[1] != 0)) rerr = 21;
   if ((gpo[2] == 0) && (gpmax[2] == 0) && (grmmax[2] != 0)) rerr = 22;
   if ((gm[1] == 0) && ((gpo[1] > 0) || (gpmax[1] > 0) || (grmmax[1] > 0))) rerr = 23;
   if ((gm[2] == 0) && ((gpo[2] > 0) || (gpmax[2] > 0) || (grmmax[2] > 0))) rerr = 24;
   if ((igdis <= 1) || (igecc < 0) || (gm[1] <= 0) || (gm[2] <= 0)) return 0;
   return 1;
}
/****************************************************************************/
int checkdisk(void)
{
   //Checks validity of particle disk parameters

   err = 0;

   //major errors
   if ((gm[1] == 0) && (gp[1] > 0)) err = 8;
   if ((gm[2] == 0) && (gp[2] > 0)) err = 9;
   if ((gring[1] > 0) && (grmin[1] == 0)) err = 10;
   if ((gring[2] > 0) && (grmin[2] == 0)) err = 11;
   if ((grmax[1] < grmin[1]) || ((gring[1] > 1) && (grmin[1] == grmax[1]))) err = 12;
   if ((grmax[2] < grmin[2]) || ((gring[2] > 1) && (grmin[2] == grmax[2]))) err = 13;
   if ((gp[1] > 0) && (gring[1] == 0)) err = 14;
   if ((gp[2] > 0) && (gring[2] == 0)) err = 15;
   if (abs(gm[1] > 0) + abs(gm[2] > 0) + gring[1] * gp[1] + gring[2] * gp[2] > maxp)
   {
      if ((gring[1] * gp[1]) > (gring[2] * gp[2])) err = 16;
      else err = 17;
   }
   if ((gm[1] > 0) && (gm[2] > 0) && (gx[1] == gx[2]) && (gy[1] == gy[2]) && (gz[1] == gz[2])) err = 18;
   if (((gpo[1] != 0) || (gpmax[1] != 0)) && (grmmax[1] == 0)) err = 19;
   if (((gpo[2] != 0) || (gpmax[2] != 0)) && (grmmax[2] == 0)) err = 20;
   if ((gpo[1] == 0) && (gpmax[1] == 0) && (grmmax[1] != 0)) err = 21;
   if ((gpo[2] == 0) && (gpmax[2] == 0) && (grmmax[2] != 0)) err = 22;
   if ((gm[1] == 0) && ((gpo[1] > 0) || (gpmax[1] > 0) || (grmmax[1] > 0))) err = 23;
   if ((gm[2] == 0) && ((gpo[2] > 0) || (gpmax[2] > 0) || (grmmax[2] > 0))) err = 24;

   //minor errors (autofix)
   if (!err)
   {
      if (gp[1] == 0) {grmin[1] = grmax[1] = gring[1] = gtilt[1] = gskew[1] = 0;}
      if (gp[2] == 0) {grmin[2] = grmax[2] = gring[2] = gtilt[2] = gskew[2] = 0;}
   }


   return err;
}
/****************************************************************************/
int checkprop(void)
{
   //Checks validity of metaproperty parameters
   err = 0;

   if (((source == 1) || (store == 2)) && (!strlen(simname))) err = 1;
   if (source && store) err = 2;
   if ( ((track == 1) && (gm[1] == 0)) || ((track == 2) && (gm[2] == 0)) || ((track == 3) && ((gm[1] == 0) || (gm[2] == 0)))) err = 3;
   if ((source) && (!err))
   {
      strcpy(filelocation, filepath);
      strcpy(filename, simname);
      strcat(filename, ".INT");
      strcat(filelocation, filename);
      fptr = fopen(filelocation, "r");
      if (!fptr) err = 4; else fclose(fptr);
   }
   if ((nump == 0) && (source == 0)) err = 5;
   if (doppler && fade) err = 6;

   return err;
}

/****************************************************************************/
void errmsg(char * errorlabel)
{
   //Copies appropriate error message into memory

   if (err == 0) strcpy(errorlabel, "");
   if (err == 1) {strcpy(errorlabel, "no file name"), errl = 1;}
   if (err == 2) {strcpy(errorlabel, "input and output cannot both be disk"), errl = 3;}
   if (err == 3) {strcpy(errorlabel, "cannot track nonexistent galaxy"), errl = 14;}
   if (err == 4) {strcpy(errorlabel, "file not found"), errl = 1;}
   if (err == 5) {strcpy(errorlabel, "simulation not ready: change particles"), errl = 0;}
   if (err == 6) {strcpy(errorlabel, "multiple color effects enabled"), errl = 16;}
   if (err == 7) {strcpy(errorlabel, "settings file not found"), errl = 4;}
   if (err == 8) {strcpy(errorlabel, "galaxy with disk is massless"), errl = 7;}
   if (err == 9) {strcpy(errorlabel, "galaxy with disk is massless"), errl = 21;}
   if (err == 10) {strcpy(errorlabel, "inner ring radius is zero"), errl = 8;}
   if (err == 11) {strcpy(errorlabel, "inner ring radius is zero"), errl = 22;}
   if (err == 12) {strcpy(errorlabel, "outer ring within inner ring"), errl = 9;}
   if (err == 13) {strcpy(errorlabel, "outer ring within inner ring"), errl = 23;}
   if (err == 14) {strcpy(errorlabel, "no rings"), errl = 10;}
   if (err == 15) {strcpy(errorlabel, "no rings"), errl = 24;}
   if (err == 16) {strcpy(errorlabel, "number of particles exceeds limit"), errl = 11;}
   if (err == 17) {strcpy(errorlabel, "number of particles exceeds limit"), errl = 25;}
   if (err == 18) {strcpy(errorlabel, "galaxies have same location"), errl = 1;}
   if (err == 19) {strcpy(errorlabel, "density not applicable under point-mass"); errl = 34;}
   if (err == 20) {strcpy(errorlabel, "density not applicable under point-mass"); errl = 37;}
   if (err == 21) {strcpy(errorlabel, "massive disk has zero density"); errl = 34;}
   if (err == 22) {strcpy(errorlabel, "massive disk has zero density"); errl = 37;}
   if (err == 23) {strcpy(errorlabel, "model not applicable to massless galaxy"); errl = 7;}
   if (err == 24) {strcpy(errorlabel, "model not applicable to massless galaxy"); errl = 21;}

   if (err == 32) {strcpy(errorlabel, "orbit decay:  increase velocity"); errl = 30;}


   return;
}
/****************************************************************************/
void printfunct(int f1, int f2, int f3, int f4, int f5)
{
   //Prints available interface functions at bottom of screen
   textbackground(0);
   gotoxy(2, 24);
   if (f1 == 0){ textcolor(8); cprintf("F1 - Check"); }
   if (f1 == 1){ textcolor(7); cprintf("F1 - Check"); }
   if (f1 == 2){ textcolor(7); cprintf("F1 - Run  "); }
   gotoxy(14, 24);
   if (f2 == 0){ textcolor(8); cprintf("F2 - Galaxies"); }
   if (f2 == 1){ textcolor(8); cprintf("F2 - Particles"); }
   if (f2 == 2){ textcolor(8); cprintf("F2 - Properties"); }
   if (f2 == 3){ textcolor(7); cprintf("F2 - Galaxies"); }
   if (f2 == 4){ textcolor(7); cprintf("F2 - Particles"); }
   if (f2 == 5){ textcolor(7); cprintf("F2 - Properties"); }
   gotoxy(31, 24);
   if (f3) textcolor(7); else textcolor(8); cprintf("F3 - Exit");
   gotoxy(43, 24);
   if (f4) textcolor(7); else textcolor(8); cprintf("F4 - Retroproject");
   gotoxy(63, 24);
   if (f5) textcolor(7); else textcolor(8); cprintf("F5 - Catalog");
   textcolor(7);

   return;
}

/*********************************************************************/
double fix(double rnum, int ndig)
{
   //Rounds to decimal place specified

   float rem;               //decimal portion of rnum

   rnum *= pow(10, ndig);
   rem = fmod(rnum, 1);
   if (((rem > 0) && (rem < 0.5)) || (rem <= -0.5)) rnum = floor(rnum); else rnum = ceil(rnum);
   rnum /= pow(10, ndig);
   if (rnum == 0) rnum = fabs(rnum);
   return rnum;
}

/*********************************************************************/
void printerrmsg(void)
{
   //Prints error message at bottom of interface screen

   char holderrmsg[48];    //string to hold error message

   gotoxy(10, 23); cprintf("                                                  ");
   errmsg(holderrmsg);
   if (err) {gotoxy(10, 23); printf("%d: %s", err, holderrmsg);}

   return;
}
/********************************************************************/
void savesettings(void)
{
   //Saves current interface parameters to disk

   fptr = fopen(filelocation, "w");
   if (!fptr) err = 7;
   else
   {
      //sdu = round(stu, 5); smu = round(smu, 5); sdu = round(smu, 5); G = round(G, 5);
      err = 0;
      if (strlen(simname) == 0) strcpy(simname, "-");//null char
      fprintf(fptr, "%s\n", simname);
      fprintf(fptr, "%d %d %d %f %f %ld %d %d %d %d %d %d %d %e %e %e %e\n", source, store, recfreq, timestep, timestep, delay, monitor, redraw, drawfreq, style, track, doppler, fade, 0, 0, 0, 0);
      fprintf(fptr, "%f %f %f %f %f %f %f %d %d %d %d %d %d %d\n", gx[1], gy[1], gz[1], gvx[1], gvy[1], gvz[1], gm[1], grmin[1], grmax[1], gring[1], gp[1], gtilt[1], gskew[1], gc[1]);
      fprintf(fptr, "%f %f %f %f %f %f %f %d %d %d %d %d %d %d\n", gx[2], gy[2], gz[2], gvx[2], gvy[2], gvz[2], gm[2], grmin[2], grmax[2], gring[2], gp[2], gtilt[2], gskew[2], gc[2]);
      fprintf(fptr, "%f %f %f %f %f %d %d %d %d %d %d %d %d %d\n", igdis, igecc, iginc, igarg, igt0, gpo[1], gpmax[1], grmmax[1], gpo[2], gpmax[2], grmmax[2], 0, 0, 0);
      fclose(fptr);
      if (simname[0] == '-') strcpy(simname, "");
   }
   return;
}
/****************************************************************************/
void loadsettings(int loadall)
{
   //Loads previously saved interface parameters from disk

   int dummy; long dumml;      //dummy variables to hold nonloaded settings

   fptr = fopen(filelocation, "r");
   if (!fptr) err = 7;
   else
   {
      err = 0;
      fscanf(fptr, "%s\n", &simname);
      if (loadall) fscanf(fptr, "%d %d %d %f %f %ld %d %d %d %d %d %d %d ", &source, &store, &recfreq, &timestep, &timestep, &delay, &monitor, &redraw, &drawfreq, &style, &track, &doppler, &fade);
      else         fscanf(fptr, "%d %d %d %f %f %ld %d %d %d %d %d %d %d ", &dummy,  &dummy, &recfreq, &timestep, &timestep, &dumml, &dummy,   &dummy,  &dummy,    &dummy, &dummy, &dummy,   &dummy);
      fscanf(fptr, "%e %e %e %e ", &sdu, &smu, &stu, &G);
      fscanf(fptr, "%f %f %f %f %f %f %f %d %d %d %d %d %d %d ", &gx[1], &gy[1], &gz[1], &gvx[1], &gvy[1], &gvz[1], &gm[1], &grmin[1], &grmax[1], &gring[1], &gp[1], &gtilt[1], &gskew[1], &gc[1]);
      fscanf(fptr, "%f %f %f %f %f %f %f %d %d %d %d %d %d %d ", &gx[2], &gy[2], &gz[2], &gvx[2], &gvy[2], &gvz[2], &gm[2], &grmin[2], &grmax[2], &gring[2], &gp[2], &gtilt[2], &gskew[2], &gc[2]);
      fscanf(fptr, "%f %f %f %f %f %d %d %d %d %d %d", &igdis, &igecc, &iginc, &igarg, &igt0, &gpo[1], &gpmax[1], &grmmax[1], &gpo[2], &gpmax[2], &grmmax[2]);
      fclose(fptr);
      if (simname[0] == '-') strcpy(simname, "");
      loadselflag = 1;
   }
   return;
}
/****************************************************************************/
void tracking(void)
{
   //Fixes camera on specified object

   if ((track == 1) || (track == 2))
   {
      ptrack = track;
      centx = x[ptrack];
      centy = y[ptrack];
      centz = z[ptrack];
   }
   if (track == 3)
   {
      centx = (x[1] * m[1] + x[2] * m[2]) / (m[1] + m[2]);
      centy = (y[1] * m[1] + y[2] * m[2]) / (m[1] + m[2]);
      centz = (z[1] * m[1] + z[2] * m[2]) / (m[1] + m[2]);
   }
   return;
}

/****************************************************************************/
int dopplercolor(int p1)
{
   //Assigns new color to particle depending on line-of-sight velocity

   float vs;                       //line-of-sight velocity

   transx0 = vx[p1];
   transy0 = vy[p1];
   transz0 = vz[p1];
   transform(rota, rotb, rotc);
   vs = transz2 - restvel / 100;

			      color = 4;   // red
   if (vs > -0.04 * contrast) color = 12;  // light red
   if (vs > -0.03 * contrast) color = 14;  // yellow
   if (vs > -0.02 * contrast) color = 10;  // light green
   if (vs > -0.01 * contrast) color = 2;   // green
   if (vs > 0.00)             color = 3;   // cyan
   if (vs > 0.01 * contrast)  color = 9;   // light blue
   if (vs > 0.02 * contrast)  color = 1;   // blue
   if (vs > 0.03 * contrast)  color = 5;   // purple
   if (vs > 0.04 * contrast)  color = 13;   // light purple

   if (c[p1] == 15) color = 15; if (c[p1] == 0) color = 0;

   return color;
}

/****************************************************************************/
int fadecolor(int p1)
{
   //Assigns new color to particle depending on line-of-sight depth

   float dz;            //line of sight distance

   transx0 = x[p1];
   transy0 = y[p1];
   transz0 = z[p1];
   transform(rota, rotb, rotc);
   dz = transz2 - restvel * 10;

			    color = 4;   // red
   if (dz > -40 * contrast) color = 12;  // light red
   if (dz > -30 * contrast) color = 14;  // yellow
   if (dz > -20 * contrast) color = 10;  // light green
   if (dz > -10 * contrast) color = 2;   // green
   if (dz > 0)              color = 3;   // cyan
   if (dz > 10 * contrast)  color = 9;   // light blue
   if (dz > 20 * contrast)  color = 1;   // blue
   if (dz > 30 * contrast)  color = 5;   // purple
   if (dz > 40 * contrast)  color = 13;  // light purple

   if (c[p1] == 15) color = 15; if (c[p1] == 0) color = 0;

   return color;
}
/****************************************************************************/
void adjcolors(int cdir)
{
   //Adjusts contrast or rest position for color functions

   if ((cdir == 0) && (contrast >= 6)) contrast -= 0.5;
   if ((cdir == 0) && (contrast > 2)) contrast -= 0.5;
   if ((cdir == 1) && (contrast < 50)) contrast += 1.0;
   if ((cdir == 1) && (contrast < 6)) contrast -= 0.5;
   if ((cdir == 2) && (restvel > -20)) restvel -= 1;
   if ((cdir == 3) && (restvel < 20)) restvel += 1;
   drawallparticles(0, 1);
   return;
}
/****************************************************************************/
void kcalc(void)
{
   //Calculates k constants to be used in acceleratory field equation

   int galaxy;        //galaxy ID

   gotoxy(1, 1);
   for (galaxy = 1; galaxy <= 2; galaxy ++)
   {
      if ((grmmax[galaxy] >= 0) && ((gpmax[galaxy] > 0) || (gpo[galaxy] > 0)))
      {
	 ki[galaxy] = ((float)gpmax[galaxy] - (float)gpo[galaxy]) / (3 * (float)grmmax[galaxy]);
	 kj[galaxy] = (float)gpo[galaxy] / 2;
	 kk[galaxy] = pow((float)grmmax[galaxy], 2) * ((float)gpmax[galaxy] / 3 + (float)gpo[galaxy] / 6);
      }
      else {ki[galaxy] = kj[galaxy] = kk[galaxy] = 0;}
   }

   return;
}
/****************************************************************************/
void printtime(void)
{
   int n;             //miscellaneous counting variable

   for (n = 0; n < 7; n ++) capsname[n] = toupper(simname[n]);
   gotoxy(11, 1); printf("%s", capsname);
   gotoxy(11, 2); printf("%.0f  ", t);
   return;
}
/****************************************************************************/