/* code released 2022-03-22 */
/*
gcc -Wall pi.c -lcairo -lX11 -lm -o pi
./pi -stroke-pi -stroke-dots -n 1 -file Z
ffmpeg -f rawvideo -s 512x512 -r 20 -pix_fmt bgra -i Z -vf palettegen palette.png
ffmpeg -f rawvideo -s 512x512 -r 25 -pix_fmt bgra -i Z -i palette.png -filter_complex paletteuse movie.gif
 */

#include <stdio.h>
#include <stdlib.h>
#include <X11/Xlib.h>
#include <X11/Xutil.h>
#include <sys/select.h>
#include <unistd.h>
#include <cairo/cairo.h>
#include <cairo/cairo-xlib.h>
#include <string.h>
#include <time.h>
#include <math.h>

/* global options */
int slow = 0;
int circle_count = 8;
int reverse = 0;
int stroke_pi      = 0;
int stroke_dots    = 0;
int stroke_circles = 1;
int stroke_path    = 1;

FILE *outfile = NULL;

/*******************************/
/* main logic                  */
/*******************************/

/*
 *      =============
 *         ||   ||
 *         ||   ||
 *         || . ||
 *         ||   ||
 *         ||   ||
 */

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

point_t pi[] = {
  { -2, -3 },
  { -2,  2 },
  { -4,  2 },
  { -4,  3 },
  {  4,  3 },
  {  4,  2 },
  {  2,  2 },
  {  2, -3 },
  {  1, -3 },
  {  1,  2 },
  { -1,  2 },
  { -1, -3 }
};

typedef struct {
  double angle;
  double phase;
  double magnitude;
} bin_t;

void dft(bin_t *out, point_t *in, int N)
{
  int k;
  int n;

  /* Xk = sum(n=0, N-1, xn * exp(-i * 2 * pi / N * k * n))
   * cs = cos(...)
   * ss = sin(...)
   * exp(...) = cs     + i * ss
   * xn       = xin[n] + i * yin[n]
   * xn * exp(...) = cs * xin[n] - ss * yin[n] + i * (ss * xin[n] + cs * yin[n])
   */
  for (k = 0; k < N; k++) {
    double cx, cy;
    cx = cy = 0;
    for (n = 0; n < N; n++) {
      double cs = cos(-2 * M_PI / N * k * n);
      double ss = sin(-2 * M_PI / N * k * n);
      cx += cs * in[n].x - ss * in[n].y;
      cy += ss * in[n].x + cs * in[n].y;
    }
    out[k].angle = 2 * M_PI / N * k;
    if (out[k].angle > M_PI)
      out[k].angle -= 2 * M_PI;
    out[k].phase = atan2(cy, cx);
    out[k].magnitude = sqrt(cx * cx + cy * cy) / N;
  }
}

double len(point_t *a, point_t *b)
{
  int x = a->x - b->x;
  int y = a->y - b->y;
  return sqrt(x * x + y * y);
}

double linear_interpolate(double x, double x0, double x1, double y0, double y1)
{
  /* (x-x0) / (x1-x0) = (y-y0) / (y1-y0)
   * y = y0 + (x-x0) * (y1-y0) / (x1-x0)
   */
  return y0 + (x-x0) * (y1-y0) / (x1-x0);
}

void lines_to_points(point_t *in, int in_len, point_t *out, int out_len)
{
  double total_len = 0;
  int i;
  for (i = 0; i < in_len - 1; i++) {
    total_len += len(&in[i], &in[i+1]);
  }
  total_len += len(&in[in_len-1], &in[0]);
  printf("total_len %g\n", total_len);

  for (i = 0; i < out_len; i++) {
    double wanted_len = total_len / out_len * i;
    int segment_start_index = 0;
    int segment_end_index = 1;
    double segment_len = len(&in[segment_start_index], &in[segment_end_index]);
    while (wanted_len >= segment_len) {
      wanted_len -= segment_len;
      segment_start_index++;
      segment_end_index++;
      segment_end_index %= in_len;
      segment_len = len(&in[segment_start_index], &in[segment_end_index]);
    }
    out[i].x = linear_interpolate(wanted_len, 0, segment_len, in[segment_start_index].x, in[segment_end_index].x);
    out[i].y = linear_interpolate(wanted_len, 0, segment_len, in[segment_start_index].y, in[segment_end_index].y);
  }
}

