/* NIGHTFALL OpenGL Interface                                              */
/* Copyright (C) 2001 Rainer Wichmann & Markus Kuster                      */
/*                                                                         */
/*  This program is free software; you can redistribute it                 */
/*  and/or modify                                                          */
/*  it under the terms of the GNU General Public License as                */
/*  published by                                                           */
/*  the Free Software Foundation; either version 2 of the License, or      */
/*  (at your option) any later version.                                    */
/*         
                                                                */
/*  This program is distributed in the hope that it will be useful,        */
/*  but WITHOUT ANY WARRANTY; without even the implied warranty of         */
/*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          */
/*  GNU General Public License for more details.                           */
/*                                                                         */
/*  You should have received a copy of the GNU General Public License      */
/*  along with this program; if not, write to the Free Software            */
/*  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.              */

/* ANSI C forbids an empty source file, so put this outside                */
/* do nothing here if we don't have OpenGL                                 */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "Light.h"

#ifdef  _WITH_OPENGL
#include "LightGLPrefs.h"


/******************************************************************
 @package   nightfall
 @author    Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
 @version   1.3
 @short     The color wedge code 
 @param     (int)  Comp the component
 @return    (void)   
 @heading   Open GL Animation
*******************************************************************/
void GLDisplayWedge( int Comp )
{
  float maxval = texture_maxval[Comp];
  float minval = texture_minval[Comp];
  char  label[128];

  /* plot borders
   */
  float xmin   =  -0.5;
  float lmin   =  -0.9;
  float xmax   =  0.5;
  float ymin   = -0.5;
  float ymax   =  0.5;

  double        ScaleMass;          /* units conversion              */
  double G_SI  = 6.6726e-11;        /* gravitational constant        */
  char   texwhat[16];

  struct GLColorSystem NTSC = {0.64, 0.33, 0.29, 0.60, 0.15, 0.06, CWP}; 
  struct GLColorSystem *cs; 
  double  x, y, z, R, G, B, max;
  
  /* R.W. 23.12.2002 set colours to BB colour
   */
  cs     = &NTSC;
  GLBBtoxyz(Binary[Comp].Temperature, &x, &y, &z);
  GLxyztoRGB(cs, x, y, z, &R, &G, &B);
  max = (R > G) ? R : G;
  max = (B > max) ? B : max;
  R = (1.0 / max) * R;
  G = (1.0 / max) * G;
  B = (1.0 / max) * B;

  if (Flags.textype == TEMP)
    {
      strcpy(texwhat, "Kelvin");
    }
  else if (Flags.textype == GRAV)
    {
      ScaleMass = Binary[Primary].Mq / ( 1.0 + Binary[Primary].Mq);
      minval = minval * Orbit.TrueMass * (1.0 - ScaleMass)
	* G_SI
	/ (Orbit.TrueDistance * Orbit.TrueDistance);
      maxval = maxval * Orbit.TrueMass * (1.0 - ScaleMass)
	* G_SI
	/ (Orbit.TrueDistance * Orbit.TrueDistance);
      strcpy(texwhat, "m kg/s^2");
    }
  else if (Flags.textype == FLUX)
    {
      minval = (1.0 / maxval) * minval;
      maxval = 1.0;
      strcpy(texwhat, "relative flux");
    }
  else if (Flags.textype == VELOCITY)
    {
      strcpy(texwhat, "m/s");
    }
  else
    {
      strcpy(texwhat, "?");
    }


  if (gtk_gl_area_make_current(GTK_GL_AREA(glArea)))
    {
      /* disable lighting, we don't want ligthing with */
      /* the axes                                      */
      glDisable(GL_LIGHTING);

      glMatrixMode(GL_PROJECTION);
      glLoadIdentity();
      glOrtho(-1, 1, -1, 1, -1, 1);
      glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);

      /* draw border for plot                                              */
      glBegin(GL_QUADS); 
      glColor3f(0.0, 0.0, 0.0);
      glVertex3f(xmin, ymin, 0.0);
      glVertex3f(xmin, ymax, 0.0);
      glColor3f(R, G, B); 
      glVertex3f(xmax, ymax, 0.0);
      glVertex3f(xmax, ymin, 0.0);
      glEnd(); 
      glPopMatrix();

#if 0
      /* draw  tickmarks                                                   */
      glBegin(GL_LINES); 
      glVertex2f(-0.5, bot_xtick);       
      glVertex2f(-0.5, top_xtick);       
      glVertex2f( 0.0, bot_xtick);       
      glVertex2f( 0.0, top_xtick);       
      glVertex2f( 0.5, bot_xtick);       
      glVertex2f( 0.5, top_xtick);
      glVertex2f( 1.0, bot_xtick);
      glVertex2f( 1.0, top_xtick);
      glEnd();
#endif

      /* select font                                   */
      GLSetFont(HELVETICA, 12);

      glDisable(GL_LIGHTING);
      
      glPushMatrix();

      /* print labels                                  */
      glColor3f(1.0, 1.0, 1.0);
      snprintf(label,  127, "%c: %10.4g",
	       (Comp == Primary) ? 'P' : 
	       (Comp == Secondary) ? 'S' : 'D', minval);
      GLDrawStr(lmin, -0.2,      0.0, label);
      snprintf(label,  127, "%10.4g (%s)", maxval, texwhat);
      GLDrawStr(xmax, -0.2,      0.0, label);

#if 0
      /* print tick marks                              */
      GLDrawStr(-0.545,bot_xtick-0.1*length_y,0.0,"0.5");
      GLDrawStr(-0.05,bot_xtick-0.1*length_y,0.0,"0.0");
      GLDrawStr( 0.455,bot_xtick-0.1*length_y,0.0,"0.5");
      GLDrawStr( 0.96,bot_xtick-0.1*length_y,0.0,"1.0");

      /* print title string                            */
      glColor3f(1.0, 1.0, 1.0);
      GLDrawStr(coord.titlepos,yaxis.max_value+0.15*length_y,0.0,coord.title);
#endif
      glPopMatrix();

      /* switch on lighing                             */
      glEnable(GL_LIGHTING);
      
      /* set line width to 1                           */
      glLineWidth(1);
    }
}


