/*   []----------------------------------------[]
      |               OpenGl.inc               |
     []----------------------------------------[]
      |                                        |
      |   AUTHOR:    MFSomers 2005.            |
      |   USE:       OpenGl routines for main  |
      |               code ready to be         |
      |               included if needed ...   |
      |                                        |
     []----------------------------------------[]
*/

// Copyright (C) 2005 M.F. Somers, Theoretical Chemistry Department, Leiden University
//
// This 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.
//
// http://www.gnu.org/licenses/gpl.txt

#if defined(__linux__)
#   include <GL/gl.h>
#   include <GL/glu.h>
#elif defined(_WIN32) || defined(WIN32)
#   include <gl\gl.h>
#   include <gl\glu.h>
#elif defined(__APPLE__)
#   include <OpenGL/gl.h>
#   include <OpenGL/glu.h>
#elif defined(FREEBSD)
#   include <GL/gl.h>
#   include <GL/glu.h>
#else
#   error No OpenGL support for this platform ...
#endif

/* -------------------------------------------------------------------------- */

#define CPU_SECONDS() ((float)(clock())/(float)(CLOCKS_PER_SEC))

#define OPENGL_MOUSE_TRIGGER              10
#define OPENGL_SPHERE_DETAIL              10

#define OPENGL_DELTA_ANGLE                1.0
#define OPENGL_DELTA_ZOOM                 0.01
#define OPENGL_MAX_ZOOM                   100.0
#define OPENGL_MIN_ZOOM                   0.01

int    OpenGL_DynamicsRunning = 1;
int    OpenGL_RefreshScreen = 1; 
double OpenGL_ModelOrientationPhi = 330.0;
double OpenGL_ModelOrientationTheta = 30.0;
double OpenGL_ModelOrientationZoomF = 1.0;
double OpenGL_LastRefreshTime = CPU_SECONDS();
double OpenGL_SinTable1[ OPENGL_SPHERE_DETAIL + 2 ];
double OpenGL_CosTable1[ OPENGL_SPHERE_DETAIL + 2 ];
double OpenGL_SinTable2[ OPENGL_SPHERE_DETAIL ];
double OpenGL_CosTable2[ OPENGL_SPHERE_DETAIL ];

/* -------------------------------------------------- 3D graphics routines - */
/* Used to return the color vector for a given ID                            */

Vector GetColorFromTable( ColorListPointer& Colors, int Id )
{
 Vector White = { 1.0, 1.0, 1.0 };
 ColorPointer TempColorPointer;
  
 // return white color by default for empty lists or non-found/defined colors
 if( Colors == NULL ) return( White );
 
 // check for fast if we request the last color...
 if( (*(Colors))[ 0 ] != NULL )
  if( (*(Colors))[ 0 ] -> Id == Id )
   return( (*(Colors))[ 0 ] -> RGB );
 
 // if not we need to lookup, if we found it, move it to the beginning of the list for speed...
 for( int n = 1; n < Colors -> size(); ++n )
  if( (*(Colors))[ n ] != NULL )
   if( (*(Colors))[ n ] -> Id == Id )
   {
    TempColorPointer = (*(Colors))[ n ];
    (*(Colors))[ n ] = (*(Colors))[ 0 ];
    (*(Colors))[ 0 ] = TempColorPointer;
    return( TempColorPointer -> RGB );
   }
   
 return( White );    
}

/* ----------------------------------------------------------------------- */
/* Called to setup the 3D OpenGL rendering scene                           */

