/*

  c console program
  It can be compiled and run under Linux, windows, Mac 
  It needs gcc

 
  immediate basin of attraction 
  of an attracting cycle  {z1,z2,...,zp}
  of period p
  = subset of components Bi in the Fatou set  
  such that zi belongs to Bi where i is in <1,p>

  1. find :
  - period
  - an attracting cycle
  2. Compute exterior ( interior is it's 
  3. and its boundary
  4. fill boundary of component containing zi
  5. mark an attracting cycle

  -----------------------------------------
  1.pgm file code is  based on the code of Claudio Rocchini
  http://en.wikipedia.org/wiki/Image:Color_complex_plot.jpg
  create 8 bit color graphic file ,  portable gray map file = pgm 
  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 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 

  to compile : 
  gcc j.c -lm -Wall
  to run ( Linux console) :
  ./a.out

 
*/
# include <stdio.h>
# include <stdlib.h>
# include <math.h>
# include <string.h>

/* iXmax/iYmax = 1 */
# define iXmax 2000 /* height of image in pixels */
# define iYmax 2000
# define iLength (iXmax*iYmax)
/* */
# define ZxMin -1.5
# define ZxMax 1.5
# define ZyMin -1.5
# define ZyMax 1.5
/* (ZxMax-ZxMin)/(ZyMax-ZyMin)==iXmax/iYmax = 1 */

# define IterationMax (iXmax/10) /* if to low then components will be not properly bordered */

# define PixelWidth  ((ZxMax-ZxMin)/iXmax)
# define PixelHeight ((ZyMax-ZyMin)/iYmax)

/* fc(z) = z*z + c */
/* internal angle of main cardioid , rotation number  = 1/3 */
# define Cx -0.12561 /* C = Cx + Cy*i center of period 3 component  */
# define Cy  0.744862

# define period 3 /* period of attracting cycle */

# define EscapeRadius 80.0 /* radius of circle around origin; its complement is a target set for escaping points */
# define ER2 (EscapeRadius*EscapeRadius)

/* colors */
# define iExterior 245 /* exterior of Julia set */
# define iJulia 0 /* border , boundary*/
# define iInterior 230

/* escape time to infinity of function fc(z) = z*z + c */
int GiveExtLastIteration(double _Zx0, double _Zy0,double C_x, double C_y, int iMax, double _ER2)
{ 
  int i; /* iteration */
  double Zx, Zy; /* Z = Zx + Zy*i */
  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; /* last iteration */
}

unsigned int f(unsigned int _iX, unsigned int _iY)
/* 
   gives position of point (iX,iY) in 1D array  ; uses also global variables 
   it does not check if index is good  so memory error is possible 
*/
{return (_iX + (iYmax-_iY-1)*iXmax );}