/******************************************************************
 @package   nightfall
 @author    Markus Kuster (kuster@astro.uni-tuebingen.de)
 @version   0.5
 @short     The main 3D window rendering code 
 @param     (GtkWidget)  *widget  discarded
 @return    (void)   
 @heading   Open GL Animation
*******************************************************************/
void GLDisplay3d( GtkWidget *widget )
{
  float  com;                         /* center of mass          */

  GLfloat m[4][4],h;
  OrbitType  *OrbPtr = &Orbit;  /* pointer to Orbit               */
  DispInfo *DispInfoPtr;

  struct GLColorSystem NTSC = {0.64, 0.33, 0.29, 0.60, 0.15, 0.06, CWP}; 
  struct GLColorSystem *cs; 
  double  x, y, z, R, G, B, max;
  
  /* R.W. 23.12.2002 set colours to BB colour
   */
  cs     = &NTSC;
  GLBBtoxyz(Binary[Primary].Temperature, &x, &y, &z);
  GLxyztoRGB(cs, x, y, z, &R, &G, &B);
  max = (R > G) ? R : G;
  max = (B > max) ? B : max;
  primary_colour[0] = (1.0 / max) * R;
  primary_colour[1] = (1.0 / max) * G;
  primary_colour[2] = (1.0 / max) * B;

  cs     = &NTSC;
  GLBBtoxyz(Binary[Secondary].Temperature, &x, &y, &z);
  GLxyztoRGB(cs, x, y, z, &R, &G, &B);
  max = (R > G) ? R : G;
  max = (B > max) ? B : max;
  secondary_colour[0] = (1.0 / max) * R;
  secondary_colour[1] = (1.0 / max) * G;
  secondary_colour[2] = (1.0 / max) * B;

#ifdef HAVE_DISK
  cs     = &NTSC;
  GLBBtoxyz(Binary[Disk].Temperature, &x, &y, &z);
  GLxyztoRGB(cs, x, y, z, &R, &G, &B);
  max = (R > G) ? R : G;
  max = (B > max) ? B : max;
  disk_colour[0] = (1.0 / max) * R;
  disk_colour[1] = (1.0 / max) * G;
  disk_colour[2] = (1.0 / max) * B;
#endif

  DispInfoPtr = (DispInfo*)gtk_object_get_data(GTK_OBJECT(glArea), "DispInfo");
  /* get the dimension of the GTK window      */

  h= (GLfloat) glArea->allocation.height / (GLfloat) glArea->allocation.width;

  if (gtk_gl_area_make_current(GTK_GL_AREA(glArea)))
    {
      glMatrixMode(GL_PROJECTION);
      glLoadIdentity();

      glOrtho((-1 - DispInfoPtr->xpos)  * DispInfoPtr->zoom,
	      ( 1 - DispInfoPtr->xpos)  * DispInfoPtr->zoom,
	       -1 * DispInfoPtr->zoom,1 * DispInfoPtr->zoom,.1,1000.);

      glMatrixMode(GL_MODELVIEW);
      
      /* the x coordinate of the center of mass                      */
      com = (((float)(Binary[Primary].Mq))/( 1.0 + (float)(Binary[Primary].Mq)));

      /* load identity transformation matrix            */
      glLoadIdentity();

      /* set the drawing mode to polygon mode           */
      if (Flags.wireframe == ON) {
	glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
      } else {
	glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
      }

      /* set view point position to (Flags.GLdistancex,0,Flags.GLdistancez)*/
      /* look to point              (Flags.GLdistancex,0,       0         )*/
      /* with head up pointing to   (        0        ,1,       0         )*/
      gluLookAt(Flags.GLdistancex,0,Flags.GLdistancez,Flags.GLdistancex,0,0,0,1,0);

      build_rotmatrix(m,DispInfoPtr->quat);
      glMultMatrixf(&m[0][0]);
      
      /* rotate orbit accoding to the inclination       */
      glRotated((270.0-((OrbPtr->Inclination-PI/2.0)*180.0/PI)),1,0,0);

      /* rotate phase steps                             */
      glRotated((OrbPtr->Phase-0.5*PI)*180.0/PI,0,0,1);
     
      glTranslatef(-(com),0.0,0.0);

      if (Flags.axes == ON) {
	/* draw x-,y- and z-axes                        */
	glPushMatrix();
	GLDrawAxes(0.4,glArea);
	glPopMatrix();
      }
      if (Flags.labels == ON) {
	glPushMatrix();
	GLDrawLabels(glArea);
	glPopMatrix();
      }

      /* save transformation matix                      */
      glPushMatrix();

      /* set the material options for primary           */
      /* primary emission; should be set according      */
      /* to the primary temperature                     */
      glMaterialfv(GL_FRONT, GL_EMISSION,  primary_colour);
      glMaterialfv(GL_FRONT, GL_DIFFUSE,   primary_mat_diffuse);
      glMaterialfv(GL_FRONT, GL_SPECULAR,  primary_mat_specular);
      glMaterialfv(GL_FRONT, GL_SHININESS, primary_mat_shininess);
      glFrontFace(GL_CCW);

      if ((Flags.texture == ON) && (Texture[Primary].IsOn == ON)) {
	/* enable texturing                                          */
	glEnable(GL_TEXTURE_2D);
	/* assign texture to the primary display list   */
	glBindTexture(GL_TEXTURE_2D, Texture[Primary].TexName);
	/* call primary lobe display list               */
	glCallList(PrimaryList);
	glBindTexture(GL_TEXTURE_2D, 0);
	/* switch off texturing                         */
	glDisable(GL_TEXTURE_2D);
      } else {
	/* call primary lobe display list               */
	glCallList(PrimaryList);
      }

      /* restore transformation matrix                  */
      glPopMatrix();
      
      /* save transformation matrix                      */
      glPushMatrix();

      /* translate to the position of the secondary                */
      glTranslatef (1.0, 0.0, 0.0);

      /* set the material options for secondary         */
      glMaterialfv(GL_FRONT, GL_EMISSION, secondary_colour);
      glFrontFace(GL_CW);

      if ((Flags.texture == ON) && (Texture[Secondary].IsOn == ON)) {
	/* enable texturing                             */
	glEnable(GL_TEXTURE_2D);
	/* assign texture to the secondary display list */
	glBindTexture(GL_TEXTURE_2D, Texture[Secondary].TexName);
	/* call secondary lobe display list             */
	glCallList(SecondaryList);
	glBindTexture(GL_TEXTURE_2D, 0);
	/* switch off texturing                         */
	glDisable(GL_TEXTURE_2D);
      } else {
	/* call secondary lobe display list             */
	glCallList(SecondaryList);
      }
      glPopMatrix();
      
#ifdef HAVE_DISK
      glPushMatrix();
      /* set the material options for the disk          */
      /* primary emission; should be set according      */
      /* to the disk temperature                        */
      glMaterialfv(GL_FRONT, GL_EMISSION,  disk_colour);
      glMaterialfv(GL_FRONT, GL_DIFFUSE,   disk_mat_diffuse);
      glMaterialfv(GL_FRONT, GL_SPECULAR,  disk_mat_specular);
      glMaterialfv(GL_FRONT, GL_SHININESS, disk_mat_shininess);

      /* translate to the position of the secondary                */
      glTranslatef (1.0, 0.0, 0.0);
      /* rotate in opposite direction according to orbital motion  */
      glRotated(-(OrbPtr->Phase-0.5*PI)*180.0/PI,0,0,1);

      if (Flags.disk == ON ) {
	if ((Flags.texture == ON) && ( Texture[Disk].IsOn == ON)){
	  /* enable texturing                           */
	  glEnable(GL_TEXTURE_2D);
	  /* assign texture to the disk display list */
	  glBindTexture(GL_TEXTURE_2D, Texture[Disk].TexName);
	  /* call accretion disk display list           */
	  glCallList(DiskList);
	  glBindTexture(GL_TEXTURE_2D, 0);
	  /* switch off texturing                       */
	  glDisable(GL_TEXTURE_2D);
	} else {
	  glCallList(DiskList);
	}
      }
      /* restore transformation matrix                  */
      glPopMatrix();
#endif 

      glPushMatrix();
      /* set light positions                            */
      glLightfv(SCENE_LIGHT, GL_POSITION      , Lights[SCENE_LIGHT].Pos);
      glLightfv(DISK_SPOT  , GL_POSITION      , Lights[DISK_SPOT].Pos);      
      glLightfv(DISK_SPOT  , GL_SPOT_DIRECTION, Lights[DISK_SPOT].Direction);
      /*glLightfv(GL_LIGHT2, GL_POSITION, lightpos);*/
      /*glLightfv(GL_LIGHT2, GL_SPOT_DIRECTION, primary_spot2_direction);*/
      glPopMatrix();
    }
}

