//                                  SSMS.C
//                       Solar System Motion Simulator

//  This is an early version of the IGMS, in the first stages of development.
//  The simulator "engine", which is now simply in main() rather than in
//  a separate core function, is essentially the same as it is in future
//  versions of the simulator.  Also remaining are the display features (zoom,
//  3D rotation capabilities, print-mode, tracking, etc.)
//
//  You will notice that not long into the simulation a 0.5 Msun star
//  (particle #10, "Wormwood" in remarks) comes careening into the Solar
//  System and captures Earth.  Modifying the parameters of this particle
//  make for some interesting variations on this theme.  Also, try setting
//  ptrack = 4 to lock on to the Earth particle to see a demonstration of
//  planetary retrograde motion.
//
//  The EGAVGA.BGI file must be included in the same directory as this file
//  for graphics to function.
//
//  Units are not the same as in IGMS.
//  1 sdu = 10^10 m
//  1 smu = 10^29 kg
//  1 stu = 10^5  s

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

#define G 0.000667           //universal gravitational constant
#define maxp 10              //maximum particles supportable
#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 * 0.0174532   //rotation increment angle

void input(void);                //finds user input and routes appropriately
void drawparticles(int gorder, int style);  //draw or erase all particles
void pause(void);                //holds simulation during user pause
void imagezoom(int zdir);        //0 = zoom in;  1 = zoom out
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

float t;                   //time
float timestep = 1;        //time increment
int nump = maxp;           //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 = 1;            //tracked particle ID
int exitflag = 0;          //flags exit command
int pauseflag = 0;         //flags pause
int displayflag = 0;       //flags display info
long delay = 0;            //steps in timing loop delay
int tracers = 1;           //flags tracers
int particlestyle = 0;	   //appearance of particles 0=point 1=cross 2=circle
int drawfreq = 10;         //redraw particles every drawfreq timesteps

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 color
int c2[maxp + 1];          //tracer color
float gx[maxp + 1];        //graphical x-coordinate
float gy[maxp + 1];        //graphical y-coordinate


int main(void)
{
   int p1, p2;                //particle IDs
   float dx, dy, dz;          //distance between particles along axes
   float D;                   //resultant distance between particles
   float A;                   //resultant acceleration
   long n, a;                 //counting variables
   float maxt = 200000;       //auto-termination time
			      //increase this to make the simulation last longer

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


   x[0] = y[0] = 0;
   if (drawfreq < timestep) drawfreq = timestep;
   if (drawfreq % (int)timestep != 0) drawfreq += abs((drawfreq % (int)timestep) - (int)timestep);



   x[1] = 0;   y[1] = 0;  vx[1] = 0; vy[1] = -0.00000015; m[1] = 19.89; c[1] = 14; c2[1] = 8;  //Sun

   x[2] = x[1] - 5.7;    m[2] = 0.00000330;  c[2] = 13;  //Mercury
   x[3] = x[1] - 10.8;   m[3] = 0.0000487;   c[3] = 14;  //Venus
   x[4] = x[1] + 14.9;   m[4] = 0.00005974;  c[4] = 9;   //Earth
   x[5] = x[1] + 22.8;   m[5] = 0.00000642;  c[5] = 4;   //Mars
   x[6] = x[1] - 77.8;   m[6] = 0.01899;     c[6] = 7;   //Jupiter
   x[7] = x[1] + 142.4;  m[7] = 0.00568;     c[7] = 14;  //Saturn
   x[8] = x[1] + 287.2;  m[8] = 0.000866;    c[8] = 3;   //Uranus
   x[9] = x[1] + 449.9;  m[9] = 0.00103;     c[9] = 3;   //Neptune

   x[10] = 57; y[10] = 760; vx[10] = 0; vy[10] = -0.022; m[10] = 10; c[10] = 15; c2[10] = 8;       //Wormwood

   for (n = 2; n < 10; n ++) //Sets up planets in circular orbits
   {
      vx[n] = 0; y[n] = 0; c2[n] = 8;
      vy[n] = (sqrt(G * m[1] / fabs(x[n] - x[1])));
      if (x[n] > x[1]) vy[n] = -vy[n];
   }

   for (n = 1; n <= 10; n ++) z[n] = vz[n] = 0;



   t = 0;
   do
   {

      if ((t == (long)t) && ((long)t % drawfreq == 0)) drawparticles(0, particlestyle);

      if (ptrack > 0) {centx = x[ptrack] + offx; centy = y[ptrack] + offy;}    //track

      if ((t == (long)t) && ((long)t % drawfreq == 0)) drawparticles(1, particlestyle);

      if ((int)t == floor(t)); gotoxy(11,1); printf("%.0f", t);   //print time


      for (p1 = 1; p1 <= nump; p1 ++)
      {
	 ax[p1] = ay[p1] = 0;
	 for (p2 = 1; p2 <= nump; p2 ++)            //determine acceleration
	 {
	    if (p1 != p2)
	    {
	       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));
	       A = G * m[p2] / pow(D, 2);
	       ax[p1] += dx * A / D;
	       ay[p1] += dy * A / D;
	       az[p1] += dz * A / D;
	    }
	 }
      }

      for (n = 1; n < delay; n ++);                            //delay


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

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

      if (kbhit()) input();

      if (!exitflag) t += timestep;
      if (t >= maxt) exitflag = 1;

   } while (!exitflag);

   gotoxy(38, 1); printf("simulation terminated by user at %.0f", t);
   getch();
   closegraph();
   return 0;
}

