/* 
 c program:
  1. draws Julia set for z^4-z
  using  boolean escape time )
 -------------------------------         
 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


(%i2) z:x+y*%i;
(%o2) %i*y+x
(%i3) realpart(z^4-z);
(%o3) y^4-6*x^2*y^2+x^4-x
(%i4) imagpart(z^4-z);
(%o4) -4*x*y^3+4*x^3*y-y



gcc z.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;
 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;
 int i; // index of array
double TwoPi=2.0*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; }


 
// save data array to pgm file 
int SaveArray2PGMFile( unsigned char data[], double 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,"%f", 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;
}

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)
{
 double angle;
 //double offset = 1.0/(2.0*6.0);


// all points tend to z=0 thru 6 petals 
// check to which sector / petal fall 
  angle=GiveTurn(x+y*I); 
  if (angle<1.0/12.0 || angle > (11.0/12.0)) return 235; // (11/12 ; 1/12)
  if (angle<3.0/12.0)     return 225; // 1/12 ; 3/12
  if (angle<5.0/12.0) return 215; // 3/12 ; 5/12
  if (angle<7.0/12.0) return 205; // 5/12 ; 7/12
  if (angle<9.0/12.0) return 195; // 7/12 ; 9/12
//if (angle<11.0/12.0) 
     return 175; // 9/12 ; 11/12
  
}

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]=245;} /* 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;
}
 
int CopyBoundaries()
{
 
  unsigned int iX,iY; /* indices of 2D virtual array (image) = integer coordinate */
  unsigned int i; /* index of 1D array  */
      

 printf("copy boundaries from edge array to data array \n");
 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;
}

 int main()
 {
 
  /* screen ( integer) coordinate */
  int iX,iY;
 
  
  iSize = iXmax*iYmax; //
 /* 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;
 //
 /* color component ( R or G or B) is coded from 0 to 255 */
 /* it is 24 bit color RGB file */
 
 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;
 const int IterationMax=2000;

  
 /* bail-out value , radius of circle ;  */
 const int EscapeRadius=2;
 int ER2=EscapeRadius*EscapeRadius;
 TargetRadiusInPixels = iXmax/100.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");


 

 /* 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=255; break;}/* exterior of Filled-in Julia set   */
             if ((x2y2)<AR2) { color=GiveColorOfInterior(x,y); 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); ;
     i=Give_i(iX,iY);
     data[i]=color;
                         
    
  
  }
 }
 SaveArray2PGMFile(data,0.101);
AddBoundaries(data);
SaveArray2PGMFile(edge,0.102);
CopyBoundaries();
SaveArray2PGMFile(data,0.103);

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

 return 0;
 }