/******************************************************************
 @package   nightfall
 @author    Markus Kuster (kuster@astro.uni-tuebingen.de)
 @version   0.5
 @short     The main light curve/rad. vel. rendering code 
 @param     (GtkWidget)  *widget  discarded
 @return    (void)   
 @heading   Open GL Animation
*******************************************************************/
void GLDisplayPlot( GtkWidget *widget, int type)
{
  AxisType xaxis;
  AxisType yaxis;
  
  xaxis = CoordSys[type].xaxis;
  yaxis = CoordSys[type].yaxis;

  if (gtk_gl_area_make_current(GTK_GL_AREA(glArea)))
    {
      glMatrixMode(GL_PROJECTION);
      glLoadIdentity();
      
      glOrtho(xaxis.min_value-0.1*xaxis.length, xaxis.max_value+0.1*xaxis.length, 
	      yaxis.min_value-0.3*yaxis.length, yaxis.max_value+0.2*yaxis.length,
	      -2.0,2.0); 
      
      glMatrixMode(GL_MODELVIEW);
      glLoadIdentity();

      /* set the drawing mode to polygon mode           */
      glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);

      /* Call plot list                                 */
      glCallList(PlotList[type]);
    }
}

/******************************************************************
 @package   nightfall
 @author    Markus Kuster (kuster@astro.uni-tuebingen.de)
 @version   1.0
 @short     The light curve drawing code 
 @param     (void)
 @return    (void)   
 @heading   Open GL Animation
*******************************************************************/
void GLMakeLCPlot( void )
{
  double    Calib;                 /* normalization for flux         */
  int       i,test;
  float     *Xplot,*Yplot,Xtick;   /* x, y position                  */
  float     MinMag, MaxMag;        /* Min,Max scaling for plot       */
  PhotoOut  *FluxPtr = FluxOut;    /* pointer to flux                */
  int    js;                       /* j minus invalid                */

  if (gtk_gl_area_make_current(GTK_GL_AREA(glArea)))
    {
      /* Initialize memory for plot                                  */
      Xplot = malloc(3 * MAXELEMENTS * sizeof(float)); 
      Yplot = malloc(3 * MAXELEMENTS * sizeof(float)); 
      
      if (Xplot == NULL || Yplot == NULL) 
#ifdef _WITH_GTK 
  	{  
	  if (Flags.interactive == ON) { 
  	    make_dialog (_(errmsg[0])); 
  	    if (Xplot != NULL) free(Xplot); 
  	    if (Yplot != NULL) free(Yplot); 
  	    return; 
  	  } else nf_error(_(errmsg[5])); 
  	} 
#else 
      nf_error(_(errmsg[5])); 
#endif 
      
      /* >>>>>>>>>>>> Light Curve <<<<<<<<<<<<<<<<<<<<<<<<<< */
      /* photometric zeropoint                               */
      Calib = FluxOut[0].Mag[Vmag];

      /* get photometric data and scale                      */
      js = 0;
      MinMag  = 10.;  
      MaxMag  = -10.; 
      FluxPtr = FluxOut; 

      /* 2002-05-03 rwichmann: the problem is that 
	 lightcurve normalization has been done already, resulting
	 in MinMag = nan  MaxMag = nan.
      */
      for (i = 0; i <= phaseind; ++i) { 
	if (i == 0 || FluxPtr->Mag[Vmag] != 0)
	  {
	    if (fabs(Calib) < DBL_EPSILON) {
	      /* already normalized                              */
	      Yplot[js] = - FluxPtr->Mag[Vmag];
	    }
	    else { 
	      /* normalize and transform flux -> magnitude       */
	      Yplot[js] = 2.5 * log10(FluxPtr->Mag[Vmag]/Calib);
	    }
	    MinMag  = MIN(MinMag,Yplot[js]); MaxMag  = MAX(MaxMag,Yplot[js]); 
	    Xplot[js] = (FluxPtr->Phase/(PI+PI)) - 0.5;
	    ++js;
	  }
  	++FluxPtr;  
      } 

      /* shift phase if necessary                            */
      test = ON;
      do { 
  	if (Xplot[0] > 0.0) { 
  	  for (i = 0; i <= js; ++i) { 
  	    Xplot[i] = Xplot[i] - 1.0; 
 	  } 
	  test = ON;  /* there was something to do           */
  	} else { 
  	  if (Xplot[0] < (-1.0)) { 
  	    for (i = 0; i <= js; ++i) { 
  	      Xplot[i] = Xplot[i] + 1.0; 
  	    } 
	    test = ON;  /* there was something to do         */ 
  	  } else { test = OFF; }
  	} 
      } while (test == ON); 
      
      if (fabs(MinMag-MaxMag) < 0.1) MinMag = -0.1; 

      /* plot box                                            */
      Xtick = (int)(100.*(MaxMag-MinMag+0.1)/3.)/100.; 

      CoordSys[LIGHTCURVE] = 
	SetupCoords(MinMag,MaxMag,"Phase","M","Light Curve");
      
      /* plot coordinate system                              */
      GLMake2DAxes(CoordSys[LIGHTCURVE], glArea);

      glPointSize(2);
      glDisable(GL_LIGHTING);
      
      glColor3f(1.0, 1.0, 0.0);
      glBegin(GL_POINTS);
      for (i = 0; i < js; ++i) { 
	glVertex2f(Xplot[i], Yplot[i]);
      }
      glEnd();

      glEnable(GL_LIGHTING);
      glPointSize(1);
     
      free(Xplot);
      free(Yplot);
    }
}

