//                                SSMS2.CPP
//                      Solar System Motion Simulator 2.0
//             Source code by Daniel Perley, last updated 16 Feb 2005

//  1 sdu = 10^10 m
//  1 smu = 10^29 kg
//  1 stu = 10^4  s

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

#include "planetps.h"

#define G 0.0006673          //universal gravitational constant
#define AUtoSDU 14.9598      //convert AU to simulation distance unit
#define DEGtoRAD 0.0174532   //convert degrees to radians
#define DAYtoSTU 8.64        //convert days to simulation time unit
#define YRtoDAY 365.242199   //convert years to days
#define maxp 8               //maximum number of particles *minus one*
#define sccx 320             //x-coordinate of center of screen
#define sccy 240             //y-coordinate of center of screen
#define zinc 1.33352         //zoom step factor
#define sinc 1.5             //shift increment
#define rinc 5               //rotation increment angle

class simulation
{
  public:
  simulation(int snump);
  ~simulation();
  void core(void);                 //applies forces to update vel., pos.
  void input(void);                //finds user input and routes appropriately
  void gtrans(int *graphx, int *graphy, float x, float y, float z);
  void drawparticles(int gorder);  //draw or erase all particles
  void draworbits(int color);      //draw theoretical orbits
  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
  double getjd();                  //get current JD
  void getdate(char *datestr, double time);   //get string with DD / MM / YYYY from "t"
  void displayhelp();              //show command keys
  void printintromsg();            //display message in top-right corner
  void printtime();                //print the time and date in top-left
  void clearbuffer();              //clear the keyboard input buffer

  void addobject(double x, double y, double z, double vx, double vy, double vz, double m, int c);
  void addsatellite(int porbit, struct orbel *elements, double m, int c);
  void addsatellitesimp(int porbit, double a, double e, double m, int c);
  void loadplanets(int jd2k);
  int opengraphicsmode(char *mode);
  void runsim();
  void flush();

  private:
  long iter;           //iteration number
  double t;            //time (STU)
  double timestep;     //time increment (STU=10^4s)
  int nump;            //number of particles *minus one*
  char com;            //user input variable
  float centx;         //x-coordinate of effective center of screen
  float centy;         //y-coordinate of effective center of screen
  float centz;         //z-coordinate of effective center of screen
  float offx;          //x-axis distance from center and tracked particle
  float offy;          //y-axis distance from center and tracked particle
  float rota;          //x-y rotation
  float rotb;          //x-z rotation
  float rotc;          //y-z rotation
  float zoom;          //zoom factor
  int ptrack;          //tracked particle ID
  int exitflag;        //flags exit command
  int pauseflag;       //flags pause
  int displayflag;     //flags display info
  long delay;          //steps in timing loop delay
  int tracers;         //flags tracers
  int drawfreq;        //redraw particles every drawfreq timesteps
  long jd2kstart;      //day since 2000/1/1 of start of simulation
  int orbits;          //draw orbits?  (value = color)
  int simnum;          //number of simulations run so far

  double px[maxp + 1];      //x-coordinate
  double py[maxp + 1];      //y-coordinate
  double pz[maxp + 1];      //z-coordinate
  double pvx[maxp + 1];     //x-axis velocity
  double pvy[maxp + 1];     //y-axis velocity
  double pvz[maxp + 1];     //z-axis velocity
  double pm[maxp + 1];      //mass
  int pc[maxp + 1];         //pixel color
  int pc2[maxp + 1];        //tracer color
  int pgx[maxp + 1];        //graphical x-coordinate
  int pgy[maxp + 1];        //graphical y-coordinate
};

simulation::simulation(int restart)
{
   nump = -1;
   timestep = 0.5;
   centx = centy = centz = 0;
   offx = offy = 0;
   rota = rotb = rotc = 0;
   zoom = 1.5;
   ptrack = 0;
   pauseflag = exitflag = displayflag = 0;
   delay = 0;
   tracers = 1;
   drawfreq = 2;

   jd2kstart = 0;
   orbits = 1;
   if (!restart) simnum = 0; else simnum++;
}

simulation::~simulation()
{

}

