/* konata0.c - Don Yang (uguu.org)

   Just like Kanata for the most part, except;
   - Pinkish color scheme instead of blue.
   - Slightly different animation parameters.

   07/06/08
                                                                           */
/*@ -realcompare @*/

#include<assert.h>
#include<sys/time.h>
#include<GL/glut.h>
#include<math.h>
#include<stdio.h>
#include<stdlib.h>

#define PI 3.14159265358979323846264338327950288419716939937510
#define WINDOW_TITLE "Konata"

#define CONTROL_POINT_OFFSET  0.4

#define MIN_TRAIL_STEP_TIME   0.06
#define MAX_TRAIL_LENGTH_TIME 30.0
#define MAX_TRAIL_POINTS      512
#define MIN_ALPHA             0.0
#define MAX_ALPHA             0.8
#define ALPHA_STEP            0.04

#define WORLD_SIZE            128.0
#define TRAIL_WIDTH           16.0
#define RANDOM_TIME() \
   (((double)(rand() & 0x7fff) / (double)0x2000) + 0.5)
#define RANDOM_POSITION() \
   (((double)((rand() & 0x7fff) - 0x4000) / (double)0x4000) * WORLD_SIZE)
#define RANDOM_ANGLE() \
   (((double)((rand() & 0x7fff) - 0x4000) / (double)0x4000) * PI)

/* Global clock, updated exactly once per frame */
static double CurrentTime;

/* Stats counters */
static double StartTime;
static int FrameCount;

/* Circular buffer of vertices */
typedef struct
{
   double t;            /* Timestamp of vertex */
   double a;            /* Roll angle */
   double x, y, z;      /* Location */
   double dx, dy, dz;   /* Forward direction */
   double sx, sy, sz;   /* Side direction */
} TrailPoint;
static TrailPoint Trail[MAX_TRAIL_POINTS];
static int HeadIndex, TailIndex;

/* Transformed vertices */
static double VertexBuffer[MAX_TRAIL_POINTS * 2][3];
static double NormalBuffer[MAX_TRAIL_POINTS * 2][3];

/* Keyframes */
typedef struct
{
   double t0, t1;           /* Keyframe times */
   double x0, y0, z0, a0;   /* Control points */
   double x1, y1, z1, a1;
   double x2, y2, z2, a2;
   double x3, y3, z3, a3;
} KeyframePoint;
static KeyframePoint Keyframe;

static double CubicPoint(double a, double b, double c, double d, double t);
static double CubicDerivative(double a, double b, double c, double d, double t);
static void CrossProduct(double ax, double ay, double az,
                         double bx, double by, double bz,
                         double *ox, double *oy, double *oz);
static void RotateVector(double ax, double ay, double az, double angle,
                         double *px, double *py, double *pz);
static void Normalize(double *x, double *y, double *z);

static /*@dependent@*/TrailPoint *AddPoint(void);
static void SetPosition(/*@dependent@*/TrailPoint *point);
static void AddKeyframe(void);
static void UpdateObjects(void);

static void Init(void);
static void Render(void);
static void Animate(void);
static void Reshape(/*@unused@*/int w, /*@unused@*/int h);
static void Quit(/*@unused@*/unsigned char c,
                 /*@unused@*/int u, /*@unused@*/int v);
static void UpdateClock(void);

int main(int argc, char **argv)
{
   glutInit(&argc, argv);
   Init();
   glutMainLoop();
   return 0;
}

/* Interpolate point on cubic curve */
static double CubicPoint(double a, double b, double c, double d, double t)
{
   double ab = a + t * (b - a);
   double bc = b + t * (c - b);
   double cd = c + t * (d - c);
   double abc = ab + t * (bc - ab);
   double bcd = bc + t * (cd - bc);
   return abc + t * (bcd - abc);
}

/* First derivative of cubic curve */
static double CubicDerivative(double a, double b, double c, double d, double t)
{
   double t2 = t * t;
   double nt = 1.0 - t;
   double nt2 = nt * nt;
   double t2nt = 2.0 * t * nt;
   return 3.0 * (-nt2 * a + (nt2 - t2nt) * b + (t2nt - t2) * c + t2 * d);
}

/* Cross product of two vectors */
static void CrossProduct(double ax, double ay, double az,
                         double bx, double by, double bz,
                         double *ox, double *oy, double *oz)
{
   *ox = ay * bz - by * az;
   *oy = bx * az - ax * bz;
   *oz = ax * by - bx * ay;
}

