/*




  based on the code    
  http://mathr.co.uk/blog/2014-11-02_practical_interior_distance_rendering.html
  and from book 


  compare with :

  http://stackoverflow.com/questions/369438/smooth-spectrum-for-mandelbrot-set-rendering



  gcc -std=c99 -Wall -Wextra -pedantic -O3 -fopenmp c.c -lm

  ./a.out  100 1000 1000 -0.75 0 1.5 s.ppm
  ./a.out  1000 1000 1000 0.3060959246356990742873 0.0237427672737983365906 2e-11 11.ppm
  ./a.out  20000 10000 10000 0.3060959246356990742873 0.0237427672737983365906 2e-10 10n.ppm


  for exp in $(seq -w 0 11)
  do
  maxiter=$(echo 500+(1000*${exp})|bc)
  echo $maxiter
  ./a.out ${maxiter} 10000 10000 0.3060959246356990742873 0.0237427672737983365906 2e-${exp} ${exp}.ppm
  done





  The escape time bands around the Mandelbrot set increase stepwise by 1. 
  We now have the final iterate radius, which we can use to smooth the colouring. 
  If the final radius is large, then the previous iterate was close to escaping. 
  If the final radius is small, then the previous iterate was further from escaping. 
  Here we assume the escape radius is large enough for the approximation ignoring |c| to hold. 
  Squaring gives a curve, which we can linearize using logarithms. 
  We have a value in [0,1], but log 0 is undefined (it heads towards negative infinity), 
  moreover we still want to end up with a value in [0,1]. 
  Adding 1 gives a value in [1,2], and taking the base-2 logarithm gives the required range, 
  as log B 1 = 0 and log B B = 1 for all B. 
  We subtract this value from the integer escape time, 
  giving a fractional escape time with continuity at the transitions between 
  bands and smoothly increasing size.



*/

#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h> // concat



// ----------- const and variables ----------------------------------


// escape time algorithm parameters 
const double escape_radius = 512.0 ;
double escape_radius2 ;
int maxiters ;
// image : integer coordinate [pixels ]
int iWidth ;
int iHeight ;
unsigned char *image ;
// view  = circle on parameter plane in floating point coordinate 
complex double center ;
double radius ;
long double pixel_spacing ; // pixel size 
int iPixelSpacingBits; // precision of computing

//
const char *filename ;
//char *comment;


// ------------ functions ----------------------------------------------


// 

static inline unsigned char *image_new(int iWidth, int iHeight) {
  return malloc(iWidth * iHeight * 3);
}

static inline void image_delete(unsigned char *image) {
  free(image);
}

static inline void image_save_ppm(unsigned char *image, int iWidth, int iHeight, const char *filename) {
  FILE *f = fopen(filename, "wb");
  if (f) { /*write header to the file*/
    fprintf(f, "P6\n");
    fprintf(f, "# Mandelbrot set, center = %.16f ; %.16f; radius = %.16f\n", creal(center), cimag(center), radius  ); 
    fprintf(f, "%d %d\n %d\n", iWidth, iHeight, 255);  
    fwrite(image, iWidth * iHeight * 3, 1, f); // data
    fclose(f);
    fprintf(stderr, "file  %s saved \n", filename);
  } else {
    fprintf(stderr, "ERROR saving `%s'\n", filename);
  }
}

static inline void image_poke(unsigned char *image, int iWidth, int i, int j, int r, int g, int b) {
  int k = (iWidth * j + i) * 3;
  image[k++] = r;
  image[k++] = g;
  image[k  ] = b;
}

static inline void hsv2rgb(double h, double s, double v, double *r, double *g, double *b) {
  double i, f, p, q, t;
  if (s == 0) { *r = *g = *b = v; } else {
    h = 6 * (h - floor(h));
    int ii = i = floor(h);
    f = h - i;
    p = v * (1 - s);
    q = v * (1 - (s * f));
    t = v * (1 - (s * (1 - f)));
    switch(ii) {
    case 0: *r = v; *g = t; *b = p; break;
    case 1: *r = q; *g = v; *b = p; break;
    case 2: *r = p; *g = v; *b = t; break;
    case 3: *r = p; *g = q; *b = v; break;
    case 4: *r = t; *g = p; *b = v; break;
    default:*r = v; *g = p; *b = q; break;
    }
  }
}

