// http://mathr.co.uk/blog/2014-11-02_practical_interior_distance_rendering.html

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

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

const double pi = 3.141592653589793;
const double infinity = 1.0 / 0.0;
const double colour_modulus = 5.7581917135421046e-2; // (1.0 + 1.0 / (phi * phi)) / 24.0;
const double escape_radius_2 = 512.0 * 512.0;

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

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

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

static inline void image_save_ppm(unsigned char *image, int width, int height, const char *filename) {
  FILE *f = fopen(filename, "wb");
  if (f) {
    fprintf(f, "P6\n%d %d\n255\n", width, height);
    fwrite(image, width * height * 3, 1, f);
    fclose(f);
  } else {
    fprintf(stderr, "ERROR saving `%s'\n", filename);
  }
}

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

static inline void colour_hsv_to_rgb(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_mandelbrot(unsigned char *image, int width, int i, int j, int period) {
  double r, g, b;
  colour_hsv_to_rgb(period * colour_modulus, 0.5, 1.0, &r, &g, &b); // changed b from tanh(distance )to 1.0 
  int ir, ig, ib;
  colour_to_bytes(r, g, b, &ir, &ig, &ib);
  image_poke(image, width, i, j, ir, ig, ib);
}

static inline void render(unsigned char *image, int maxiters, int width, int height, complex double center, double radius) {

  double pixel_spacing = radius / (height / 2.0);

   #pragma omp parallel for schedule(dynamic, 1)
    for (int j = 0; j < height; ++j) {
     for (int i = 0; i < width; ++i) {
      double x = i + 0.5 - width / 2.0;
      double y = height / 2.0 - j - 0.5;
      complex double c = center + pixel_spacing * (x + I * y);
      complex double z = 0;
      complex double dc = 0;
      double minimum_z2 = infinity; // atom domain
      int period = 0;

      // iteration 
      for (int n = 1; n <= maxiters; ++n) {
        dc = 2 * z * dc + 1;
        z = z * z + c;
        double z2 = cabs2(z);

        if (z2 < minimum_z2) {
          minimum_z2 = z2;
          period = n;}}

         colour_mandelbrot(image, width, i, j, period);
 
    }
  }
}

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

  if (argc != 8) {

     
    fprintf(stderr,
	    "usage: %s maxiters width height creal cimag radius filename\n example :\n ./mandelbrot  100 1024 1024 -0.75 0 1.5 1.ppm\n now argc = %d \n", argv[0], argc);
    return 1;
  }

  
  int maxiters = atoi(argv[1]);
  int width = atoi(argv[2]);
  int height = atoi(argv[3]);
  complex double center = atof(argv[4]) + I * atof(argv[5]);
  double radius = atof(argv[6]);
  const char *filename = argv[7];
 
  unsigned char *image = image_new(width, height);

  render(image, maxiters, width, height, center, radius);
  image_save_ppm(image, width, height, filename);
  image_delete(image);

  return 0;
}