#include <complex>
extern "C" {
#include "ppm.h"
}
#include <stdio.h>
#include <math.h>

//****************************************************************
// the following switches are available at compile time:
// HIREZ (1600x1200), LOREZ (800x600), DIRTYREZ (fast drawing)
// SUPERPOSE (draw the dynamical julia set on top of the moduli one)
// FORMAT_AVI, FORMAT_MPG (produce a movie)
//****************************************************************

using namespace std;

//****************************************************************
// determine quality of image
//****************************************************************
#if defined(HIREZ)
#define CSIZE 1600
#define RSIZE 1200
#define JSIZE 1200
#define JFRAME 2.5
#define JROFFSET 200
#define JCOFFSET 0

#ifndef SUPERPOSE
#define MSIZE 500
#define MFRAME 2.5
#define MROFFSET (CSIZE-MSIZE)
#define MCOFFSET -90
#endif

#define MJULIA_ITER 20
#define JULIA_ITER 48
#elif defined(LOREZ)
#define CSIZE 800
#define RSIZE 600
#define JSIZE 600
#define JFRAME 2.5
#define JROFFSET 100
#define JCOFFSET 0

#ifndef SUPERPOSE
#define MSIZE 250
#define MROFFSET (CSIZE-MSIZE)
#define MCOFFSET -45
#define MFRAME 2.5
#endif

#define MJULIA_ITER 19
#define JULIA_ITER 47
#elif defined(DIRTYREZ)
#define CSIZE 500
#define RSIZE 400
#define JSIZE 400
#define JFRAME 2.5
#define JROFFSET 50
#define JCOFFSET 0

#ifndef SUPERPOSE
#define MSIZE 400
#define MFRAME 2.5
#define MROFFSET 50
#define MCOFFSET 0
#endif

#define MJULIA_ITER 20
#define JULIA_ITER 17
#endif

#ifdef SUPERPOSE
#define MSIZE RSIZE
#define MROFFSET JROFFSET
#define MCOFFSET JCOFFSET
#define MFRAME JFRAME
#endif

//****************************************************************
// determine type of image
//****************************************************************
#ifdef MAKE_GEN_S
#define PATH_COLOR PPM_MAXMAXVAL,0,0
#define GEN_S
#define PATH0_ITER 18
#define PATH1_ITER 17
#define PATH_DELTA 0.015
#define SKIP_FRAMES 50
#define CONSTANT_SPEED
#elif defined(MAKE_GEN_T)
#define PATH_COLOR 0,0,PPM_MAXMAXVAL
#define GEN_T
#define PATH1_ITER 17
#define PATH2_ITER 16
#define PATH_DELTA 0.015
#define SKIP_FRAMES 50
#define CONSTANT_SPEED
#elif defined(MAKE_NBHD_A)
#define PATH_COLOR PPM_MAXMAXVAL,0,0
#define NBHD_A
#define PATH0_ITER 19
#define PATH_DELTA 0.001
#define SKIP_FRAMES 50
#define CONSTANT_SPEED
#define IMAG_THOLD 0.1
#elif defined(MAKE_NBHD_R)
#define PATH_COLOR 0,PPM_MAXMAXVAL,0
#define NBHD_R
#define PATH0_ITER 23
#define PATH1_ITER 22
#define PATH_DELTA 0.0002
#define SKIP_FRAMES 50
#define CONSTANT_SPEED
#define ABS_THOLD 0.1
#endif

//****************************************************************
// variables
//****************************************************************
typedef complex<double> cdouble;

const cdouble Rabbit = cdouble(0.87743883312334638002,0.74486176661974423659),
  Corabbit = cdouble(0.87743883312334638002,-0.74486176661974423659),
  Airplane = -0.75487766624669276005;

#define MJULIA_COLOR PPM_MAXMAXVAL/2,PPM_MAXMAXVAL/2,PPM_MAXMAXVAL/2
#define JULIA_COLOR 0,0,0
#define JULIA_DOT_COLOR PATH_COLOR
#define PATH_LENGTH 1000000

pixel _array[RSIZE][CSIZE];
pixel backup_array[RSIZE][CSIZE];
pixel *array[RSIZE];

//****************************************************************
// hash coding of images. we remember in a bitmask which points
// were already visited at which iteration number
//****************************************************************
#define HSIZE 8192
#define HFRAME 2.5
typedef unsigned long hmask;

hmask harray[HSIZE][HSIZE];
inline hmask *hash(cdouble z) {
  int i = (int) round(HSIZE*(-imag(z) + HFRAME)/(2.0*HFRAME)),
    j = (int) round(HSIZE*(real(z) + HFRAME)/(2.0*HFRAME));
  if (i >= 0 && i < HSIZE && j >= 0 && j < HSIZE)
    return harray[i]+j;
  else
    return harray[0];
}