/*******************************************************************/
void input(void)
{
   com = getch();
   //gotoxy(1,2); printf("%c = %d", com, com);
   if (com == 27) exitflag = 1; else exitflag = 0;
   if ((com == '-') || (com == '_')) changedelay(0);
   if ((com == '+') || (com == '=')) changedelay(1);
   if ((com == 'i') || (com == 'I')) imagezoom(0);
   if ((com == 'o') || (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 == ',') || (com == '<')) imagerotate(0);
   if ((com == '.') || (com == '>')) imagerotate(1);
   if ((com == 'l') || (com == 'L')) imagerotate(2);
   if ((com == ';') || (com == ':')) imagerotate(3);
   if  (com == 'p')                  imagerotate(4);
   if ((com == '[') || (com == '{')) imagerotate(5);
   if ((com == 'c') || (com == 'C')) {clearviewport(); drawparticles(1, particlestyle);}
   if ((com == ' ') && (!pauseflag)) pause();
   if ((com == '/') || (com == '?')) {displayflag = 1; pause();}
   if ((com == 't') || (com == 'T')) tracers = !tracers;
   if ((com == 'a') || (com == 'A')) printmode();
   return;
}

/*******************************************************************/
void drawparticles(int gorder, int style)
{

   int p1;
   int color;
   //float offax, offbx, offcx;
   //float offay, offby, offcy;
   //float offaz, offbz, offcz;
   float transx0, transy0, transz0;
   float transx1, transy1, transz1;
   float transx2, transy2, transz2;

   for (p1 = 1; p1 <= nump; p1++)
   {
     if (gorder == 1)
     {
	transx0 = x[p1] - centx;
	transy0 = y[p1] - centy;
	transz0 = z[p1] - centz;

	transx1 = transx0 * cos(rota) + transy0 * sin(rota);
	transy1 = transy0 * cos(rota) - transx0 * sin(rota);

	transx2 = transx1 * cos(rotb) + transz0 * sin(rotb);
	transz1 = transz0 * cos(rotb) - transx1 * sin(rotb);

	transy2 = transy1 * cos(rotc) + transz1 * sin(rotc);
      //transz2 = transz1 * cos(rotc) + transy1 * sin(rotc); //not used


	gx[p1] = zoom * (transx2) + sccx;
	gy[p1] = zoom * (transy2) + sccy;

	color = c[p1];
	if (getbkcolor() == 15)
	{
	  if (color == 15) color = 0;
	  //if (color == 14)) color = 6;
	  if (color == 7) color = 0;
	}
      }
      if (gorder == 0) color = c2[p1];
      if ((gorder == -1) || (gorder == 0) && (tracers == 0)) color = 0;
      if (gorder == -2) color = 15;

      //draw
      if (style == 0) putpixel (gx[p1], gy[p1], color);
      if (style == 1)
      {
	 putpixel (gx[p1], gy[p1], color);
	 putpixel (gx[p1] + 1, gy[p1], color);
	 putpixel (gx[p1] - 1, gy[p1], color);
	 putpixel (gx[p1], gy[p1] + 1, color);
	 putpixel (gx[p1], gy[p1] - 1, color);
      }
      if ((style == 2) || (style == 3))
      {
	 setcolor(color);
	 setfillstyle(1, color);
	 circle(gx[p1], gy[p1], 2);
	 if (style == 3) floodfill(gx[p1], gy[p1], color);
      }
   }
   return;
}