int compare_dots_dft(const void *p1, const void *p2)
{
  const bin_t *a = p1;
  const bin_t *b = p2;
  if (a->magnitude > b->magnitude) return -1;
  if (a->magnitude < b->magnitude) return 1;
  return 0;
}

/*******************************/
/* GUI code and main function  */
/*******************************/

int fps = 25;

#define WIDTH (512)
#define HEIGHT (512)
#define SCALE (50.)

cairo_surface_t *cs;
cairo_surface_t *csout;

void plot_circle_0(cairo_t *c, double center_x, double center_y, double x, double y, double magnitude)
{
  cairo_set_line_width(c, 5);
  cairo_set_source_rgb(c, 0.4, 0.2, 0.2);
  cairo_arc(c, center_x * SCALE + 256, -center_y * SCALE + 256, magnitude * SCALE, 0, M_PI * 2);
  cairo_fill(c);
}

void plot_circle_1(cairo_t *c, double center_x, double center_y, double x, double y, double magnitude)
{
  cairo_set_line_width(c, 5);
  cairo_set_source_rgb(c, 0.6, 0.4, 0.4);
  cairo_arc(c, center_x * SCALE + 256, -center_y * SCALE + 256, magnitude * SCALE, 0, M_PI * 2);
  cairo_stroke(c);
  cairo_arc(c, center_x * SCALE + 256, -center_y * SCALE + 256, 5, 0, M_PI * 2);
  cairo_fill(c);
}

void plot_circle_2(cairo_t *c, double center_x, double center_y, double x, double y, double magnitude)
{
  cairo_set_line_width(c, 3);
  cairo_set_source_rgb(c, 0.8, 0.9, 0.4);
  cairo_move_to(c, center_x * SCALE + 256, -center_y * SCALE + 256);
  cairo_line_to(c, (center_x + x) * SCALE + 256, -(center_y + y) * SCALE + 256);
  cairo_stroke(c);
  cairo_arc(c, (center_x + x) * SCALE + 256, -(center_y + y) * SCALE + 256, 5, 0, M_PI * 2);
  cairo_fill(c);
}

