/*
  c console program
  -----------------------------------------
  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 color values of pixels,
  fills tha array with data and after that writes the data from array to pgm file.
  It allows free ( non sequential) acces 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 
 
*/
#include <stdio.h>
#include <math.h>
#include <complex.h>

/* 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;
}

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;
}

/* attracting time to finite attractor ZA */
int GiveIntLastIteration(double _Zx0, double _Zy0,double C_x, double C_y, int iMax, double _AR2, double _ZAx, double _ZAy )
{ 
  int i;
  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 + C_y;
      Zx=Zx2-Zy2 +C_x;
      Zx2=Zx*Zx;
      Zy2=Zy*Zy;
      dX=Zx-_ZAx;
      dY=Zy-_ZAy;
      d=dX*dX+dY*dY;
    };
  return i;
}

int main(){
  const double Cx=-0.12565651,
    Cy=0.65720;
  unsigned char period=3; /* period of finite attractor */
  const int iXmax = 1000; 
  const int iYmax = 1000;
  int iX,iY,
    iWidth = iXmax + 1,
    iHeight= iYmax + 1;
  /* world ( double) coordinate = parameter plane*/
  const double ZxMin=-1.4;
  const double ZxMax=1.4;
  const double ZyMin=-1.4;
  const double ZyMax=1.4;
  double Zx, Zy;    /* Z=Zx+Zy*i   */
  double complex ZA;  /* atractor ZA = ZAx + ZAy*i */
  /* */
  double PixelWidth=(ZxMax-ZxMin)/iWidth;
  double PixelHeight=(ZyMax-ZyMin)/iHeight;
  const double EscapeRadius=80.0; /* radius of circle around origin; its complement is a target set for escaping points */
  double ER2=EscapeRadius*EscapeRadius;
  const double AR = PixelWidth/10; /* radius of circle around attractor ZA = target set for attracting points */
  double       AR2 = AR* AR;
  const int IterationMax=60,
    IterationMaxBig= 1000001;
  int eLastIteration, iLastIteration;
  unsigned char G, Gh, Gv; /* sobel filter */
  unsigned char data[iYmax][iXmax]; /* 2D array for colors ( shades of gray ) */
  unsigned char edge[iYmax][iXmax]; /* 2D array for colors ( shades of gray ) */
  unsigned char color[]={255,230,180};
  const int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
  FILE * fp;
  char *filename ="ct02.pgm";
  char *comment="# this is binary pgm  file";/* comment should start with # */
 
 ZA = GiveAttractor( Cx, Cy, ER2, IterationMaxBig); /* find attractor ZA  using forward iteration of critical point Z = 0  */
 /* fill the data array */
  for(iY=0;iY<iYmax;++iY){ 
    Zy=ZyMax - iY*PixelHeight; /* reverse Y  axis */
    if (fabs(Zy)<PixelHeight/2) Zy=0.0; /*  */    
    for(iX=0;iX<iXmax;++iX){ 
      Zx=ZxMin + iX*PixelWidth;
      eLastIteration = GiveExtLastIteration(Zx, Zy, Cx, Cy, IterationMax, ER2 );
      if ( IterationMax != eLastIteration ) 
	{data[iY][iX]=245;} /* exterior */
      else /* interior */
	{ iLastIteration = GiveIntLastIteration(Zx, Zy, Cx, Cy, IterationMaxBig, AR2, creal(ZA), cimag(ZA));
          data[iY][iX]=color[iLastIteration % period];} /* 255 - 30*(iLastIteration % period);}*/
      /* if (Zx>0 && Zy>0) data[iY][iX]=255-data[iY][iX];    check the orientation of Z-plane by marking first quadrant */
   }
 }
  /* find boundaries - Sobel filter
    for points not the borders of data array 
 */
  for(iY=1;iY<iYmax-1;++iY){ 
    for(iX=1;iX<iXmax-1;++iX){ 
      Gv= data[iY+1][iX-1] + 2*data[iY+1][iX] + data[iY+1][iX-1] - data[iY-1][iX-1] - 2*data[iY][iX-1] - data[iY-1][iX+1];
      Gh= data[iY+1][iX+1] + 2*data[iY][iX+1] + data[iY-1][iX-1] - data[iY-1][iX+1] - 2*data[iY][iX-1] - data[iY-1][iX-1];
      G = sqrt(Gh*Gh + Gv*Gv);
      if (G==0) {edge[iY][iX]=255;} /* background */
      else {edge[iY][iX]=0;}  /* boundary */
    }
  }

  for(iY=1;iY<iYmax-1;++iY){ 
    for(iX=1;iX<iXmax-1;++iX){ if (edge[iY][iX]==0) data[iY][iX]=0;}}

  /* write the whole data array to ppm file in one step */      
  fp= fopen(filename,"wb"); /*create new file,give it a name and open it in binary mode  */
  fprintf(fp,"P5\n %s\n %d\n %d\n %d\n",comment,iXmax,iYmax,MaxColorComponentValue);  /*write header to the file*/
  fwrite(data,sizeof data,1,fp);  /*write image data bytes to the file*/
  fclose(fp);
  printf("OK - file %s saved\n", filename);
  return 0;
}