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

   07/18/08
                                                                           */

#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 "Kanata"

#define CONTROL_POINT_OFFSET  0.4

#define MIN_TRAIL_STEP_TIME   0.05
#define MAX_TRAIL_LENGTH_TIME 20.0
#define MAX_TRAIL_POINTS      512
#define MIN_ALPHA             0.0
#define MAX_ALPHA             0.8
#define ALPHA_STEP            0.08

#define WORLD_SIZE            128.0
#define TRAIL_WIDTH           16.0

double CurrentTime;

double StartTime;
int FrameCount;

typedef struct
{
   double t;
   double a;
   double x, y, z;
   double dx, dy, dz;
   double sx, sy, sz;
} TrailPoint;
TrailPoint Trail[MAX_TRAIL_POINTS];
int HeadIndex, TailIndex;

double VertexBuffer[MAX_TRAIL_POINTS * 2][3];
double NormalBuffer[MAX_TRAIL_POINTS * 2][3];

GLfloat light_diffuse[4] = {1.0f, 1.0f, 1.0f, 1.0f};
GLfloat light_specular[4] = {1.0f, 1.0f, 1.0f, 1.0f};
GLfloat light_position[4] = {0.0f, 1.0f, -1.0f, 0.0f};
GLfloat light_direction[4] = {0.0f, -1.0f, 1.0f, 0.0f};
GLfloat material_specular[4] = {0.9f, 0.93f, 1.0f, 0.8f};

typedef struct
{
   double t0, t1;
   double x0, y0, z0, a0;
   double x1, y1, z1, a1;
   double x2, y2, z2, a2;
   double x3, y3, z3, a3;
} KeyframePoint;
KeyframePoint Keyframe;

double NextRandom(double range)
{
   return range * ((rand() & 0x7fff) - 0x4000) / (double)0x4000;
}

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);
}

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);
}

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;
}

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;
}

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;
   }
}

TrailPoint *AddPoint(void)
{
   for(; TailIndex != HeadIndex &&
         CurrentTime - Trail[TailIndex].t >= MAX_TRAIL_LENGTH_TIME;
       TailIndex = (TailIndex + 1) % MAX_TRAIL_POINTS);

   if( CurrentTime - Trail[HeadIndex].t >= MIN_TRAIL_STEP_TIME )
   {
      HeadIndex = (HeadIndex + 1) % MAX_TRAIL_POINTS;
      Trail[HeadIndex].t = CurrentTime;
   }
   return &Trail[HeadIndex];
}

void AddKeyframe(void)
{
   double dx, dy, dz;
   int i;

   Keyframe.t0 = Keyframe.t1;
   Keyframe.t1 = CurrentTime + 0.5 + (rand() & 0x7fff) / (double)0x2000;
   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);
   for(i = 0; i < 8; i++)
   {
      Keyframe.x3 = NextRandom(WORLD_SIZE);
      Keyframe.y3 = NextRandom(WORLD_SIZE);
      Keyframe.z3 = NextRandom(WORLD_SIZE);
      dx = Keyframe.x3 - Keyframe.x0;
      dy = Keyframe.y3 - Keyframe.y0;
      dz = Keyframe.z3 - Keyframe.z0;
      if( dx * dx + dy * dy + dz * dz > WORLD_SIZE * WORLD_SIZE / 9.0 )
         break;
   }
   for(i = 0; i < 8; i++)
   {
      Keyframe.x2 = Keyframe.x3 + CONTROL_POINT_OFFSET * NextRandom(WORLD_SIZE);
      Keyframe.y2 = Keyframe.y3 + CONTROL_POINT_OFFSET * NextRandom(WORLD_SIZE);
      Keyframe.z2 = Keyframe.z3 + CONTROL_POINT_OFFSET * NextRandom(WORLD_SIZE);
      dx = Keyframe.x2 - Keyframe.x3;
      dy = Keyframe.y2 - Keyframe.y3;
      dz = Keyframe.z2 - Keyframe.z3;
      if( dx * dx + dy * dy + dz * dz > WORLD_SIZE * WORLD_SIZE / 9.0 )
         break;
   }

   Keyframe.a3 = NextRandom(PI);
   Keyframe.a2 = Keyframe.a3 + CONTROL_POINT_OFFSET * NextRandom(PI);
}

void SetPosition(TrailPoint *point)
{
   double dt, d;

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

   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);

   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));

   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;
   }

   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));
}

void UpdateObjects(void)
{
   int a, b;

   SetPosition(AddPoint());

   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;

   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];
}

void Animate(void)
{
   glutPostRedisplay();
}

void Reshape(int w, int h)
{
   glutPostRedisplay();
}

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

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

void Render(void)
{
   int width, height, i;
   double dw, dh, alpha;

   UpdateClock();
   UpdateObjects();

   glDrawBuffer(GL_BACK);
   glClear(GL_COLOR_BUFFER_BIT);

   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);

   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);

   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(0.4, 0.5, 1.0, alpha);
         }
      }
   glEnd();

   glutSwapBuffers();
   glFlush();
}

void Init(void)
{
   glutInitDisplayMode(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);

   glBlendFunc(GL_SRC_ALPHA, GL_ONE);

   UpdateClock();
   srand((unsigned)CurrentTime);
   StartTime = CurrentTime;
   FrameCount = 0;

   HeadIndex = TailIndex = 0;
   Trail[0].t = CurrentTime - MIN_TRAIL_STEP_TIME * 2.0;
   Keyframe.t1 = Trail[0].t;
   Keyframe.x2 = NextRandom(WORLD_SIZE);
   Keyframe.y2 = NextRandom(WORLD_SIZE);
   Keyframe.z2 = NextRandom(WORLD_SIZE);
   Keyframe.a2 = NextRandom(PI);
   Keyframe.x3 = 0.0;
   Keyframe.y3 = 0.0;
   Keyframe.z3 = 0.0;
   Keyframe.a0 = 0.0;
}

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