/******************************************************************
 @package   nightfall
 @author    Markus Kuster (kuster@astro.uni-tuebingen.de)
 @version   1.0
 @short     The radial velocity drawing code 
 @param     (void)
 @return    (void)   
 @heading   Open GL Animation
*******************************************************************/
void GLMakeRVPlot ( void )
{
  double    Calib;                 /* normalization for flux         */
  int       i,test;
  float     *Xplot,*Yplot,*Zplot;  /* x, y, z position                  */
  float     MinVel, MaxVel;        /* plot scaling           */
  PhotoOut  *FluxPtr = FluxOut;    /* pointer to flux                */
  int    js;                       /* j minus invalid                */

  if (gtk_gl_area_make_current(GTK_GL_AREA(glArea)))
    {
      /* Initialize memory for plot                               */
      /* Initialize memory for plot                               */
      Xplot = malloc(3 * MAXELEMENTS * sizeof(float)); 
      Yplot = malloc(3 * MAXELEMENTS * sizeof(float)); 
      Zplot = malloc(3 * MAXELEMENTS * sizeof(float)); 

      if (Xplot == NULL || Yplot == NULL) 
#ifdef _WITH_GTK 
  	{  
	  if (Flags.interactive == ON) { 
  	    make_dialog (_(errmsg[0])); 
  	    if (Xplot != NULL) free(Xplot); 
  	    if (Yplot != NULL) free(Yplot); 
  	    return; 
  	  } else nf_error(_(errmsg[5])); 
  	} 
  #else 
        nf_error(_(errmsg[5])); 
#endif 
      /* >>>>>>>>>>>> Light Curve <<<<<<<<<<<<<<<<<<<<<<<<<< */
      /* photometric zeropoint                               */
      Calib = FluxOut[0].Mag[Vmag];
      
      /* get photometric data and scale                      */
      js = 0;
      MinVel  = (FluxPtr->RV[Primary])  /1000.;
      MaxVel  = (FluxPtr->RV[Secondary])/1000.; 

      if (Flags.elliptic == ON) {
	for (i = 0; i <= phaseind; ++i) {
	  if (i == 0 || FluxPtr->Mag[Vmag] != 0)
	    {
	      Xplot[js] = (FluxPtr->Phase/(PI+PI)) - 0.5;  
	      Yplot[js] =   (FluxPtr->RV[Primary])  /1000.;
	      Zplot[js] =   (FluxPtr->RV[Secondary])/1000.;
	      MinVel  = MIN(MinVel,Yplot[js]); MaxVel = MAX(MaxVel,Yplot[js]);
	      MinVel  = MIN(MinVel,Zplot[js]); MaxVel = MAX(MaxVel,Zplot[js]);
	      ++js;
	    }
	  ++FluxPtr;
	}
      } else {
	for (i = 0; i <= phaseind; ++i) {
	  if (i == 0 || FluxPtr->Mag[Vmag] != 0)
	    {
	      Xplot[js] = (FluxPtr->Phase/(PI+PI)) - 0.5;  
	      Yplot[js] = - (FluxPtr->RV[Primary])  /1000.;
	      Zplot[js] = - (FluxPtr->RV[Secondary])/1000.;
	      MinVel  = MIN(MinVel,Yplot[js]); MaxVel  = MAX(MaxVel,Yplot[js]);
	      MinVel  = MIN(MinVel,Zplot[js]); MaxVel  = MAX(MaxVel,Zplot[js]);
	      ++js;
	    }
	  ++FluxPtr;
	}
      }

      /* shift phase if necessary                            */
      test = ON;
      do { 
  	if (Xplot[0] > 0.0) { 
  	  for (i = 0; i <= js; ++i) { 
  	    Xplot[i] = Xplot[i] - 1.0; 
 	  } 
	  test = ON;  /* there was something to do           */
  	} else { 
  	  if (Xplot[0] < (-1.0)) { 
  	    for (i = 0; i <= js; ++i) { 
  	      Xplot[i] = Xplot[i] + 1.0; 
  	    } 
	    test = ON;  /* there was something to do         */ 
  	  } else { test = OFF; }
  	} 
      } while (test == ON); 
     
      if (fabs(MinVel-MaxVel) < 5.) { 
	MinVel = -5.01;
	MaxVel = 5.01; 
      }

      /* plot box                                            */
      CoordSys[RADVELOCITY] = 
	SetupCoords(MinVel, MaxVel,"Phase","v","Radial Velocity");

      /* plot coordinate system                         */
      GLMake2DAxes(CoordSys[RADVELOCITY], glArea);

      glPointSize(2);
      glDisable(GL_LIGHTING);
      
      glBegin(GL_POINTS);
      for (i = 0; i < js; ++i) { 
	/* Primary data (red)               */
	glColor3f(1.0, 0.37, 0.0);
	glVertex2f(Xplot[i], Yplot[i]);
	/* Secondary data (blue)            */
	glColor3f(0.6, 0.86, 1.0);
	glVertex2f(Xplot[i], Zplot[i]);
      }
      glEnd();

      glEnable(GL_LIGHTING);
      glPointSize(1);
     
      free(Xplot);
      free(Yplot);
    }
}