int FillContour(double dXseed, double dYseed,  unsigned char color, unsigned char _data[])
{ 
  /* 
     fills contour with black border ( color = iJulia)  using seed point inside contour 
     and horizontal lines 
     it starts from seed point, saves max right( iXmaxLocal) and max left ( iXminLocal) interior points of horizontal line,
     in new line ( iY+1 or iY-1) it computes new interior point  : iXmidLocal=iXminLocal + (iXmaxLocal-iXminLocal)/2;
     result is stored in _data array : 1D array of 1-bit colors ( shades of gray)
     it does not check if index of _data array is good  so memory error is possible 
  */

  int iXseed = (int)((dXseed - ZxMin)/PixelWidth);
  int iYseed = (int)((dYseed - ZyMin)/PixelHeight);
  
  int iX, /* seed integer coordinate */
    iY=iYseed,
    /* most interior point of line iY */
    iXmidLocal=iXseed, 
    /* min and max of interior points of horizontal line iY */
    iXminLocal, 
    iXmaxLocal; 
  int i ; /* index of _data array */;

  
  /* ---------  move up --------------- */ 
  do{
    iX=iXmidLocal;
    i =f(iX,iY); /* index of _data array */;
  
    /* move to right */
    while (_data[i]==iInterior) 
      { _data[i]=color;
	iX+=1; 
	i=f(iX,iY);  
      }
    iXmaxLocal=iX-1;

    /* move to left */
    iX=iXmidLocal-1; 
    i=f(iX,iY);
    while (_data[i]==iInterior) 
      { _data[i]=color;
	iX-=1; 
	i=f(iX,iY); 
      }
    iXminLocal=iX+1; 

    iY+=1; /* move up */
    iXmidLocal=iXminLocal + (iXmaxLocal-iXminLocal)/2; /* new iX inside contour */
    i=f(iXmidLocal,iY); /* index of _data array */;
    if ( _data[i]==iJulia)  break; /*  it should not cross the border */
 
  } while  (iY<iYmax); 
  
  
  /* ------  move down ----------------- */
  iXmidLocal=iXseed;
  iY=iYseed-1;
  
  
  do{
    iX=iXmidLocal;
    i =f(iX,iY); /* index of _data array */;
  
    /* move to right */
    while (_data[i]==iInterior) /*  */
      { _data[i]=color;
	iX+=1;
	i=f(iX,iY);  
      }
    iXmaxLocal=iX-1;

    /* move to left */
    iX=iXmidLocal-1; 
    i=f(iX,iY);
    while (_data[i]==iInterior) /*  */
      { _data[i]=color;
	iX-=1; /* move to right */
	i=f(iX,iY);  
      }
    iXminLocal=iX+1; 
  
    iY-=1; /* move down */
    iXmidLocal=iXminLocal + (iXmaxLocal-iXminLocal)/2; /* new iX inside contour */
    i=f(iXmidLocal,iY); /* index of _data array */;
    if ( _data[i]==iJulia)  break; /*  it should not cross the border */
  } while  (0<iY); 

  /* mark seed point by big pixel */
  const int iSide =iXmax/500; /* half of width or height of big pixel */
  for(iY=iYseed-iSide;iY<=iYseed+iSide;++iY){ 
    for(iX=iXseed-iSide;iX<=iXseed+iSide;++iX){ 
      i= f(iX,iY); /* index of _data array */
      _data[i]=10;}}

  return 0;
}

/* --------------------------------------------------------------------------------------------------------- */

int main(){

  
  /* attracting cycle : Z1=0, Z2=f(0)=c, z3=f(f(0))=c*c + c */
  const double Z1x = 0.0,
    Z1y = 0.0,
    Z2x = Cx,
    Z2y = Cy,
    Z3x = (Cx*Cx)-(Cy*Cy) +Cx,
    Z3y = 2*Cx*Cy + Cy;

  unsigned int iX,iY, /* indices of 2D virtual array (image) = integer coordinate */
    i; /* index of 1D array  */
   
  
  /* */
  double Zx,Zy;
  
    
  int LastIteration;
  /* sobel filter */
  unsigned char G, Gh, Gv; 
  /* color */
  //unsigned char ColorList[]={255,230,180};
  const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
  

  /* 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");
      return 1;
    }
  else printf(" memory is OK\n");

   
 
  
  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);    
    for(iX=0;iX<iXmax;++iX){ 
      Zx=ZxMin + iX*PixelWidth;
      LastIteration = GiveExtLastIteration(Zx, Zy, Cx, Cy, IterationMax, ER2 );
      i= f(iX,iY); /* compute index of 1D array from indices of 2D array */
      if ( IterationMax != LastIteration ) 
	{data[i]=iExterior;} /* exterior */
      else { data[i]=iInterior;}/* interior */
	 
      /* 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]=iJulia;}  /* 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]==iJulia) data[i]=iJulia;}}

  printf(" mark immediate basin of attracting cycle \n");
  FillContour(Z1x, Z1y, 110, data);
  FillContour(Z2x, Z2y, 140, data);
  FillContour(Z3x, Z3y, 170, data);
   

  /* ---------- file  -------------------------------------*/
  printf(" save  data array to the pgm file \n");
  FILE * fp;
  char name [10]; /* name of file */
  i = sprintf(name,"jPeriod%uiXmax%u",period,iXmax); /* result (is saved in i) but is not used */
  char *filename =strcat(name,".pgm");
  char *comment="# C= ";/* 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\n %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);

  /* --------------free memory ---------------------*/
  free(data);
  free(edge);
  
  

  return 0;
}