void paint(int _curtime, point_t *dots, bin_t *dots_dft, int dots_size)
{
  cairo_t *c;
  int i;
  point_t center;
  double t, t_start, t_end;
  double curtime;

  if (slow) {
    /* go 5 times slower */
    _curtime %= dots_size * 2 * 5;
    curtime = _curtime / 5.;
  } else {
    _curtime %= dots_size * 2;
    curtime = _curtime;
  }

  c=cairo_create(cs);

  /* clear screen */
  cairo_rectangle(c, 0, 0, WIDTH, HEIGHT);
  cairo_set_source_rgb(c, 0.2, 0.2, 0.2);
  cairo_fill(c);

  if (stroke_pi) {
    /* stroke PI */
    cairo_set_source_rgb(c, 0.9, 0.8, 0.8);
    cairo_move_to(c, pi[0].x * SCALE + 256, -pi[0].y * SCALE + 256);
    for (i = 0; i < sizeof(pi) / sizeof(pi[0]); i++) {
      cairo_line_to(c, pi[i].x * SCALE + 256, -pi[i].y * SCALE + 256);
    }
    cairo_close_path(c);
    cairo_stroke(c);
  }

  if (stroke_dots) {
    /* stroke dots */
    cairo_set_source_rgb(c, 0.9, 0.5, 0.5);
    for (i = 0; i < dots_size; i++) {
      cairo_arc(c, dots[i].x * 50 + 256, -dots[i].y * 50 + 256, 5, 0, M_PI * 2);
      cairo_fill(c);
    }
  }

  if (stroke_circles) {
    center = (point_t){ 0 , 0 };
    for (i = 0; i < circle_count; i++) {
      int j = reverse ? circle_count - 1 - i : i;
      double x = dots_dft[j].magnitude * cos(dots_dft[j].angle * curtime + dots_dft[j].phase);
      double y = dots_dft[j].magnitude * sin(dots_dft[j].angle * curtime + dots_dft[j].phase);
      point_t new_center = center;
      new_center.x += x;
      new_center.y += y;
      plot_circle_0(c, center.x, center.y, x, y, dots_dft[j].magnitude);
      center = new_center;
    }
    center = (point_t){ 0 , 0 };
    for (i = 0; i < circle_count; i++) {
      int j = reverse ? circle_count - 1 - i : i;
      double x = dots_dft[j].magnitude * cos(dots_dft[j].angle * curtime + dots_dft[j].phase);
      double y = dots_dft[j].magnitude * sin(dots_dft[j].angle * curtime + dots_dft[j].phase);
      point_t new_center = center;
      new_center.x += x;
      new_center.y += y;
      plot_circle_1(c, center.x, center.y, x, y, dots_dft[j].magnitude);
      center = new_center;
    }
    center = (point_t){ 0 , 0 };
    for (i = 0; i < circle_count; i++) {
      int j = reverse ? circle_count - 1 - i : i;
      double x = dots_dft[j].magnitude * cos(dots_dft[j].angle * curtime + dots_dft[j].phase);
      double y = dots_dft[j].magnitude * sin(dots_dft[j].angle * curtime + dots_dft[j].phase);
      point_t new_center = center;
      new_center.x += x;
      new_center.y += y;
      plot_circle_2(c, center.x, center.y, x, y, dots_dft[j].magnitude);
      center = new_center;
    }
  }

  if (stroke_path) {
    cairo_set_source_rgb(c, 0.7, 1, 0.8);
    if (curtime <= dots_size) {
      t_start = 0;
      t_end = curtime;
    } else {
      t_start = curtime - dots_size;
      t_end = dots_size;
    }
    for (t = t_start; t <= t_end + 0.05; t+=.1) {
      center = (point_t){ 0 , 0 };
      for (i = 0; i < circle_count; i++) {
        double x = dots_dft[i].magnitude * cos(dots_dft[i].angle * t + dots_dft[i].phase);
        double y = dots_dft[i].magnitude * sin(dots_dft[i].angle * t + dots_dft[i].phase);
        point_t new_center = center;
        new_center.x += x;
        new_center.y += y;
        center = new_center;
      }
      if (t == 0)
        cairo_move_to(c, center.x * SCALE + 256, -center.y * SCALE + 256);
      else
        cairo_line_to(c, center.x * SCALE + 256, -center.y * SCALE + 256);
    }
    if (curtime == dots_size)
      cairo_close_path(c);
    cairo_stroke(c);
  }

  cairo_destroy(c);

  if (outfile == NULL) {
    /* and put it all to the xpixmap */
    c=cairo_create(csout);
    cairo_set_source_surface(c, cs, 0, 0);
    cairo_rectangle(c, 0, 0, WIDTH, HEIGHT);
    cairo_paint(c);
    cairo_destroy(c);
  } else {
    cairo_surface_t *img = cairo_image_surface_create(CAIRO_FORMAT_ARGB32, WIDTH, HEIGHT);
    c=cairo_create(img);
    cairo_set_source_surface(c, cs, 0, 0);
    cairo_rectangle(c, 0, 0, WIDTH, HEIGHT);
    cairo_paint(c);
    cairo_destroy(c);
    cairo_surface_flush(img);
    unsigned char *d = cairo_image_surface_get_data(img);
    fwrite(d, 512*512*4, 1, outfile);
    cairo_surface_destroy(img);
  }
}

void expose(Display *d, Window w, Pixmap p)
{
  XCopyArea(d, p, w, DefaultGC(d, DefaultScreen(d)), 0,0,
            WIDTH, HEIGHT, 0,0);
}

void usage(void)
{
  printf("options:\n");
  printf("    -slow\n");
  printf("        move 5 times slower (plots 5 times more images)\n");
  printf("    -n <count>\n");
  printf("        use <count> circles, default 8 max 64\n");
  printf("    -stroke-pi\n");
  printf("    -stroke-dots\n");
  printf("    -stroke-circles\n");
  printf("    -stroke-path\n");
  printf("    -no-stroke\n");
  printf("        controls what is plotted, default is circles+path\n");
  printf("        options are processed in order, so to stroke only PI,\n");
  printf("        pass '-no-stroke -stroke-pi'\n");
  printf("    -frame <number>\n");
  printf("        only plot the given frame\n");
  printf("    -file <filename>\n");
  printf("        dump to file instead of screen\n");
  printf("    -reverse\n");
  printf("        use circles from smaller to bigger rather than the opposite\n");
  exit(1);
}

