/****************************************************************/
/****   illiSol: A Modern Orrery - William Davis Spring '08  ****/
/****    (C) 1994 Board of Trustees University of Illinois   ****/ 
/****    A Model Real-Time Interactive OpenGL Application    ****/
/****   Produced by William Davis, built from the skeleton   ****/
/****            Syzygy Application developed by             ****/
/****    George Francis, Glenn Chappell, Chris Hartman       ****/
/****                 wldavis2@uiuc.edu                      ****/
/******* revised May 2, 2008 by W. Davis ************************/
/****************************************************************/

#include "arPrecompiled.h"  //must be the first line
#include "arMasterSlaveFramework.h"

#define GOVERNOR 5
//#define UNIX
#define WIN32 
//#define DEBUG
//#define VOYAGER_DEBUG
/****************************************************************/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#ifdef WIN32
#include <sys\timeb.h>    /* for windows */ 
#define random rand
#endif
#ifdef UNIX
#include <sys/time.h>    /* for unix */ 
#endif
//#include <glut.h>
#include "KeplerianElements.h"

/*Most of the following was left untouched from the skeleton
 *file, just reformatted to be more readable */
#define  MAX(x,y)         (((x)<(y))?(y):(x))
#define  MIN(x,y)         (((x)<(y))?(x):(y))
#define  ABS(u)           ((u)<0 ? -(u): (u))
#define  FOR(a,b,c)       for((a)=(b);(a)<(c);(a)++)
#define  DOT(p,q)         ((p)[0]*(q)[0]+(p)[1]*(q)[1]+(p)[2]*(q)[2])
#define  NRM(p)           sqrt(DOT((p),(p)))
#define  DG               M_PI/180
#define  S(u)             sin(u*DG)
#define  C(u)             cos(u*DG)
#define  CLAMP(x,u,v)     (x<u? u : (x>v ? v: x))
#define  PI               3.14159265
#define  FLYMODE          (0)
#define  TURNMODE         (1) 
/* GLOBAL VARIABLES */

  /*arguments() should later set default gap0*/
  float gap, gap0=.8; 

  /*World light direction vector*/
  float lux[3]={1.,2.,3.};

  /*Object space direction vector*/
  float luxx[3];

  /*Ambient fraction, pseudo-specular power*/
  float amb, pwr ;

  /*Console Navigation Variables*/
  float mysiz,speed, torq, focal, wfar;
  /*Window code: 2 for full screen, 0 for demand sized*/
  int win = 1;

  /*These are used in chaptrack gluttery*/
  unsigned int BUT, XX, YY,SHIF;

  /*Used for describing the viewport, used to be xt, yt, xm, ym*/
  int xt,yt;
  
  /*More variables that can be toggled*/
  int mode,morph,msg,binoc;

  /*From siz, and caveyes in CAVEskels*/
  int sizz;

  /*Enumerate the different kinds of eyes*/
  enum {console, cave, cube} eyes;

  /*More globals for computations*/  
  int gesture; float eps; float del;

  /*Another variable*/
  int resethand;

  /*Apparently opengl_stars.c requires this*/
  int caveyes=0;

  /*These were part of the old skeleton, and shouldn't be
   *used anymore (they were parameters for the torus). Once
   *everything works, these should be removed just to clean
   *up*/
  int th0, th1, dth, ta0, ta1, dta;

  /*Counter and Temporary Variables - SHOULD NOT BE GLOBAL,
   *but I don't have time to go through and fix that right now*/
  int ii, jj, kk;  float tmp, temp;

  /*Different transformation matrices for use with OpenGL*/
  float affmat[16], starmat[16], mat[16];

  /*Nose-to-eye distance in console*/ 
  float nose;

  
/*My Globals - these are all of my own creation. Originally,
 *a structure was going to be used for each planet and an
 *array of these structures was to be initialized in main,
 *but this led to headaches when trying to draw each planet -
 *namely, I think my pointer logic was iffy. To save time and
 *energy, I just converted everything to a global 2D array
 *to hold the necessary characteristics of each planet. Messy,
 *but it works. */

/*A little structure to hold our Gregorian Date*/
typedef struct timestamp Timestamp;

struct timestamp{
    int year;
    int month;
    int day;
    int hour;
    int minute;
    int second;
};

/*This will hold the texture maps for each planet (i.e. the surface
 *of each planet*/
GLubyte* textures[9];

/*This is the basic template for each planet in GLUT*/
GLUquadricObj* sph;

/*The global 2D double array for planet characteristics.
 * planets[n][0] = holds the radius (in astronomical units)
 *                 of the nth planet
 * planets[n][1] = holds the X coordinate of the nth planet
 *                 in its orbital plane
 * planets[n][2] = holds the Y coordinate of the nth planet
 *                 in its orbital plane
 * planets[n][3] = holds the X coordinate of the nth planet
 *                 in the reference plane
 * planets[n][4] = holds the Y coordinate of the nth planet
 *                 in the reference plane
 * planets[n][5] = holds the Z coordinate of the nth planet
 *                 above the reference plane
 */
double planets[9][6];

/*This is used for the Voyager II's position*/
double voyager[1][3];

/*May be used, maybe not, but will be needed to draw the orbits
 *of each planet*/
double planetsOrbit[8][3];

/*A global for our current Julian Date, and a variable to change it*/
double julianDay, timeChange, tempTimeChange, moonAngle, IoAngle,
       EuropaAngle, GanymedeAngle, CallistoAngle, TitanAngle;

Timestamp* startDate;

/*A couple parameters for the magnification factors of the
 *RADII of the planets and the Sun, respectively */
double solMag, planetMag;
int magOn, pause;
/*The equations used in this program are divided into two time
 *periods: era = 1 is 1850 AD - 2050 AD
 *         era = 2 is 3000 BC - 3000 AD
 */
int era;

/* These are all of our function declarations: to get rid of 
 * warnings and whatnot */
void compute_Position(double julianDate, int planet, int era, int orbit);
void compute_Position_AD(double julianDate, int planet, int orbit);
void compute_Position_BCandAD(double julianDate, int planet, int orbit);
void compute_Voyager_Position();
double convertGregToJDN(Timestamp* startDate);
void updateGregorian();
int getDate(Timestamp* startDate);
void drawSol();
void drawPlanet(int planetNum);
void drawVoyager();
void drawSaturnRings();
void drawMoon();
void drawJupiterMoons();
void drawSaturnMoons();
void drawOrbit(int planet);
void initPlanets();
GLubyte* loadTexture (char *fileName, int width, int height, int depth);
/**********************************************************************
 **********************************************************************
 ****************       Function Definitions        *******************
 **********************************************************************
 **********************************************************************/
bool cubeinit(arMasterSlaveFramework& fw, arSZGClient& )
{
   /*Set initial position and tell Syzygy what we're dealing with*/
   /*First set the initial position*/
   ar_setNavMatrix(ar_translationMatrix(0,-5,0));

   /*Then add all of the fields to the Syzygy framework*/
   fw.addTransferField("affmat", affmat, AR_FLOAT, 16);
   fw.addTransferField("starmat", starmat, AR_FLOAT, 16);
   fw.addTransferField("sizz", (&sizz), AR_FLOAT, 1);
   fw.addTransferField("del", (&del), AR_FLOAT, 1);
   fw.addTransferField("eps", (&eps), AR_FLOAT, 1);
   fw.addTransferField("gesture", (&gesture), AR_INT, 1);
   fw.addTransferField("resethand", (&resethand), AR_INT, 1);
   fw.addTransferField("mode", (&mode), AR_INT, 1);
   fw.addTransferField("luxx", luxx, AR_FLOAT, 3);
   fw.addTransferField("th0", (&th0), AR_INT, 1);
   fw.addTransferField("th1", (&th1), AR_INT, 1);
   fw.addTransferField("ta0", (&ta0), AR_INT, 1);
   fw.addTransferField("ta1", (&ta1), AR_INT, 1);
   fw.addTransferField("gap", (&gap), AR_FLOAT, 1);
   fw.addTransferField("magOn", (&magOn), AR_INT, 1);

   return true;
}

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

/*This was for the torus, but just to keep things working I'm
 *not going to touch it for now*/
void autotymer(int reset)
{ 
  /*This was unreadable before*/
  #define  TYME(cnt,max,act) {static int cnt;if(first) cnt=max; else\
                              if(cnt?cnt--:0){act; goto Break;}}
  static int first = 1;  /* the first time autymer is called */
  if( reset)
    first=1;             /* or if it is reset to start over  */
  TYME(shrink, 150, th0++; 
                    th1--;
                    ta0++;
                    ta1--)
  
  TYME(pause, 20,        )
  
  TYME(grow, 150,   th0--;
                    th1++;
                    ta0--;
                    ta1++)
  
  TYME(dwell, 30,        )
  
  TYME(finish, 1, first=1)
  
  first = 0;
  
  Break:   ;
}

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

void deFault(void)
{ 
  float tmp; 

  th0=5; th1=355;  ta0=5; ta1=355; gap = gap0; 
  msg=1; binoc=0; nose=.06; mode=TURNMODE;  
  speed=.01; torq=.002; focal = .0001; wfar =1000000.; mysiz=10./*.01*/; morph=0; sizz=5; 
  FOR(ii, 0, 16)
    starmat[ii]= affmat[ii] = (ii/4==ii%4);   /* identities */ 
  tmp=NRM(lux); FOR(ii,0,3)lux[ii]/=tmp; amb = 0.3; pwr = 10.0;

  /*Some globals I have added in*/
  era = 0;
  planetMag = 1.0;
  solMag = 1.0;
  timeChange = 2.0;
  tempTimeChange = 0;
  moonAngle = 0.0;
  IoAngle = 25.0;
  EuropaAngle = 56.0;
  GanymedeAngle = 0.0;
  CallistoAngle = 184.0;
  TitanAngle = 0.0;
  pause = 0; 
  
  if(eyes==cave)
  {
    affmat[12]=0;
    affmat[13]= 5;
    affmat[14]= -5;
  }

  if(eyes==cube)
  {
    affmat[12]=0; 
    affmat[13]= 0; 
    affmat[14]=  0; 
  }

  if(eyes==console)
  {
    affmat[12]=0;
    affmat[13]= 0;
    affmat[14]= -4.2;
  }

  /* place where we can see it */
  autotymer(1); /* reset autotymer to start at the beginning */
  gesture=0;
  eps=0.002;    /*hand rotation gesture attenuation*/
  del=0.02;     /*hand translation gesture attenuation*/
  resethand=0;
}

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

