#include "GraphicsForGravity.h"
#include "gravity_support.h"
#include "GraphicsSupport.h"
#include <stdio.h>
#include <ctype.h>
#include <math.h>
/* plus define pi */
#define pi		3.14159265358979323846
#include <X11/Xlib.h>	/* Xgraf unfortunately no longer does this for you! */
/*
 *	There are local copies of the "Xgraf and allied" files
 *	here in this directory
 *	(modified to use "filename.h" syntax
 *	rather than <filename.h> syntax as needed)
 *	because it would be far too complicated to try to
 *	resurrect and re-install
 *	the various mixed Modula-2 / C {compiling, linking, library-generating}
 *	mechanisms of old.
 *	Hence the "filename.h" syntax rather than <filename.h> syntax below.
 */
#include "Xgraf.h"
#include "Xcolours.h"

#define bool	int
#define TRUE	1
#define FALSE	0

#define WindowEventMask		\
		(MaskKeypress | MaskAllMotion | MaskMouseDown | MaskResize)
/* note: MaskAllMotion is for mouse movements */
#define WindowLeftX		(-1)
#define WindowTopY		(-1)
#define InitWindowWidth		500
#define InitWindowHeight	500
#define WindowDistance		500	/* eye-to-window, in pixels */

#define SMALL_FACTOR		0.1
/* for re-initializing the user translation distance and boost speed */

#define COARSE_ADJ_FACTOR	2.0
/* for changing the user translation distance and boost speed */

#define RotationAngle		(pi/128.0)


static Window		Win;
static int		WindowWidth, WindowHeight;
#define WindowShortSideLength	\
			(WindowWidth<WindowHeight ? WindowWidth : WindowHeight)
static GRAPHICS_STATE	gstate;



/* private forward references */
static void SetUpSensibleGraphicsState (STATE *s, GRAPHICS_STATE *gs);





void InitializeGraphics (STATE *AState)
{
  /* get window ready */
  InitialiseGraf();
  AllowNoEvent();	/* for clean semantics of EventPending() & GetEvent() */
  SetEventMask(WindowEventMask);
  SetWBackingStore();	/* asks server to try to provide off-screen pixmap */

  WindowWidth  = InitWindowWidth;
  WindowHeight = InitWindowHeight;

  Win = MakeWindow(RootWin(), "gravity", WindowLeftX, WindowTopY,
                   WindowWidth, WindowHeight);
  /*
  now set this window to be inverse video.
  NOTE: must do this here, after the call to MakeWindow; any earlier and the
        current window is still RootWin()!
  */
  SetColour(White);
  SetBackColour(Black);

  /* now get graphics state ready. */
  SetUpSensibleGraphicsState(AState, &gstate);
}





