/* "Experimental Power Stack Mandelbulb for Owen"
 * By: Chris M. Thomasson
 * Source: comp.lang.c digest, Dec 5th 2019.
 * Topic: A 3d volumentric renderer [sic] of my PowerBulb...
 */
 
#include <assert.h>
#include <complex.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <tgmath.h>
 
#define CT_WIDTH  256 // width of plane
#define CT_HEIGHT 256 // height of plane
#define CT_FRAMES 256 // slices of the mbulb
#define CT_ITERS   64 // iterations per pixel
 
/* the canvas */
struct ct_rgb { unsigned char r, g, b; };
 
struct ct_canvas {
  unsigned long width;
  unsigned long height;
  struct ct_rgb* buf;
};
 
static bool ct_canvas_create(struct ct_canvas* const self,
                             unsigned long width,
                             unsigned long height) {
  size_t size = width * height * sizeof(*self->buf);

  self->buf = calloc(1, size);

  if (self->buf) {
    self->width = width;
    self->height = height;
    return true;
  }

  return false;
}
 
static void ct_canvas_destroy( struct ct_canvas const* const self) {
  free(self->buf);
}

static bool ct_canvas_save_ppm(struct ct_canvas const* const self,
                               char const* fname) {
  FILE* fout = fopen(fname, "wb");

  if (fout) {
    char const ppm_head[] = "P6\n"
    "# Chris M. Thomasson Simple 2d Plane ver:0.0.0.0 (pre-alpha)";

    fprintf(fout, "%s\n%lu %lu\n%u\n",
            ppm_head,
            self->width, self->height,
            255U);

    size_t size = self->width * self->height;

    for (size_t i = 0; i < size; ++i) {
//    unsigned int c = self->buf[i];
      struct ct_rgb* c = self->buf + i;
      fprintf(fout, "%c%c%c", c->r, c->g, c->b);
    }

    if (!fclose(fout))
      return true;
  }

  return false;
}
 
 
 
/* the axes */
struct ct_axes { double xmin, xmax, ymin, ymax; };
 
static struct ct_axes ct_axes_from_point(double complex z, double radius) {
  struct ct_axes axes = {
    creal(z) - radius, creal(z) + radius,
    cimag(z) - radius, cimag(z) + radius
  };
  return axes;
}
 
 
 
/* The Plane */
struct ct_plane {
  struct ct_axes axes;
  double xstep;
  double ystep;
};
 
static void ct_plane_init(struct ct_plane* const self,
                          struct ct_axes const* axes,
                          unsigned long width,
                          unsigned long height) {
  self->axes = *axes;
 
  double awidth = self->axes.xmax - self->axes.xmin;
  double aheight = self->axes.ymax - self->axes.ymin;
 
  assert(width > 0 && height > 0 && awidth > 0.0);
 
  double daspect = fabs((double)height / (double)width);
  double waspect = fabs(aheight / awidth);
 
  if (daspect > waspect) {
    double excess = aheight * (daspect / waspect - 1.0);
    self->axes.ymax += excess / 2.0;
    self->axes.ymin -= excess / 2.0;
  }
  else if (daspect < waspect) {
    double excess = awidth * (waspect / daspect - 1.0);
    self->axes.xmax += excess / 2.0;
    self->axes.xmin -= excess / 2.0;
  }
 
  self->xstep = (double)(self->axes.xmax - self->axes.xmin) / (double)width;
  self->ystep = (double)(self->axes.ymax - self->axes.ymin) / (double)height;

  return;
}
 
 
 
/* the plot */
struct ct_plot {
  struct ct_plane plane;
  struct ct_canvas* canvas;
};
 
static void ct_plot_init(struct ct_plot* const self,
                         struct ct_axes const* axes,
                         struct ct_canvas* canvas) {
  ct_plane_init(&self->plane, axes, canvas->width - 1, canvas->height - 1);
  self->canvas = canvas;
}
 