/* Rotate vector about an axis */
static void RotateVector(double ax, double ay, double az, double angle,
                         double *px, double *py, double *pz)
{
   double sa = sin(angle);
   double ca = cos(angle);
   double nca = 1.0 - ca;
   double xy = ax * ay;
   double xz = ax * az;
   double yz = ay * az;
   double x0 = *px;
   double y0 = *py;
   double z0 = *pz;
   *px = (nca * ax * ax + ca) * x0
       + (nca * xy - az * sa) * y0
       + (nca * xz + ay * sa) * z0;
   *py = (nca * xy + az * sa) * x0
       + (nca * ay * ay + ca) * y0
       + (nca * yz - ax * sa) * z0;
   *pz = (nca * xz - ay * sa) * x0
       + (nca * yz + ax * sa) * y0
       + (nca * az * az + ca) * z0;
}

/* Normalize a vector */
static void Normalize(double *x, double *y, double *z)
{
   double d = *x * *x + *y * *y + *z * *z;
   if( d > 0.000001 )
   {
      d = sqrt(d);
      *x /= d;
      *y /= d;
      *z /= d;
   }
}

/* Add/expire points from trail, return next point to update */
static /*@dependent@*/TrailPoint *AddPoint(void)
{
   /* Expire oldest point */
   for(; TailIndex != HeadIndex &&
         CurrentTime - Trail[TailIndex].t >= MAX_TRAIL_LENGTH_TIME;
       TailIndex = (TailIndex + 1) % MAX_TRAIL_POINTS);

   /* Add new point, or replace latest point */
   if( CurrentTime - Trail[HeadIndex].t >= MIN_TRAIL_STEP_TIME )
   {
      HeadIndex = (HeadIndex + 1) % MAX_TRAIL_POINTS;
      Trail[HeadIndex].t = CurrentTime;
   }
   return &Trail[HeadIndex];
}

/* Set position of next trail point */
static void SetPosition(/*@dependent@*/TrailPoint *point)
{
   double dt, d;

   if( CurrentTime >= Keyframe.t1 )
      AddKeyframe();
   dt = (CurrentTime - Keyframe.t0) / (Keyframe.t1 - Keyframe.t0);

   /* Position */
   point->x = CubicPoint(
      Keyframe.x0, Keyframe.x1, Keyframe.x2, Keyframe.x3, dt);
   point->y = CubicPoint(
      Keyframe.y0, Keyframe.y1, Keyframe.y2, Keyframe.y3, dt);
   point->z = CubicPoint(
      Keyframe.z0, Keyframe.z1, Keyframe.z2, Keyframe.z3, dt);

   /* Direction */
   point->dx = CubicDerivative(
      Keyframe.x0, Keyframe.x1, Keyframe.x2, Keyframe.x3, dt);
   point->dy = CubicDerivative(
      Keyframe.y0, Keyframe.y1, Keyframe.y2, Keyframe.y3, dt);
   point->dz = CubicDerivative(
      Keyframe.z0, Keyframe.z1, Keyframe.z2, Keyframe.z3, dt);
   Normalize(&(point->dx), &(point->dy), &(point->dz));

   /* Find vector perpendicular to direction vector and positive Z axis.
      This is the "side direction" vector.                               */
   CrossProduct(point->dx, point->dy, point->dz, 0.0, 0.0, 1.0,
                &(point->sx), &(point->sy), &(point->sz));
   d = point->sx * point->sx + point->sy * point->sy + point->sz * point->sz;
   if( d < 0.000001 )
   {
      point->sx = 1.0;
      point->sy = point->sz = 0.0;
   }
   else
   {
      d = sqrt(d);
      point->sx /= d;
      point->sy /= d;
      point->sz /= d;
   }

   /* Rotate side direction vector about direction vector */
   point->a = CubicPoint(
      Keyframe.a0, Keyframe.a1, Keyframe.a2, Keyframe.a3, dt);
   RotateVector(point->dx, point->dy, point->dz, point->a,
                &(point->sx), &(point->sy), &(point->sz));
}