void simulation::core(void)
{
   //Use Euler's method to update the particle positions after a duration
   // equal to one timestep.

   int p1, p2;                //particle IDs
   double dx, dy, dz;         //distance between particles along axes
   double D;                  //resultant distance between particles
   double A;                  //resultant acceleration
   double pax, pay, paz;      //aceleration along each axis


   for (p1 = 0; p1 <= nump; p1 ++)
   {
      pax = pay = paz = 0.0;
      for (p2 = 0; p2 <= nump; p2 ++)            //determine acceleration
      {
	 if (p1 != p2)
	 {
	    dx = px[p2] - px[p1];
	    dy = py[p2] - py[p1];
	    dz = pz[p2] - pz[p1];
	    D = sqrt(dx*dx + dy*dy + dz*dz);
	    A = G * pm[p2] / (D*D);
	    pax += dx * A / D;
	    pay += dy * A / D;
	    paz += dz * A / D;
	 }
      }
      pvx[p1] += pax * timestep;            //accelerate
      pvy[p1] += pay * timestep;
      pvz[p1] += paz * timestep;
   }

   for (p1 = 0; p1 <= nump; p1++)
   {
      px[p1] += pvx[p1] * timestep;          //move
      py[p1] += pvy[p1] * timestep;
      pz[p1] += pvz[p1] * timestep;
   }
   return;
}

/*******************************************************************/
void simulation::input(void)
{
   com = getch();
   //gotoxy(1,2); printf("%c = %d", com, com);
   if (com == 27 || com == 'q') exitflag = 1; else exitflag = 0;
   if ((com == '-') || (com == '_') || com == 'f') changedelay(0);
   if ((com == '+') || (com == '=') || com == 's') changedelay(1);
   if ((com == 'i') || (com == 'I')) imagezoom(0);
   if ((com == 'o') || (com == 'O')) imagezoom(1);
   if (com == 72) {imagerotate(4);}
   if (com == 75) {imagerotate(0);}
   if (com == 77) {imagerotate(1);}
   if (com == 80) {imagerotate(5);}
   //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();
      if (orbits) draworbits(orbits);
      drawparticles(2);
      printtime();
   }
   if ((com == ' ') && (!pauseflag)) pause();
   if ((com == '?') || ((com == '/' || com == 'h') && displayflag == 1))
   {
      clearinfo();
      displayflag = 2; pause();
   }
   else if ((com == '/' || com == 'h'))
   {
     if (displayflag == 2) clearinfo();
     displayflag = 1; pause();
   }
   if ((com == 't') || (com == 'T')) tracers = !tracers;
   if ((com == 'a') || (com == 'A')) printmode();
   if (com == 'r') {timestep = -timestep;}
   if (com >= '0' && com <= '0'+nump && com <= '9')
   {  orbits = 0;
      clearviewport();
      ptrack = com - '0';
      drawparticles(2);
      printtime();
   }
   if (com == '`' || com == 'm') {clearviewport(); drawparticles(2); ptrack = -1;}
   if (((com == 'b') || (com == 'B')) && (ptrack == 0)) {orbits = 1; draworbits(orbits); drawparticles(2);}
   if ((com == 'v') || (com == 'V')) {orbits = 0; clearviewport(); drawparticles(2);}
   //if ((com == '/') || (com == '?')) displayhelp();

   if (com) clearbuffer();  // this prevents the annoying beeping when you
			    // hold down a key


   return;
}
/*******************************************************************/

void simulation::clearbuffer()
{
  while (kbhit()) {getch();}
  return;
}

/*******************************************************************/
void simulation::gtrans(int *graphx, int *graphy, float x, float y, float z)
{
   float transx0, transy0, transz0;
   float transx1, transy1, transz1;
   float transx2, transy2, transz2;

   float rotarad = rota * DEGtoRAD;
   float rotcrad = rotc * DEGtoRAD;

   transx0 = x - centx;
   transy0 = y - centy;
   transz0 = z - centz;

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

   transy2 = transy1 * cos(rotcrad) + transz0 * sin(rotcrad);
   transz1 = transz0 * cos(rotcrad) - transy1 * sin(rotcrad);

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

   transx2 = transx1; //transy2 = transy1;

   *graphx = (int)(zoom * (transx2) + sccx + offx);
   *graphy = (int)(-zoom * (transy2) + sccy - offy);

   return;
}

