/* 2022-11-13
 *
 * a variation on http://sedcore.eu.org/small/koch.c
 * check it for basic instructions
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/****************************************************************************/
/* bezier curve begin                                                       */
/****************************************************************************/

/* bezier curves are super nice for trajectories */

typedef struct {
  float x;
  float y;
} point_t;

typedef struct {
  point_t p0;
  point_t p1;
  point_t p2;
  point_t p3;
} bezier_t;

/* cubic bezier (4 control points) */
point_t bezier(bezier_t *b, float t)
{
  float _1_t = 1 - t;
  float _1_t_2 = _1_t * _1_t;
  float _1_t_3 = _1_t_2 * _1_t;
  float t_2 = t * t;
  float t_3 = t * t_2;
  return (point_t){
    x: _1_t_3 * b->p0.x + 3 * _1_t_2 * t * b->p1.x + 3 * _1_t * t_2 * b->p2.x + t_3 * b->p3.x,
    y: _1_t_3 * b->p0.y + 3 * _1_t_2 * t * b->p1.y + 3 * _1_t * t_2 * b->p2.y + t_3 * b->p3.y
  };
}

void position_bezier(bezier_t *b, float angle, float x1, float y1)
{
  float c = cos(angle);
  float s = sin(angle);
  float x;
  float y;
  x = x1 + c * b->p0.x - s * b->p0.y;
  y = y1 + s * b->p0.x + c * b->p0.y;
  b->p0.x = x;
  b->p0.y = y;
  x = x1 + c * b->p1.x - s * b->p1.y;
  y = y1 + s * b->p1.x + c * b->p1.y;
  b->p1.x = x;
  b->p1.y = y;
  x = x1 + c * b->p2.x - s * b->p2.y;
  y = y1 + s * b->p2.x + c * b->p2.y;
  b->p2.x = x;
  b->p2.y = y;
  x = x1 + c * b->p3.x - s * b->p3.y;
  y = y1 + s * b->p3.x + c * b->p3.y;
  b->p3.x = x;
  b->p3.y = y;
}

/****************************************************************************/
/* bezier curve end                                                         */
/****************************************************************************/

typedef struct {
  unsigned char r, g, b;
} color_t;

float screen[160*90];

void point(int x, int y, float color)
{
  if (x < 0 || x > 159 || y < 0 || y > 89) return;
  y = 89 - y;

  if (color < screen[y * 160 + x])
    screen[y * 160 + x] = color;
}

int is_right(int x, int y, float x1, float y1, float x2, float y2)
{
  float v1x = x2 - x1;
  float v1y = y2 - y1;
  float v2x = x - x1;
  float v2y = y - y1;
  float det = v1x * v2y - v2x * v1y;
  return det <= 0;
}

color_t palette(float col)
{
  if (col == 10000) return (color_t){140, 160, 255};
  color_t a = (color_t) { 0, 70, 255 };
  color_t b = (color_t) { 255, 255, 255 };
  if (col < 0) col = 0;
  if (col > 1) col = 1;
  return (color_t){
    r : b.r * col + a.r * (1 - col),
    g : b.g * col + a.g * (1 - col),
    b : b.b * col + a.b * (1 - col),
  };
}

float weight_color(float x, float y,
                   float x1, float y1, float x2, float y2, float x3, float y3,
                   float color_a, float color_b, float color_c)
{
  /* https://codeplea.com/triangular-interpolation
   * basically Barycentric Coordinates
   * x = w1 * x1 + w2 * x2 + w3 * x3
   * y = w1 * y1 + w2 * y2 + w3 * y3
   * w1 + w2 + w3 = 1
   * and solve
   */
  float w1, w2, w3;
  w1 = ((y2 - y3) * (x - x3) + (x3 - x2) * (y - y3))
     / ((y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3));
  w2 = ((y3 - y1) * (x - x3) + (x1 - x3) * (y - y3))
     / ((y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3));
  w3 = 1 - w1 - w2;
  return w1 * color_a + w2 * color_b + w3 * color_c;
}

void triangle(float x1, float y1, float x2, float y2, float x3, float y3,
              float color_a, float color_b, float color_c)
{
  int x, y;
  float xmin = x1-1;
  float ymin = y1-1;
  float xmax = x1+1;
  float ymax = y1+1;
  if (x2-1 < xmin) xmin = x2-1;
  if (x3-1 < xmin) xmin = x3-1;
  if (y2-1 < ymin) ymin = y2-1;
  if (y3-1 < ymin) ymin = y3-1;
  if (x2+1 > xmax) xmax = x2+1;
  if (x3+1 > xmax) xmax = x3+1;
  if (y2+1 > ymax) ymax = y2+1;
  if (y3+1 > ymax) ymax = y3+1;
  for (x = xmin; x <= xmax; x++) for (y = ymin; y <= ymax; y++)
    if (is_right(x, y, x1, y1, x2, y2) &&
        is_right(x, y, x2, y2, x3, y3) &&
        is_right(x, y, x3, y3, x1, y1)) {
      point(x, y, weight_color(x, y, x1, y1, x2, y2, x3, y3,
                               color_a, color_b, color_c));
  }
}