void DrawState (STATE *AState)
{
  /* MouseX and MouseY are for squirreling away the mouse position. */
  static int MouseX = 0.5*InitWindowWidth;
  static int MouseY = 0.5*InitWindowHeight;

  /*
  UserDimension is for storing a user dimension number from 1 to NDIMS, and is
  initialized to 0 to show it hasn't yet been obtained from the keyboard.
  */
  static int UserDimension = 0;

  /* TrailsWanted stores whether or not the bodies are to leave trails. */
  static bool TrailsWanted = TRUE;

  Window     dummywin;
  EventType  event;
  EventData  eventdata;

  GRAPHICS_VECTOR apparent_user_velocity;	/* in user coordinates */

  GRAPHICS_VECTOR boost_velocity;		/* in actual coordinates */
  GRAPHICS_VECTOR apparent_boost_velocity;	/* in user coordinates */

  VECTOR P;		/* a body in user coordinates */
  double     Rx, Ry;	/* 2-D retinal coordinates */
  double     Rlength;	/* used for moving towards a given retinal direction */
  int        i,j,k,d;
  GRAPHICS_MATRIX temp_matrix;


  /*
  first update the graphics state to take into account the drift of
  user position since the last time we updated it.
  */
  MultiplyGraphicsVectorByMatrix
   (gstate.user_velocity, gstate.M, apparent_user_velocity);
  for (d=0; d<NDIMS+1; d++) {
    gstate.M[ NDIMS ][ d ] -=
     apparent_user_velocity[d] * (AState->Time - gstate.TimeLastUpdated);
  }/*endfor*/
  gstate.TimeLastUpdated = AState->Time;


  /*
  now see if any events need special attention. Note that EventExpose,
  perhaps surprisingly, does not need special attention: we're going to
  draw the latest state anyway.
  (albeit without clearing the window, if trails are wanted)
  */
  while (EventPending()) {
    GetEvent(&dummywin, &event, &eventdata);
    switch (event) {
      case EventMove :	/* mouse movement, not moving of the whole window! */

        /* squirrel away the mouse position for future use by "m" or "M" */
        MouseX = eventdata.x;
        MouseY = eventdata.y;
        break;

      case EventResize :

        WindowWidth  = eventdata.w;
        WindowHeight = eventdata.h;
        ClearWindow();
        break;

      case EventKeypress :

        #define FLIP_FACTOR	( islower(eventdata.ch) ? 1.0 : -1.0 )

        switch (eventdata.ch) {
          case '0' : case '1' : case '2' : case '3' : case '4' :
          case '5' : case '6' : case '7' : case '8' : case '9' :
            UserDimension = eventdata.ch == '0' ? 10 : eventdata.ch - '0';
            if (UserDimension > NDIMS) {
              beep();
              UserDimension = 0;
            }/*endif*/
            break;

          case 'b' : case 'B' :
            if (UserDimension > 0) {
              for (d=0; d<NDIMS+1; d++) {
                apparent_boost_velocity[d] =
                 d==UserDimension-1 ? FLIP_FACTOR * gstate.BoostSpeed : 0.0;
              }/*endfor*/
              InvertOrthonormal(gstate.M, temp_matrix);
              MultiplyGraphicsVectorByMatrix
               (apparent_boost_velocity, temp_matrix, boost_velocity);
              AddGraphicsVectors
               (gstate.user_velocity, boost_velocity, gstate.user_velocity);
              ClearWindow();
            }else{
              beep();
            }/*endif*/
            break;

          case 'c'            : ClearWindow();
                                break;

          case 'C'            : TrailsWanted = ! TrailsWanted;
                                break;

          case 'g' : case 'G' : /* go to a good vantage point */
                                SetUpSensibleGraphicsState(AState, &gstate);
                                ClearWindow();
                                break;

          case 'h' : case 'H' :
            printf("*****HELP MENU*****\n");
            printf("h,H	display this help menu\n");
            printf("q,Q	quit\n\n");
            printf("In the examples below any valid dimension-numbers can be used.\n");
            if (NDIMS == 10) {
              printf("[For dimension 10 press 0.]\n");
            }/*endif*/
            printf("Also, t,T,r,R,b,B will auto-repeat without repeating the dimension number.\n\n");
            printf("m	translate in direction indicated by mouse cursor\n");
            printf("M	translate in opposite direction to mouse cursor\n");
            printf("3t	translate along positive 3-axis\n");
            printf("3T	translate along negative 3-axis\n");
            printf("LeftMouseButton: rotate to face direction clicked on\n");
            printf("2r	rotate towards positive 2-axis\n");
            printf("2R	rotate towards negative 2-axis\n");
            printf("3b	boost along positive 3-axis\n");
            printf("3B	boost along negative 3-axis\n");
            printf("c  	clear display window\n");
            printf("C  	toggle between clearing vs. leaving trails\n");
            printf("x	reduce translation distance\n");
            printf("X	increase translation distance\n");
            printf("v	reduce boost speed\n");
            printf("V	increase boost speed\n");
            printf("g,G	go to a good vantage point\n\n\n");
            break;

          case 'm' : case 'M' :
            Rx =   (MouseX + 0.5 - 0.5*WindowWidth ) / WindowDistance;
            Ry = - (MouseY + 0.5 - 0.5*WindowHeight) / WindowDistance;
            Rlength = sqrt(Rx*Rx+Ry*Ry+1.0);

            /* now move the correct distance in the required direction */
            gstate.M[ NDIMS ][0] -=
             FLIP_FACTOR *  Rx / Rlength * gstate.TranslationDistance;
            gstate.M[ NDIMS ][1] -=
             FLIP_FACTOR *  Ry / Rlength * gstate.TranslationDistance;
            gstate.M[ NDIMS ][ NDIMS-1 ] -=
             FLIP_FACTOR * 1.0 / Rlength * gstate.TranslationDistance;

            ClearWindow();
            break;

          case 'q' : case 'Q' : CloseDownGraf();
                                exit(0);
                                break;

          case 'r' : case 'R' :
            if (UserDimension > 0 && UserDimension < NDIMS) {
              CreateRotationMatrix
               (NDIMS-1, UserDimension-1, - FLIP_FACTOR * RotationAngle,
                temp_matrix);
              MultiplyMatrices(gstate.M, temp_matrix, gstate.M);
              ClearWindow();
            }else{
              beep();
            }/*endif*/
            break;

          case 't' : case 'T' :
            if (UserDimension > 0) {
              CreateTranslationMatrix
               (UserDimension-1, - FLIP_FACTOR * gstate.TranslationDistance,
                temp_matrix);
              MultiplyMatrices(gstate.M, temp_matrix, gstate.M);
              ClearWindow();
            }else{
              beep();
            }/*endif*/
            break;

          case 'x'            : gstate.TranslationDistance /= COARSE_ADJ_FACTOR;
                                break;

          case 'X'            : gstate.TranslationDistance *= COARSE_ADJ_FACTOR;
                                break;

          case 'v'            : gstate.BoostSpeed /= COARSE_ADJ_FACTOR;
                                break;

          case 'V'            : gstate.BoostSpeed *= COARSE_ADJ_FACTOR;
                                break;

          default             : break;
        }/*endswitch*/
        break;

      case EventMouseDown :

        switch (eventdata.button) {
          case 1 :
            Rx =   (eventdata.x + 0.5 - 0.5*WindowWidth ) / WindowDistance;
            Ry = - (eventdata.y + 0.5 - 0.5*WindowHeight) / WindowDistance;
            if (Rx != 0.0 || Ry != 0.0) {
              /* first twist so retinal point is on positive 1-axis */
              CreateRotationMatrix
               (1, 0, atan2(Ry,Rx), temp_matrix);
              MultiplyMatrices(gstate.M, temp_matrix, gstate.M);

              /* now rotate point to exact centre of view */
              CreateRotationMatrix
               (0, NDIMS-1, atan2(sqrt(Rx*Rx+Ry*Ry), 1.0), temp_matrix);
              MultiplyMatrices(gstate.M, temp_matrix, gstate.M);

              /* now undo twist above. */
              CreateRotationMatrix
               (0, 1, atan2(Ry,Rx), temp_matrix);
              MultiplyMatrices(gstate.M, temp_matrix, gstate.M);

              ClearWindow();
            }/*endif*/
            break;

          default :
            break;
        }/*endswitch*/
        break;

      default :

        break;

    }/*endswitch*/
  }/*endwhile*/


  if (! TrailsWanted) {
    ClearWindow();
  }/*endif*/


  /* now, just draw all the bodies! */
  for (i=0; i<AState->n; i++) {
    /*
    turn body i into user coordinates. We use a tight little loop here as
    it's actually clearer to do this than to turn everything into graphics
    types and then use the graphics matrix/vector functions.
    (not to mention faster, in this "fairly inner" position.....)
    */
    CreateZeroVector(P);
    for (j=0; j<NDIMS; j++) {
      for (k=0; k<NDIMS+1; k++) {
        P[j] +=
         (k<NDIMS ? AState->bodydata[i].position[k] : 1.0) * gstate.M[k][j];
      }/*endfor*/
    }/*endfor*/

    /* check if point is in front of user */
    if (P[NDIMS-1] > 0) {
      /* find retinal coordinates */
      Rx = P[0] / P[NDIMS-1];
      Ry = P[1] / P[NDIMS-1];

      /* now plot in Xgraf coordinates */
      DrawPoint( (int) (0.5*WindowWidth  + Rx*WindowDistance),
                 (int) (0.5*WindowHeight - Ry*WindowDistance)
               );
    }/*endif*/
  }/*endfor*/
}