void drawvert(int th, int ta)
{  
  float lmb, spec, nn[3], dog, cat;
  nn[0] = C(th)*C(ta); /* unit sphere radius vector */  
  nn[1] = S(th)*C(ta);
  nn[2] =       S(ta);
  lmb = DOT(nn,luxx);lmb =(lmb<0 ? lmb : lmb);  lmb = MAX(amb, lmb); 
  spec = CLAMP((1.1 - pwr+pwr*lmb), 0., 1.);  
  dog = (ta-ta0)/(float)(ta1-ta0); cat = (th-th0)/(float)(th1-th0);
  glColor3f(
      MAX(spec, lmb*dog),           /* map R^2(dog,cat)->R^3(RGBspace */
      MAX(spec, lmb*(.25 + ABS(cat -.5))), /* dog cat model of Hartman*/
      MAX(spec, lmb*(1 - cat)));    /* illiLight by Ray Idaszak 1989  */
  glVertex3f(
      C(th) + .5*nn[0],
      S(th) + .5*nn[1],
      .5*nn[2]);
} 

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

void drawall(void)
{
    glPushMatrix(); 
        int p;
        drawSol();
        for(p = 0; p < 8; p++)
        {
            //if(julianDay > 2443375.5)
            //    drawVoyager();
            drawPlanet(p);

            /*We stil need to get drawOrbit working correctly*/
            /*Unfortunately, the computations for drawing the
             *orbits are far too stressful on the machine...
             *performance was reduced to 3.5 fps */
            //drawOrbit(p);
        }
    glPopMatrix();

    /*Advance our point in time so we can actually see the solar sytem
     *in action - this should be come a changable global */
    if(pause == 1)
    {
        tempTimeChange = timeChange;
        timeChange = 0.0;
    }else if (timeChange == 0.0)
        timeChange = 2.0;

    julianDay += timeChange;
    updateGregorian();

    /*Recompute the positions of every planet for the new updated date*/
    for(p = 0; p < 8; p++)
    {
        #ifdef DEBUG
          cout << "Computing position for planets..." << endl;
        #endif
        compute_Position(julianDay, p + 1, era, 0);
        //compute_Voyager_Position();
    }
}

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

#include "opengl_stars.c"

/**********************************************************************/
/*Developed by Pat Hanrahan, 1989... wow that was the year I was born */
void arguments(int argc,char **argv)
{ 
  //eyes = console;    
  eyes=cube; 
  while(--argc)
  {
    ++argv; 
    if(argv[0][0]=='-')
    {
      switch(argv[0][1])
      {
        case 'w': 
          win =atoi(argv[1]); 
          argv++;
          argc--;
          break; 
        case 'c':
          eyes=console;
          break;
        case 'g':
          gap0 =atof(argv[1]);
          argv++;
          argc--;
          break;
        case 'L':
          lux[0]=atof(argv[1]);
          lux[1]=atof(argv[2]);
          lux[2]=atof(argv[3]);
          argv +=3;
          argc -=3;
          break;
      };
    }
  }
}

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

void keyboard(unsigned char key, int x, int y)
{
    #define  IF(K)            if(key==K)
    #define  PRESS(K,A,b)     IF(K){b;}IF((K-32)){A;}  
    #define  TOGGLE(K,flg)    IF(K){(flg) = 1-(flg); }
    #define  CYCLE(K,f,m)     PRESS((K), (f)=(((f)+(m)-1)%(m)),(f)=(++(f)%(m)))
    #define  SLIDI(K,f,m,M)   PRESS(K,(--f<m?m:f), (++f>M?M:f))
    #define  SLIDF(K,f,m,M,d) PRESS(K,((f -= d)<m?m:f), ((f += d)>M?M:f))
   /* Only ASCII characters can be processes by this GLUT callback function */ 
    IF(27) { exit(0); }              /* ESC exit            */
    PRESS('v', timeChange *= 1.1, timeChange /= 1.1); 
    CYCLE('m', mode,TURNMODE+1);                 /* fly/turn modes      */
    PRESS('h',planetMag += 40.0; solMag += 2.0,
              planetMag -= 40.0; solMag -= 2.0); /* alter magnification */
    PRESS('k',planetMag += 5.0; solMag += 5.0,
              planetMag -= 5.0; solMag -= 5.0);  /* Alter with equal ds */
    TOGGLE('w',msg);                             /* writing on/off      */
    TOGGLE('f', pause);                          /* freeze animation    */
    PRESS('n', nose -= .001 , nose += .001 );    /* for binoculars      */
    PRESS('s',speed /= 1.02, speed *= 1.02);     /* flying speed        */
    PRESS('q',torq /= 1.02, torq *= 1.02);       /* turning speed       */
    PRESS('o', focal *= 1.1 , focal /= 1.1)      /* telephoto           */
    PRESS('i', mysiz /= 1.1, mysiz *= 1.1)       /* rescale the world   */
    PRESS('p', wfar *= 1.01 , wfar   /= 1.01)    /* rear clipping plane */
    PRESS('z', julianDay += 365, julianDay -= 365);            /* zap changes         */
    PRESS('g',gap /= .9, gap *= .9);             /* gap parameter       */
    PRESS('a',amb /= .9, amb *= .9);             /* ambient fraction    */
    PRESS('r',pwr /= .9, pwr *= .9);             /* pseudo-spec power   */
}

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

void clefs(arMasterSlaveFramework& , unsigned char key, int x, int y)
{ 
    keyboard(key, x, y);
}

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

void special_keybo(int key, int x, int y)
{
    /*non-ASCII keypresses go here, if you're lucky enough to know their names */
    fprintf(stderr," non-ASCII character was pressed.\n");
    fprintf(stderr," use special_keybo() to process it\n");
}

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

#ifdef UNIX
    float speedometer(void)
    {
        double dbl;
        static double rate;
        static int ii=0;
        static struct timezone notused; 
        static struct timeval now, then;
        if(++ii % 8 == 0)
        {  
            /* 8 times around measure time */
            gettimeofday(&now, &notused); /* elapsed time */
            dbl =  (double)(now.tv_sec - then.tv_sec) +
                   (double)(now.tv_usec - then.tv_usec)/1000000;
            then = now;  rate = 8/dbl;
        }

        return ((float)rate);
    }
#endif

#ifdef WIN32
    float speedometer(void)
    { 
        double dbl;
        static double rate;
        static int ii=0;
        static struct _timeb lnow, lthen;
        if(++ii % 8 == 0)
        {
            /* 8 times around measure time */
	          _ftime(&lnow);
		        dbl =  (double)(lnow.time - lthen.time) +
			             (double)(lnow.millitm - lthen.millitm)/1000;
        		lthen = lnow;  rate = 8/dbl;
        }

        return ((float)rate);
    }
#endif

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

void char2wall(float x,float y,float z, char buf[])
{
     char *p;
     glRasterPos3f(x,y,z);
     for(p=buf;*p;p++)glutBitmapCharacter(GLUT_BITMAP_9_BY_15,*p);
}

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

void graffiti(void)
{     
      int myYear = startDate->year;
      int myMonth = startDate->month;
      int myDay = startDate->day;
      char* monthName = NULL;

      switch(myMonth)
      {
          case 1:
              monthName = "January";
              break;
          case 2:
              monthName = "February";
              break;
          case 3:
              monthName = "March";
              break;
          case 4:
              monthName = "April";
              break;
          case 5:
              monthName = "May";
              break;
          case 6:
              monthName = "June";
              break;
          case 7:
              monthName = "July";
              break;
          case 8:
              monthName = "August";
              break;
          case 9:
              monthName = "September";
              break;
          case 10:
              monthName = "October";
              break;
          case 11:
              monthName = "November";
              break;
          case 12:
              monthName = "December";
              break;
      };
      char buf[256]; /* console graffiti are done differently from cave */
      #define LABEL2(x,y,W,u) {sprintf(buf,(W),(u));char2wall(x,y,0.,buf);}
      glMatrixMode(GL_PROJECTION); 
      glPushMatrix(); 
          glLoadIdentity();
          gluOrtho2D(0,3000,0,3000);
          glMatrixMode(GL_MODELVIEW);
          glPushMatrix();
              glLoadIdentity();
              /* writings */
              if(mode==TURNMODE) 
                glColor3f(50./255.,205./255.,50./255.);
              else
                glColor3f(50./255.,205./255.,50./255.);
              
           LABEL2(80, 80, "%4.1f fps", speedometer());
           /*LABEL2(80, 2840,\
           "(ESC)ape {wandbutn} (V)Binoc (MAUS2)Fore Fly(M)ode %s {A}(H)omotopy (W)riting",
                  mode?"FLYING":"CONTROL");*/
           LABEL2(80, 2840, "(ESC) to Exit ***************************", ""); 
           LABEL2(80, 2770, "(H) to Scale the System Non-Uniformly ***", "");
           LABEL2(80, 2700, "(h) to Descale Non-Uniformly ************", "");
           LABEL2(80, 2630, "(K) to Scale the System Uniformly *******", "");
           LABEL2(80, 2560, "(k) to Descale Uniformly ****************", "");
           LABEL2(80, 2490, "(V) to Speed Up Time ********************", "");
           LABEL2(80, 2420, "(v) to Slow Down Time *******************", "");
           LABEL2(80, 2350, "(f) to Freeze the System ****************", "");
           LABEL2(80, 2280, "(z) to Go Back a Year *******************", "");
           LABEL2(80, 2210, "(Z) to Go Forward a Year ****************", "");
           LABEL2(80, 2140, "*****************************************", ""); 
           LABEL2(10,10,"illiSol: A Modern Orrery Application in Syzygy"
               /* April 27, 2008 by William Davis \
                  Anchored in the Skeleton by aszgard14mar08 Francis, Schaeffer, Blank, Woodruff, Bourd, Hartman & Chappell, U Illinois, 1995..2004 %s"*/,"");
           //LABEL2(80,2770,"(N)ose   %0.3f",nose);
           //LABEL2(80,2700,"(S)peed  %0.4f",speed);
           //LABEL2(80,2630," tor(Q) %0.4f",torq);
           //LABEL2(80,2560,"near clipper %g", mysiz*focal);
           //LABEL2(80,2490,"f(O)cal factor %g",focal);
           //LABEL2(80,2420,"my s(I)ze %.2g",mysiz);
           //LABEL2(80,2350,"wfar cli(P)er= %.2g",wfar);
           //LABEL2(80,2280,"(Z)ap{Z} %s","");
           //LABEL2(80,2210,"(G)ap %.2g",gap);
           //LABEL2(80,2140,"(A)mb %.2g",amb);
           //LABEL2(80,2070,"pw(R) %.2g",pwr);
           LABEL2(80, 2040,"Planet Magnification: %3.1f", planetMag);
           LABEL2(80, 1960,"Solar Magnification : %1.1f", solMag);

           LABEL2(80, 810, "Julian Day Number: %10.1f", julianDay);
           LABEL2(80, 720, "Gregorian Date:  ", "");
           LABEL2(160, 630, "Month: %s", monthName); 
           LABEL2(160, 540, "Day  : %d", myDay);
           LABEL2(160, 450, "Year : %d", myYear);       
           //LABEL2(80,1790,"head {B} %d",resethand);
         glPopMatrix();
         glMatrixMode(GL_PROJECTION); 
      glPopMatrix();
} 

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