void hclear (void) {
  memset (harray, 0, sizeof harray);
}

//****************************************************************
// square roots with appropriate sign of real / imaginary part
//****************************************************************
inline cdouble sqrt_re (cdouble z, bool pos)
{
  cdouble w = sqrt(z);
  if (pos ^ (real(w) > 0))
    return -w;
  else return w;
}

inline cdouble sqrt_im (cdouble z, bool pos)
{
  cdouble w = sqrt(z);
  if (pos ^ (imag(w) > 0))
    return -w;
  else return w;
}

// bug in C++ 3.2.3 libraries. should be abs
double cabs (cdouble z) { return sqrt(real(z)*real(z)+imag(z)*imag(z)); }

//****************************************************************
// plot a pixel in the "mother" julia set or the standard one
//****************************************************************
inline void mjulia_plot (cdouble z, int r, int g, int b, int radius)
{
  int i = (int) round(MCOFFSET+MSIZE*(-imag(z) + MFRAME)/(2.0*MFRAME)),
    j = (int) round(MROFFSET+MSIZE*(real(z) + MFRAME)/(2.0*MFRAME));
  if (i >= radius && i < RSIZE-radius && j >= radius && j < CSIZE-radius) {
    for (int ii = i-radius; ii < i+radius; ii++)
      for (int jj = j-radius; jj < j+radius; jj++)
	PPM_ASSIGN(array[ii][jj],r,g,b);
  }
}

inline void julia_plot (cdouble z, int r, int g, int b, int radius)
{
  double y = JCOFFSET+JSIZE*(-imag(z) + JFRAME)/(2.0*JFRAME),
    x = JROFFSET+JSIZE*(real(z) + JFRAME)/(2.0*JFRAME);
  int i = (int) round(y), j = (int) round(x);

  if (i < radius || i >= RSIZE-radius || j < radius || j >= CSIZE-radius)
    return;
  if (radius > 1)
    for (int ii = i-radius; ii < i+radius; ii++)
      for (int jj = j-radius; jj < j+radius; jj++)
	PPM_ASSIGN(array[ii][jj],r,g,b);
  else
    PPM_ASSIGN(array[i][j],r,g,b);
}

//****************************************************************
// construct the list of points on a path in the "mother" julia set
//****************************************************************
int path (int n, int type, bool dir, cdouble list[])
{
  int len, t;

  if (n > 0)
    switch (type) {
    case 0:
      len = path (n-1, 1, dir, list);
      for (int i = 0; i < len; i++)
	list[i] = 1.0 / sqrt_re(1.0-list[i], true);
      return len;
    case 1:
      len = path (n-1, 2, dir, list);
      for (int i = 0; i < len; i++)
	list[i] = 1.0 / sqrt_re(1.0-list[i], true);
      return len;
    case 2:
      len = path (n-1, 0, dir, list);
      for (int i = 0; i < len; i++)
	list[i] = 1.0 / sqrt_im(1.0-list[i], !dir);
      t = path (n-1, 2, !dir, list+len);
      for (int i = len; i < len+t; i++)
	list[i] = 1.0 / sqrt_re(1.0-list[i], false);
      for (int i = 0; i < len; i++)
	list[i+len+t] = -list[i];
      return 2*len+t;
    }
  else {
    len = 0;
    list[0] = Rabbit;
    return 1;
    if (dir)
      list[len++] = Rabbit;
    else
      list[len++] = Corabbit;
    if (type == 2)
      list[len++] = Airplane;
    return len;
  }
}

//****************************************************************
// draw the "mother" julia set
//****************************************************************
void mjulia_recur (cdouble z, int n)
{
  if (n > 0) {
    hmask *h = hash (z);
    if (*h & 1 << (n-1))
      return;
    *h |= 1 << (n-1);
    
    cdouble p = 1.0/sqrt(1.0-z);

    mjulia_recur (p, n-1);
    mjulia_recur (-p, n-1);
  } else
    mjulia_plot (z, MJULIA_COLOR, 1);
}

void mjulia (void)
{
  hclear ();
  mjulia_recur (Rabbit, MJULIA_ITER);
}

//****************************************************************
// draw a julia set
//****************************************************************
void julia_recur (cdouble z, cdouble w[], int n)
{
  if (n > 0) {
    hmask *h = hash (z);
    if (*h & 1 << (n-1))
      return;
    *h |= 1 << (n-1);
    
    cdouble p = w[n-1]*sqrt(1.0-z);
    n--;
    julia_recur (p, w, n);
    julia_recur (-p, w, n);
  } else
    julia_plot (z, JULIA_COLOR, 1);
}