#if (0)
static bool ct_plot_add(struct ct_plot* const self,
                        double complex z,
                        struct ct_rgb const* color) {
  long x = (creal(z) - self->plane.axes.xmin) / self->plane.xstep;
  long y = (self->plane.axes.ymax - cimag(z)) / self->plane.ystep;
 
  if (x > -1 &&
      x < (long)self->canvas->width &&
      y > -1 &&
      y < (long)self->canvas->height) {
    // Now, we can convert to index.
    size_t i = x + y * self->canvas->width;
 
    assert(i < self->canvas->height * self->canvas->width);

    struct ct_rgb exist = self->canvas->buf[i];

    exist.r += 1;

    self->canvas->buf[i] = exist;
    return true;
  }

  return true;
}
#endif

static bool ct_plot_pixel(struct ct_plot* const self,
                          long x,
                          long y,
                          struct ct_rgb const* color) {
  if (x > -1 &&
      x < (long)self->canvas->width &&
      y > -1 &&
      y < (long)self->canvas->height) {
    // Now, we can convert to index.
    size_t i = x + y * self->canvas->width;

    assert(i < self->canvas->height * self->canvas->width);

    self->canvas->buf[i] = *color;
    return true;
  }
 
return false;
}
 
static void ct_plot_clear(struct ct_plot* const self,
                          struct ct_rgb const* color) {
  size_t size = self->canvas->height * self->canvas->width;
  for (size_t i = 0; i < size; ++i)
    self->canvas->buf[i] = *color;
  return;
}
 
static double complex ct_project_to_xy(struct ct_plot* const self,
                                       long x,
                                       long y) {
  double complex p =
          (self->plane.axes.xmin + x * self->plane.xstep) +
          (self->plane.axes.ymax - y * self->plane.ystep) * I;
  return p;
}

static bool ct_plot_point(struct ct_plot* const self,
                          double complex z,
                          struct ct_rgb const* color) {
  long x = (creal(z) - self->plane.axes.xmin) / self->plane.xstep;
  long y = (self->plane.axes.ymax - cimag(z)) / self->plane.ystep;
 
  if (x > -1 &&
      x < (long)self->canvas->width &&
      y > -1 &&
      y < (long)self->canvas->height) {
    // Now, we can convert to index.
    size_t i = x + y * self->canvas->width;

    assert(i < self->canvas->height * self->canvas->width);

    self->canvas->buf[i] = *color;
    return true;
  }
 
  return false;
}
 
/* slow, so what for now... ;^) */
static void ct_circle(struct ct_plot* const plot,
                      double complex c,
                      double radius,
                      unsigned int n) {
  double abase = 6.2831853071 / n;
 
  for (unsigned int i = 0; i < n; ++i) {
    double angle = abase * i;
    double complex z =
            (creal(c) + cos(angle) * radius) +
            (cimag(c) + sin(angle) * radius) * I;
 
    struct ct_rgb rgb = { 255, 255, 255 };
    ct_plot_point(plot, z, &rgb);
  }
  return;
}
 
/* simple vec3 for just what I need */
struct ct_vec3 { double x, y, z; };
 
static double ct_vev3_length(struct ct_vec3 const p) {
  double sum = p.x * p.x + p.y * p.y + p.z * p.z;
  return sqrt(sum);
}
 
static struct ct_vec3 ct_vec3_add(struct ct_vec3 const p,
                                  struct ct_vec3 const addend) {
  struct ct_vec3 sum = {
    p.x + addend.x,
    p.y + addend.y,
    p.z + addend.z
  };
  return sum;
}
 
/* the Power Stack Mandelbulb */
static void ct_iterate_mbulb_pixel(struct ct_plot* const plot,
                                   struct ct_vec3 const z0,
                                   struct ct_vec3 const c0,
                                   double power,
                                   long x,
                                   long y,
                                   unsigned int n) {
  struct ct_vec3 zs = z0;
  struct ct_vec3 cs = c0;
 
  for (unsigned long i = 0; i < n; ++i) {
    double r = ct_vev3_length(zs);
    double rpower = pow(r, power);
 
    double angle0 = atan2(zs.z, sqrt(zs.x * zs.x + zs.y * zs.y));
    double angle1 = atan2(zs.y, zs.x);
 
    struct ct_vec3 znew;
    znew.x = rpower * cos(power * angle1) * cos(power * angle0);
    znew.y = rpower * sin(power * angle1) * cos(power * angle0);
    znew.z = rpower * sin(power * angle0);
 
    zs = ct_vec3_add(znew, cs);

    if (r > 2.0)
      return;
  }
 
  struct ct_rgb rgb = { 0, 255, 0 };
  ct_plot_pixel(plot, x, y, &rgb);
  return;
}
 