void chaptrack(int but,int xx,int yy,int shif)
{  
   //was long
   long dx,dy; 
   dx = xx -(int)(.5*xt);
   dx = abs(dx)>5?dx:0;
   dy = yy -(int)(.5*yt);
   dy = abs(dy)>5?dy:0;

   glMatrixMode(GL_MODELVIEW);
   glPushMatrix();
      glLoadIdentity();

      if(mode==TURNMODE) 
        glTranslatef(affmat[12],affmat[13],affmat[14]);
      glRotatef(dx*torq,0.,1.,0.); 
      glRotatef(dy*torq,1.,0.,0.);
      if(but&(1<<GLUT_RIGHT_BUTTON ))
        glRotatef(shif?-10:-1,0.,0.,1.);
      if(but&(1<<GLUT_LEFT_BUTTON  ))
        glRotatef(shif?10:1,0.,0.,1.);
      if(mode==FLYMODE)
      {
        glPushMatrix();
          glMultMatrixf(starmat);
          glGetFloatv(GL_MODELVIEW_MATRIX,starmat);
        glPopMatrix();
      }

      if(but&(1<<GLUT_MIDDLE_BUTTON)) 
        glTranslatef(0.,0.,shif ? -speed : speed);
      if(mode==TURNMODE) 
        glTranslatef(-affmat[12],-affmat[13],-affmat[14]);
   
      glMultMatrixf(affmat); 
      glGetFloatv(GL_MODELVIEW_MATRIX,affmat);
      FOR(ii,0,3)
      {
        luxx[ii]=0; 
        FOR(jj,0,3)
          luxx[ii] += affmat[ii*4 + jj]*lux[jj];
      }
   glPopMatrix();

}

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

void reshaped(int xx, int yy)
{
  /*Origin of the moved window*/
  xt=xx; 
  yt=yy;
} 

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

void drawcons(void)
{  
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  if(binoc) 
      glViewport(0,yt/4,xt/2,yt/2);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  glFrustum(-mysiz*xt/yt, mysiz*xt/yt, -mysiz, mysiz, mysiz*focal, wfar); 
  glMatrixMode(GL_MODELVIEW); 
  glLoadIdentity();
  drawstars();
  glTranslatef(-binoc*nose,0.0,0.0);
  glMultMatrixf(affmat);
  drawall();
  if(binoc){
    glViewport(xt/2, yt/4, xt/2, yt/2); 
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    drawstars();
    glTranslatef(binoc*nose,0.0,0.0);
    glMultMatrixf(affmat);
    drawall();
  }
  glViewport(0,0,xt,yt);
  if(msg)
  {
    glDisable(GL_TEXTURE_2D); 
      graffiti();
    glEnable(GL_TEXTURE_2D);
  }
  glutSwapBuffers();
}

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

void idle(void)
{ 
   /*do this when nothing else is happening*/
   if(morph)
      autotymer(0);      /* advance autotymer */ 
   glutPostRedisplay();  /* redraw the window */
   chaptrack(BUT,XX,YY,SHIF);
}

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

void mousepushed(int but,int stat,int x,int y)
{
  if(stat==GLUT_DOWN) 
    BUT |= (1<<but);
  else 
    BUT &= (-1 ^ (1<<but));  
  XX=x;
  YY=y;
  SHIF=(glutGetModifiers()==GLUT_ACTIVE_SHIFT)?1:0;
}

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

void mousemoved(int x,int y)
{  
  XX=x; 
  YY=y; 
}

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

void dataprep(void)
{
    deFault();
    glDisable(GL_TEXTURE_2D);
    initstars();
    glEnable(GL_TEXTURE_2D);

    /*Initialize the Sun's Characteristics*/
    planets[8][0] = 0.0464913;
    planets[8][1] = planets[8][2] = planets[8][3] = planets[8][4] = planets[8][5] = 0.0;

    initPlanets();

    /*In a desperate attempt to get the texture mapping working I placed
     *the following block of code in many places, we'll eliminate them
     *by trial and error */
    glEnable(GL_COLOR_MATERIAL);
    glColorMaterial(GL_FRONT, GL_AMBIENT_AND_DIFFUSE);
    glEnable(GL_TEXTURE_2D);
    glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
    glEnable(GL_CULL_FACE);
    glEnable(GL_DEPTH_TEST);
    glDepthFunc(GL_LESS);
    glCullFace(GL_BACK);
    sph = gluNewQuadric();
    gluQuadricDrawStyle(sph, GLU_FILL);
    gluQuadricNormals(sph, GLU_SMOOTH);
    gluQuadricTexture(sph, GL_TRUE);
    gluQuadricOrientation(sph, GLU_OUTSIDE);
}
/**********************************************************************/

void overlay(arMasterSlaveFramework&)
{
  if(msg){
    glDisable(GL_TEXTURE_2D);
      graffiti();
    glEnable(GL_TEXTURE_2D);
  }
}

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

//gesturetrack gkf: 2jun04 
/***************************************************************/
/* Given a rotation matrix convert it into a quaternion    */
/* Given a quaternion convert it into a rotation matrix    */
/* By Robert Shuttleworth, illiMath2001, July 2001         */
/*                                                         */
/***************************************************************/

//gkf float mat[16]={1,0,0,0,   0,1,0,0,   0,0,1,0,   0,0,0,1};
//gkf float quat[4] ={1,0,0,0}; 

void normalize(float *vec)
{ 
  int ii;
  float temp = vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2] + vec[3]*vec[3];
  if (temp == 1)
  {
    /*printf(" The vector is of unit length.\n");*/
  }
  else
  { 
    temp = sqrt(temp);
    for(ii=0;ii<4;ii++)
        vec[ii] /= temp;
  }
}

void printquat(float *quat)
{
  printf("\n Quaternion: \n");
	printf(" w = %f ", quat[0]);
	printf(" x = %f ", quat[1]);
	printf(" y = %f ", quat[2]);
	printf(" z = %f ", quat[3]);
}

void printmat(float *mat)
{ 
  printf("\n Matrix: \n");
	printf(" %f %f %f %f \n", mat[0], mat[4], mat[8],  mat[12]);
	printf(" %f %f %f %f \n", mat[1], mat[5], mat[9],  mat[13]);
	printf(" %f %f %f %f \n", mat[2], mat[6], mat[10], mat[14]);
	printf(" %f %f %f %f \n", mat[3], mat[7], mat[11], mat[15]);
}

void q2m(float *quat, float *mat)
{ 
  /* quaternion yields a rotation */
  float w = quat[0], x = quat[1], y = quat[2], z = quat[3]; 
  float ww, xx, yy, zz, xy, xz, yz, wz, wy, wx;
  ww = w*w; xx = x*x; yy = y*y; zz = z*z; 
  xy = x*y; xz = x*z; yz = y*z; wz = w*z; wy = w*y; wx = w*x; xz = x*z;
  mat[0]=ww+xx-yy-zz;  mat[4]=2*xy - 2*wz;  mat[8]=2*wy + 2*xz;
  mat[1]=2*wz + 2*xy;  mat[5]=ww-xx+yy-zz;  mat[9]=2*yz - 2*wx;
  mat[2]=2*xz - 2*wy;  mat[6]=2*wx + 2*yz; mat[10]=ww-xx-yy+zz;
  mat[3] = mat[7] = mat[11] = mat[12] = mat[13] = mat[14] = 0;
  mat[15] = ww + xx + yy + zz; /* this is a check */
}

void m2q(float *mat, float *quat)
{  
    /* isometry yields a quaternion */ 
    float tmp; int ii, jj;

    /*pseudo traces*/ 
    float tr[4] ={mat[0] + mat[5] + mat[10],
                  mat[0] - mat[5] - mat[10],
                 -mat[0] + mat[5] - mat[10],
                 -mat[0] - mat[5] + mat[10]};

    //printf("trace  %f %f %f %f \n" , tr[0],tr[1],tr[2],tr[3]);  
    ii=1; jj=2;  /* find the largest tr */
    if(tr[0] > tr[1])
      ii = 0; 
    if(tr[3] > tr[2])
      jj = 3; 
    if(tr[jj] > tr[ii])
      ii = jj;
    //printf("\n %i is the largest trace. \n", ii);
    tmp = 4*sqrt(tr[ii]+1)/2; 
    switch(ii)
    {
        case 0: 
          quat[0] = tmp/4;
          quat[1] = (mat[6] - mat[9]) / tmp; 
          quat[2] = (mat[8] - mat[2]) / tmp; 
          quat[3] = (mat[1] - mat[4]) / tmp; 
          break; 
        case 1: 
          quat[0] = (mat[6] - mat[9]) / tmp;
          quat[1] = tmp/4; 
          quat[2] = (mat[1] + mat[4]) / tmp;
          quat[3] = (mat[8] + mat[2]) / tmp;
          break;
        case 2: 
          quat[0] = (mat[8] - mat[2]) / tmp;
          quat[1] = (mat[4] + mat[1]) / tmp; 
          quat[2] = tmp/4; 
          quat[3] = (mat[9] + mat[6]) / tmp; 
          break; 
        case 3:
          quat[0] = (mat[1] - mat[4]) / tmp;
          quat[1] = (mat[8] + mat[2]) / tmp; 
          quat[2] = (mat[9] + mat[6]) / tmp; 
          quat[3] = tmp/4;
          break; 
    }
}

void loadIdentity(float* mat)    
{
    int ii; 
    for(ii=0;ii<16;ii++) 
        mat[ii]=float(ii/4 == ii%4);
}

void getmat(float *target, float *source)
{ 
     memcpy(target,source, 16*sizeof(float));
}

