#include <assert.h>
#include <math.h>
#include <stdio.h>
#include "mandelbrot.h"

static RGBColor get_color_in_scale(double value,
      RGBColor* colors, unsigned int nof_colors,
      double color_shift, double color_progression) {
   double t = fmod(color_shift + value *
      color_progression, (double) nof_colors);
   unsigned int cp1 = floor(t); t = t - cp1;
   unsigned int cp2 = (cp1 + 1) % nof_colors;
   RGBColor color = interpolate(colors[cp1], colors[cp2], t);
   return color;
}

static RGBColor get_color(unsigned int iterations, double t,
      RGBColor* colors, unsigned int nof_colors,
      double color_shift, double color_progression,
      double maxiterations) {
   return get_color_in_scale(
      log(iterations + t) / log(maxiterations),
      colors, nof_colors, color_shift, color_progression);
}

/* Algorithm based on pseudo-code at
   https://en.wikipedia.org/wiki/Mandelbrot_set#Continuous_.28smooth.29_coloring
*/
RGBColor mandelbrot(double c_x, double c_y, unsigned int maxiterations,
      RGBColor* colors, unsigned int nof_colors,
      double color_shift, double color_progression) {
   const double bailout = 2 << 8;
   const double bailout2 = bailout * bailout;
   const double e2 = M_E * M_E;
   double log_bailout = log(bailout);
   double x = 0;
   double y = 0;
   for (unsigned int iteration = 0; iteration < maxiterations; ++iteration) {
      double x_square = x * x;
      double y_square = y * y;
      double distance2 = x_square + y_square;
      if (distance2 >= bailout2) {
	 double hue; double value; double saturation;
	 double distance = sqrt(distance2);
	 double nu = log(log(distance) / log_bailout) / M_LN2;
	 if (nu > 1) nu = 1;
	 return get_color(iteration, 1 - nu,
	    colors, nof_colors,
	    color_shift, color_progression, maxiterations);
      }
      double next_x = x_square - y_square + c_x;
      double next_y = 2 * x * y + c_y;
      x = next_x; y = next_y;
   }
   return (RGBColor){0, 0, 0};
}