/******************************************************************
 @package   nightfall
 @author    Markus Kuster (kuster@astro.uni-tuebingen.de)
 @version   1.0
 @short     Draw x-, y- and z-axes of the binary systems reference 
            coordinate system according to Djurasevic, 1991
	    Note: this is not the OpenGL Viewing coordinate 
	          system !!
 @param     (GLfloat) length  length of the axes
 @return    (void)   
 @heading   Open GL Animation
*******************************************************************/
void GLDrawAxes(GLfloat length, GtkWidget *widget)
{
  float com;

  /* the x coordinate of the center of mass                      */
  com = (((float)(Binary[Primary].Mq))/( 1.0 + (float)(Binary[Primary].Mq)));

  if (gtk_gl_area_make_current(GTK_GL_AREA(glArea)))
    {
      /* draw axes with thick lines                    */
      glLineWidth(2);
      
      /* disable lighting, we don't want ligthing with */
      /* the axes                                      */
      glDisable(GL_LIGHTING);
      
      /* draw axes as lines                            */
      glBegin(GL_LINES);
      /* blue x-axis                                   */
      glColor3f(0.0, 0.0, 1.0);
      glVertex3f(com, 0.0, 0.0);
      glVertex3f(com+length, 0.0, 0.0);

      /* green y-axis                                  */
      glColor3f(0.0, 1.0, 0.0);
      glVertex3f(com, 0.0, 0.0);
      glVertex3f(com, -length, 0.0);
      
      /* red z-axis                                    */
      glColor3f(1.0, 0.0, 0.0);
      glVertex3f(com, 0.0, 0.0);
      glVertex3f(com, 0.0, length);
      glEnd();
      
      /* switch on lighing                             */
      glEnable(GL_LIGHTING);
      
      /* set line width to 1                           */
      glLineWidth(1);
    }
}