void koch(float x1, float y1, float x2, float y2, int frame, int depth,
          float color_a, float color_b)
{
  float x3, y3;
  float x4, y4;
  float x5, y5;
  float color_c = depth * .3;

  frame %= 100;

  if (fabs(x1-x2) < .5 && fabs(y1-y2) < .5)
    return;

  x3 = x1 + (x2 - x1) / 3;
  y3 = y1 + (y2 - y1) / 3;

  x5 = x2 - (x2 - x1) / 3;
  y5 = y2 - (y2 - y1) / 3;

  if (frame < 60) {
    float angle = atan2(y2-y1, x2-x1);
    float length = sqrt((y2-y1)*(y2-y1) + (x2-x1)*(x2-x1));
    bezier_t b = (bezier_t){
      (point_t){x: length / 2, y: 0 },
      (point_t){x: length    , y: length /3 },
      (point_t){x: length,     y: 0 },
      (point_t){x: length,     y: 0 }
    };
    position_bezier(&b, angle, x1, y1);
    point_t p = bezier(&b, frame / 60.);
    x4 = p.x;
    y4 = p.y;
  } else {
    x4 = (x5 + x3) / 2;
    y4 = (y5 + y3) / 2;
  }

  float len1 = sqrt((y4-y1)*(y4-y1) + (x4-x1)*(x4-x1));
  float len2 = sqrt((y4-y2)*(y4-y2) + (x4-x2)*(x4-x2));
  if (len2 > len1) len1 = len2;
  float col_c = (len1 - 90) / 90;
  if (col_c < 0) col_c = 0;
  if (col_c > 1) col_c = 1;
  if (col_c > color_c) color_c = col_c;
  if (color_c < color_a) color_c = color_a;
  if (color_c < color_b) color_c = color_b;

  float c1 = color_a;
  float c2 = color_b;
  float c4 = color_c;
  float c3 = weight_color(x3, y3, x1, y1, x2, y2, x4, y4, c1, c2, c4);
  float c5 = weight_color(x5, y5, x1, y1, x2, y2, x4, y4, c1, c2, c4);

  triangle(x3, y3, x4, y4, x5, y5, c3, c4, c5);
  koch(x1, y1, x3, y3, frame, depth + 1, c1, c3);
  koch(x3, y3, x4, y4, frame, depth + 1, c3, c4);
  koch(x4, y4, x5, y5, frame, depth + 1, c4, c5);
  koch(x5, y5, x2, y2, frame, depth + 1, c5, c2);
}

void clear(void)
{
  float *s = screen;
  int i;
  for (i = 0; i < 160*79; i++, s++) s[0] = 10000;
  for (i = 160*79; i < 160*90; i++, s++) s[0] = 0;
}

void dump(void)
{
  unsigned char out[160*90*4];
  int i;
  for (i = 0; i < 160*90; i++) {
    color_t c = palette(screen[i]);
    out[i * 4] = c.b;
    out[i * 4 + 1] = c.g;
    out[i * 4 + 2] = c.r;
    out[i * 4 + 3] = 255;
  }
  fwrite(out, 160*90*4, 1, stdout);
}

void up(int x, int p)
{
  int y;
  for (y = 0; y < 90; y++) {
    if (y + p > 89)
      screen[y * 160 + x]   = 0;
    else
      screen[y * 160 + x] = screen[(y+p) * 160 + x];
  }
}

void wave(int frame)
{
  int x;
  frame %= 100;
  for (x = 0; x < 160; x++) {
    float c = (cos(2 * M_PI * x / 360. + 2 * M_PI * 1 * frame/100. + 3) + 1) *4 
            + (cos(2 * M_PI * -x / 90. + 2 * M_PI * 1 * frame/100. + 3) + 1) * 2;
    up(x, c);
  }
}

int main(void)
{
  int frame = 60;
  //while (1) {
  while (frame < 100 + 60) {
fprintf(stderr, "frame %d\n", frame);
    clear();
    koch(0, 10, 160, 10, frame, 0, 0, 0);
    wave(frame);
    dump();
    frame++;
  }
  fflush(stdout);
  return 0;
}