int main(int n, char **v)
{
  Display *d;
  Window w;
  Pixmap p;
  fd_set r;
  int repaint;
  int redraw;
  int i;
  int ret;
  struct timeval time_out;
  int curtime = 0;
  point_t dots[64];
  bin_t dots_dft[64];
  int frame = -1;
  char *filename = NULL;
  int first_frame, last_frame;

  for (i = 1; i < n; i++) {
    if (!strcmp(v[i], "-slow")) { slow = 1; continue; }
    if (!strcmp(v[i], "-n")) {if (i>n-2)usage(); circle_count = atoi(v[++i]); continue; }
    if (!strcmp(v[i], "-no-stroke")) { stroke_pi = stroke_dots = stroke_circles = stroke_path = 0; continue; }
    if (!strcmp(v[i], "-stroke-pi")) { stroke_pi = 1; continue; }
    if (!strcmp(v[i], "-stroke-dots")) { stroke_dots = 1; continue; }
    if (!strcmp(v[i], "-stroke-circles")) { stroke_circles = 1; continue; }
    if (!strcmp(v[i], "-stroke-path")) { stroke_path = 1; continue; }
    if (!strcmp(v[i], "-frame")) {if (i>n-2)usage(); frame = atoi(v[++i]); continue; }
    if (!strcmp(v[i], "-file")) {if (i>n-2)usage(); filename = v[++i]; continue; }
    if (!strcmp(v[i], "-reverse")) { reverse = 1; continue; }
    usage();
  }

  if (circle_count < 0 || circle_count > 64) {
    printf("error: -n takes values between 0 and 64\n");
    return 1;
  }

  lines_to_points(pi, sizeof(pi) / sizeof(pi[0]), dots, sizeof(dots) / sizeof(dots[0]));
  dft(dots_dft, dots, sizeof(dots) / sizeof(dots[0]));
  qsort(dots_dft, sizeof(dots_dft) / sizeof(dots_dft[0]), sizeof(dots_dft[0]), compare_dots_dft);
  for (i = 0; i < sizeof(dots_dft) / sizeof(dots_dft[0]); i++)
    printf("%d angle %g phase %g magnitude %g\n", i,
           dots_dft[i].angle,
           dots_dft[i].phase,
           dots_dft[i].magnitude);

  cs = cairo_image_surface_create(CAIRO_FORMAT_ARGB32,
                                  WIDTH, HEIGHT);

  if (frame == -1) {
    first_frame = 0;
    last_frame = 64*2;
    if (slow) last_frame *= 5;
  } else {
    first_frame = frame;
    last_frame = first_frame + 1;
  }

  /* output to file */
  if (filename) {
    outfile = fopen(filename, "w");
    if (outfile == NULL) { perror(filename); return 1; }
    for (i = first_frame; i < last_frame; i++)
      paint(i, dots, dots_dft, sizeof(dots)/ sizeof(dots[0]));
    if (fclose(outfile)) { perror(filename); return 1; }
    return 0;
  }

  /* x window mode */

  d = XOpenDisplay(0); if (!d) abort();
  w = XCreateSimpleWindow(d, DefaultRootWindow(d), 0, 0,
                          WIDTH, HEIGHT, 0, 0, 0);
  XSelectInput(d, w, KeyPressMask | ExposureMask);
  XMapWindow(d, w);

  p = XCreatePixmap(d, w, WIDTH,HEIGHT,
                    DefaultDepth(d, DefaultScreen(d)));
  csout = cairo_xlib_surface_create(d, p, DefaultVisual(d, DefaultScreen(d)),
                                    WIDTH, HEIGHT);

  repaint = 1;
  redraw = 0;

  while (1) {
    XEvent ev;
    while (XPending(d)) {
      XNextEvent(d, &ev);
      switch (ev.type) {
      case Expose: redraw = 1; break;
      case KeyPress: {
        char buf[65];
        KeySym ks;
        XLookupString(&ev.xkey, buf, 64, &ks, NULL);
        switch (ks) {
          case XK_Right: case 'n': printf("n"); fflush(stdout); break;
          case XK_Left:  case 'p': printf("p"); fflush(stdout); break;
          case XK_Escape: case 'q': exit(1);
        }
        break;
      }
      default: break;
      }
    }
    if (repaint) { paint(frame == -1 ? curtime : frame, dots, dots_dft, sizeof(dots)/ sizeof(dots[0])); redraw = 1; repaint = 0; }
    if (redraw) { expose(d, w, p); redraw = 0; }
    FD_ZERO(&r);
    FD_SET(ConnectionNumber(d), &r);
    XFlush(d);
    /* this time keeping is wrong, we sleep too much because we don't
     * take into account computation time of previous frame.
     */
    time_out.tv_sec = 0;
    time_out.tv_usec = 1000000/fps;
    ret = select(ConnectionNumber(d)+1, &r, 0, 0, &time_out);
    if (ret == -1) abort();
    if (ret == 0) { repaint = 1; curtime++; }
  }

  return 0;
}