/* 2022-09-07 */

/* inspired by https://twitter.com/Marina_Costant/status/1560219380482539522 */

/* gcc -Wall -o norm norm.c -lm -O3 -ffast-math -fomit-frame-pointer -march=native
 * ./norm > XX
 * ffmpeg -f rawvideo -s 640x360 -r 25 -pix_fmt bgra -i XX -vf palettegen palette.png
 * ffmpeg -f rawvideo -s 640x360 -r 25 -pix_fmt bgra -i XX -i palette.png -filter_complex paletteuse movie.gif
 */

/* not realtime at all, set delta = 0.05 or higher for faster but uglier */
/* or improve dig_t(): we look for the smallest t>1 for which norm() = 1,
 * if it exists of course
 */

/* not super good comments, see norm_infinity.c for better */

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

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

unsigned char data[160*90*4];

double norm(double p, double t, double x0, double y0, double z0, double x1, double y1, double z1)
{
  double x = x0 + t * (x1 - x0);
  double y = y0 + t * (y1 - y0);
  double z = z0 + t * (z1 - z0);
  return   pow(fabs(x), p)
         + pow(fabs(y), p)
         + pow(fabs(z), p);
}

double euclid(double t, double x0, double y0, double z0, double x1, double y1, double z1)
{
  double x = x0 + t * (x1 - x0);
  double y = y0 + t * (y1 - y0);
  double z = z0 + t * (z1 - z0);
  return sqrt(x*x + y*y + z*z);
}

double dig_t(double p, double x0, double y0, double z0, double x1, double y1, double z1)
{
  double t0, t1;
  double n0, n1;
  double delta = .01;

  /* norm is abs(x)^p + abs(y)^p + abs(z)^p
   * we want norm = 1 (we want to be on the "unit sphere" of the norm)
   * current ray is from (x0, y0, z0) to (x1,y1,z1)
   * a point (x,y,z) is on the ray if there is t such that:
   * (x,y,z) = t * ray
   * for t=0 we decide to be at (x0,y0,z0)
   * for t=1 we decide to be at (x1,y1,z1)
   * so for generic point (x,y,z), we have:
   * 0 x0  y0  z0
   * 1 x1  y1  z1
   * t x   y   z
   * so:
   * t = (x - x0) / (x1 - x0)
   * same for y and z
   * and so:
   * x = x0 + t * (x1 - x0)
   * y = y0 + t * (y1 - y0)
   * z = z0 + t * (z1 - z0)
   */

  /* start at t = 1 and end at t = 1 + delta */
  t0 = 1;
  n0 = norm(p, t0, x0, y0, z0, x1, y1, z1);
  if (n0 <= 1) { fprintf(stderr, "bad start\n"); exit(1); }
  t1 = t0 + delta;
  n1 = norm(p, t1, x0, y0, z0, x1, y1, z1);
  while ((   (n0 > 1 && n1 > 1)
          || (n0 < 1 && n1 < 1))
         && t1 < 5) {
    n0 = n1;
    t0 = t1;
    t1 += delta;
    n1 = norm(p, t1, x0, y0, z0, x1, y1, z1);
  }
  if (t1 >= 5) return 1000;
  return t1;
}

#define rotate_x(x, y, z, angle) do {          \
  double y2 = y * cos(angle) - z * sin(angle); \
  double z2 = z * cos(angle) + y * sin(angle); \
  y = y2;                                      \
  z = z2;                                      \
} while (0)

#define rotate_y(x, y, z, angle) do {          \
  double x2 = x * cos(angle) - z * sin(angle); \
  double z2 = z * cos(angle) + x * sin(angle); \
  x = x2;                                      \
  z = z2;                                      \
} while (0)

#define clip(v, a, b) ((v) < (a) ? (a) : (v) > (b) ? (b) : (v))
#define max(a, b) ((a) > (b) ? (a) : (b))

