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

   07/19/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], *point;
int HeadIndex, TailIndex, i, j;

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

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

double Keyframe_t0, Keyframe_t1;
double Keyframe_x0, Keyframe_y0, Keyframe_z0, Keyframe_a0;
double Keyframe_x1, Keyframe_y1, Keyframe_z1, Keyframe_a1;
double Keyframe_x2, Keyframe_y2, Keyframe_z2, Keyframe_a2;
double Keyframe_x3, Keyframe_y3, Keyframe_z3, Keyframe_a3;

double dd, ab, bc, cd, abc, bcd, dt;
double xy, xz, yz, dx, dy, dz;

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

double CubicPoint(double a, double b, double c, double d)
{
   ab = a + dt * (b - a);
   bc = b + dt * (c - b);
   cd = c + dt * (d - c);
   abc = ab + dt * (bc - ab);
   bcd = bc + dt * (cd - bc);
   return abc + dt * (bcd - abc);
}

double CubicDerivative(double a, double b, double c, double d)
{
   ab = dt * dt;
   bc = 1.0 - dt;
   cd = bc * bc;
   abc = 2.0 * dt * bc;
   return 3.0 * (-cd * a + (cd - abc) * b + (abc - ab) * c + ab * 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 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;
   double dw, dh, alpha;

   UpdateClock();

   if( CurrentTime >= Keyframe_t1 )
   {
      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);
   }

   dt = (CurrentTime - Keyframe_t0) / (Keyframe_t1 - Keyframe_t0);

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

   point = &Trail[HeadIndex];
   point->x = CubicPoint(
      Keyframe_x0, Keyframe_x1, Keyframe_x2, Keyframe_x3);
   point->y = CubicPoint(
      Keyframe_y0, Keyframe_y1, Keyframe_y2, Keyframe_y3);
   point->z = CubicPoint(
      Keyframe_z0, Keyframe_z1, Keyframe_z2, Keyframe_z3);

   point->dx = CubicDerivative(
      Keyframe_x0, Keyframe_x1, Keyframe_x2, Keyframe_x3);
   point->dy = CubicDerivative(
      Keyframe_y0, Keyframe_y1, Keyframe_y2, Keyframe_y3);
   point->dz = CubicDerivative(
      Keyframe_z0, Keyframe_z1, Keyframe_z2, Keyframe_z3);
   dd = (point->dx) * (point->dx) +
        (point->dy) * (point->dy) +
        (point->dz) * (point->dz);
   if( dd > 0.000001 )
   {
      dd = sqrt(dd);
      point->dx /= dd;
      point->dy /= dd;
      point->dz /= dd;
   }

   CrossProduct(point->dx, point->dy, point->dz, 0.0, 0.0, 1.0,
                &(point->sx), &(point->sy), &(point->sz));
   dd = point->sx * point->sx + point->sy * point->sy + point->sz * point->sz;
   if( dd < 0.000001 )
   {
      point->sx = 1.0;
      point->sy = point->sz = 0.0;
   }
   else
   {
      dd = sqrt(dd);
      point->sx /= dd;
      point->sy /= dd;
      point->sz /= dd;
   }

   point->a = CubicPoint(
      Keyframe_a0, Keyframe_a1, Keyframe_a2, Keyframe_a3);

   dx = point->dx;
   dy = point->dy;
   dz = point->dz;

   ab = sin(point->a);
   bc = cos(point->a);
   cd = 1.0 - bc;
   xy = dx * dy;
   xz = dx * dz;
   yz = dy * dz;
   abc = point->sx;
   bcd = point->sy;
   dd = point->sz;
   point->sx = (cd * dx * dx + bc) * abc
             + (cd * xy - dz * ab) * bcd
             + (cd * xz + dy * ab) * dd;
   point->sy = (cd * xy + dz * ab) * abc
             + (cd * dy * dy + bc) * bcd
             + (cd * yz - dx * ab) * dd;
   point->sz = (cd * xz - dy * ab) * abc
             + (cd * yz + dx * ab) * bcd
             + (cd * dz * dz + bc) * dd;

   i = HeadIndex * 2;
   j = i + 1;
   VertexBuffer[i][0] = Trail[HeadIndex].x - TRAIL_WIDTH * Trail[HeadIndex].sx;
   VertexBuffer[i][1] = Trail[HeadIndex].y - TRAIL_WIDTH * Trail[HeadIndex].sy;
   VertexBuffer[i][2] = Trail[HeadIndex].z - TRAIL_WIDTH * Trail[HeadIndex].sz;
   VertexBuffer[j][0] = Trail[HeadIndex].x + TRAIL_WIDTH * Trail[HeadIndex].sx;
   VertexBuffer[j][1] = Trail[HeadIndex].y + TRAIL_WIDTH * Trail[HeadIndex].sy;
   VertexBuffer[j][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[i][0]),
                &(NormalBuffer[i][1]),
                &(NormalBuffer[i][2]));
   NormalBuffer[j][0] = NormalBuffer[i][0];
   NormalBuffer[j][1] = NormalBuffer[i][1];
   NormalBuffer[j][2] = NormalBuffer[i][2];

   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_specular);
   glLightfv(GL_LIGHT0, GL_SPECULAR, light_diffuse_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);

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