void SetupOpenGLScene( int Width, int Height, const Vector& TheBox )
{
 GLfloat specular[]  = { 1.0, 1.0, 1.0, 1.0 };
 GLfloat shininess[] = { 75.0 };
 GLfloat position[]  = { 1.75 * TheBox.X, 1.75 * TheBox.Y, 1.75 * TheBox.Z, 0.0 };
 double Ratio;
 
 // Check sizes...
 if( Width < 10 ) Width = 10;
 if( Height < 10 ) Height = 10;
 Ratio = ( double )( Width ) / ( double )( Height );
 
 // Setup vars
 OpenGL_LastRefreshTime = CPU_SECONDS();
 
 // setup the OpenGL default global state
 glClearColor( 0.0, 0.0, 0.0, 0.0 );
 glEnable( GL_DEPTH_TEST );
 glShadeModel( GL_SMOOTH );
 glMaterialfv( GL_FRONT_AND_BACK, GL_SPECULAR, specular );
 glMaterialfv( GL_FRONT_AND_BACK, GL_SHININESS, shininess );
 glColorMaterial( GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE );
 glEnable( GL_COLOR_MATERIAL );
 glLightfv( GL_LIGHT0, GL_POSITION, position );
 glEnable( GL_LIGHTING );
 glEnable( GL_LIGHT0 );
  
 // setup the OpenGL default projection state
 glMatrixMode( GL_PROJECTION );
 glLoadIdentity();
 glPushMatrix();
 glViewport( 0, 0, Width, Height ); 
 gluPerspective( 90, Ratio, 0.01, 1.75 * Abs( TheBox ) );
 gluLookAt( 0.75 * Abs( TheBox ), 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 );
 
 // and setup default model state
 glMatrixMode( GL_MODELVIEW );
 glLoadIdentity();
 glPushMatrix();
 glRotatef( OpenGL_ModelOrientationTheta, 0.0, 1.0, 0.0 );
 glRotatef( OpenGL_ModelOrientationPhi, 0.0, 0.0, 1.0 ); 
 glScalef( OpenGL_ModelOrientationZoomF, OpenGL_ModelOrientationZoomF, OpenGL_ModelOrientationZoomF );
 
 // now setup sin and cos tables for rendering/drawing the spheres
 OpenGL_SinTable1[ 0 ] = 0.0;
 OpenGL_CosTable1[ 0 ] = 1.0;
 for( int i = 1 ; i < OPENGL_SPHERE_DETAIL+1; ++i )
 {
  Ratio = ( double )( i ) * ( 2.0 * CONST_PI / (OPENGL_SPHERE_DETAIL+1) );
  OpenGL_SinTable1[ i ] = sin( Ratio );
  OpenGL_CosTable1[ i ] = cos( Ratio );
 }
 OpenGL_SinTable1[ OPENGL_SPHERE_DETAIL+1 ] = 0.0;
 OpenGL_CosTable1[ OPENGL_SPHERE_DETAIL+1 ] = 1.0;
 
 OpenGL_SinTable2[ 0 ] = 0.0;
 OpenGL_CosTable2[ 0 ] = 1.0;
 for( int i = 1 ; i < OPENGL_SPHERE_DETAIL; ++i )
 {
  Ratio = ( double )( i ) * ( CONST_PI / OPENGL_SPHERE_DETAIL );
  OpenGL_SinTable2[ i ] = sin( Ratio );
  OpenGL_CosTable2[ i ] = cos( Ratio );
 }
 
} 

/* ----------------------------------------------------------------------- */
/* Resizes the 3D OpenGL rendering scene                                   */

void ResizeOpenGLScene( int Width, int Height, const Vector& TheBox )
{
 double Ratio;
 
 if( Height == 0 ) Height = 1;

 Ratio = ( double )( Width ) / ( double )( Height );

 // build new projection matrix
 glMatrixMode( GL_PROJECTION );
 glPopMatrix();
 glPushMatrix();
        
 // set the viewport to be the entire window
 glViewport( 0, 0, Width, Height );

 // set the correct perspective
 gluPerspective( 90, Ratio, 0.01, 1.75 * Abs( TheBox ) );
 gluLookAt( 0.75 * Abs( TheBox ), 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 );
 
 glMatrixMode( GL_MODELVIEW );
 
 OpenGL_RefreshScreen = 1;
}

/* ----------------------------------------------------------------------- */
/* Changes the way you look at the rendering scene                         */