/*******************************************************************/
void simulation::drawparticles(int gorder)
{
   int p1;
   int color;
   int style;

   int oldgx, oldgy;
   int newgx, newgy;

   for (p1 = 0; p1 <= nump; p1++)
   {
     if (zoom < 0.3)
     {
	if (pm[p1] > 1) style = 1; else style = 0;
     }
     else if (zoom < 0.7)
     {
	//if (pm[p1] > 1) style = 1;
	if (pm[p1] > 0.0002) style = 1;
	else style = 0;
     }
     else if (zoom < 8)
     {
       if (pm[p1] > 1) style = 2;
       else if (pm[p1]  > 0.000001) style = 1;
       else style = 0;
     }
     else if (zoom < 11)
     {
       if (pm[p1] > 1) style = (0.005*AUtoSDU*zoom > 3 ? 0.005*AUtoSDU*zoom : 3);
       else if (pm[p1] > 0.00001) style = 2;
       else if (pm[p1] > 0.000001) style = 1;
       else style = 0;
     }
     else
     {
       if (pm[p1] > 1) style = (0.005*AUtoSDU*zoom > 3 ? 0.005*AUtoSDU*zoom : 3);
       else if (pm[p1] > 0.000001) style = 2;
       else if (pm[p1] > 0.0000001) style = 1;
       else style = 0;
     }


     oldgx = pgx[p1]; oldgy = pgy[p1];
     if (ptrack >= 0) {centx = px[ptrack]; centy = py[ptrack]; centz = pz[ptrack];}    //track

     gtrans(&newgx,&newgy, px[p1], py[p1], pz[p1]);

     //if ((t < 5*timestep*drawfreq) && (p1 == 0))
     //  {gotoxy(1,++q); printf("%.0f: %d %d %d %d (%d) %.5e", t, oldgx, oldgy, newgx, newgy, gorder, px[p1]);}

     if (oldgx == newgx && oldgy == newgy && gorder < 2) continue;


     if (gorder >= 1)
     {
	//gtrans(&(pgx[p1]),&(pgy[p1]), px[p1], py[p1], pz[p1]);

	pgx[p1] = newgx;  pgy[p1] = newgy;

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

      //draw
      if (style == 0) putpixel (pgx[p1], pgy[p1], color);
      if (style == 1)
      {
	 putpixel (pgx[p1], pgy[p1], color);
	 putpixel (pgx[p1] + 1, pgy[p1], color);
	 putpixel (pgx[p1] - 1, pgy[p1], color);
	 putpixel (pgx[p1], pgy[p1] + 1, color);
	 putpixel (pgx[p1], pgy[p1] - 1, color);
      }
      if (style >= 2)
      {
	 setcolor(color);
	 setfillstyle(1, color);
	 fillellipse(pgx[p1], pgy[p1], style, style);
	 //floodfill(pgx[p1], pgy[p1], color);
      }
   }
   return;
}
/*******************************************************************/
void simulation::draworbits(int color)
{
  int i;
  struct orbel el[9];
  double x, y, z;
  int gx, gy;
  double d;
  double period;
  double pinc;
  double jd2know;
  int prevx, prevy=-100;

  getplanetels(el);

  if (ptrack >= 0) {centx = px[ptrack]; centy = py[ptrack]; centz = pz[ptrack];}    //track
  jd2know = jd2kstart + t / DAYtoSTU;

  setcolor(color);
  for (i = 1; i <= 8; i++)
  {
    period = 364.25 * pow(el[i].a,1.5); //Use Kepler's law to get planet's period
    pinc = 6*(period/80.0) / pow(el[i].a,0.5) / sqrt(zoom);

    for (d = jd2know; d < jd2know + period /*+ 2*pinc*/; d += pinc)
    {
      calcposel(&x, &y, &z, &el[i], d); //inefficient; should store in memory
      gtrans(&gx, &gy, x*AUtoSDU + px[0], y*AUtoSDU + py[0], z*AUtoSDU + pz[0]);
      if (prevy != -100) line(gx,gy,prevx,prevy);            //this is fairly inefficient, too
      prevx = gx; prevy = gy;
    }
    calcposel(&x, &y, &z, &el[i], jd2know);
    gtrans(&gx, &gy, x*AUtoSDU + px[0], y*AUtoSDU + py[0], z*AUtoSDU + pz[0]);
    line(gx,gy,prevx,prevy);
    prevy = -100;
  }

  return;
}


/*******************************************************************/
void simulation::pause(void)
{
   pauseflag = 1;
   do
   {
      if (displayflag == 1) displayinfo();
      if (displayflag == 2) displayhelp();
      input();
   } while ((com != ' ') && (exitflag == 0));

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

/*******************************************************************/
void simulation::imagezoom(int zdir)
{
   if ((!zdir) && (zoom < 1000)) zoom *= zinc;
   if ((zdir) && (zoom > 0.01)) zoom /= zinc;
   clearviewport();
   if (orbits) draworbits(orbits);
   drawparticles(2);
   printtime();
   return;
}

/*******************************************************************/
void simulation::imageshift(int sdir)
{
   int p;

   //if (orbits) draworbits(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();}
   if (tracers || orbits) clearviewport(); else drawparticles(-1);

   if (orbits) draworbits(orbits);
   drawparticles(2);
   printtime();

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

   return;
}
/*******************************************************************/
void simulation::changedelay(int ddir)
{
   if (ddir)
   {
      if (delay < 1000000) delay += 100000;
      else delay += 500000;
   }
   else
   {
     if (delay > 0 && delay <= 1000000) delay -= 100000;
     else if (delay > 1000000) delay -= 500000;
   }
   return;
}
/*******************************************************************/
void simulation::printtime()
{
   char date[11];
   gotoxy(11,1); printf("%.0f", t);
   getdate(date, t);
   gotoxy(25,1); printf("%s", date);
}

/*******************************************************************/
void simulation::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 / DAYtoSTU / YRtoDAY);
   gotoxy(1, 4); printf("   delay: %ld   ", delay/1000);
   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(" x-y rot: %.0f  ", rota);
   //gotoxy(1, 9); printf("b rotate: %.0f  ", rotb / 0.0174532);
   gotoxy(1, 9); printf("xy-z rot: %.0f  ", rotc);
   gotoxy(1, 10); printf("tracking: %d ", ptrack);
   //gotoxy(1, 11); printf("drawfreq: %d        ", drawfreq);
   return;
}
/*******************************************************************/
void simulation::clearinfo(void)
{
   int n;
   gotoxy(1,1);
   if (displayflag == 1)
     for (n = 1; n <= 12; n++) printf("                   \n");
   if (displayflag == 2)
     for (n = 1; n <= 17; n++) printf("                        \n");

   if (orbits) draworbits(1);
   drawparticles(2);

   return;
}
/*******************************************************************/
void simulation::printmode(void)
{
   //int printmodestyle;
   //if (!particlestyle) printmodestyle = 1;
   //else printmodestyle = particlestyle;

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

   drawparticles(-1);

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

   if (tracers == 1) {tracers = 2; clearviewport();}
   else if (!pauseflag || orbits) 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 ) 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 ;

   drawparticles(-1);
   if (orbits) draworbits(orbits);
   drawparticles(2);
   printtime();

   return;
}