void multmatinv(float *T)
{  
   /* of an Euclidean isometry */ 
   float N[16];  getmat(N,T);
   N[1] = T[4]; N[2] = T[8]; N[6] = T[9];
   N[4] = T[1]; N[8] = T[2]; N[9] = T[6];
   N[12]= -T[0]*T[12]-T[1]*T[13]-T[2]*T[14];
   N[13]= -T[4]*T[12]-T[5]*T[13]-T[6]*T[14];
   N[14]= -T[8]*T[12]-T[9]*T[13]-T[10]*T[14];
   glMultMatrixf(N);
}

/* uses eps to shrink rotation and del to shring translation */

void updateaffmat(float *oldhand, float *newhand)
{ 
  float dx=0, dy=0, dz=0 ;
  float dhand[16];
  float quat[4];
  glMatrixMode(GL_TEXTURE);  //do some arithmetic on an unused stack 
  glPushMatrix();
      glLoadIdentity();
      multmatinv(newhand);
      glMultMatrixf(oldhand); 
      glGetFloatv(GL_TEXTURE_MATRIX, dhand);
  glPopMatrix(); //dhand= newhand^-1 * oldhand 
  dx=del*dhand[12]; dhand[12]=0;  //attenuation is about right
  dy=del*dhand[13]; dhand[13]=0; 
  dz=del*dhand[14]; dhand[14]=0; 
  //use a bit of the translation in dhand and then make it the rotation
  if(gesture)
  { //shrink dhand
    m2q(dhand,quat);
    quat[0]= 1-eps + eps*quat[0];
    quat[1] *=  eps; quat[2] *=  eps; quat[3] *=  eps;
    normalize(quat); 
    q2m(quat,dhand); 
  }


 if(mode==FLYMODE)
 {
   glMatrixMode(GL_TEXTURE); glLoadIdentity();
   glMultMatrixf(dhand); //rotate the stars
   glMultMatrixf(starmat);
   glGetFloatv(GL_TEXTURE_MATRIX,starmat);
 }

  glMatrixMode(GL_TEXTURE); 
  glPushMatrix();
      glLoadIdentity();
      if(mode==TURNMODE) 
        glTranslatef(affmat[12],affmat[13],affmat[14]);
      glMultMatrixf(dhand);  //just the rotation of it
      glTranslatef(dx,dy,dz);
      if(mode==TURNMODE)
        glTranslatef(-affmat[12],-affmat[13],-affmat[14]);
      glMultMatrixf(affmat);
      glGetFloatv(GL_TEXTURE_MATRIX,affmat);
}


void readbuttons(arMasterSlaveFramework& fw)
{
  const int butn[6]= { fw.getButton(0), fw.getButton(1), fw.getButton(2), 
                      fw.getButton(3), fw.getButton(4), fw.getButton(5) };
  static int obutn[6]={0,0,0,0,0,0};

  if(obutn[3]==0 & butn[3] == 1)
  {
    //Should increase the scale uniformly
    keyboard('H', 0,0); //{X} button 
  }
  if(obutn[2]==0 & butn[2] == 1)
  {
    //Should decrease the speed of time
    keyboard('v', 0,0);
    //gesture=1-gesture; //{C}button
  }
  if(obutn[1]==0 & butn[1] == 1)
  {
    //Should increase the speed of time
    keyboard('V', 0,0);
    //resethand=1;       //{B}button
  }
  if(obutn[4]==0 & butn[4] == 1)
  {
    //Should decrease the scale non-uniformly
    keyboard('h', 0,0);  
    mode=1-mode;       // {Y}button
  } 
  if(obutn[5]==0 & butn[5] == 1)
  {
    //Should go back one year in time
    keyboard('z', 0,0); //{Z} button
    //pause = 1 - pause; 
  }
  if(obutn[0]==0 & butn[0] == 1)
  {
    //Should pause time
    keyboard('f', 0,0); //{A} button?  
  }

  {
    static int click=GOVERNOR;
    if(morph&&click==0)
        autotymer(0);
    click=click>0?click-1:GOVERNOR;
  }
  
  memcpy(obutn,butn,6*sizeof(int));
}

void gesturetrack(arMasterSlaveFramework& fw)
{
    //the present logic is legacy stuff from zyspace, rationalize it!
    if (!fw.getMaster())
        return; //only the master should listen
    //this is the gesture navigator
    static float  oldhand[16];
    float handmat[16];

    arMatrix4 arhand=fw.getMatrix(1);
    getmat(handmat,arhand.v);
    //geotrack(arhand.v) gets arhand.v and used as float* handmat
    static bool  reshand=true;
    if(resethand)
    {
      reshand=true;
      resethand=false;
    }
    
    //do we really need two flags?
    if(gesture) 
      updateaffmat(oldhand,handmat);

    if(reshand)
    {
      reshand=false;
      getmat(oldhand, handmat);
    }



    //this is the Pape navigator
    //this should really be done on the arithmetic (Texture) stack 
    glMatrixMode(GL_MODELVIEW); 
    glPushMatrix();
      glLoadIdentity();
      const float joy0 = fw.getAxis(0);  //joystick right-left 
      const float joy1 = fw.getAxis(1); //joystick fore-aft
      const float xx = arhand.v[8];
      const float yy = arhand.v[9];
      const float zz = arhand.v[10];
      // printf(" xx %g yy %g zz %g \n",xx,yy,zz);

      const float da = 1.2*joy0; /*.5*joy0;*/ 
      const float ds = .65*joy1; /*.1*joy1;*/
      //joystick translates in the wand direction
      //if(ds < 0)
      //  glTranslatef(-1*xx*ds, -1*yy*ds, -1*zz*ds);
      //else
          glTranslatef(xx*ds, yy*ds, zz*ds);
      //if(ds < 0)
      //  glRotatef(da, -1*xx*ds,-1*yy*ds,-1*zz*ds);
      //else
          glRotatef(da, 1.,1.,0.); //was da,xx*ds,yy*ds,xx*ds
      glMultMatrixf(affmat);
      glGetFloatv(GL_MODELVIEW_MATRIX,affmat);
    glPopMatrix();

    //this is an illogical place for this but maybe it works again

    readbuttons(fw);

}

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

void diddle(arMasterSlaveFramework&)
{
  //now we want to do any calculations that would rely on the shared stuff
}

void cubedraw(arMasterSlaveFramework& fw)
{
  glPushMatrix();
    glTranslatef(0,5,0);
    glScalef(sizz,sizz,sizz);
    drawstars();
  glPopMatrix();

  glMultMatrixf(affmat);
  glScalef(sizz,sizz,sizz);
  drawall();
}

int getInitDate(Timestamp* starting)
{
    int first, i;
    first = 1;
    while(era == 0)
    {
        #ifdef DEBUG
          cout << "In date while loop..." << endl;
        #endif
        /*Check to see what "era" we want to compute dates for*/
        era = getDate(starting);

        if(era == 0)
        {
            printf("Error: The date you have chosen is outside");
            printf(" the domain of our modeling equations.\n");
            return 0;
        }else if(first == 1)
        {    
          /*Here we convert the Gregorian Calendar date to
           *the Julian Day Number */
          julianDay = convertGregToJDN(starting);
          first = 0;
        }
    }

   for(i = 0; i < 8; i++)
   {
       #ifdef DEBUG
         cout << "Computing position for planets..." << endl;
       #endif
       compute_Position(julianDay, i + 1, era, 0);
   }
}

int main(int argc, char **argv)
{ 
   arguments(argc, argv);
   dataprep();
   startDate = (Timestamp *)malloc(sizeof(Timestamp));
   if(startDate == NULL)
   {
       return 0;
   }
   getInitDate(startDate);

   if(eyes==console) /*Console Main*/
   { 
       glutInit(&argc, argv);
       glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB |GLUT_DEPTH);
       
       /*This is another set of that block of code described in
        *dataprep()... */
       glEnable(GL_COLOR_MATERIAL);
       glColorMaterial(GL_FRONT, GL_AMBIENT_AND_DIFFUSE);
       glEnable(GL_TEXTURE_2D);
       glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
       glEnable(GL_CULL_FACE);
       glEnable(GL_DEPTH_TEST);
       glDepthFunc(GL_LESS);
       glCullFace(GL_BACK); 
       switch(win)
       { 
         case 0:
           break;
         case 1: 
           glutInitWindowSize(640, 480);
           glutInitWindowPosition(0,1024-480);
           break;
         case 2: 
           glutInitWindowPosition(0,0);
	         break;
       }

       glutCreateWindow("illiSol: A Modern Orrery for Virtual Environments");
       if(glutGet(GLUT_WINDOW_COLORMAP_SIZE) != 0)
       {
         exit(1);
       }
       if(win==2) 
         glutFullScreen();
       /* Here's that same block of code again... silly */
       glEnable(GL_COLOR_MATERIAL);
       glColorMaterial(GL_FRONT, GL_AMBIENT_AND_DIFFUSE);
       glEnable(GL_TEXTURE_2D);
       glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
       glEnable(GL_CULL_FACE);
       glEnable(GL_DEPTH_TEST); 
       glutDisplayFunc(drawcons);       
       glutKeyboardFunc(keyboard);
       glutSpecialFunc(special_keybo);
       glutMouseFunc(mousepushed);
       glutMotionFunc(mousemoved);       
       glutPassiveMotionFunc(mousemoved); 
       glutReshapeFunc(reshaped);
       glutIdleFunc(idle);             
       glutMainLoop();
   }
   else if(eyes==cube) /*Syzygy Main*/
   { 
       arMasterSlaveFramework* fff = new arMasterSlaveFramework();
       if(!fff->init(argc,argv)) 
         return 1;
       fff->setStartCallback(cubeinit); 
       fff->setPreExchangeCallback(gesturetrack);
       fff->setDrawCallback(cubedraw);
       fff->setPostExchangeCallback(diddle);
       fff->setKeyboardCallback(clefs);
       fff->setOverlayCallback(overlay);

       /*Same block of code for the Syzygy Main */
       glEnable(GL_COLOR_MATERIAL);
       glColorMaterial(GL_FRONT, GL_AMBIENT_AND_DIFFUSE);
       glEnable(GL_TEXTURE_2D);
       glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
       glEnable(GL_CULL_FACE);
       glEnable(GL_DEPTH_TEST);
       
       return fff->start() ? 0:1;
   }

   /*Once we're done with everything, delete the template we
    *made for the planets*/
   gluDeleteQuadric(sph);
}

void compute_Position(double julianDate, int planet, int era, int orbit)
{
    #ifdef DEBUG
      cout << "In compute_Position()..." << endl;
    #endif

    if(era == 1)
        compute_Position_AD(julianDate, planet, orbit);
    else if(era == 2)
        compute_Position_BCandAD(julianDate, planet, orbit);
}

