/* 
 c console program:
  


https://arxiv.org/abs/cs/0609129
One approach to the digital visualization of hedgehogs in holomorphic dynamics
by Alessandro Rosa
see Tuning the number of branches ( page 27)

 -------------------------------         
 2. technic of creating ppm file is  based on the code of Claudio Rocchini
 http://en.wikipedia.org/wiki/Image:Color_complex_plot.jpg
 create 24 bit color graphic file ,  portable pixmap file = PPM 
 see http://en.wikipedia.org/wiki/Portable_pixmap
 to see the file use external application ( graphic viewer)
 ---------------------------------
 I think that creating graphic can't be simpler

Adam Majewski
fraktal.republika.pl

gcc b.c -Wall -lm
 
 */
 #include <stdio.h>
#include <stdlib.h> // malloc
#include <math.h>
#include <complex.h>
#include <string.h> // strcat
 
  const int iXmax = 1000; 
 const int iYmax = 1000;

 /* world ( double) coordinate = parameter plane*/
 const double xMin=-1.0;
 const double xMax=1.4;
 const double yMin=-1.2;
 const double yMax=1.2;
 /* */
 double PixelWidth; //=(xMax-xMin)/iXmax;
 double PixelHeight; //=(yMax-yMin)/iYmax;

 
 /*  */
 
 const int IterationMax=2000;

  
 /* bail-out value , radius of circle ;  */
 const int EscapeRadius=2;
 int ER2; //=EscapeRadius*EscapeRadius;

 int iSize ; //= iXmax*iYmax; //





// target set for points falling into alfa fixed point
// is a circle around alfa fixed point
// with radius = AR
double AR ; // radius of target set around alfa fixed point in world coordinate AR = PixelWidth*TargetWidth; 
double AR2; // =AR*AR;
double TargetRadiusInPixels; // radius of target set in pixels ; 

// memmory 1D arrays 
 unsigned char *data;
 unsigned char *edge;
 
 
//
#define  p 6 // number of petals= denominator of angle in turn 
// petals are from angle 0 to 1/p
// 1/p to 2/p
// and so on 
// if you will increase subdivision to angles 1/(n*p) then !!!
#define  nMax 21 // number of branches of star
unsigned char Colors[(nMax*p)];
unsigned char iExterior = 245;
// angles 
double TwoPi=2*M_PI;


/* gives position of 2D point (iX,iY) in 1D array  ; uses also global variable iWidth */
unsigned int Give_i(unsigned int ix, unsigned int iy)
{ return ix + iy*iYmax; }


 

// colors of components interior = shades of gray
// uses global vars Colors
int InitColors(int iMax)
{
  int i;
  // int iMax = p*nMax; 
  unsigned int iStep;
 
  iStep=150/iMax;
 
  for (i = 1; i <= iMax; ++i)
    {Colors[i-1] = iExterior -i*iStep; 
     //printf("i= %d color = %i  \n",i-1, Colors[i-1]); // info
     }
  return 0;
}

double GiveTurn(double complex z)
{
  double argument;
 
  argument = carg(z); //   argument in radians from -pi to pi
  if (argument<0) argument=argument + TwoPi; //   argument in radians from 0 to 2*pi
  return argument/TwoPi ; // argument in turns from 0.0 to 1.0
}




unsigned char GiveColorOfInterior(double x, double y, int iMax)
{
 
double angle;
unsigned char sector=0;
int i ; // number of the sector

  angle=GiveTurn(x+y*I);
   for (i = 0; i < iMax; ++i)
      if (angle<(i+1.0)/iMax) {sector=i; return Colors[sector] ;} 
  return Colors[sector];

}



// fill array 
// uses global var :  ...
// scanning complex plane 
int FillArray(unsigned char data[], int b )
{
  unsigned int iX, iY; // pixel coordinate /* screen ( integer) coordinate */
  unsigned char color;
  double x, y, tempx,    /* Z=x+y*i   */
  x0, y0,  /* Z0 = x0 + y0*i */
  x2, y2; /* x2=x*x;  y2=y*y  */
  double x2y2; //, y4, x4, x2my2;
  int Iteration;
  int i; // index of array
 
   /* compute and write image data bytes to the file*/
 for(iY=0;iY<iYmax;++iY)
    {
     y0=yMax - iY*PixelHeight; /* reverse Y  axis */
     if (fabs(y0)<PixelHeight/2) y0=0.0; /*  */  
     printf(" iy = %d from %d\r", iY, iYmax); //info   
     for(iX=0;iX<iXmax;++iX)
       {    /* initial value of orbit Z0 */
        x0=xMin + iX*PixelWidth;
        /* Z = Z0 */
       x=x0;
       y=y0;
            
       x2=x*x;
       y2=y*y;
       x2y2= x2+y2;
       /* */
       for (Iteration=0;Iteration<IterationMax ;Iteration++)
          { 
             if ((x2y2)>ER2) { color=iExterior; break;}/* exterior of Filled-in Julia set   */
             if ((x2y2)<AR2) { color=GiveColorOfInterior(x,y,b); break;}/* interior of Filled-in Julia set   */
             // z^4- z
             tempx = y*y*y*y -6*y2*x2 + x*x*x*x - x ; // y^4-6*x^2*y^2+x^4-x
             y = -4*x*y*y*y +4*x*x*x*y - y ; // -4*x*y^3+4*x^3*y-y
	     x=tempx;
             x2=x*x;
             y2=y*y;
	     x2y2= x2+y2;
             
             
          };

     if (Iteration==IterationMax) color=GiveColorOfInterior(x,y,b); ;
     i=Give_i(iX,iY);
     data[i]=color;
   //printf("color = %d\n", color);
                         
    
  
  }
 }
 
  return 0;
}


