/*
 
  c console program
  
  https://commons.wikimedia.org/wiki/File:Quadratic_Julia_set_with_Internal_tile_for_internal_ray_0.ogv
  -----------------------------------------
  1.ppm file code 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
  ---------------------------
  2. first it creates data array which is used to store rgb color values of pixels,
  fills tha array with data and after that writes the data from array to pgm file.
  It alows free ( non sequential) access to "pixels"
 
  -------------------------------------------
  Adam Majewski   fraktal.republika.pl 
 
  Sobel filter 
  Gh = sum of six values ( 3 values of matrix are equal to 0 ). Each value is = pixel_color * filter_coefficients 
 
 
 
 
  gcc i.c -lm -Wall -o2
  gcc i.c -lm -Wall -march=native
  time ./a.out

 
 
 
*/
# include <stdio.h>
# include <stdlib.h>
# include <math.h>
# include <complex.h>
# include <string.h>
 
/* iXmax/iYmax = 1 */
unsigned int iXmax = 20000; /* height of image in pixels */
unsigned int iYmax = 20000;
unsigned int iLength; 
/* fc(z) = z*z + c */
# define denominator 1 /* denominator of internal angle */

double  AR = 0.0014998955;  /* PixelWidth*1.5   radius of circle around attractor ZA = target set for attracting points */

//#define alfa (1-sqrt(1-4*Cx))/2 /* attracting or parabolic fixed point z = alfa */
//#define beta (1+sqrt(1-4*Cx))/2 /* repelling or parabolic fixed point z = beta */

/* color */
unsigned char color;
// Arrays are 0-indexed, so the first array element is at index = 0, and the highest is =(size_of_array – 1) 
  unsigned char colorArray[2][2]={{255,231},
				   {123,99}}; /* shades of gray used in image */
const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */

 
/* escape time to infinity */
int GiveExtLastIteration(double _Zx0, double _Zy0,double C_x, double C_y, int iMax, double _ER2)
{ 
  int i;
  double Zx, Zy;
  double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
  Zx=_Zx0; /* initial value of orbit  */
  Zy=_Zy0;
  Zx2=Zx*Zx;
  Zy2=Zy*Zy;
  for (i=0;i<iMax && ((Zx2+Zy2)<_ER2);i++)
    {
      Zy=2*Zx*Zy + C_y;
      Zx=Zx2-Zy2 +C_x;
      Zx2=Zx*Zx;
      Zy2=Zy*Zy;
    };
  return i;
}
 
 
/* find attractor ZA  using forward iteration of critical point Z = 0  */
/* if period is >1 gives one point from attracting cycle */
double complex GiveAttractor(double _Cx, double _Cy, double ER2, int _IterationMax)
{
  int Iteration;
  double Zx, Zy; /* z = zx+zy*i */
  double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
  /* -- find attractor ZA  using forward iteration of critical point Z = 0  */
  Zx=0.0;
  Zy=0.0;
  Zx2=Zx*Zx;
  Zy2=Zy*Zy;
  for (Iteration=0;Iteration<_IterationMax && ((Zx2+Zy2)<ER2);Iteration++)
    {
      Zy=2*Zx*Zy + _Cy;
      Zx=Zx2-Zy2 + _Cx;
      Zx2=Zx*Zx;
      Zy2=Zy*Zy;
    };
  return Zx+Zy*I;
}

// color is related to 2 measures : iLastIteration and part of internal tile 
unsigned char GiveIntColor(double Zx0, double Zy0, double Cx, double Cy, int iMax, double AR2, double ZAx, double ZAy, unsigned char colorArray[2][2])
{  int i, m, n;
  double Zx, Zy; /* z = zx+zy*i */
  double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
  double d, dX, dY; /* distance from z to Alpha  */
  Zx= Zx0; /* initial value of orbit  */
  Zy= Zy0;
  Zx2=Zx*Zx;
  Zy2=Zy*Zy;
  dX=Zx- ZAx;
  dY=Zy- ZAy;
  d=dX*dX+dY*dY;
  for (i=0;i<iMax && (d> AR2);i++)
    {
      Zy=2*Zx*Zy + Cy;
      Zx=Zx2-Zy2 +Cx;
      Zx2=Zx*Zx;
      Zy2=Zy*Zy;
      dX=Zx- ZAx;
      dY=Zy- ZAy;
      d=dX*dX+dY*dY;
    }
   m = (Zy > 0 ? 0 : 1); // petal part
   n = (i % 2); // attraction time 
   return colorArray[m][n]; //iColor
}

 
/* gives position of point (iX,iY) in 1D array  ; uses also global variables */
unsigned int f(unsigned int _iX, unsigned int _iY)
{return (_iX + (iYmax-_iY-1)*iXmax );}

// save data array to pgm file 
int SavePGMFile(int iWidth, unsigned char data[])
{
FILE * fp;
  char name [15]; /* name of file */
  sprintf(name,"%d", iWidth); /*  basename = file name without extension*/
  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,iLength,1,fp);  /*write image data bytes to the file in one step */
  printf("File %s saved. \n", filename);
  fclose(fp);
  return 0;
}

 
/* --------------------------------------------------------------------------------------------------------- */
 
