/* kanata6.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>

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

double CurrentTime,

       VertexBuffer[1024][3],
       NormalBuffer[1024][3],

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

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

GLfloat materia[4][4] =
{
   {1.0f, 1.0f, 1.0f, 1.0f},
   {0.0f, 1.0f, -1.0f, 0.0f},
   {0.0f, -1.0f, 1.0f, 0.0f},
   {0.9f, 0.93f, 1.0f, 0.8f}
};

double NextRandom(double range)
{
   return range * ((rand() & 32767) - 16384) / 16384.;
}

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 - dt;
   cd = bc * bc;
   abc = 2 * dt * bc;
   return 3 * (-cd * a + (cd - abc) * b + (abc - ab) * c + ab * d);
}

void CrossProduct(double ax, double ay, double az,
                  double bx, double by, double bz)
{
   dx = ay * bz - by * az;
   dy = bx * az - ax * bz;
   dz = ax * by - bx * ay;
}

void Animate(void)
{
   glutPostRedisplay();
}

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

void Quit(unsigned char c, int u, int v)
{
   glFlush();
   exit(0);
}

void UpdateClock()
{
   struct timeval t;
   gettimeofday(&t, NULL);
   CurrentTime = t.tv_sec + t.tv_usec / 1000000.;
}

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() & 32767) / 8192.;
      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;)
      {
         Keyframe_x3 = NextRandom(128);
         Keyframe_y3 = NextRandom(128);
         Keyframe_z3 = NextRandom(128);
         dx = Keyframe_x3 - Keyframe_x0;
         dy = Keyframe_y3 - Keyframe_y0;
         dz = Keyframe_z3 - Keyframe_z0;
         if( dx * dx + dy * dy + dz * dz > 1820 )
            break;
      }
      for(i = 0; i++ < 8;)
      {
         Keyframe_x2 = Keyframe_x3 + 0.4 * NextRandom(128);
         Keyframe_y2 = Keyframe_y3 + 0.4 * NextRandom(128);
         Keyframe_z2 = Keyframe_z3 + 0.4 * NextRandom(128);
         dx = Keyframe_x2 - Keyframe_x3;
         dy = Keyframe_y2 - Keyframe_y3;
         dz = Keyframe_z2 - Keyframe_z3;
         if( dx * dx + dy * dy + dz * dz > 1820 )
            break;
      }

      Keyframe_a3 = NextRandom(pi);
      Keyframe_a2 = Keyframe_a3 + 0.4 * NextRandom(pi);
   }

   for(dt = (CurrentTime - Keyframe_t0) / (Keyframe_t1 - Keyframe_t0);
       TailIndex - HeadIndex &&
         CurrentTime - Trail[TailIndex].t >= 20; /*MAX_TRAIL_LENGTH_TIME*/
       TailIndex = (TailIndex + 1) & 511);

   if( CurrentTime - Trail[HeadIndex].t >= 0.05 ) /*MIN_TRAIL_STEP_TIME*/
   {
      HeadIndex = (HeadIndex + 1) & 511;
      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);
   if( (dd = point->dx * point->dx +
             point->dy * point->dy +
             point->dz * point->dz) > epsilon )
   {
      dd = sqrt(dd);
      point->dx /= dd;
      point->dy /= dd;
      point->dz /= dd;
   }

   CrossProduct(point->dx, point->dy, point->dz, 0, 0, 1);
   point->sx = 1;
   point->sy = point->sz = 0;
   if( (dd = dx * dx + dy * dy + dz * dz) > epsilon )
   {
      dd = sqrt(dd);
      point->sx = dx / dd;
      point->sy = dy / dd;
      point->sz = dz / dd;
   }

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

   ab = sin(point->a =
            CubicPoint(Keyframe_a0, Keyframe_a1, Keyframe_a2, Keyframe_a3));
   bc = cos(point->a);
   cd = 1 - 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;

   VertexBuffer[i = HeadIndex * 2][0] =
      Trail[HeadIndex].x - 16 * Trail[HeadIndex].sx;
   VertexBuffer[i][1] = Trail[HeadIndex].y - 16 * Trail[HeadIndex].sy;
   VertexBuffer[i][2] = Trail[HeadIndex].z - 16 * Trail[HeadIndex].sz;
   VertexBuffer[j = i + 1][0] = Trail[HeadIndex].x + 16 * Trail[HeadIndex].sx;
   VertexBuffer[j][1] = Trail[HeadIndex].y + 16 * Trail[HeadIndex].sy;
   VertexBuffer[j][2] = Trail[HeadIndex].z + 16 * Trail[HeadIndex].sz;

   CrossProduct(Trail[HeadIndex].dx, Trail[HeadIndex].dy, Trail[HeadIndex].dz,
                Trail[HeadIndex].sx, Trail[HeadIndex].sy, Trail[HeadIndex].sz);
   NormalBuffer[j][0] = NormalBuffer[i][0] = dx;
   NormalBuffer[j][1] = NormalBuffer[i][1] = dy;
   NormalBuffer[j][2] = NormalBuffer[i][2] = dz;

   glDrawBuffer(GL_BACK);
   glClear(GL_COLOR_BUFFER_BIT);

   width = glutGet(GLUT_WINDOW_WIDTH);
   height = glutGet(GLUT_WINDOW_HEIGHT);
   dw = dh = 1;
   dd = (width > height)
      ? (dw = width / (double)height)
      : (dh = height / (double)width);
   glViewport(0, 0, width, height);
   glMatrixMode(GL_PROJECTION);
   glLoadIdentity();
   glFrustum(-dw, dw, -dh, dh, 10, 12800);
   gluLookAt(0, 0, 1280,  0, 0, 0,  0, 1, 0);

   glLightfv(GL_LIGHT0, GL_DIFFUSE, *materia);
   glLightfv(GL_LIGHT0, GL_SPECULAR, *materia);
   glLightfv(GL_LIGHT0, GL_POSITION, materia[1]);
   glLightfv(GL_LIGHT0, GL_SPOT_DIRECTION, materia[2]);
   glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, GL_TRUE);
   glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, materia[3]);
   glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 4);

   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, VertexBuffer);
   glNormalPointer(GL_DOUBLE, 0, NormalBuffer);

   glBegin(GL_QUAD_STRIP);
      glColor4d(alpha = 0, 0, 0, 0);
      for(i = TailIndex; i != HeadIndex; i = (i + 1) & 511)
      {
         glArrayElement(i * 2);
         glArrayElement(i * 2 + 1);
         if( alpha < 0.8 )
            glColor4d(0.4, 0.5, 1, alpha += 0.08); /*ALPHA_STEP*/
      }
   glEnd();

   glutSwapBuffers();
   glFlush();
}

void Init()
{
   glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE);
   glutSetWindow(glutCreateWindow("Kanata"));
   pi = atan2(0, -1);
   epsilon = 0.000001;

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

   Keyframe_t1 = Trail[0].t = CurrentTime - 0.05 * 2.0; /*MIN_TRAIL_STEP_TIME*/
   Keyframe_x2 = NextRandom(128);
   Keyframe_y2 = NextRandom(128);
   Keyframe_z2 = NextRandom(128);
   Keyframe_a2 = NextRandom(pi);
   Keyframe_x3 =
   Keyframe_y3 =
   Keyframe_z3 =
   Keyframe_a0 =
   HeadIndex = TailIndex = 0;
}

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