int AddBoundaries(unsigned char data[])
{
 
  unsigned int iX,iY; /* indices of 2D virtual array (image) = integer coordinate */
  unsigned int i; /* index of 1D array  */
  /* sobel filter */
  unsigned char G, Gh, Gv; 
 
 
 
 
 // printf(" find boundaries in data array using  Sobel filter\n");   

  for(iY=1;iY<iYmax-1;++iY){ 
    for(iX=1;iX<iXmax-1;++iX){ 
      Gv= data[Give_i(iX-1,iY+1)] + 2*data[Give_i(iX,iY+1)] + data[Give_i(iX-1,iY+1)] - data[Give_i(iX-1,iY-1)] - 2*data[Give_i(iX-1,iY)] - data[Give_i(iX+1,iY-1)];
      Gh= data[Give_i(iX+1,iY+1)] + 2*data[Give_i(iX+1,iY)] + data[Give_i(iX-1,iY-1)] - data[Give_i(iX+1,iY-1)] - 2*data[Give_i(iX-1,iY)] - data[Give_i(iX-1,iY-1)];
      G = sqrt(Gh*Gh + Gv*Gv);
      i= Give_i(iX,iY); /* compute index of 1D array from indices of 2D array */
      if (G==0) {edge[i]=255;} /* background */
      else {edge[i]=0;}  /* boundary */
    }
  }
 
  // copy boundaries from edge array to data array 
  //for(iY=1;iY<iyMax-1;++iY){ 
  //   for(iX=1;iX<ixMax-1;++iX){i= Give_i(iX,iY); if (edge[i]==0) data[i]=0;}}
 
 
 
  return 0;
}


// save data array to pgm file 
int SaveArray2PGMFile( unsigned char data[], int t)
{
 
  FILE * fp;
  const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
  char name [10]; /* name of file */
  sprintf(name,"%d", t); /*  */
  char *filename =strcat(name,".pgm");
  char *comment="# ";/* comment should start with # */
 
  /* save image to the pgm file  */      
  fp= fopen(filename,"wb"); /*create new file,give it a name and open it in binary mode  */
  fprintf(fp,"P5\n %s\n %u %u\n %u\n",comment,iXmax,iYmax,MaxColorComponentValue);  /*write header to the file*/
  fwrite(data,iSize,1,fp);  /*write image data bytes to the file in one step */
  printf("File %s saved. \n", filename);
  fclose(fp);
 
  return 0;
}
 
/* ------------------------------- main ----------------------------- */

 int main()
 {
 
  int n;

  //int iX,iY;
 PixelWidth=(xMax-xMin)/iXmax;
 PixelHeight=(yMax-yMin)/iYmax;
  
  iSize = iXmax*iYmax; //
ER2=EscapeRadius*EscapeRadius;
 TargetRadiusInPixels = iXmax/50.0; // = 15.0; // radius of target set in pixels ; Maybe increase to 20 = 1000/50
 AR = PixelWidth*TargetRadiusInPixels; // !!!! important value  ( and precision and time of creation of the pgm image ) 
 AR2= AR*AR;

/* create dynamic 1D arrays for colors ( shades of gray ) */
  data = malloc( iSize * sizeof(unsigned char) );
  edge = malloc( iSize * sizeof(unsigned char) );
  if (edge == NULL || edge == NULL )
    {
      fprintf(stderr," Could not allocate memory\n");
      return 1;
    }
  else fprintf(stderr," memory is OK \n");


  



for(n=1; n<nMax; ++n){ 
  InitColors(n*p);
  FillArray( data , n*p); // no symmetry
  //SaveArray2PGMFile(data,1);
  AddBoundaries(data);
  SaveArray2PGMFile(edge,n);
}

//
  free(data);
  free(edge);

 return 0;
 }