void compute_Position_AD(double julianDate, int planet, int orbit)
{
    #ifdef DEBUG
      cout << "In compute_Position_AD()" << endl;
    #endif

    /*Yeah, a lot of temporary variables to hold onto
     *things while we compute the planetary positions */
    double sm_axis, sm_error, eccentricity, inclination, 
           meanLong, longPeri, longAscNode, relativeDate,
           argPeri, meanAnom, eccDegrees, oldEccAnom,
           dEccAnom, newEccAnom, dMeanAnom, absDAnom;
    double ecc_error, inc_error, mLong_error, lPeri_error,
           lAscNode_error;
    double helioCX, helioCY, helioCZ;

    sm_axis = sm_error = eccentricity = inclination = meanLong = 0.0;
    longPeri = longAscNode = relativeDate = dEccAnom = newEccAnom = 0.0;
    dMeanAnom = absDAnom = ecc_error = inc_error = mLong_error = 0.0;
    lPeri_error = lAscNode_error = helioCX = helioCY = helioCZ = 0.0;
    
    /*This switch statement simply retrieves the Keplerian
     *elements that we are going to be computing with. Nothing
     *interesting here.*/

    switch(planet)
    {
        case 1: 
          sm_axis = MERCURY_SM_AXIS_A;
          sm_error = MERCURY_SM_ERROR_A;
          eccentricity = MERCURY_ECC_A;
          ecc_error = MERCURY_ECC_ERROR_A;
          inclination = MERCURY_INC_A;
          inc_error = MERCURY_INC_ERROR_A;
          meanLong = MERCURY_MLONG_A;
          mLong_error = MERCURY_MLONG_ERROR_A;
          longPeri = MERCURY_PERI_A;
          lPeri_error = MERCURY_PERI_ERROR_A;
          longAscNode = MERCURY_ASC_A;
          lAscNode_error = MERCURY_ASC_ERROR_A;
          break;
        case 2:
          sm_axis = VENUS_SM_AXIS_A;
          sm_error = VENUS_SM_ERROR_A;
          eccentricity = VENUS_ECC_A;
          ecc_error = VENUS_ECC_ERROR_A;
          inclination = VENUS_INC_A;
          inc_error = VENUS_INC_ERROR_A; 
          meanLong = VENUS_MLONG_A;
          mLong_error = VENUS_MLONG_ERROR_A; 
          longPeri = VENUS_PERI_A;
          lPeri_error = VENUS_PERI_ERROR_A; 
          longAscNode = VENUS_ASC_A;
          lAscNode_error = VENUS_ASC_ERROR_A; 
          break;
        case 3:
          sm_axis = EARTH_SM_AXIS_A;
          sm_error = EARTH_SM_ERROR_A;
          eccentricity = EARTH_ECC_A;
          ecc_error = EARTH_ECC_ERROR_A;
          inclination = EARTH_INC_A;
          inc_error = EARTH_INC_ERROR_A;
          meanLong = EARTH_MLONG_A;
          mLong_error = EARTH_MLONG_ERROR_A; 
          longPeri = EARTH_PERI_A;
          lPeri_error = EARTH_PERI_ERROR_A; 
          longAscNode = EARTH_ASC_A;
          lAscNode_error = EARTH_ASC_ERROR_A;
          break;
        case 4: 
          sm_axis = MARS_SM_AXIS_A;
          sm_error = MARS_SM_ERROR_A;
          eccentricity = MARS_ECC_A;
          ecc_error = MARS_ECC_ERROR_A;
          inclination = MARS_INC_A;
          inc_error = MARS_INC_ERROR_A;
          meanLong = MARS_MLONG_A;
          mLong_error = MARS_MLONG_ERROR_A; 
          longPeri = MARS_PERI_A;
          lPeri_error = MARS_PERI_ERROR_A; 
          longAscNode = MARS_ASC_A;
          lAscNode_error = MARS_ASC_ERROR_A;
          break;
        case 5: 
          sm_axis = JUPITER_SM_AXIS_A;
          sm_error = JUPITER_SM_ERROR_A;
          eccentricity = JUPITER_ECC_A;
          ecc_error = JUPITER_ECC_ERROR_A;
          inclination = JUPITER_INC_A;
          inc_error = JUPITER_INC_ERROR_A;
          meanLong = JUPITER_MLONG_A;
          mLong_error = JUPITER_MLONG_ERROR_A; 
          longPeri = JUPITER_PERI_A;
          lPeri_error = JUPITER_PERI_ERROR_A; 
          longAscNode = JUPITER_ASC_A;
          lAscNode_error = JUPITER_ASC_ERROR_A;
          break;
        case 6: 
          sm_axis = SATURN_SM_AXIS_A;
          sm_error = SATURN_SM_ERROR_A;
          eccentricity = SATURN_ECC_A;
          ecc_error = SATURN_ECC_ERROR_A;
          inclination = SATURN_INC_A;
          inc_error = SATURN_INC_ERROR_A;
          mLong_error = SATURN_MLONG_ERROR_A; 
          meanLong = SATURN_MLONG_A;
          longPeri = SATURN_PERI_A;
          lPeri_error = SATURN_PERI_ERROR_A; 
          longAscNode = SATURN_ASC_A;
          lAscNode_error = SATURN_ASC_ERROR_A;
          break;
        case 7: 
          sm_axis = URANUS_SM_AXIS_A;
          sm_error = URANUS_SM_ERROR_A;
          eccentricity = URANUS_ECC_A;
          ecc_error = URANUS_ECC_ERROR_A;
          inclination = URANUS_INC_A;
          inc_error = URANUS_INC_ERROR_A;
          meanLong = URANUS_MLONG_A;
          mLong_error = URANUS_MLONG_ERROR_A; 
          longPeri = URANUS_PERI_A;
          lPeri_error = URANUS_PERI_ERROR_A; 
          longAscNode = URANUS_ASC_A;
          lAscNode_error = URANUS_ASC_ERROR_A; 
          break;
        case 8: 
          sm_axis = NEPTUNE_SM_AXIS_A;
          sm_error = NEPTUNE_SM_ERROR_A;
          eccentricity = NEPTUNE_ECC_A;
          ecc_error = NEPTUNE_ECC_ERROR_A;
          inclination = NEPTUNE_INC_A;
          inc_error = NEPTUNE_INC_ERROR_A;
          meanLong = NEPTUNE_MLONG_A;
          mLong_error = NEPTUNE_MLONG_ERROR_A; 
          longPeri = NEPTUNE_PERI_A;
          lPeri_error = NEPTUNE_PERI_ERROR_A; 
          longAscNode = NEPTUNE_ASC_A;
          lAscNode_error = NEPTUNE_ASC_ERROR_A; 
          break;
    }
    

    /* Okay, now we have our elements, let's start the computing. */
    /* First we need to find the elements at this point in time,
     * so we offset according to the date */
    relativeDate = (julianDate - 2451545.0)/36525.0;
    sm_axis = sm_axis + (sm_error * relativeDate);
    eccentricity = eccentricity + (ecc_error * relativeDate);
    inclination = inclination + (inc_error * relativeDate);
    meanLong = meanLong + (mLong_error * relativeDate);
    longPeri = longPeri + (lPeri_error * relativeDate);
    longAscNode = longAscNode + (lAscNode_error * relativeDate);

    argPeri = longPeri - longAscNode;
    meanAnom = meanLong - longPeri;
    
    #ifdef DEBUG
      cout << "Current Mean Anomaly: ";
      cout << meanAnom << endl;
    #endif

    double half = 360.0; //Was 180, but only half of the orbit drew
    int temp;
    temp = (int)(meanAnom / half);
    meanAnom = meanAnom - ((float)(temp))*half;
    
    /*This part was added to make it go from -180 <= M <= 180 */
    if(meanAnom > 180)
      meanAnom -= 360;

    eccDegrees = eccentricity * (180 / PI);
    newEccAnom = meanAnom + eccDegrees * (S(meanAnom));
    do{
        #ifdef DEBUG
          cout << "Attempting to get more accurate" << endl;
        #endif
        oldEccAnom = newEccAnom;
        dMeanAnom = meanAnom - 
              (oldEccAnom - eccDegrees * (S(oldEccAnom)));
        dEccAnom = dMeanAnom/(1 - eccentricity * C(oldEccAnom));
        //dEccAnom = dMeanAnom/(1-eccDegrees*C(oldEccAnom));
        newEccAnom = oldEccAnom + dEccAnom;
        if(dEccAnom < 0)
            absDAnom = -1*dEccAnom;
        else
            absDAnom = dEccAnom;
    }
    while(absDAnom > (pow(10.0, -6.0)));
    
    #ifdef DEBUG
      cout << "Returned from Mean Anomaly Calculation..." << endl;
    #endif

    /* We now have our Eccentric Anomaly, and we can get the coordinates*/
   
   if(orbit == 0)
   { 
      planets[planet-1][1] = sm_axis * (cos((PI/180)*newEccAnom) - eccentricity);
      planets[planet-1][2] = sm_axis * (sqrt(1 - (eccentricity)*(eccentricity)))*
                             S(newEccAnom);
   }
      
   helioCX = (C(argPeri)*C(longAscNode) -
              S(argPeri)*S(longAscNode)*C(inclination))*(planets[planet-1][1]);
   helioCX = helioCX - (S(argPeri)*C(longAscNode) +
             C(argPeri)*S(longAscNode)*C(inclination))*(planets[planet-1][2]);

   helioCY = (C(argPeri)*S(longAscNode) +
              S(argPeri)*C(longAscNode)*C(inclination))*(planets[planet-1][1]);
   helioCY = helioCY - (S(argPeri)*S(longAscNode) - 
             C(argPeri)*C(longAscNode)*C(inclination))*(planets[planet-1][2]);

   helioCZ = (S(argPeri)*S(inclination))*(planets[planet-1][1]) +
             (C(argPeri)*S(inclination))*(planets[planet-1][2]);

   if(orbit == 0)
   {
      planets[planet-1][3] = helioCX;
      planets[planet-1][4] = helioCY;
      planets[planet-1][5] = helioCZ;
   }else
   {
      planetsOrbit[planet-1][0] = helioCX;
      planetsOrbit[planet-1][1] = helioCY;
      planetsOrbit[planet-1][2] = helioCZ;
   } 
}