color_t raytrace_norm(double x, double y, int frame)
{
  double x0 = 0, y0 = 0, z0 = 3, x1 = x, y1 = y, z1 = 2;
  double p;
  double t;

  /* rotate (x0,y0,z0) and (x1,y1,z1) */
  double angle = 2 * M_PI * frame / 100;

  rotate_x(x0, y0, z0, -M_PI/7 * cos(2 * M_PI * frame / 80));
  rotate_x(x1, y1, z1, -M_PI/7 * cos(2 * M_PI * frame / 80));

  rotate_y(x0, y0, z0, angle);
  rotate_y(x1, y1, z1, angle);

  /* x0   y0
   * x    y
   * x1   y1
   * (x-x0) / (x1-x0) = (y-y0) / (y1-y0)
   * y = y0 + (x-x0) / (x1-x0) * (y1-y0)
   */
#define P3(x, x0, x1, y0, y1) ((double)(y0) + ((double)(x)-(double)(x0)) / ((double)(x1)-(double)(x0)) * ((double)(y1)-(double)(y0)))
  /* p in [1/2 .. 4] */
  frame %= 200;
  /* 1/2   1   2    4 */
  if (frame < 50)       p = P3(frame, 0, 50, 1./2, 1);
  else if (frame < 75)  p = P3(frame, 50, 75, 1, 2);
  else if (frame < 100) p = P3(frame, 75, 100, 2, 4);
  else if (frame < 125) p = P3(frame, 100, 125, 4, 2);
  else if (frame < 150) p = P3(frame, 125, 150, 2, 1);
  else                  p = P3(frame, 150, 200, 1, 1./2);

  color_t blue = { r : 15, g : 17, b : 23 };
  t = dig_t(p, x0, y0, z0, x1, y1, z1);
  if (t == 1000) return blue;

  /* a palette based on distance from (0,0,0) to the point */
  double d = euclid(t, x0, y0, z0, x1, y1, z1);
  double g = P3(d, 0, 0.6, 125, 10);
  double rb = P3(d, 0, 0.6, 125, 10);
  double ra = P3(d, 0, 2, 0, 255);
  g = clip(g, 0, 255);
  rb = clip(rb, 0, 255);
  double w = clip(P3(d, 1, 1.4, 0, 120), 0, 255);
  return (color_t){
    r: max(clip(ra+rb, 0, 255), w),
    g : max(g, w),
    b : max(20, w)
  };
}

color_t raytrace(int x, int y, int frame)
{
  x -= 80;
  y = 45 - y;

  /* 100 pixels mean 1 */
  double scale = 1. / 100.;
  double x0 = x * scale;
  double y0 = y * scale;

  return raytrace_norm(x0, y0, frame);
}

void point(int x, int y, color_t color)
{
  data[(y * 160 + x) * 4] = color.b;
  data[(y * 160 + x) * 4 + 1] = color.g;
  data[(y * 160 + x) * 4 + 2] = color.r;
  data[(y * 160 + x) * 4 + 3] = 255;
}

void draw(int frame)
{
  int x, y;
  for (y = 0; y < 90; y++) {
#pragma omp parallel for
    for (x = 0; x < 160; x++) {
      point(x, y, raytrace(x, y, frame));
    }
  }
}

unsigned char bigscreen[160*90*4*4*4];

void dump(void)
{
  unsigned char *b = bigscreen;
  unsigned char *screen = data;
  int x, y;
  int i;
  for (y = 0; y < 90; y++) for (i = 0; i < 4; i++) for (x = 0; x < 160; x++) {
    memcpy(b, screen+(y*160+x)*4, 4); b+=4;
    memcpy(b, screen+(y*160+x)*4, 4); b+=4;
    memcpy(b, screen+(y*160+x)*4, 4); b+=4;
    memcpy(b, screen+(y*160+x)*4, 4); b+=4;
  }
  fwrite(bigscreen, 160*90*4*4*4, 1, stdout);
}

int main(void)
{
  int i;

//  while (1)
  for (i = 0; i < 400; i++) {
    fprintf(stderr, "%d\n", i);
    draw(i);
    dump();
  }
  fflush(stdout);
  return 0;
}