int main(){
 
 unsigned int nMax ; /* number of steps = number of images */
 unsigned int n;
//double CxMin =  0.0837;
//double CxMax =  0.0839; /* C = Cx + Cy*i */
double carray[]={
0.0, 0.01,  0.02, 0.03, 0.035, 0.03722, 0.04, 0.045, 0.05, 0.06, 0.07, 0.077, 0.08378,  0.09, 0.10,
0.11, 0.12214, 0.13, 0.14, 0.150656, 0.155, 0.16, 0.165, 0.171528, 0.175, 0.18, 0.183, 0.186948,  0.19, 0.193, 0.197,
0.207388, 0.203, 0.206, 0.209, 0.2142728, 0.215, 0.2165, 0.2175, 0.21971144, 0.221, 0.222, 0.223, 0.22406680, 0.226,
0.228, 0.23050176, 0.232, 0.2335, 0.23492, 0.2353, 0.2356, 0.2366192000, 0.24, 0.245, 0.25};

nMax=sizeof(carray)/sizeof(double);

double Cx;
//double stepCx;
double Cy =  0.0;
// stepCx = (CxMax - CxMin)/ nMax;
 
  unsigned int iX,iY, /* indices of 2D virtual array (image) = integer coordinate */
    i; /* index of 1D array  */
    iLength = iXmax*iYmax;/* length of array in bytes = number of bytes = number of pixels of image * number of bytes of color */
  /* world ( double) coordinate = dynamic plane = z-plane */
  const double dSide = 1.3;
  const double ZxMin=-dSide;
  const double ZxMax=dSide;
  const double ZyMin=-dSide;
  const double ZyMax=dSide;
  double PixelWidth=(ZxMax-ZxMin)/iXmax;
  double PixelHeight=(ZyMax-ZyMin)/iYmax;
  /* */
  double Zx, Zy;    /* Z=Zx+Zy*i   */
 // double alfa; // define alfa (1-sqrt(1-4*Cx))/2 /* attracting or parabolic fixed point z = alfa */
  double complex ZA;  /* atractor ZA = ZAx + ZAy*i */
  double AR2 =  AR*AR;
  /* */
 
  const double EscapeRadius=2.0; /* radius of circle around origin; its complement is a target set for escaping points */
  double ER2=EscapeRadius*EscapeRadius;
 
  const int IterationMax=60,
    IterationMaxBig= 1000001;
  int eLastIteration; // iLastIteration;
  //int InternalTile;
  /* sobel filter */
  unsigned char G, Gh, Gv; 
  
 
 
  /* dynamic 1D arrays for colors ( shades of gray ) */
  unsigned char *data, *edge;
  data = malloc( iLength * sizeof(unsigned char) );
  edge = malloc( iLength * sizeof(unsigned char) );
  if (data == NULL || edge==NULL)
    {
      fprintf(stderr," Could not allocate memory");
      getchar(); 
      return 1;
    }
   
 
  //for(n=0;n<=nMax;++n)
  // {

   Cx = 0.25; //carray[n];
   //alfa = (1-sqrt(1-4*Cx))/2 ; /* attracting or parabolic fixed point z = alfa */
 
  ZA = GiveAttractor( Cx, Cy, ER2, IterationMaxBig); /* find attractor ZA  using forward iteration of critical point Z = 0  */
 
 
 
 // printf(" fill the data array \n");
  for(iY=0;iY<iYmax;++iY){ 
    Zy=ZyMin + iY*PixelHeight; /*  */
    if (fabs(Zy)<PixelHeight/2) Zy=0.0; /*  */
   //printf(" row %u from %u \n",iY, iYmax); /* info */   
    for(iX=0;iX<iXmax;++iX){ 
      Zx=ZxMin + iX*PixelWidth;
      eLastIteration = GiveExtLastIteration(Zx, Zy, Cx, Cy, IterationMax, ER2 );
      i= f(iX,iY); /* compute index of 1D array from indices of 2D array */
       if ( IterationMax != eLastIteration ) 
        color = 245; /* exterior */
        /* interior */
        else color = GiveIntColor(Zx, Zy, Cx,  Cy, IterationMaxBig, AR2, creal(ZA), cimag(ZA), colorArray); /*   */
     data[i]=color;

        //if (Zx>=0 && Zx <= 0.5 && (Zy > 0 ? Zy : -Zy) <= 0.5 - Zx) data[i]=255-data[i]; // show petal

         
      /*  if (Zx>0 && Zy>0) data[i]=255-data[i];    check the orientation of Z-plane by marking first quadrant */

    }
  }
 
 
 // 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[f(iX-1,iY+1)] + 2*data[f(iX,iY+1)] + data[f(iX-1,iY+1)] - data[f(iX-1,iY-1)] - 2*data[f(iX-1,iY)] - data[f(iX+1,iY-1)];
      Gh= data[f(iX+1,iY+1)] + 2*data[f(iX+1,iY)] + data[f(iX-1,iY-1)] - data[f(iX+1,iY-1)] - 2*data[f(iX-1,iY)] - data[f(iX-1,iY-1)];
      G = sqrt(Gh*Gh + Gv*Gv);
      i= f(iX,iY); /* compute index of 1D array from indices of 2D array */
      if (G==0) {edge[i]=255;} /* background */
      else {edge[i]=0;}  /* boundary */
    }
  }
 
    //printf(" copy boundaries from edge to data array \n");
    for(iY=1;iY<iYmax-1;++iY){ 
     for(iX=1;iX<iXmax-1;++iX)
      {i= f(iX,iY); /* compute index of 1D array from indices of 2D array */
    if (edge[i]==0) data[i]=0;}}
 
 
  /* ---------- file  -------------------------------------*/
  //printf(" save  data array to the file \n");
  SavePGMFile( iXmax, data);

 // } // for n ....
 
  /* --------------free memory ---------------------*/
  free(data);
  free(edge);
 
 
 
  return 0;
}