void ChangeViewPointOfOpenGLScene( double Phi, double Theta, double ZoomF )
{
 Phi = Phi >= 360 ? Phi - 360 : Phi;
 Phi = Phi <= 0 ? 360 - Phi : Phi;
 OpenGL_ModelOrientationPhi = Phi;
  
 Theta = Theta >= 360 ? Theta - 360 : Theta;
 Theta = Theta <= 0 ? 360 - Theta : Theta;
 OpenGL_ModelOrientationTheta = Theta;
 
 ZoomF = ZoomF >= OPENGL_MAX_ZOOM ? OPENGL_MAX_ZOOM : ZoomF;
 ZoomF = ZoomF <= OPENGL_MIN_ZOOM ? OPENGL_MIN_ZOOM : ZoomF;
 OpenGL_ModelOrientationZoomF = ZoomF;
 
 glMatrixMode( GL_MODELVIEW );
 glPopMatrix();
 glPushMatrix();
 glRotatef( Theta, 0.0, 1.0, 0.0 );
 glRotatef( Phi, 0.0, 0.0, 1.0 ); 
 glScalef( ZoomF, ZoomF, ZoomF );
 
 OpenGL_RefreshScreen = 1;
}

/* ----------------------------------------------------------------------- */
/* Outputs a solid sphere with given parameters                            */
/* This code was taken from freeglut and speed-up for the sake of it       */

void OpenGLSolidSphere( double Radius )
{
 double z0, z1;
 double r0, r1;
 int    i,  j;
 
  // The top stack is covered with a triangle fan
  z0 = 1.0; z1 = OpenGL_CosTable2[ 1 ];
  r0 = 0.0; r1 = OpenGL_SinTable2[ 1 ];

  glBegin( GL_TRIANGLE_FAN );
   glNormal3d( 0.0, 0.0, 1.0 );
   glVertex3d( 0.0, 0.0, Radius );
   for( j = 0; j <= OPENGL_SPHERE_DETAIL+1; ++j )
   {
    glNormal3d( OpenGL_CosTable1[ j ] * r1, OpenGL_SinTable1[ j ] * r1, z1 );
    glVertex3d( OpenGL_CosTable1[ j ] * r1 * Radius, OpenGL_SinTable1[ j ] * r1 * Radius, z1 * Radius );
   }
  glEnd();
    
  // Cover each stack with a quad strip, except the top and bottom stacks
  for( i = 1 ; i < OPENGL_SPHERE_DETAIL - 1; ++i )
  {
   z0 = z1; z1 = OpenGL_CosTable2[ i + 1 ];
   r0 = r1; r1 = OpenGL_SinTable2[ i + 1 ];
   
   glBegin( GL_QUAD_STRIP );
   for( j = 0; j <= OPENGL_SPHERE_DETAIL+1; ++j )
   {
    glNormal3d( OpenGL_CosTable1[ j ] * r1, OpenGL_SinTable1[ j ] * r1, z1 );
    glVertex3d( OpenGL_CosTable1[ j ] * r1 * Radius, OpenGL_SinTable1[ j ] * r1 * Radius, z1 * Radius );
    glNormal3d( OpenGL_CosTable1[ j ] * r0, OpenGL_SinTable1[ j ] * r0, z0 );
    glVertex3d( OpenGL_CosTable1[ j ] * r0 * Radius, OpenGL_SinTable1[ j ] * r0 * Radius, z0 * Radius );
   }
   glEnd();
  }

  // The bottom stack is covered with a triangle fan
  z0 = z1;
  r0 = r1;
  
  glBegin( GL_TRIANGLE_FAN );
   glNormal3d( 0.0, 0.0, -1.0 );
   glVertex3d( 0.0, 0.0, -Radius );
   for( j = 0; j <= OPENGL_SPHERE_DETAIL+1; ++j )
   {
    glNormal3d( OpenGL_CosTable1[ j ] * r0, OpenGL_SinTable1[ j ] * r0, z0 );
    glVertex3d( OpenGL_CosTable1[ j ] * r0 * Radius, OpenGL_SinTable1[ j ] * r0 * Radius, z0 * Radius );
   }
  glEnd();
    
}

