// Mandelbrot set, adapted from:
// https://progur.com/2017/02/create-mandelbrot-fractal-javascript.html
//
// Usage: prog [scale [x0 [y0]]]
//
// Mandelbrot:
// 800 1.8 0.4 -> fractal_003
// 1000 1.65 0.3 > fractal_004
//
// Julia:
// 600 1.5 0.5 -> fractal_06, _07 (#if 0)
//
#include <stdio.h>
#include <stdlib.h>
#include <tgmath.h>
#include "ppm.h"

#define ITER 50

// RGB values for Roy G. Biv
// from https://simple.wikipedia.org/wiki/Rainbow
//
struct { int r, g, b; } rgb[] =
{
  {   0,   0,   0 },  // black
  { 255,   0,   0 },  // red
  { 255, 127,   0 },  // orange
  { 255, 255,   0 },  // yellow
  {   0, 255,   0 },  // green
  {   0,   0, 255 },  // blue
  {  75,   0, 130 },  // indigo
  { 139,   0, 255 },  // violet
  { 255, 255, 255 }   // white
};

// interpolate, a + alpha*(b-a), 0.0 <= alpha <= 1.0
//
int interp( double alpha, int a, int b)
{
  return 0.5 + a + alpha*(b-a); // round to int
}

// iterate complex Z = Z*Z + (x,y), starting with Z=(0,0)
// if it converges it is in the set
//
int Mandelbrot( double x, double y)
{
  double complex XY = CMPLX(x,y), Z = XY;

  for( int i = 1; i <= ITER; ++i) { Z = Z*Z + XY; if( fabs(Z) > 2) return i; }

  return 0; // in the set
}

// iterate complex Z = Z*Z + C, starting with Z=(x,y)
// if it converges it is in the set
//
int Julia( double x, double y, double cx, double cy)
{
  double complex C = CMPLX(cx,cy), Z = CMPLX(x,y);

  for( int i = 1; i <= ITER; ++i) { Z = Z*Z + C; if( fabs(Z) > 2) return i; }

  return 0; // in the set
}

int main( int argc, char *argv[])
{
#ifdef MANDELBROT
  ppm a = ppm_new( 600, 600);
#else
  ppm a = ppm_new( 600, 800);
#endif

  double scale = 200.0, x0 = 2, y0 = 1.5;

  if( argc > 1) scale = atof(argv[1]);
  if( argc > 2) x0 = atof(argv[2]);
  if( argc > 3) y0 = atof(argv[3]);

  for( int i = 0; i < a.rows; ++i)
    for( int j = 0; j < a.cols; ++j)
    {
#ifdef MANDELBROT
      int k = Mandelbrot( j/scale - x0, i/scale - y0);
#else
      int k = Julia( j/scale - x0, i/scale - y0, -0.8, 0.156);
#endif

#if 0
      if( k == 0)
      {
        a.r[i][j] = a.g[i][j] = a.b[i][j] = 0;
      }
      else
#endif
      {
        // a.r[i][j] = 255*(v/ITER);  a.g[i][j] = a.b[i][j] = 0;
        //
	// int k = 8*(2*v/ITER + 0.5);  if( k > 8) k = 8;
	// if( k > 8) k = 8;
        // k = 8 - k;
        // a.r[i][j] = rgb[k].r;
        // a.g[i][j] = rgb[k].g;
        // a.b[i][j] = rgb[k].b;

#if 0
        double v = k/2.0;

        if( v >= 8) // too big, make it white
        {
	  a.r[i][j] = a.g[i][j] = a.b[i][j] = 255;
	  continue;
        }
#else
        double v = k/4.0;

        if( v >= 8) v = 8;
        v = 8-v;
#endif

        k = v; // RBG lower index

        double alpha = v - k; // interpolation factor

        a.r[i][j] = interp( alpha, rgb[k].r, rgb[k+1].r);
        a.g[i][j] = interp( alpha, rgb[k].g, rgb[k+1].g);
        a.b[i][j] = interp( alpha, rgb[k].b, rgb[k+1].b);
      }
    }

  ppm_write(a);
}