/* Add keyframe point */
static void AddKeyframe(void)
{
   Keyframe.t0 = Keyframe.t1;
   Keyframe.t1 = CurrentTime + RANDOM_TIME();
   Keyframe.x0 = Keyframe.x3;
   Keyframe.y0 = Keyframe.y3;
   Keyframe.z0 = Keyframe.z3;
   Keyframe.a0 = Keyframe.a3;
   Keyframe.x1 = Keyframe.x3 + (Keyframe.x3 - Keyframe.x2);
   Keyframe.y1 = Keyframe.y3 + (Keyframe.y3 - Keyframe.y2);
   Keyframe.z1 = Keyframe.z3 + (Keyframe.z3 - Keyframe.z2);
   Keyframe.a1 = Keyframe.a3 + (Keyframe.a3 - Keyframe.a2);
   Keyframe.x3 = RANDOM_POSITION();
   Keyframe.y3 = RANDOM_POSITION();
   Keyframe.z3 = RANDOM_POSITION();
   Keyframe.x2 = Keyframe.x3 + CONTROL_POINT_OFFSET * RANDOM_POSITION();
   Keyframe.y2 = Keyframe.y3 + CONTROL_POINT_OFFSET * RANDOM_POSITION();
   Keyframe.z2 = Keyframe.z3 + CONTROL_POINT_OFFSET * RANDOM_POSITION();
   Keyframe.a3 = RANDOM_ANGLE();
   Keyframe.a2 = Keyframe.a3 + CONTROL_POINT_OFFSET * RANDOM_ANGLE();
}

/* Update animation */
static void UpdateObjects(void)
{
   int a, b;

   /* Compute point */
   SetPosition(AddPoint());

   /* Compute vertices */
   a = HeadIndex * 2;
   b = a + 1;
   VertexBuffer[a][0] = Trail[HeadIndex].x - TRAIL_WIDTH * Trail[HeadIndex].sx;
   VertexBuffer[a][1] = Trail[HeadIndex].y - TRAIL_WIDTH * Trail[HeadIndex].sy;
   VertexBuffer[a][2] = Trail[HeadIndex].z - TRAIL_WIDTH * Trail[HeadIndex].sz;
   VertexBuffer[b][0] = Trail[HeadIndex].x + TRAIL_WIDTH * Trail[HeadIndex].sx;
   VertexBuffer[b][1] = Trail[HeadIndex].y + TRAIL_WIDTH * Trail[HeadIndex].sy;
   VertexBuffer[b][2] = Trail[HeadIndex].z + TRAIL_WIDTH * Trail[HeadIndex].sz;

   /* Compute normal */
   CrossProduct(Trail[HeadIndex].dx, Trail[HeadIndex].dy, Trail[HeadIndex].dz,
                Trail[HeadIndex].sx, Trail[HeadIndex].sy, Trail[HeadIndex].sz,
                &(NormalBuffer[a][0]),
                &(NormalBuffer[a][1]),
                &(NormalBuffer[a][2]));
   NormalBuffer[b][0] = NormalBuffer[a][0];
   NormalBuffer[b][1] = NormalBuffer[a][1];
   NormalBuffer[b][2] = NormalBuffer[a][2];
}

/* Initialize animation loop and objects */
static void Init(void)
{
   /* Initialize OpenGL + GLUT */
   glutInitDisplayMode((unsigned)(GLUT_RGBA | GLUT_DOUBLE));
   glutSetWindow(glutCreateWindow(WINDOW_TITLE));

   glutDisplayFunc(Render);
   glutIdleFunc(Animate);
   glutReshapeFunc(Reshape);
   glutKeyboardFunc(Quit);

   glEnable(GL_BLEND);
   glEnable(GL_LIGHTING);
   glEnable(GL_LIGHT0);
   glEnable(GL_NORMALIZE);
   glShadeModel(GL_SMOOTH);

   /* Note: because this blend function is additive, and addition of
      is commutative, we don't need to sort the vertices before
      drawing them since the result will be independent of Z order.  */
   glBlendFunc(GL_SRC_ALPHA, GL_ONE);

   /* Initialize objects */
   UpdateClock();
   srand((unsigned int)CurrentTime);
   StartTime = CurrentTime;
   FrameCount = 0;

   /* Set initial trail point to be in the past, so that it's replaced
      in the first frame.                                              */
   HeadIndex = TailIndex = 0;
   Trail[0].t = CurrentTime - MIN_TRAIL_STEP_TIME * 2.0;
   Keyframe.t1 = Trail[0].t;
   Keyframe.x2 = RANDOM_POSITION();
   Keyframe.y2 = RANDOM_POSITION();
   Keyframe.z2 = RANDOM_POSITION();
   Keyframe.a2 = RANDOM_ANGLE();
   Keyframe.x3 = 0.0;
   Keyframe.y3 = 0.0;
   Keyframe.z3 = 0.0;
   Keyframe.a0 = 0.0;
}