double simulation::getjd()
{
   return jd2kstart + t/DAYtoSTU;
}

void simulation::getdate(char *datestr, double time)
{
   long days = 366;
   int day;
   int mon;
   long y = 2000;
   long day2k = (long)getjd();

   if (day2k == 0) {sprintf(datestr,"1999/12/31"); return;}
   if (day2k > 0)
     while (day2k > days)
     {
       //printf("%ld %d\n", day2k, days);
       day2k -= days;
       y++;
       if (y % 4 == 0  && (y % 400 != 0 || y % 2000 == 0)) days = 366;
						    else days = 365;
     }
   if (day2k < 0)
     while (day2k < 1)
     {
       //printf("%ld %d\n", day2k, days);
       y--;
       if (y % 4 == 0  && (y % 400 != 0 || y % 2000 == 0)) days = 366;
						    else days = 365;
       day2k += days;

     }

   if (days == 365)
   {
      if (day2k <= 31) {mon = 1; day = day2k;}
      if (day2k >= 32 && day2k <= 59) {mon = 2; day = day2k - 31;}
      if (day2k >= 60 && day2k <= 90) {mon = 3; day = day2k - 59;}
      if (day2k >= 91 && day2k <= 120) {mon = 4; day = day2k - 90;}
      if (day2k >= 121 && day2k <= 151) {mon = 5; day = day2k - 120;}
      if (day2k >= 152 && day2k <= 181) {mon = 6; day = day2k - 151;}
      if (day2k >= 182 && day2k <= 212) {mon = 7; day = day2k - 181;}
      if (day2k >= 213 && day2k <= 243) {mon = 8; day = day2k - 212;}
      if (day2k >= 244 && day2k <= 273) {mon = 9; day = day2k - 243;}
      if (day2k >= 274 && day2k <= 304) {mon = 10; day = day2k - 273;}
      if (day2k >= 305 && day2k <= 334) {mon = 11; day = day2k - 304;}
      if (day2k >= 335 && day2k <= 365) {mon = 12; day = day2k - 334;}
   }
   if (days == 366)
   {
      if (day2k <= 31) {mon = 1; day = day2k;}
      if (day2k >= 32 && day2k <= 60) {mon = 2; day = day2k - 31;}
      if (day2k >= 61 && day2k <= 91) {mon = 3; day = day2k - 60;}
      if (day2k >= 92 && day2k <= 121) {mon = 4; day = day2k - 91;}
      if (day2k >= 122 && day2k <= 152) {mon = 5; day = day2k - 121;}
      if (day2k >= 153 && day2k <= 182) {mon = 6; day = day2k - 152;}
      if (day2k >= 183 && day2k <= 213) {mon = 7; day = day2k - 182;}
      if (day2k >= 214 && day2k <= 244) {mon = 8; day = day2k - 213;}
      if (day2k >= 245 && day2k <= 274) {mon = 9; day = day2k - 244;}
      if (day2k >= 275 && day2k <= 305) {mon = 10; day = day2k - 274;}
      if (day2k >= 306 && day2k <= 335) {mon = 11; day = day2k - 305;}
      if (day2k >= 336 && day2k <= 366) {mon = 12; day = day2k - 335;}
   }

   sprintf(datestr, "%4ld/%2d/%2d", y, mon, day);
   return;
}