/******************************************************************
 @package   nightfall
 @author    Markus Kuster (kuster@astro.uni-tuebingen.de)
 @version   1.1
 @short     Build a display list for x- and y-axes of the 
            coordinate system
 @param     (GLfloat) length  length of the axes
 @return    (void)   
 @heading   Open GL Animation
*******************************************************************/
void GLMake2DAxes(CoordType coord, GtkWidget *widget)
{
  float length_x, length_y, eye_x;
  float top_xtick, bot_xtick;
  float xmin_border,xmax_border,ymin_border,ymax_border;

  AxisType xaxis;
  AxisType yaxis;
  
  xaxis = coord.xaxis;
  yaxis = coord.yaxis;

  /* length of x- and y-axes                           */
  length_x = xaxis.length;
  length_y = yaxis.length;

  if (gtk_gl_area_make_current(GTK_GL_AREA(glArea)))
    {
      ymin_border = yaxis.min_value-coord.yborder_size;
      ymax_border = yaxis.max_value+coord.yborder_size;

      xmin_border = xaxis.min_value-coord.xborder_size;
      xmax_border = xaxis.max_value+coord.xborder_size;

      /* viewers position                              */
      eye_x     = xaxis.min_value + 0.5*length_x;
      /* tickmarks on x-axis                           */
      top_xtick = ymin_border     + 0.02*length_y;
      bot_xtick = ymin_border     - 0.02*length_y;

      /* draw axes with thick lines                    */
      glLineWidth(2);
      
      /* disable lighting, we don't want ligthing with */
      /* the axes                                      */
      glDisable(GL_LIGHTING);

      glColor3f(1.0, 1.0, 1.0); 
      /* draw border for plot                                              */
      glBegin(GL_LINES); 
      glVertex2f(xmin_border, ymin_border);
      glVertex2f(xmin_border, ymax_border);      

      glVertex2f(xmin_border, ymax_border);
      glVertex2f(xmax_border, ymax_border);

      glVertex2f(xmax_border, ymax_border);
      glVertex2f(xmax_border, ymin_border);

      glVertex2f(xmax_border, ymin_border);
      glVertex2f(xmin_border, ymin_border);
      glEnd(); 

      /* draw  tickmarks                                                   */
      glBegin(GL_LINES); 
      glVertex2f(-0.5, bot_xtick);       
      glVertex2f(-0.5, top_xtick);       
      glVertex2f( 0.0, bot_xtick);       
      glVertex2f( 0.0, top_xtick);       
      glVertex2f( 0.5, bot_xtick);       
      glVertex2f( 0.5, top_xtick);
      glVertex2f( 1.0, bot_xtick);
      glVertex2f( 1.0, top_xtick);
      glEnd();

      /* select font                                   */
      GLSetFont(HELVETICA, 12);

      glDisable(GL_LIGHTING);
      
      glPushMatrix();

      /* print labels for axes                         */
      glColor3f(1.0, 1.0, 1.0);
      GLDrawStr(eye_x-0.1,yaxis.min_value-0.25*length_y,0.0,xaxis.label);
      GLDrawStr(-0.6,yaxis.min_value+0.5*length_y,0.0,yaxis.label);

      /* print tick marks                              */
      GLDrawStr(-0.545,bot_xtick-0.1*length_y,0.0,"0.5");
      GLDrawStr(-0.05,bot_xtick-0.1*length_y,0.0,"0.0");
      GLDrawStr( 0.455,bot_xtick-0.1*length_y,0.0,"0.5");
      GLDrawStr( 0.96,bot_xtick-0.1*length_y,0.0,"1.0");

      /* print title string                            */
      glColor3f(1.0, 1.0, 1.0);
      GLDrawStr(coord.titlepos,yaxis.max_value+0.15*length_y,0.0,coord.title);

      glPopMatrix();

      /* switch on lighing                             */
      glEnable(GL_LIGHTING);
      
      /* set line width to 1                           */
      glLineWidth(1);
    }
}