void compute_Position_BCandAD(double julianDate, int planet, int orbit)
{
    
    double sm_axis, sm_error, eccentricity, inclination, 
           meanLong, longPeri, longAscNode, relativeDate,
           argPeri, meanAnom, eccDegrees, oldEccAnom,
           dEccAnom, newEccAnom, dMeanAnom, b, c, s, f,
           temp, absDAnom;
    double ecc_error, inc_error, mLong_error, lPeri_error,
           lAscNode_error;
    double helioCX, helioCY, helioCZ;

    sm_axis = sm_error = eccentricity = inclination = meanLong = 0.0;
    longPeri = longAscNode = relativeDate = dEccAnom = newEccAnom = 0.0;
    dMeanAnom = absDAnom = ecc_error = inc_error = mLong_error = 0.0;
    lPeri_error = lAscNode_error = helioCX = helioCY = helioCZ = 0.0;
    
    /*This switch statement simply retrieves the Keplerian
     *elements that we are going to be computing with. Nothing
     *interesting here.*/

    switch(planet)
    {
        case 1:
          sm_axis = MERCURY_SM_AXIS_B;
          sm_error = MERCURY_SM_ERROR_B;
          eccentricity = MERCURY_ECC_B;
          ecc_error = MERCURY_ECC_ERROR_B;
          inclination = MERCURY_INC_B;
          inc_error = MERCURY_INC_ERROR_B;
          meanLong = MERCURY_MLONG_B;
          mLong_error = MERCURY_MLONG_ERROR_B;
          longPeri = MERCURY_PERI_B;
          lPeri_error = MERCURY_PERI_ERROR_B;
          longAscNode = MERCURY_ASC_B;
          lAscNode_error = MERCURY_ASC_ERROR_B;
          break;
        case 2:
          sm_axis = VENUS_SM_AXIS_B;
          sm_error = VENUS_SM_ERROR_B;
          eccentricity = VENUS_ECC_B;
          ecc_error = VENUS_ECC_ERROR_B;
          inclination = VENUS_INC_B;
          inc_error = VENUS_INC_ERROR_B; 
          meanLong = VENUS_MLONG_B;
          mLong_error = VENUS_MLONG_ERROR_B; 
          longPeri = VENUS_PERI_B;
          lPeri_error = VENUS_PERI_ERROR_B; 
          longAscNode = VENUS_ASC_B;
          lAscNode_error = VENUS_ASC_ERROR_B; 
          break;
        case 3:
          sm_axis = EARTH_SM_AXIS_B;
          sm_error = EARTH_SM_ERROR_B;
          eccentricity = EARTH_ECC_B;
          ecc_error = EARTH_ECC_ERROR_B;
          inclination = EARTH_INC_B;
          inc_error = EARTH_INC_ERROR_B;
          meanLong = EARTH_MLONG_B;
          mLong_error = EARTH_MLONG_ERROR_B; 
          longPeri = EARTH_PERI_B;
          lPeri_error = EARTH_PERI_ERROR_B; 
          longAscNode = EARTH_ASC_B;
          lAscNode_error = EARTH_ASC_ERROR_B;
          break;
        case 4: 
          sm_axis = MARS_SM_AXIS_B;
          sm_error = MARS_SM_ERROR_B;
          eccentricity = MARS_ECC_B;
          ecc_error = MARS_ECC_ERROR_B;
          inclination = MARS_INC_B;
          inc_error = MARS_INC_ERROR_B;
          meanLong = MARS_MLONG_B;
          mLong_error = MARS_MLONG_ERROR_B; 
          longPeri = MARS_PERI_B;
          lPeri_error = MARS_PERI_ERROR_B; 
          longAscNode = MARS_ASC_B;
          lAscNode_error = MARS_ASC_ERROR_B;
          break;
        case 5: 
          sm_axis = JUPITER_SM_AXIS_B;
          sm_error = JUPITER_SM_ERROR_B;
          eccentricity = JUPITER_ECC_B;
          ecc_error = JUPITER_ECC_ERROR_B;
          inclination = JUPITER_INC_B;
          inc_error = JUPITER_INC_ERROR_B;
          meanLong = JUPITER_MLONG_B;
          mLong_error = JUPITER_MLONG_ERROR_B; 
          longPeri = JUPITER_PERI_B;
          lPeri_error = JUPITER_PERI_ERROR_B; 
          longAscNode = JUPITER_ASC_B;
          lAscNode_error = JUPITER_ASC_ERROR_B;
          b = JUPITER_b;
          c = JUPITER_c;
          s = JUPITER_s;
          f = JUPITER_f;
          break;
        case 6: 
          sm_axis = SATURN_SM_AXIS_B;
          sm_error = SATURN_SM_ERROR_B;
          eccentricity = SATURN_ECC_B;
          ecc_error = SATURN_ECC_ERROR_B;
          inclination = SATURN_INC_B;
          inc_error = SATURN_INC_ERROR_B;
          mLong_error = SATURN_MLONG_ERROR_B; 
          meanLong = SATURN_MLONG_B;
          longPeri = SATURN_PERI_B;
          lPeri_error = SATURN_PERI_ERROR_B; 
          longAscNode = SATURN_ASC_B;
          lAscNode_error = SATURN_ASC_ERROR_B;
          b = SATURN_b;
          c = SATURN_c;
          s = SATURN_s;
          f = SATURN_f;
          break;
        case 7: 
          sm_axis = URANUS_SM_AXIS_B;
          sm_error = URANUS_SM_ERROR_B;
          eccentricity = URANUS_ECC_B;
          ecc_error = URANUS_ECC_ERROR_B;
          inclination = URANUS_INC_B;
          inc_error = URANUS_INC_ERROR_B;
          meanLong = URANUS_MLONG_B;
          mLong_error = URANUS_MLONG_ERROR_B; 
          longPeri = URANUS_PERI_B;
          lPeri_error = URANUS_PERI_ERROR_B; 
          longAscNode = URANUS_ASC_B;
          lAscNode_error = URANUS_ASC_ERROR_B; 
          b = URANUS_b;
          c = URANUS_c;
          s = URANUS_s;
          f = URANUS_f;
          break;
        case 8: 
          sm_axis = NEPTUNE_SM_AXIS_B;
          sm_error = NEPTUNE_SM_ERROR_B;
          eccentricity = NEPTUNE_ECC_B;
          ecc_error = NEPTUNE_ECC_ERROR_B;
          inclination = NEPTUNE_INC_B;
          inc_error = NEPTUNE_INC_ERROR_B;
          meanLong = NEPTUNE_MLONG_B;
          mLong_error = NEPTUNE_MLONG_ERROR_B; 
          longPeri = NEPTUNE_PERI_B;
          lPeri_error = NEPTUNE_PERI_ERROR_B; 
          longAscNode = NEPTUNE_ASC_B;
          lAscNode_error = NEPTUNE_ASC_ERROR_B; 
          b = NEPTUNE_b;
          c = NEPTUNE_c;
          s = NEPTUNE_s;
          f = NEPTUNE_f;
          break;
    }

    /* Okay, now we have our elements, let's start the computing. */
    /* First we need to find the semimajor axis at this point in time */
    relativeDate = (julianDate - 2451545.0)/36525.0;
    sm_axis = sm_axis + (sm_error * relativeDate);
    eccentricity = eccentricity + (ecc_error * relativeDate);
    inclination = inclination + (inc_error * relativeDate);
    meanLong = meanLong + (mLong_error * relativeDate);
    longPeri = longPeri + (lPeri_error * relativeDate);
    longAscNode = longAscNode + (lAscNode_error * relativeDate);
    
    argPeri = longPeri - longAscNode;
    if(planet >= 5)
    {
        temp = b * (relativeDate)*(relativeDate);
        temp = temp + c * cos(relativeDate * f);
        temp = temp + s * sin(relativeDate * f);
        meanAnom = meanLong - longPeri + temp;
    }
    else
        meanAnom = meanLong - longPeri;

    /*Everything that follows is the exact same as the function above
     *but since it's only used twice, copy and past was easier
     *than writing a function and passing all the orbital parameters
     *to it */

    #ifdef DEBUG
      cout << "Current Mean Anomaly: ";
      cout << meanAnom << endl;
    #endif

    double half = 360.0; //Was 180, but only half of the orbit drew
    int temp1;
    temp1 = (int)(meanAnom / half);
    meanAnom = meanAnom - ((float)(temp1))*half;
    if(meanAnom > 180)
      meanAnom -= 360;

    eccDegrees = eccentricity * (180 / PI);
    newEccAnom = meanAnom + eccDegrees * (S(meanAnom));
    do{
        #ifdef DEBUG
          cout << "Attempting to get more accurate" << endl;
        #endif
        oldEccAnom = newEccAnom;
        dMeanAnom = meanAnom - 
              (oldEccAnom - eccDegrees * (S(oldEccAnom)));
        dEccAnom = dMeanAnom/(1 - eccentricity * C(oldEccAnom));
        newEccAnom = oldEccAnom + dEccAnom;
        if(dEccAnom < 0)
            absDAnom = -1*dEccAnom;
        else
            absDAnom = dEccAnom;
    }
    while(absDAnom > (pow(10.0, -6.0)));
    
    #ifdef DEBUG
      cout << "Returned from Mean Anomaly Calculation..." << endl;
    #endif

    /* We now have our Eccentric Anomaly, and we can get the coordinates*/
   
   if(orbit == 0)
   { 
      planets[planet-1][1] = sm_axis * (cos((PI/180)*newEccAnom) - eccentricity);
      planets[planet-1][2] = sm_axis * (sqrt(1 - (eccentricity)*(eccentricity)))*
                             S(newEccAnom);
   }
      
   helioCX = (C(argPeri)*C(longAscNode) -
              S(argPeri)*S(longAscNode)*C(inclination))*(planets[planet-1][1]);
   helioCX = helioCX - (S(argPeri)*C(longAscNode) +
             C(argPeri)*S(longAscNode)*C(inclination))*(planets[planet-1][2]);

   helioCY = (C(argPeri)*S(longAscNode) +
              S(argPeri)*C(longAscNode)*C(inclination))*(planets[planet-1][1]);
   helioCY = helioCY - (S(argPeri)*S(longAscNode) - 
             C(argPeri)*C(longAscNode)*C(inclination))*(planets[planet-1][2]);

   helioCZ = (S(argPeri)*S(inclination))*(planets[planet-1][1]) +
             (C(argPeri)*S(inclination))*(planets[planet-1][2]);

   if(orbit == 0)
   {
      planets[planet-1][3] = helioCX;
      planets[planet-1][4] = helioCY;
      planets[planet-1][5] = helioCZ;
   }else
   {
      planetsOrbit[planet-1][0] = helioCX;
      planetsOrbit[planet-1][1] = helioCY;
      planetsOrbit[planet-1][2] = helioCZ;
   } 
}