/* ----------------------------------------------------------------------- */
/* Called whenever the 3D scene needs to be rendered                       */

void RenderOpenGLScene( ParticleListPointer TheState, ColorListPointer& TheColors, const Vector& TheBox, int Periodicity )
{
 ParticlePointer TheParticle;
 Vector          ThePosition;
 Vector          HalfBox = 0.5 * TheBox;
 Vector          TheColor;
 double          TheRadius;
 int             n;
 
 // store model matrix and clear window
 glPushMatrix();
 glClear( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT );

 // draw a red box, got the code from freeglut
 glColor3f( 1.0, 0.0, 0.0 );
 glBegin( GL_LINE_LOOP );
      glNormal3f( 1.0, 0.0, 0.0 );
      glVertex3f(+HalfBox.X,-HalfBox.Y,+HalfBox.Z);
      glVertex3f(+HalfBox.X,-HalfBox.Y,-HalfBox.Z);
      glVertex3f(+HalfBox.X,+HalfBox.Y,-HalfBox.Z);
      glVertex3f(+HalfBox.X,+HalfBox.Y,+HalfBox.Z);
 glEnd();
 glBegin( GL_LINE_LOOP );
      glNormal3f( 0.0, 1.0, 0.0 );
      glVertex3f(+HalfBox.X,+HalfBox.Y,+HalfBox.Z);
      glVertex3f(+HalfBox.X,+HalfBox.Y,-HalfBox.Z);
      glVertex3f(-HalfBox.X,+HalfBox.Y,-HalfBox.Z);
      glVertex3f(-HalfBox.X,+HalfBox.Y,+HalfBox.Z);
 glEnd();
 glBegin( GL_LINE_LOOP );
      glNormal3f( 0.0, 0.0, 1.0 );
      glVertex3f(+HalfBox.X,+HalfBox.Y,+HalfBox.Z);
      glVertex3f(-HalfBox.X,+HalfBox.Y,+HalfBox.Z);
      glVertex3f(-HalfBox.X,-HalfBox.Y,+HalfBox.Z);
      glVertex3f(+HalfBox.X,-HalfBox.Y,+HalfBox.Z);
 glEnd();
 glBegin( GL_LINE_LOOP );
      glNormal3f( -1.0, 0.0, 0.0 );
      glVertex3f(-HalfBox.X,-HalfBox.Y,+HalfBox.Z);
      glVertex3f(-HalfBox.X,+HalfBox.Y,+HalfBox.Z);
      glVertex3f(-HalfBox.X,+HalfBox.Y,-HalfBox.Z);
      glVertex3f(-HalfBox.X,-HalfBox.Y,-HalfBox.Z);
 glEnd();
 glBegin( GL_LINE_LOOP );
      glNormal3f( 0.0, -1.0, 0.0 );
      glVertex3f(-HalfBox.X,-HalfBox.Y,+HalfBox.Z);
      glVertex3f(-HalfBox.X,-HalfBox.Y,-HalfBox.Z);
      glVertex3f(+HalfBox.X,-HalfBox.Y,-HalfBox.Z);
      glVertex3f(+HalfBox.X,-HalfBox.Y,+HalfBox.Z);
 glEnd();
 glBegin( GL_LINE_LOOP );
      glNormal3f( 0.0, 0.0, -1.0 );
      glVertex3f(-HalfBox.X,-HalfBox.Y,-HalfBox.Z);
      glVertex3f(-HalfBox.X,+HalfBox.Y,-HalfBox.Z);
      glVertex3f(+HalfBox.X,+HalfBox.Y,-HalfBox.Z);
      glVertex3f(+HalfBox.X,-HalfBox.Y,-HalfBox.Z);
 glEnd();

 // now draw the X, Y and Z axis in green
 glColor3f( 0.0, 1.0, 0.0 );
 glBegin( GL_LINES );
   glNormal3f( 0.25 * HalfBox.X, 0.25 * HalfBox.Y, 0.25 * HalfBox.Y );
   glVertex3f( 0.0, 0.0, 0.0 );
   glVertex3f( 0.25 * HalfBox.X, 0.0, 0.0 );
   glVertex3f( 0.0, 0.0, 0.0 );
   glVertex3f( 0.0, 0.25 * HalfBox.Y, 0.0 );
   glVertex3f( 0.0, 0.0, 0.0 );
   glVertex3f( 0.0, 0.0, 0.25 * HalfBox.Z );
 glEnd();  

 // now start drawing the scene using spheres...
 for( n = 0; n < TheState -> size(); ++n )
  if( (*(TheState))[ n ] != NULL )
  {
   TheParticle = (*(TheState))[ n ];
   TheRadius   = TheParticle -> Radius;
   ThePosition = TheParticle -> Position;
   
   TheColor = GetColorFromTable( TheColors, TheParticle -> ColorIndex );       
    
   glPushMatrix();
   glColor3f( TheColor.X, TheColor.Y, TheColor.Z );
         
   if( Periodicity )                  // if we use periodicity, make sure the position is within the cell
   {
    ThePosition = MimimumSystemImageConventionVector( ThePosition, TheBox );
    glTranslatef( ThePosition.X, ThePosition.Y, ThePosition.Z );
    OpenGLSolidSphere( TheRadius );
    
    if( ThePosition.X >= HalfBox.X - 1.25 * TheRadius )
    {
     glPopMatrix();
     glPushMatrix();
     glTranslatef( ThePosition.X - TheBox.X, ThePosition.Y, ThePosition.Z );
     OpenGLSolidSphere( TheRadius );
    }
    if( ThePosition.X <= 1.25 * TheRadius - HalfBox.X )
    {
     glPopMatrix();
     glPushMatrix();    
     glTranslatef( ThePosition.X + TheBox.X, ThePosition.Y, ThePosition.Z );
     OpenGLSolidSphere( TheRadius );
    }    
    if( ThePosition.Y >= HalfBox.Y - 1.25 * TheRadius )
    {
     glPopMatrix();
     glPushMatrix();    
     glTranslatef( ThePosition.X, ThePosition.Y - TheBox.Y, ThePosition.Z );
     OpenGLSolidSphere( TheRadius );
    }
    if( ThePosition.Y <= 1.25 * TheRadius - HalfBox.Y )
    {
     glPopMatrix();
     glPushMatrix();    
     glTranslatef( ThePosition.X, ThePosition.Y + TheBox.Y, ThePosition.Z );
     OpenGLSolidSphere( TheRadius );
    }       
    if( ThePosition.Z >= HalfBox.Z - 1.25 * TheParticle -> Radius )
    {
     glPopMatrix();
     glPushMatrix();    
     glTranslatef( ThePosition.X, ThePosition.Y, ThePosition.Z - TheBox.Z );
     OpenGLSolidSphere( TheRadius );
    }
    if( ThePosition.Z <= 1.25 * TheRadius - HalfBox.Z )
    {
     glPopMatrix();
     glPushMatrix();    
     glTranslatef( ThePosition.X, ThePosition.Y, ThePosition.Z + TheBox.Z );
     OpenGLSolidSphere( TheRadius );
    }      
        
   }
   else
   {
    glTranslatef( ThePosition.X, ThePosition.Y, ThePosition.Z );
    OpenGLSolidSphere( TheRadius );
   }
    
   glPopMatrix();
  }
	 
 // finally restore model matrix
 glPopMatrix();
 
 OpenGL_RefreshScreen = 0;
 OpenGL_LastRefreshTime = CPU_SECONDS(); 
}

/* ----------------------------------------------------------------------- */