/* gain the field... */
static void ct_iterate_mbulb_plane(struct ct_plot* const plot,
                                   double zaxis,
                                   double power,
                                   unsigned int n) {
  for (long y = 0; y < plot->canvas->height; ++y) {
    for (long x = 0; x < plot->canvas->width; ++x) {
      double complex cz = ct_project_to_xy(plot, x, y);
      // Well, our per slice coords...
      struct ct_vec3 z;
      z.x = creal(cz);
      z.y = cimag(cz);
      z.z = zaxis;
      // Iterate the slice
      ct_iterate_mbulb_pixel(plot, z, z, power, x, y, n);
    }
 
    if (! (y % (plot->canvas->height / 10))) {
      printf("ct_iterate_mbulb_plane: %lu of %lu\r",
              y + 1, plot->canvas->height);
      fflush(stdout);
    }
  }
 
  printf("ct_iterate_mbulb_plane: %lu of %lu\n\n",
  plot->canvas->height, plot->canvas->height);
  fflush(stdout);
 
  return;
}
 
/* The Animation Driver */
static void ct_anime(struct ct_plot* const plot,
                     unsigned long frame_n,
                     unsigned long n) {
  assert(frame_n > 0);

  double zmin = 0.;
  double zmax = 1.;
  double zdif = zmax - zmin;
  double zbase = zdif / ((double)frame_n - 1.);
 
  double pmin = 2.;
  double pmax = 10.;
  double pdif = pmax - pmin;
  double pbase = pdif / ((double)frame_n - 1.);
 
  char fname[256] = { '\0' };
 
  for (unsigned long frame_i = 0; frame_i < frame_n; ++frame_i) {
    double zaxis = zmin + zbase * (double)frame_i;
    double power = pmin + pbase * (double)frame_i;
 
    sprintf(fname, "ct_anime_%03lu.ppm", frame_i);

    printf("ct_anime: %lu of %lu, zaxis(%lf), power(%lf), (%s)\n",
            frame_i, frame_n, zaxis, power, fname);
 
    fflush(stdout);
 
    struct ct_rgb black = { 0, 0, 0 };
    ct_plot_clear(plot, &black);
 
    ct_iterate_mbulb_plane(plot, zaxis, power, n);
 
    ct_canvas_save_ppm(plot->canvas, fname);
  }
 
  printf("\n\nct_anime complete! %lu frames: \n", frame_n);
  fflush(stdout);

  return;
}
 
static void ct_test_plane(struct ct_plot* const plot) {
  printf("Generating...\n\n");
  fflush(stdout);
  ct_anime(plot, CT_FRAMES, CT_ITERS);
  return;
}
 
int main(void) {
  struct ct_canvas canvas;
 
  if (ct_canvas_create(&canvas, CT_WIDTH, CT_HEIGHT)) {
    double complex plane_origin = 0+0*I;
    double plane_radius = 1.5;
 
    struct ct_axes axes = ct_axes_from_point(plane_origin, plane_radius);
 
    struct ct_plot plot;
    ct_plot_init(&plot, &axes, &canvas);
 
    ct_test_plane(&plot);
 
    // Our unit circle
    ct_circle(&plot, 0+0*I,  1.0, 2048);
    ct_circle(&plot, 2+0*I,  2.0, 2048);
    ct_circle(&plot, -2+0*I, 2.0, 2048);
    ct_circle(&plot, 0+2*I,  2.0, 2048);
    ct_circle(&plot, 0-2*I,  2.0, 2048);
 
    ct_canvas_save_ppm(&canvas, "ct_canvas.ppm");
 
    ct_canvas_destroy(&canvas);
 
    return EXIT_SUCCESS;
  }

  return EXIT_FAILURE;
}