void simulation::displayhelp()
{
   int fl = 1;
   gotoxy(1,fl+0); printf("Esc End simulation");
   gotoxy(1,fl+1); printf("  + Increase delay");
   gotoxy(1,fl+2); printf("  - Decrese delay");
   gotoxy(1,fl+3); printf("  i Zoom in");
   gotoxy(1,fl+4); printf("  o Zoom out");
   gotoxy(1,fl+5); printf("%c%c%c Rotate", 24, 25, 29);
   //gotoxy(1,fl+6); printf(" <> Rotate (ecliptic)");
   //gotoxy(1,fl+7); printf(" p[ Rotate (polar)");
   gotoxy(1,fl+6); printf("  c Clear screen");
   gotoxy(1,fl+7); printf("spc Pause/resume");
   gotoxy(1,fl+8); printf("  / Display info");
   gotoxy(1,fl+9); printf("  t Toggle tracers");
   gotoxy(1,fl+10); printf("  a Use Print Mode");
   gotoxy(1,fl+11); printf("  r Reverse direction");
   gotoxy(1,fl+12); printf("0-9 Track Sun/planet");
   gotoxy(1,fl+13); printf("  ~ Track CM");
   gotoxy(1,fl+14); printf("b/v Display/hide orbits");
}
/********************************************************/

void simulation::printintromsg()
{
   outtextxy(360,2, "Solar System Motion Simulator v2.0");
   outtextxy(340,15,"Press [space] to begin or ? for help.");
   return;
}

/********************************************************/

void simulation::addobject(double x, double y, double z, double vx, double vy,
			   double vz, double m, int c)
{
  if (nump < maxp) nump += 1;
  px[nump] = x;  py[nump] = y;  pz[nump] = z;
  pvx[nump] = vx; pvy[nump] = vy; pvz[nump] = vz;
  pm[nump] = m; pc[nump] = c; pc2[nump] = 8;
  return;
}

void simulation::addsatellite(int porbit, struct orbel *elements, double m, int c)
{
  if (porbit > nump || porbit < 0) return;

  double x, y, z, vx, vy, vz;
  calcphaseel(&x,&y,&z,&vx,&vy,&vz,elements,getjd(),0.05);

  x += px[porbit];  y += py[porbit];  z += pz[porbit];
  vx += pvx[porbit];  vy += pvy[porbit];  vz += pvz[porbit];

  addobject(x,y,z,vx,vy,vz,m,c);
}

void simulation::addsatellitesimp(int porbit, double a, double e, double m, int c)
{
  if (porbit > nump || porbit < 0) return;

  double x, y, z, vx, vy, vz;

  x = a*AUtoSDU / (1+e);
  y = z = 0;
  vy = -sqrt(G * pm[porbit] * (1+e) / ((m/pm[porbit] + 1) * a*AUtoSDU));
  vx = vz = 0;

  x += px[porbit];  y += py[porbit];  z += pz[porbit];
  vx += pvx[porbit];  vy += pvy[porbit];  vz += pvz[porbit];

  addobject(x,y,z,vx,vy,vz,m,c);
}