void updateGregorian()
{
    int j, g, dg, c, dc, b, db, a, da, y, m, d, Y, M, D;
    j = julianDay + 32044;
    g = (int)(j / 146097);
    dg = j % 146097;
    c = ((dg / 36524 + 1) * 3) / 4;
    dc = dg - c * 36524;
    b = dc / 1461;
    db = dc % 1461;
    a = ((db / 365 + 1) * 3) / 4;
    da = db - a * 365;
    y = g*400 + c*100 + b*4 + a;
    m = (da*5 + 308) / 153 - 2;
    d = da - (m+4)*153 / 5 + 122;
    Y = y - 4800 + (m + 2) / 12;
    M = ((m + 2) % 12) + 1;
    D = d + 1;
    startDate->year = Y;
    startDate->month = M;
    startDate->day = D;
}

double convertGregToJDN(Timestamp* startDate)
{   
    double julDate;
    int a, y, m, myyear, mymonth, myday, myhour, myminute, mysecond;
    
    myyear = startDate->year;
    mymonth = startDate->month;
    myday = startDate->day;
    myhour = startDate->hour;
    myminute = startDate->minute;
    mysecond = startDate->second;
    
    a = (14 - mymonth) / 12; //Note the truncation to type int
    y = myyear + 4800 - a;
    m = mymonth + 12 * a - 3;

    //Again, truncation will be used in the next line, since
    //we are only dealing with integer days

    julDate = myday + (153 * m + 2) / 5 + 365 * y + y / 4 - y / 100;
    julDate = julDate + y / 400 - 32045;

    //Now, since we want the fraction of the day, we cast
    //our integers to type double
    julDate = julDate + ((double)myhour - 12.0) / 24.0;
    julDate = julDate + ((double)myminute) / 1440.0;
    julDate = julDate + ((double)mysecond) / 86400.0;
    return julDate;
}

int getDate(Timestamp* startDate)
{
    int tempY, tempM, tempD, tempH, tempMin, tempS;
    
    tempY = 1977;
    tempM = 8;
    tempD = 20;
    tempH = 12;
    tempMin = 0;
    tempS = 0;
    
    startDate->year = tempY;
    startDate->month = tempM;
    startDate->day = tempD;
    startDate->hour = tempH;
    startDate->minute = tempMin;
    startDate->second = tempS;

    if((tempY >= 1850) && (tempY <= 2050))
        return 1;
    
    if((tempY >= -2999) && (tempY <= 3000))
        return 2;

    return 0; //This is a flag that the date given is not a valid
              //start date for our equations
}

GLuint sun_texture;

void drawSol()
{
   double radius_au, helioCX, helioCY, helioCZ;
   GLubyte* tex;
   radius_au = planets[8][0];
   helioCX = planets[8][3];
   helioCY = planets[8][4];
   helioCZ = planets[8][5];

   if(solMag < 1)
     solMag = 1.0;

   glPushMatrix();
      glEnable(GL_TEXTURE_2D);
      glTranslatef(helioCX, helioCY, helioCZ - 15.0 );
      #ifdef DEBUG
        cout << "Now making the sphere..." << endl;
      #endif
      glColor4f(1.0,1.0,1.0,1.0);
      tex = textures[8];  /*Get the texture map for the Sun*/
      if(tex != NULL)
      {
            if(sun_texture == 0) {
			  	    glGenTextures( 1, &sun_texture );
				      glBindTexture(GL_TEXTURE_2D, sun_texture );
              glTexEnvf( GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_DECAL ); 
              gluBuild2DMipmaps(GL_TEXTURE_2D, GL_RGB, 512, 256,
     		                        GL_RGB, GL_UNSIGNED_BYTE, tex);
			      } else {
			         glBindTexture(GL_TEXTURE_2D, sun_texture);
			      }         
      }
      gluSphere(sph, radius_au*solMag, 25, 25);
   glPopMatrix();

}

GLuint planet_texture[8];

void drawPlanet(int planetNum)
{
    glEnable(GL_TEXTURE_2D);
    #ifdef DEBUG
      cout << "Attempting to draw planet ";
      cout << planetNum << endl;
    #endif
    GLubyte* tex;
    double radius_au, helioCX, helioCY, helioCZ;
    radius_au = planets[planetNum][0];
    helioCX = planets[planetNum][3];
    helioCY = planets[planetNum][4];
    helioCZ = planets[planetNum][5];
    
    if(planetMag < 1.0)
      planetMag = 1.0;

    #ifdef DEBUG
      cout << "HelioX = ";
      cout << helioCX << endl;
      cout << "HelioY = ";
      cout << helioCY << endl;
      cout << "HelioZ = ";
      cout << helioCZ << endl;
    #endif

    #ifdef VOYAGER_DEBUG
      if((julianDay > 2443390) && (julianDay < 2443393) && (planetNum == 2))
      {
        cout << "Earth Coordinates: X = ";
        cout << helioCX;
        cout << "; Y = ";
        cout << helioCY;
        cout << "; Z = ";
        cout << helioCZ;
        cout << ";" << endl;
      }

      if((julianDay > 2443937) && (julianDay < 2443940) && (planetNum == 4))
      {
        cout << "Jupiter Coordinates: X = ";
        cout << helioCX;
        cout << "; Y = ";
        cout << helioCY;
        cout << "; Z = ";
        cout << helioCZ;
        cout << ";" << endl;
      }

      if((julianDay > 2444554) && (julianDay < 2444557) && (planetNum == 5))
      {
        cout << "Saturn Coordinates: X = ";
        cout << helioCX;
        cout << "; Y = ";
        cout << helioCY;
        cout << "; Z = ";
        cout << helioCZ;
        cout << ";" << endl;
      }
    #endif

    glPushMatrix();
        glTranslatef(helioCX, helioCY, helioCZ - 15.0);
        glColor4f(1.0,1.0,1.0,1.0);
        tex = textures[planetNum];
        if(tex != NULL)
        {
           if(planet_texture[planetNum] == 0) {
				    glGenTextures( 1, &planet_texture[planetNum] );
				    glBindTexture(GL_TEXTURE_2D, planet_texture[planetNum] );
            glTexEnvf( GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_DECAL );    
            gluBuild2DMipmaps(GL_TEXTURE_2D, GL_RGB, 512, 256,
     		                      GL_RGB, GL_UNSIGNED_BYTE, tex);
			      } else {
			         glBindTexture(GL_TEXTURE_2D, planet_texture[planetNum]);
			      }

        }
        gluSphere(sph, radius_au*planetMag,20, 20);
        if(planetNum == 5)
        {
            drawSaturnRings();
            drawSaturnMoons();
        }else if(planetNum == 4)
            drawJupiterMoons();
        else if(planetNum == 2)
            drawMoon();
    glPopMatrix();
}

void drawSaturnRings()
{
    int theta;
    double radius_au = planets[5][0];
    glEnable(GL_LINE_SMOOTH);
    glDisable(GL_LINE_STIPPLE);

    glLineWidth(0.2);
    glBegin(GL_LINE_LOOP);
    for(theta = 0; theta < 360; theta += 10)
    {
      glVertex2f((1.7*radius_au*planetMag)*S(theta),
                 (1.7*radius_au*planetMag)*C(theta));
    }
    glEnd();
    
    glLineWidth(1.5);
    glBegin(GL_LINE_LOOP);
    for(theta = 0; theta < 360; theta += 10)
    {
      glVertex2f((2.0*radius_au*planetMag)*S(theta),
                 (2.0*radius_au*planetMag)*C(theta));
    }
    glEnd();
    
    glLineWidth(0.05);
    glBegin(GL_LINE_LOOP);
    for(theta = 0; theta < 360; theta += 10)
    {
      glVertex2f((1.9*radius_au*planetMag)*S(theta),
                 (1.9*radius_au*planetMag)*C(theta));
    }
    glEnd();

    
    glLineWidth(0.05);
    glBegin(GL_LINE_LOOP);
    for(theta = 0; theta < 360; theta += 10)
    {
      glVertex2f((1.4*radius_au*planetMag)*S(theta),
                 (1.4*radius_au*planetMag)*C(theta));
    }
    glEnd();

    glLineWidth(0.05);
    glBegin(GL_LINE_LOOP);
    for(theta = 0; theta < 360; theta += 10)
    {
      glVertex2f((1.5*radius_au*planetMag)*S(theta),
                 (1.5*radius_au*planetMag)*C(theta));
    }
    glEnd();

    glLineWidth(0.05);
    glBegin(GL_LINE_LOOP);
    for(theta = 0; theta < 360; theta += 10)
    {
      glVertex2f((1.65*radius_au*planetMag)*S(theta),
                 (1.65*radius_au*planetMag)*C(theta));
    }
    glEnd();

    glLineWidth(0.05);
    glBegin(GL_LINE_LOOP);
    for(theta = 0; theta < 360; theta += 10)
    {
      glVertex2f((1.95*radius_au*planetMag)*S(theta),
                 (1.95*radius_au*planetMag)*C(theta));
    }
    glEnd();
}