static inline void colour_to_bytes(double r, double g, double b, int *r_out, int *g_out, int *b_out) {
  *r_out = fmin(fmax(255 * r, 0), 255);
  *g_out = fmin(fmax(255 * g, 0), 255);
  *b_out = fmin(fmax(255 * b, 0), 255);
}




static inline void colour_pixel(unsigned char *image, int iWidth, int i, int j, int final_n, complex double final_z) {

  double hue, sat, val;
  double dR, dG, dB;
  int iR, iG, iB;


  
  if (final_n == 0) // interior 
    { // black is defined as V=0, independently of H and S
      hue = 0.0;
      sat = 0.0;
      val = 0.0;
    }
  else  // exterior 
    {
      double final_z_abs= log(cabs(final_z)) / log(escape_radius) - 1.0; 
      // hsv color model       
      hue = (final_n - log2(final_z_abs + 1.0)) / 64.0;
      sat = 0.7;
      val = 1.0; 
    }
  
  hsv2rgb(hue, sat, val, &dR, &dG, &dB);
  colour_to_bytes(dR, dG, dB, &iR, &iG, &iB);
  image_poke(image, iWidth, i, j, iR, iG, iB);
}


static inline double cabs2(complex double z) {
  return creal(z) * creal(z) + cimag(z) * cimag(z);
}



static inline void render(unsigned char *image, int n_max, int iWidth, int iHeight, complex double center, double pixel_spacing) {

  

#pragma omp parallel for schedule(dynamic, 1)
  for (int j = 0; j < iHeight; ++j) {
    for (int i = 0; i < iWidth; ++i) {
 
      double x = i + 0.5 - iWidth / 2.0;
      double y = iHeight / 2.0 - j - 0.5;
 
      complex double c = center + pixel_spacing * (x + I * y);
      complex double z = 0;
      int final_n = 0;
      complex double final_z = 0;
 
      // iteration 
      for (int n = 1; n < n_max; ++n) {
	z = z * z + c;
	if (cabs2(z) > escape_radius2) {
	  final_n = n;
	  final_z = z;
	  break;  
	}
      }
      //    
      colour_pixel(image, iWidth, i, j, final_n, final_z);
    }
  }
}



int PrintInfo(long double lPixelSpacing)
{
 
  int iPixelSpacingBits = -log2l( lPixelSpacing); // 
  
  if (iPixelSpacingBits <= 40) { printf("use float \n"); return iPixelSpacingBits;}
  if (iPixelSpacingBits > 60)  
    {printf("use arbitrary precision number, like MPFR \n");
      return iPixelSpacingBits; }
 
  if (iPixelSpacingBits > 50)   
    {printf("use long double \n");
      return iPixelSpacingBits; }
 
  if (iPixelSpacingBits > 40)  
    printf("use double \n");
 
  return iPixelSpacingBits;
}
 




int main(int argc, char **argv) {

  if (argc != 8) {

    fprintf(stderr,
	    "usage: %s maxiters iWidth iHeight creal(center) cimag(center) radius filename\n"
            "example :\n" 
            "%s  100 1024 1024 -0.75 0 1.5 1.ppm\n "
            "number of arguments (argc) should be 7 but now is %d \n",
            argv[0], argv[0], argc-1);
    return 1;
  }

  // read the arguments 
  maxiters = atoi(argv[1]);
  iWidth = atoi(argv[2]);
  iHeight = atoi(argv[3]);
  center = atof(argv[4]) + I * atof(argv[5]);
  radius = atof(argv[6]);
  filename = argv[7];
  // setup 
  pixel_spacing = radius / (iHeight / 2.0); 
  escape_radius2 = escape_radius * escape_radius;
  //
  iPixelSpacingBits =  PrintInfo(pixel_spacing);
  // 
  if ( iPixelSpacingBits <50 )
    {
      image = image_new(iWidth, iHeight);
      render(image, maxiters, iWidth, iHeight, center, pixel_spacing);
      image_save_ppm(image, iWidth, iHeight, filename);
      image_delete(image);
    }
  

  



  return 0;
}