/*******************************************************************/
void pause(void)
{
   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)
{
   if ((!zdir) && (zoom < 1000)) zoom *= zinc;
   if ((zdir) && (zoom > 0.01)) zoom /= zinc;
   clearviewport();
   drawparticles(1, particlestyle);
   return;
}

/*******************************************************************/
void imageshift(int sdir)
{
   if (sdir == 1) offy -= sinc / zoom;
   if (sdir == 2) offx -= sinc / zoom;
   if (sdir == 3) offx += sinc / zoom;
   if (sdir == 4) offy += sinc / zoom;

   centx = x[ptrack] + offx; centy = y[ptrack] + offy;

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

   drawparticles(-1, particlestyle);
   drawparticles(1, particlestyle);

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

   return;
}
/*******************************************************************/
void changedelay(int ddir)
{
   if (ddir)
   {
      if ((delay > 0) && (delay < 4000000)) delay *= 2;
      if (delay == 0) delay = 100;
   }
   else
   {
      if (delay > 100) delay /= 2;
      if (delay == 100) delay = 0;
   }
   return;
}
/*******************************************************************/
void displayinfo(void)
{
   if (timestep >= 1)
   {
      gotoxy(1, 1); printf("    time: %.0f      ", t);
      gotoxy(1, 3); printf("timestep: %.0f      ", timestep);
   }
   else
   {
      gotoxy(1, 1); printf("    time: %.3f      ", t);
      gotoxy(1, 3); printf("timestep: %.3f      ", timestep);
   }
   gotoxy(1, 2); printf("          %.3f yr   ", t * 0.00031688);
   gotoxy(1, 4); printf("   delay: %ld       ", delay);
   gotoxy(1, 5); printf("    zoom: %.0f\%    ", zoom * 100);
   gotoxy(1, 6); printf("x center: %.2f      ", centx);
   gotoxy(1, 7); printf("y center: %.2f      ", centy);
   gotoxy(1, 8); printf("a rotate: %.0f      ", rota / 0.0174532);
   gotoxy(1, 9); printf("b rotate: %.0f      ", rotb / 0.0174532);
   gotoxy(1, 10); printf("c rotate: %.0f      ", rotc / 0.0174532);
   gotoxy(1, 11); printf("tracking: %d        ", ptrack);
   gotoxy(1, 12); printf("drawfreq: %d        ", drawfreq);
   return;
}
/*******************************************************************/
void clearinfo(void)
{
   int n;
   gotoxy(1,1);
   for (n = 1; n <= 12; n++) printf("                   \n");
   return;
}
/*******************************************************************/
void printmode(void)
{
   int printmodestyle;

   if (!particlestyle) printmodestyle = 1;
   else printmodestyle = particlestyle;

   drawparticles(-1, particlestyle);
   setbkcolor(15);
   drawparticles(1, printmodestyle);
   setpalette(15, 0);
   setpalette(8, 7);
   setcolor(15);
   drawparticles(-2, printmodestyle);
   getch();
   setbkcolor(0);
   setpalette(8, 56);
   setpalette(15, 63);

   drawparticles(-1, printmodestyle);

   return;
}
/*******************************************************************/
void imagerotate(int rdir)
{

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

   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 * 0.0174532) rota = 180 * 0.0174532;
   if (rotb < -179 * 0.0174532) rotb = 180 * 0.0174532;
   if (rotc < -179 * 0.0174532) rotc = 180 * 0.0174532;

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

   drawparticles(-1, particlestyle);
   drawparticles(1, particlestyle);

   return;
}