static void SetUpSensibleGraphicsState
             (STATE *AState, GRAPHICS_STATE *AGraphicsState)
/*
Treats its first arg as read-only.
Sets up sensible values in the indicated graphics state, with the user
co-moving with the system of bodies at a well-chosen vantage point.
*/
{
  double total_mass;
  VECTOR total_moment;
  VECTOR total_momentum;
  VECTOR centre_of_mass_position;
  VECTOR centre_of_mass_velocity;
  VECTOR position_from_CM, velocity_from_CM;
  double max_CM_disp, max_CM_speed;
  int    i, d;

  /*
  find centre-of-mass position and velocity, and maximum deviations therefrom,
  in order to set up user to be co-moving with the system
  at a well-chosen vantage point.
  */
  total_mass   = 0.0;
  CreateZeroVector(total_moment);
  CreateZeroVector(total_momentum);
  max_CM_disp  = 0.0;
  max_CM_speed = 0.0;

  for (i=0; i<AState->n; i++) {
    total_mass += AState->bodydata[i].mass;
    for (d=0; d<NDIMS; d++) {
      total_moment[d]   += AState->bodydata[i].mass * AState->bodydata[i].position[d];
      total_momentum[d] += AState->bodydata[i].mass * AState->bodydata[i].velocity[d];
    }/*endfor*/;
  }/*endfor*/

  for (d=0; d<NDIMS; d++) {
    centre_of_mass_position[d] = total_moment[d] / total_mass;
    centre_of_mass_velocity[d] = total_momentum[d] / total_mass;
  }/*endfor*/;

  for (i=0; i<AState->n; i++) {
    for (d=0; d<NDIMS; d++) {
      position_from_CM[d] =
       AState->bodydata[i].position[d] - centre_of_mass_position[d];
      velocity_from_CM[d] =
       AState->bodydata[i].velocity[d] - centre_of_mass_velocity[d];
    }/*endfor*/;

    if (sqrt(lengthsquared(position_from_CM)) > max_CM_disp) {
      max_CM_disp = sqrt(lengthsquared(position_from_CM));
    }/*endif*/

    if (sqrt(lengthsquared(velocity_from_CM)) > max_CM_speed) {
      max_CM_speed = sqrt(lengthsquared(velocity_from_CM));
    }/*endif*/
  }/*endfor*/;


  /* now use this hard-won data to set a sensible initial user position etc. */

  /* first put user exactly *on* the centre-of-mass..... */
  CreateIdentityMatrix(AGraphicsState->M);
  for(d=0; d<NDIMS+1; d++) {
    AGraphicsState->M[ NDIMS ][ d ] =
     d<NDIMS ? (- centre_of_mass_position[d]) : 1.0;
  }/*endfor*/

  /* .....and now back away a well-chosen length. */
  AGraphicsState->M[ NDIMS ][ NDIMS-1 ] +=
   sqrt(16.0 * pow(WindowDistance, 2.0) / pow(WindowShortSideLength, 2.0) + 4.0)
    * max_CM_disp;

  /* start user off co-moving with the system. */
  for (d=0; d<NDIMS+1; d++) {
    AGraphicsState->user_velocity[d] =
     d<NDIMS ? centre_of_mass_velocity[d] : 0.0;
  }/*endfor*/

  /* set up sensible user translation distance and boost speed */
  AGraphicsState->TranslationDistance = SMALL_FACTOR * max_CM_disp;
  AGraphicsState->BoostSpeed          = SMALL_FACTOR * max_CM_speed;

  /*
  finally, record the time, to enable future updates of the graphics state.
  */
  AGraphicsState->TimeLastUpdated = AState->Time;
}