const cdouble Rfixed = cdouble(0.68650690670529541302,0.25803248958615534467);

void julia (cdouble w0)
{
  cdouble w[JULIA_ITER];

  w[0] = w0;
  for (int i = 1; i < JULIA_ITER; i++)
    w[i] = 1.0 - 1.0 / (w[i-1]*w[i-1]);

  hclear ();
  julia_recur (Rfixed, w, JULIA_ITER);
}

cdouble list[PATH_LENGTH];

int main(int argc, char *argv[]) {
  char tempdir[] = "frames_XXXXXX";
  mkdtemp (tempdir);

  ppm_init(&argc,argv);

  if (argc != 2) {
    fprintf(stderr, "Use: %s <filename>\n", argv[0]);
    return -1;
  }

  for (int i = 0; i < RSIZE; array[i++] = _array[i]);
  for (int i = 0; i < RSIZE; i++)
    for (int j = 0; j < CSIZE; j++)
      PPM_ASSIGN(array[i][j],PPM_MAXMAXVAL,PPM_MAXMAXVAL,PPM_MAXMAXVAL);

  if (0) {
    julia (Rabbit);
    FILE *f = popen("ppmtojpeg > test.jpeg","w");
    ppm_writeppm(f,array,CSIZE,RSIZE,PPM_MAXMAXVAL,0);
    pclose(f);
    return 0;
  }

  mjulia ();

  int len;
#ifdef GEN_S
  len = path (PATH0_ITER, 0, true, list);
  len += path (PATH1_ITER, 1, false, list+len);
#elif defined(GEN_T)
  len = path (PATH1_ITER, 1, true, list);
  len += path (PATH2_ITER, 2, false, list+len);
#elif defined(NBHD_A)
  len = path (PATH0_ITER, 0, true, list);
  int top, bottom;
  for (top = 0; imag(list[top]) > IMAG_THOLD; top++);
  for (bottom = len-1; imag(list[bottom]) < -IMAG_THOLD; bottom--);
  len = bottom - top;
  for (int i = 0; i < len; i++)
    list[i] = list[i+top];
#elif defined(NBHD_R)
  len = path (PATH1_ITER, 1, false, list);
  int sublen;
  for (sublen = len-1; cabs(list[sublen]-Rabbit) < ABS_THOLD; sublen--);
  for (int i = 0; i < len-sublen; i++)
    list[i] = list[sublen+i];
  len = len-sublen;
  sublen = path(PATH0_ITER, 0, true, list+len);
  for (; cabs(list[len]-Rabbit) < ABS_THOLD; len++);
#else
#error Should define GEN_S or GEN_T
#endif
  fprintf(stderr, "# Path with %d points\n", len);

  cdouble z = Rabbit+1.0; // out of the way
  int frame = 0;
  for (int i = 0; i < len; i++) {
    mjulia_plot (list[i], PATH_COLOR, 2);
#ifdef CONSTANT_SPEED
    if (cabs(z-list[i]) > PATH_DELTA)
#else
    if (i % SKIP_FRAMES == 0)
#endif
      { z = list[i];

      fprintf(stderr,"# Computing frame %d\n", frame);

      memcpy (backup_array, _array, sizeof _array); // save clean page
      julia (list[i]);
      julia_plot(list[i], JULIA_DOT_COLOR, 3);

      char s[100];
#ifdef FORMAT_AVI
      sprintf (s, "ppmtojpeg > %s/frame_%04d.jpeg", tempdir, frame++);
#elif defined(FORMAT_MPG__)
      sprintf (s, "pnmtopng > %s/frame_%04d.png", tempdir, frame++);
#else
      sprintf (s, "cat > %s/frame_%04d.ppm", tempdir, frame++);
#endif      
      FILE *f = popen(s,"w");
      ppm_writeppm(f,array,CSIZE,RSIZE,PPM_MAXMAXVAL,0);
      pclose(f);
      memcpy (_array, backup_array, sizeof _array); // restore page w/o julia
    }
  }
  char encodestr[200];
#ifdef FORMAT_AVI
  sprintf(encodestr, "mencoder \"mf://%s/frame_*.jpeg\" -mf fps=25 -o %s -ovc lavc -lavcopts vcodec=mpeg4", tempdir, argv[1]);
#elif defined(FORMAT_MPG)
  sprintf(encodestr, "ffmpeg/ffmpeg -i \"%s/frame_%%04d.ppm\" %s", tempdir, argv[1]);
#endif
  system(encodestr);

  return 0;
}