/* Animate and render a single frame */
static void Render(void)
{
   /* Light settings */
   static GLfloat light_diffuse[4] = {1.0f, 1.0f, 1.0f, 1.0f};
   static GLfloat light_specular[4] = {1.0f, 1.0f, 1.0f, 1.0f};
   static GLfloat light_position[4] = {0.0f, 1.0f, -1.0f, 0.0f};
   static GLfloat light_direction[4] = {0.0f, -1.0f, 1.0f, 0.0f};
   static GLfloat material_specular[4] = {1.0f, 0.9f, 0.9f, 0.8f};

   int width, height, i;
   double dw, dh, alpha;

   UpdateClock();
   UpdateObjects();

   glDrawBuffer(GL_BACK);
   glClear(GL_COLOR_BUFFER_BIT);

   /* Set projection */
   width = glutGet(GLUT_WINDOW_WIDTH);
   height = glutGet(GLUT_WINDOW_HEIGHT);
   if( width > height )
   {
      dw = (double)width / (double)height;
      dh = 1.0;
   }
   else
   {
      dw = 1.0;
      dh = (double)height / (double)width;
   }
   glViewport(0, 0, width, height);
   glMatrixMode(GL_PROJECTION);
   glLoadIdentity();
   glFrustum(-dw, dw, -dh, dh, 10.0, WORLD_SIZE * 100.0);
   gluLookAt(0.0, 0.0, 10.0 * WORLD_SIZE,  0.0, 0.0, 0.0,  0.0, 1.0, 0.0);

   /* Set lighting */
   glLightfv(GL_LIGHT0, GL_DIFFUSE, light_diffuse);
   glLightfv(GL_LIGHT0, GL_SPECULAR, light_specular);
   glLightfv(GL_LIGHT0, GL_POSITION, light_position);
   glLightfv(GL_LIGHT0, GL_SPOT_DIRECTION, light_direction);
   glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, GL_TRUE);
   glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, material_specular);
   glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 4.0f);

   glColorMaterial(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE);
   glEnable(GL_COLOR_MATERIAL);

   /* Draw trail */
   glMatrixMode(GL_MODELVIEW);
   glLoadIdentity();

   glEnableClientState(GL_VERTEX_ARRAY);
   glEnableClientState(GL_NORMAL_ARRAY);
   glVertexPointer(3, GL_DOUBLE, 0, (GLvoid*)VertexBuffer);
   glNormalPointer(GL_DOUBLE, 0, (GLvoid*)NormalBuffer);

   glBegin(GL_QUAD_STRIP);
      alpha = MIN_ALPHA;
      glColor4d(0.0, 0.0, 0.0, 0.0);
      for(i = TailIndex; i != HeadIndex; i = (i + 1) % MAX_TRAIL_POINTS)
      {
         glArrayElement(i * 2);
         glArrayElement(i * 2 + 1);
         if( alpha < MAX_ALPHA )
         {
            alpha += ALPHA_STEP;
            glColor4d(1.0, 0.53, 0.67, alpha);
         }
      }
   glEnd();

   glutSwapBuffers();
   glFlush();
}

/* Force redraw on clock */
static void Animate(void)
{
   glutPostRedisplay();
}

/* Force redraw on reshape */
static void Reshape(/*@unused@*/int w, /*@unused@*/int h)
{
   glutPostRedisplay();
}

/* Exit */
static void Quit(/*@unused@*/unsigned char c,
                 /*@unused@*/int u, /*@unused@*/int v)
{
   glFlush();
   if( CurrentTime - StartTime > 1.0 )
      printf("fps = %f\n", FrameCount / (CurrentTime - StartTime));
   exit(EXIT_SUCCESS);
}

/* Update animation time */
static void UpdateClock(void)
{
   struct timeval t;
   int gettimeofday_status = gettimeofday(&t, NULL);
   assert(gettimeofday_status == 0);
   CurrentTime = (double)t.tv_sec + (double)(t.tv_usec / 1000000.0);
   FrameCount++;
}