void drawJupiterMoons()
{
    double current_sm_axis, current_inc, current_radius;
    double x, y, dAngle;

    /*No texture mapping here - it slows down the program*/
    /*First do Io*/
    current_sm_axis = 0.00281889;
    current_inc =     0.05;
    current_radius =  0.000012232784;
    glPushMatrix();
      glRotatef(current_inc, 0.0, 1.0, 0.0);
      if(timeChange != 0)
        dAngle = (timeChange/1.769137786)*360;
      IoAngle += dAngle;
      if(IoAngle > 360)
        IoAngle -= 360;
      
      glDisable(GL_TEXTURE_2D);
      x = current_sm_axis*C(IoAngle);
      y = current_sm_axis*S(IoAngle);
      glTranslatef(x*10*planetMag/5.0, y*10*planetMag/5.0, 0.0);
      glColor3f(205./255., 205./255, 0.0);
      gluSphere(sph, 10*planetMag*current_radius, 14, 14);
    glPopMatrix();

    /*Now for Europa*/
    current_sm_axis = 0.00448558136;
    current_inc =     0.471;
    current_radius =  0.0000104279469;
    glPushMatrix();
      glRotatef(current_inc, 0.0, 1.0, 0.0);
      if(timeChange != 0)
        dAngle = (timeChange/3.551181041)*360;
      EuropaAngle += dAngle;
      if(EuropaAngle > 360)
        EuropaAngle -= 360;
      
      glDisable(GL_TEXTURE_2D);
      x = current_sm_axis*C(EuropaAngle);
      y = current_sm_axis*S(EuropaAngle);
      glTranslatef(x*10*planetMag/5.0, y*10*planetMag/5.0, 0.0);
      glColor3f(139./255., 37./255, 0.0);
      gluSphere(sph, 10*planetMag*current_radius, 14, 14);
    glPopMatrix();
 
    /*And Ganymede*/
    current_sm_axis = 0.00715525609;
    current_inc =     0.204;
    current_radius =  0.0000175871335;
    glPushMatrix();
      glRotatef(current_inc, 0.0, 1.0, 0.0);
      if(timeChange != 0)
        dAngle = (timeChange/7.15455296)*360;
      GanymedeAngle += dAngle;
      if(GanymedeAngle > 360)
        GanymedeAngle -= 360;
      
      glDisable(GL_TEXTURE_2D);
      x = current_sm_axis*C(GanymedeAngle);
      y = current_sm_axis*S(GanymedeAngle);
      glTranslatef(x*10*planetMag/7.0, y*10*planetMag/7.0, 0.0);
      glColor3f(209./255., 209./255., 209./255.);
      gluSphere(sph, 10*planetMag*current_radius, 14, 14);
    glPopMatrix();

    /*Callisto*/
    current_sm_axis = 0.0125851215;
    current_inc =     0.205;
    current_radius =  0.000016109841;
    glPushMatrix();
      glRotatef(current_inc, 0.0, 1.0, 0.0);
      if(timeChange != 0)
        dAngle = (timeChange/16.6890184)*360;
      CallistoAngle += dAngle;
      if(CallistoAngle > 360)
        CallistoAngle -= 360;
      
      glDisable(GL_TEXTURE_2D);
      x = current_sm_axis*C(CallistoAngle);
      y = current_sm_axis*S(CallistoAngle);
      glTranslatef(x*10*planetMag/8.0, y*10*planetMag/8.0, 0.0);
      glColor3f(133./255., 133./255., 133./255.);
      gluSphere(sph, 10*planetMag*current_radius, 14, 14);
    glPopMatrix();
}

void drawSaturnMoons()
{
    double current_sm_axis, current_inc, current_radius;
    double x, y, dAngle;

    /*No texture mapping here - it slows down the program*/
    /*First do Io*/
    current_sm_axis = 0.0081685584;
    current_inc =     0.33;
    current_radius =  0.000017212797;
    glPushMatrix();
      glRotatef(current_inc, 0.0, 1.0, 0.0);
      if(timeChange != 0)
        dAngle = (timeChange/1.769137786)*360;
      TitanAngle += dAngle;
      if(TitanAngle > 360)
        TitanAngle -= 360;
      
      glDisable(GL_TEXTURE_2D);
      x = current_sm_axis*C(TitanAngle);
      y = current_sm_axis*S(TitanAngle);
      glTranslatef(x*10*planetMag/5.0, y*10*planetMag/5.0, 0.0);
      glColor3f(205./255., 205./255, 0.0);
      gluSphere(sph, 10*planetMag*current_radius, 14, 14);
    glPopMatrix();
}
void drawMoon()
{
  /*Note: This is NOT an accurate representation of the moon.
   *      It is simply a start to adding satellites to the system
   */
    double dMoonAngle =      0.0;
    double MOON_RADIUS =     0.0000116117863;
    double MOON_SM_AXIS =    0.00257;
    double MOON_INC =        5.14500;
    glPushMatrix();
      glRotatef(MOON_INC, 0.0, 1.0, 0.0);
      if(timeChange != 0)
       dMoonAngle = (timeChange/ 27.321582)*360;
      moonAngle += dMoonAngle;
      if(moonAngle > 360)
       moonAngle -= 360;

      glDisable(GL_TEXTURE_2D);
      double x, y;
      x = MOON_SM_AXIS*C(moonAngle);
      y = MOON_SM_AXIS*S(moonAngle);
      glTranslatef(x*planetMag, y*planetMag, 0.0);
      glColor3f(0.47, 0.53, 0.6);
      gluSphere(sph, 10*planetMag*MOON_RADIUS, 14, 14);
    glPopMatrix(); 
}



void drawOrbit(int planet)
{
   double fastJulian = julianDay;
   double period = 0.0;
   switch(planet){
      case 0: 
        period = 90.0;
        break;
      case 1:
        period = 225.0;
        break;
      case 2:
        period = 366.0;
        break;
      case 3:
        period = 688.0;
        break;
      case 4:
        period = 4332.0;
        break;
      case 5:
        period = 10833.0;
        break;
      case 6:
        period = 30800.0;
        break;
      case 7:
        period = 60191.0;
        break; 
   };

   /*int temp;
   temp = (int)(julianDay / period);
   fastJulian = julianDay - ((float)(temp)) * period;
   if(fastJulian > 10)
     return;*/
 
   glEnable(GL_LINE_SMOOTH);
   glBegin(GL_LINE_STRIP);
       glColor3f(1.0, 1.0, 1.0);
       glLineWidth(5.0);   
       for(fastJulian = 0.0; fastJulian <= period; fastJulian += 50)
       {
         cout << "Placing an orbital vertex..." << endl;
         compute_Position(julianDay + fastJulian, planet+1, era, 1); 
         glVertex3d(planetsOrbit[planet][0],
                    planetsOrbit[planet][1],
                    planetsOrbit[planet][2]);
       }
   glEnd();
}

void initPlanets()
{
      /*Initialize the radius of all the planets
       *(The Sun was already initialized */
      planets[0][0] = MERCURY_RADIUS;
      planets[1][0] = VENUS_RADIUS;
      planets[2][0] = EARTH_RADIUS;
      planets[3][0] = MARS_RADIUS;
      planets[4][0] = JUPITER_RADIUS;
      planets[5][0] = SATURN_RADIUS;
      planets[6][0] = URANUS_RADIUS;
      planets[7][0] = NEPTUNE_RADIUS;

      /*Load in the images for texture mapping to the
       *planets*/
      textures[0] = loadTexture("iSImages/mercury.raw", 512, 256, 3); 
      textures[1] = loadTexture("iSImages/venus.raw", 512, 256, 3); 
      textures[2] = loadTexture("iSImages/earth.raw", 512, 256, 3); 
      textures[3] = loadTexture("iSImages/mars.raw", 512, 256, 3); 
      textures[4] = loadTexture("iSImages/jupiter.raw", 512, 256, 3); 
      textures[5] = loadTexture("iSImages/saturn.raw", 512, 256, 3); 
      textures[6] = loadTexture("iSImages/uranus.raw", 512, 256, 3); 
      textures[7] = loadTexture("iSImages/neptune.raw", 512, 256, 3); 
      textures[8] = loadTexture("iSImages/sun.raw", 512, 256, 3); 
}


/*Here is a texture mapping function I got from an online tutorial - 
 *Thanks to UI Chicago for posting the shell program for sphere
 *mapping to the gluSphere and the raw image files for the job
 *http://www.evl.uic.edu/julian/cs488/2005-11-10/index.html
 */ 
GLubyte* loadTexture (char *fileName, int width, int height, int depth)
{
  GLubyte *raw_bitmap, *reversed_bitmap;
  FILE *texFile;

  int byteSize, textureSize;
  byteSize = sizeof(GLubyte);
  textureSize = width * height * depth * byteSize;

  texFile = fopen(fileName, "rb");

  if (texFile == NULL)
  {
      printf ("Texture %s not found\n", fileName);
      //exit(1);
      return NULL;
  }

  raw_bitmap = (GLubyte *) malloc (textureSize);
  reversed_bitmap = (GLubyte *) malloc (textureSize);

  if ((raw_bitmap == NULL) || (reversed_bitmap == NULL))
  {
      printf("Cannot allocate memory for texture %s\n", fileName);
      fclose(texFile);
      exit(1);
  }

  fread (raw_bitmap , width * height * depth, 1 , texFile );
  fclose (texFile);

  /* need to reverse the texture horizontally here */
  /*************************************************/

  int i, j, k;

  for (i=0; i<height;i++)
    for (j=0; j<width*depth;j+=depth)
      for (k=0;k<depth;k++)
	reversed_bitmap[i*width*depth + (j+k)] =
	                  raw_bitmap[(i+1)*width*depth-(j+depth)+k];


  glTexParameteri (GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
  glTexParameteri (GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
  glTexEnvf (GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE );

  free (raw_bitmap);

  return (reversed_bitmap);
}

void compute_Voyager_Position()
{
    double relativeDate, tempor;
    relativeDate = (julianDay - 2451545.0)/36525.0;
    
    /*If it hasn't launched yet, put it in the middle
     *of the Earth... just make sure it never gets
     *bigger than the Earth. */
    if(julianDay < 2443375.0)
    {
      voyager[0][0] = planets[2][3];
      voyager[0][1] = planets[2][4];
      voyager[0][2] = planets[2][5];
      return;
    }

    /* Otherwise, compute the Voyager II's position */

    tempor = 36239*(relativeDate)*(relativeDate);
    tempor = tempor + 15319*relativeDate + 1613.5;
    voyager[0][0] = tempor;

    tempor = -6120.8*(relativeDate)*(relativeDate);
    tempor = tempor - 2531.9*relativeDate - 260.59;
    voyager[0][2] = tempor;

    tempor = 6504.5*(relativeDate)*(relativeDate);
    tempor += 2824.9*relativeDate + 306.48;
    voyager[0][1] = tempor;
}

void drawVoyager()
{
    glEnable(GL_TEXTURE_2D);
    #ifdef DEBUG
      cout << "Attempting to draw planet ";
      cout << planetNum << endl;
    #endif
    //GLubyte* tex;
    double helioCX, helioCY, helioCZ;
    helioCX = voyager[0][0];
    helioCY = voyager[0][1];
    helioCZ = voyager[0][2];
    
    #ifdef DEBUG
      cout << "HelioX = ";
      cout << helioCX << endl;
      cout << "HelioY = ";
      cout << helioCY << endl;
      cout << "HelioZ = ";
      cout << helioCZ << endl;
    #endif

    glPushMatrix();
        glTranslatef(helioCX, helioCY, helioCZ - 15.0);
        glColor4f(1.0,1.0,1.0,1.0);

        /*If you wanted to do a texture for the voyager,
         *it would go here.... use the kind of format I
         *used for the spheres, but I don't know how to
         *do it for any other kind of shape */

        gluSphere(sph, .0005*planetMag,20, 20);
    glPopMatrix();
}