void simulation::loadplanets(int jd2k)
{
   int p1;
   double vxcm, vycm, vzcm, mt;

   //if (maxp < 8) return;
   if (nump < 8) nump = 8;
   if (nump > maxp) return;

   px[0] = 0; py[0] = 0; pz[0] = 0;
   pvx[0] = 0; pvy[0] = 0; pvz[0] = 0;
   pm[0] = 19.8892; pc[0] = 14; pc2[0] = 8;  //Sun

   pm[1] = 0.00000330;  pc[1] = 13;  //Mercury
   pm[2] = 0.0000487;   pc[2] = 14;  //Venus
   pm[3] = 0.00005974;  pc[3] = 9;   //Earth
   pm[4] = 0.00000642;  pc[4] = 4;   //Mars
   pm[5] = 0.01899;     pc[5] = 7;   //Jupiter
   pm[6] = 0.00568;     pc[6] = 14;  //Saturn
   pm[7] = 0.000866;    pc[7] = 3;   //Uranus
   pm[8] = 0.00103;     pc[8] = 3;   //Neptune

   getplanetphases(px,py,pz,pvx,pvy,pvz,jd2k,0.05);
   jd2kstart = jd2k;

   for (p1 = 1; p1 <= 8; p1++)
   {
     //convert AU to SDU
     px[p1] *= AUtoSDU;
     py[p1] *= AUtoSDU;
     pz[p1] *= AUtoSDU;
     pvx[p1] *= AUtoSDU;
     pvy[p1] *= AUtoSDU;
     pvz[p1] *= AUtoSDU;

     //convert per day to per STU
     pvx[p1] /= DAYtoSTU;
     pvy[p1] /= DAYtoSTU;
     pvz[p1] /= DAYtoSTU;

     pc2[p1] = 8;
   }

   // Change reference frames so the center of mass is stationary
   vxcm = vycm = vzcm = mt = 0;
   for (p1 = 0; p1 <= 8; p1++)
   {
     vxcm += pvx[p1] * pm[p1];
     vycm += pvy[p1] * pm[p1];
     vzcm += pvz[p1] * pm[p1];
     mt += pm[p1];
   }
   vxcm /= mt;
   vycm /= mt;
   vzcm /= mt;
   for (p1 = 0; p1 <= 8; p1++)
   {
     pvx[p1] -= vxcm;
     pvy[p1] -= vycm;
     pvz[p1] -= vzcm;
   }
   //printf("%d: %.2f %.2f %.2f %f %f %f\n", p1, px[p1], py[p1], pz[p1], pvx[p1], pvy[p1], pvz[p1]);

  return;
}

int simulation::opengraphicsmode(char *mode)
{
   int gdriver, gmode, errorcode;
   gdriver = DETECT;
   initgraph(&gdriver, &gmode, mode);
   errorcode = graphresult();

  if (errorcode != grOk)
  {
     printf("Graphics error: %s\n", grapherrormsg(errorcode));
     return errorcode;
  }

  return 0;
}

void simulation::runsim()
{
   t = 0.0;
   iter = 0;
   long n;

   clearviewport();
   if (orbits) draworbits(orbits);
   drawparticles(2);

   do
   {
      if (iter % drawfreq == 0)
      {
	 drawparticles(0);
	 drawparticles(1);
	 if (iter % (drawfreq*10)) drawparticles(2);
	 printtime();
      }

      if (iter == 0)
      {
	 if (!simnum) {setcolor(7); printintromsg();}
	 pause();
	 if (com == 27) com = 'q';
	 if (!simnum) {setcolor(0); printintromsg();}
      }

      core();

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

      if (kbhit()) input();

      if (!exitflag) {iter++; t += timestep;}

   } while (!exitflag);

   gotoxy(42, 1); printf("simulation terminated at %.0f", t);
}

void simulation::flush()
{
   simulation::simulation(1);
}

int main(int argc, char **argv)
{
   char com2;
   int gerr;
   int startjd2k;
   simulation ssystem(0);

   gerr = ssystem.opengraphicsmode("EGAVGA.BGI");
   if (gerr) {getch(); exit(1);}

   while (com2 != 27)
   {
     startjd2k = (argc > 1) ? atoi(argv[1]) : 1;
     ssystem.loadplanets(startjd2k);

     ssystem.runsim();

     gotoxy(37, 2); printf("(press any key to restart or [esc] to quit)");
     ssystem.flush();
     com2 = getch();
   }

   closegraph();

   return 0;
}