/******************************************************************
 @package   nightfall
 @author    Markus Kuster (kuster@astro.uni-tuebingen.de)
 @version   0.9
 @short     Set default parameters for coordinate system
 @param     (float)  ymin    minimum y value
 @param     (float)  ymax    maximum y value
 @param     (char *) xlabel  label of the x-axis
 @param     (char *) ylabel  label of the y-axis
 @param     (char *) title   title of the plot
 @return    (CoordType) structure containing information of the 
                        coordinate system 
 @heading   Open GL Animation
*******************************************************************/
CoordType SetupCoords(float ymin, float ymax, 
		      char *xlabel, char *ylabel, char *title) 
{
  AxisType x,y;
  CoordType coord;
  size_t str_length;

  x = coord.xaxis;
  y = coord.yaxis;
  
  strncpy(coord.title, title , MAX_STR_LEN-1);
  str_length = strlen(title);
  if (str_length <= 12) {
    coord.titlepos = 0.08;
  } else {
    coord.titlepos = 0.02;
  }    

  strncpy(x.label, xlabel, MAX_STR_LEN-1);
  strncpy(y.label, ylabel, MAX_STR_LEN-1);
  
  /* min/max values for x/y axis                                           */
  x.max_value = 1.0;
  x.min_value = -0.5;
  y.max_value = ymax;
  y.min_value = ymin;

  /* set size of border between plot and window  */
  /* to 0.1 axes length                          */
  x.length           = x.max_value - x.min_value;
  coord.xborder_size = 0.0166 * x.length;

  y.length           = y.max_value - y.min_value;
  coord.yborder_size = 0.03 * y.length;

  /* number of tickmarks                        */
  x.numticks = 4;
  y.numticks = 4;

  coord.xaxis = x; 
  coord.yaxis = y;
  return(coord);
}

#endif /* OpenGL end */