/*

  Adam Majewski
  adammaj1 aaattt o2 dot pl  // o like oxygen not 0 like zero 
  fraktal.republika.pl

  c console progam 

  gcc c.c -lm -Wall -march=native 
  time ./a.out


  gcc c.c -lm -Wall -march=native -fopenmp

  method of filling gaps in rays turns 
  should be the same as for drawing rays !!!!!!!!!!!!!

  !!!!!---------- algorithm ---------!!!!!!!!
  trap ( target set) for iteration :
  * of interior points = interior of circle centered at alfa and radius = dMaxDistance2Alfa, with iPeriodChild sectors described by periodic external rays landing on alfa
  * of exterior points = exterior of circle centered at origin ( z=0) and radius = ER
  Iterate point until it fails to one of 2 traps.
  If it it fail into interior trap then check in which sector is it. Color = number of sector

  Interior target set should be all inside Julia set 
  ( of course not counting external rays that cross it and exterior parts which are very thin = means width < pixelwidth)
  if it is a circle around alfa then it shrinks as iPeriodOfChild grows
  ( = time of creating image grows)
  -----------------------------------------
  Better is to take as a interior target set :
  part of this circle bordered by 2 external rays 
  (One can choose the longest one !!!)
  and color the interior proportionaly to 
  (i % iPeriodOfChild)


  gcc c.c -lm -Wall -march=native

  ./a.out



  without openmp
  real	7m54.710s
  user	7m54.861s
  sys	0m0.008s

  with openmp : 
  real	1m47.034s
  user	14m8.420s
  sys	0m0.280s







*/

#include <stdio.h>
#include <stdlib.h> // malloc
#include <string.h> // strcat
#include <math.h> // M_PI; needs -lm also 
#include <complex.h>
#include <omp.h>

/* --------------------------------- global variables and consts ------------------------------------------------------------ */

#define iPeriodChild 4 // Period of secondary component joined by root point with the parent component 
int iPeriodParent = 2; // main cardioid of Mandelbrot set



// virtual 2D array and integer ( screen) coordinate
// Indexes of array starts from 0 not 1 
//unsigned int ix, iy; // var
static unsigned int ixMin = 0; // Indexes of array starts from 0 not 1
static unsigned int ixMax ; //
static unsigned int iWidth ; // horizontal dimension of array

static unsigned int iyMin = 0; // Indexes of array starts from 0 not 1
static unsigned int iyMax ; //

static unsigned int iHeight = 5000; //  
// The size of array has to be a positive constant integer 
static unsigned int iSize ; // = iWidth*iHeight; 

// memmory 1D array 

unsigned char *data; // fatou components
unsigned char *julia; // julia set = common edge between fatou components
unsigned char *edge; //  pther edges 
unsigned char *zero; // zero algorithm

// unsigned int i; // var = index of 1D array
//static unsigned int iMin = 0; // Indexes of array starts from 0 not 1
static unsigned int iMax ; // = i2Dsize-1  = 
// The size of array has to be a positive constant integer 
// unsigned int i1Dsize ; // = i2Dsize  = (iMax -iMin + 1) =  ;  1D array with the same size as 2D array

/* world ( double) coordinate = dynamic plane */
static   const double ZxMin=-1.8;
static  const double ZxMax=1.8;
static  const double ZyMin=-0.9;
static  const double ZyMax=0.9;
static  double PixelWidth; // =(ZxMax-ZxMin)/ixMax;
static  double PixelHeight; // =(ZyMax-ZyMin)/iyMax;
static  double ratio ;
 

// complex numbers of parametr plane 
double Cx; // c =Cx +Cy * i
double Cy;
double complex c; // parameter of function fc(z)=z^2 + c



static unsigned long int iterMax  = 10000; //iHeight*100;

static double ER = 2.0; // Escape Radius for bailout test 
static double ER2;

/* colors = shades of gray from 0 to 255 */
// 8 bit color = int number from 0 to 255


// Arrays are 0-indexed, so the first array element is at index = 0, and the highest is =(size_of_array – 1) 
//unsigned char iColorInterior[2][iPeriodChild]={{255,231}, {123,99}}; /* shades of gray used in image */
unsigned char iColors[4]={240,200,160, 120}; //
static unsigned char iColorOfExterior = 245;
// static unsigned char iColorOfInterior = 200;
unsigned char iColorOfUnknown = 100;
int iNumberOfUknknown = 0;

// period 2 cycle = {Zx2a, Zx2b} , Zy = 0 
double Zx2a = -1.207106781186548;
double Zx2b = 0.207106781186547;

// limits : Zx2aMinus < Zx2a < Zx2aPlus < Zx2bMinus < Zx2b < Zx2bPlus
double Zx2aMinus; 
double Zx2bMinus;
double Zx2aPlus; 
double Zx2bPlus;


//static double TwoPi=2*M_PI;


/* ------------------------------------------ functions -------------------------------------------------------------*/





//------------------complex numbers -----------------------------------------------------


// from screen to world coordinate ; linear mapping
// uses global cons
double GiveZx(unsigned int ix)
{ return (ZxMin + ix*PixelWidth );}

// uses globaal cons
double GiveZy(unsigned int iy)
{ return (ZyMax - iy*PixelHeight);} // reverse y axis

/* -----------  array functions = drawing -------------- */

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

// plots raster point (ix,iy) 
int iDrawPoint(unsigned char A[], unsigned int ix, unsigned int iy, unsigned char iColor)
{ 

  /* i =  Give_i(ix,iy) compute index of 1D array from indices of 2D array */
  A[Give_i(ix,iy)] = iColor;

  return 0;
}

// draws point to memmory array data
// uses complex type so #include <complex.h> and -lm 
int dDrawPoint(unsigned char A[], complex double point,unsigned char iColor )
{

  unsigned int ix, iy; // screen coordinate = indices of virtual 2D array
  //unsigned int i; // index of 1D array
  
  ix = (creal(point)- ZxMin)/PixelWidth; 
  iy = (ZyMax - cimag(point))/PixelHeight; // inverse Y axis 
  iDrawPoint(A, ix, iy, iColor);
  return 0;
}



//;;;;;;;;;;;;;;;;;;;;;;  setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

int setup(int ParentPeriod, int ChildPeriod)
{

  
  
  
  
  
  printf("setup\n");

  ratio = (ZxMax-ZxMin)/(ZyMax-ZyMin);

  
  Cx=-1.25;
  Cy=0.0;
  c=Cx+Cy*I;
  
  

  /* 2D array ranges */
  if (!(iHeight % 2)) iHeight+=1; // it sholud be even number (variable % 2) or (variable & 1)
  iWidth = ratio*iHeight;
  iSize = iWidth*iHeight; // size = number of points in array 
  // iy
  iyMax = iHeight - 1 ; // Indexes of array starts from 0 not 1 so the highest elements of an array is = array_name[size-1].
  //ix
  
  ixMax = iWidth - 1;

  /* 1D array ranges */
  // i1Dsize = i2Dsize; // 1D array with the same size as 2D array
  iMax = iSize-1; // Indexes of array starts from 0 not 1 so the highest elements of an array is = array_name[size-1].

  /* Pixel sizes */
  PixelWidth = (ZxMax-ZxMin)/ixMax; //  ixMax = (iWidth-1)  step between pixels in world coordinate 
  PixelHeight = (ZyMax-ZyMin)/iyMax;
  ratio = ((ZxMax-ZxMin)/(ZyMax-ZyMin))/((float)iWidth/(float)iHeight); // it should be 1.000 ...
  
  

  // for numerical optimisation in iteration
  ER2 = ER * ER;
  
  //
  Zx2bMinus = Zx2b - 0.2;
  Zx2aMinus = Zx2a - 0.2;  
  Zx2bPlus = Zx2b + 0.2;
  Zx2aPlus = Zx2a + 0.2;  

  /* create dynamic 1D arrays for colors ( shades of gray ) */
  data = malloc( iSize * sizeof(unsigned char) );
  edge = malloc( iSize * sizeof(unsigned char) );
  julia = malloc( iSize * sizeof(unsigned char) );
  zero = malloc( iSize * sizeof(unsigned char) );

  if (edge== NULL || data == NULL || zero==NULL || julia ==NULL)
    {
      fprintf(stderr," Could not allocate memory");
      getchar(); 
      return 1;
    }

  
  

 
  

  
   
  
  printf(" end of setup \n");
  
  return 0;

} // ;;;;;;;;;;;;;;;;;;;;;;;;; end of the setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

unsigned char ComputeColor(unsigned int ix, unsigned int iy, int IterationMax)
{ 
  // check behavour of z under fc(z)=z^2+c
  // using 1 target set:
  // 1. exterior or circle (center at origin and radius ER ) 
  // as a target set containing infinity = for escaping points ( bailout test)
  // for points of exterior of julia set
  

  double Zx2, Zy2;
  int i=0;  // number of the iteration = fc(z)
  
 
  double Zx, Zy;
  
 
 
  
  
  // from screen to world coordinate 
  Zx = GiveZx(ix);
  Zy = GiveZy(iy);
  

  // if not inside target set around attractor ( alfa fixed point )
  while (1 )
    { // then iterate 
      
      Zx2 = Zx*Zx; 
      Zy2 = Zy*Zy;
       
      // bailout test 
      if (Zx2 + Zy2 > ER2) return iColorOfExterior; // if escaping stop iteration
       
      // if not escaping or not attracting then iterate = check behaviour
      // new z : Z(n+1) = Zn * Zn  + C
      Zy = 2*Zx*Zy + Cy; 
      Zx = Zx2 - Zy2 + Cx; 
      //
      i+=1;
      if (i > IterationMax) break; 
      
     
       
                  
	  
    }
    

  if (Zx2aMinus < Zx && Zx < Zx2a )     return iColors[0];
  if (Zx2a      < Zx && Zx < Zx2aPlus ) return iColors[1];
  if (Zx2bMinus < Zx && Zx < Zx2b )     return iColors[2];
  if (Zx2b      < Zx && Zx < Zx2bPlus ) return iColors[3];
      
    
  iNumberOfUknknown +=1; 
  return   iColorOfUnknown;   //
}

// plots raster point (ix,iy) 
int PlotPoint(unsigned char A[] , unsigned int ix, unsigned int iy, int IterationMax)
{
  unsigned i; /* index of 1D array */
  unsigned char iColor;
  

  i = Give_i(ix,iy); /* compute index of 1D array from indices of 2D array */
  iColor = ComputeColor(ix, iy, IterationMax);
  A[i] = iColor;

  return 0;
}



int FillArrayWithColor(unsigned char A[], unsigned char color )
{
  int i;
  for(i = 0; i <= iSize; ++i) A[i]=color; 
  return 0;


}

// fill array 
// uses global var :  ...
// scanning complex plane 
int ComputeFatouComponents(unsigned char A[], int IterationMax )
{
  unsigned int ix, iy; // pixel coordinate 

  //printf("compute image \n");
  // for all pixels of image 
#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(ixMax , iyMax, IterationMax)
  for(iy = iyMin; iy<=iyMax; ++iy) 
    { printf(" %d z %d \r", iy, iyMax); //info 
      for(ix= ixMin; ix<=ixMax; ++ix) PlotPoint(A, ix, iy, IterationMax ) ; //  
    } 
   
  return 0;
}




// from Source to Destination
int ComputeBoundaries(unsigned char S[], unsigned char D[])
{
 
  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; 
  // boundaries are in D  array ( global var )
 
  // clear D array
  FillArrayWithColor(D , iColorOfExterior);
 
  // printf(" find boundaries in S array using  Sobel filter\n");   
#pragma omp parallel for schedule(dynamic) private(i,iY,iX,Gv,Gh,G) shared(iyMax,ixMax, ER2)
  for(iY=1;iY<iyMax-1;++iY){ 
    for(iX=1;iX<ixMax-1;++iX){ 
      Gv= S[Give_i(iX-1,iY+1)] + 2*S[Give_i(iX,iY+1)] + S[Give_i(iX-1,iY+1)] - S[Give_i(iX-1,iY-1)] - 2*S[Give_i(iX-1,iY)] - S[Give_i(iX+1,iY-1)];
      Gh= S[Give_i(iX+1,iY+1)] + 2*S[Give_i(iX+1,iY)] + S[Give_i(iX-1,iY-1)] - S[Give_i(iX+1,iY-1)] - 2*S[Give_i(iX-1,iY)] - S[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) {D[i]=255;} /* background */
      else {D[i]=0;}  /* boundary */
    }
  }
 
   
 
  return 0;
}



// copy from Source to Destination
int CopyBoundaries(unsigned char S[],  unsigned char D[])
{
 
  unsigned int iX,iY; /* indices of 2D virtual array (image) = integer coordinate */
  unsigned int i; /* index of 1D array  */
 
 
  //printf("copy boundaries from S array to D array \n");
  for(iY=1;iY<iyMax-1;++iY)
    for(iX=1;iX<ixMax-1;++iX)
      {i= Give_i(iX,iY); if (S[i]==0) D[i]=0;}
 
 
 
  return 0;
}


 
 
int DrawCriticalOrbit(unsigned char A[], unsigned int IterMax)
{
  // integer = pixel coordinate
  unsigned int ix, iy; 
  // double = world coordinate
  // critical point Z= Zx+ZY*i;
  double Zx=0.0;
  double Zy=0.0; 
  double Zx2=0.0;
  double Zy2=0.0;
  
  unsigned int i; /* index of 1D array */
  unsigned int j; // number of iteration
 
 
  // draw critical point  
  ix = (int)((Zx-ZxMin)/PixelWidth);
  iy = (int)((ZyMax-Zy)/PixelHeight); // reverse y axis
  i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */
  A[i]=255-A[i]; 
 
  // iterate
  for (j = 1; j <= IterMax; j++) //larg number of iteration s
    {  Zx2 = Zx*Zx; 
      Zy2 = Zy*Zy;
 
      // bailout test 
      if (Zx2 + Zy2 > ER2) return iColorOfExterior; // if escaping stop iteration
 
      // if not escaping iterate
      // Z(n+1) = Zn * Zn  + C
      Zy = 2*Zx*Zy + Cy; 
      Zx = Zx2 - Zy2 + Cx;
      //compute integer coordinate  
      ix = (int)((Zx-ZxMin)/PixelWidth);
      iy = (int)((ZyMax-Zy)/PixelHeight); // reverse y axis
      i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */
      A[i]=0; //255-A[i];   // mark the critical orbit
 
    }
  return 0;
}

// Check Orientation of image : mark first quadrant 
// it should be in the upper right position
// uses global var :  ...
int CheckOrientation(unsigned char A[] )
{
  unsigned int ix, iy; // pixel coordinate 
  double Zx, Zy; //  Z= Zx+ZY*i;
  unsigned i; /* index of 1D array */
  for(iy=iyMin;iy<=iyMax;++iy) 
    {
      Zy = GiveZy(iy);
      for(ix=ixMin;ix<=ixMax;++ix) 
	{

	  // from screen to world coordinate 
	  Zx = GiveZx(ix);
	  i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */
	  if (Zx>0 && Zy>0) A[i]=255-A[i];   // check the orientation of Z-plane by marking first quadrant */

	}
    }
   
  return 0;
}

unsigned char ComputeColorZero(unsigned int ix, unsigned int iy, int iMax)
{ 
  // check behavour of z under fc(z)=z^2+c
  // using 1 target set:
  // 1. exterior or circle (center at origin and radius ER ) 
  // as a target set containing infinity = for escaping points ( bailout test)
  // for points of exterior of julia set
  

  double Zx2, Zy2;
  int i;  // number of the iteration = fc(z)
  
  unsigned char iColor;
  double Zx, Zy;
  
 
 
  i = 0;
  
  // from screen to world coordinate 
  Zx = GiveZx(ix);
  Zy = GiveZy(iy);
  

  
  while (i <iMax) 
    { // then iterate 
      
      Zx2 = Zx*Zx; 
      Zy2 = Zy*Zy;
       
             
      // if not escaping or not attracting then iterate = check behaviour
      // new z : Z(n+1) = Zn * Zn  + C
      Zy = 2*Zx*Zy + Cy; 
      Zx = Zx2 - Zy2 + Cx; 
      //

      i+=1;
    }
  //
  if (  Zy>0) iColor = 150; //re(z_n) > 0 and im(z_n) > 0 (first quadrant)
  // if ( Zx<0 && Zy>0) iColor = 170; //re(z_n) < 0 and im(z_n) > 0 (second)
  //if ( Zx<0 && Zy<0) iColor = 190; //re(z_n) < 0 and im(z_n) < 0 (third)
  if (  Zy<0) iColor = 200; //re(z_n) > 0 and im(z_n) < 0 (fourth).
  //
  return iColor;  
}


// plots raster point (ix,iy) 
int PlotPointZero(unsigned char A[] , unsigned int ix, unsigned int iy, int n)
{
  unsigned i; /* index of 1D array */
  unsigned char iColor;
  

  i = Give_i(ix,iy); /* compute index of 1D array from indices of 2D array */
  iColor = ComputeColorZero(ix, iy, n);
  // mark Fatou components, first interior 
  if (data[i]==iColors[0] || data[i]==iColors[1] || data[i]==iColors[2] || data[i]==iColors[3] )
    A[i] =  iColor+15; // !!! 
  else // exterior
    {if (n<10) A[i]= iColor; else A[i]= iColorOfExterior;} // exterior , only for low n 

  return 0;
}

// fill array 
// uses global var :  ...
// scanning complex plane 
int ComputeZerosOfQnc(unsigned char A[], int n)
{
  unsigned int ix, iy; // pixel coordinate 

  //printf("compute image \n");
  // for all pixels of image 
#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(ixMax , iyMax, n)
  for(iy = iyMin; iy<=iyMax; ++iy) 
    { printf(" %d z %d \r", iy, iyMax); //info 
      for(ix= ixMin; ix<=ixMax; ++ix) PlotPointZero(A, ix, iy, n ) ; //  
    } 
   
  return 0;
}
 

// save "A" array to pgm file 
int SaveArray2PGMFile( unsigned char A[], double k)
{
  
  FILE * fp;
  const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
  char name [30]; /* name of file */
  sprintf(name,"%.0f", k); /*  */
  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,iWidth,iHeight,MaxColorComponentValue);  /*write header to the file*/
  fwrite(A,iSize,1,fp);  /*write A array to the file in one step */
  printf("File %s saved. \n", filename);
  fclose(fp);

  return 0;
}

int info()
{
  // diplay info messages
  printf("Numerical approximation of parabolic Julia set for fc(z)= z^2 + c \n");
  printf("iPeriodParent = %d \n", iPeriodParent);
  printf("iPeriodOfChild  = %d \n", iPeriodChild);
  printf("parameter c = ( %f ; %f ) \n", Cx, Cy); 
  
  printf("Image Width = %f \n", ZxMax-ZxMin);
  printf("PixelWidth = %f \n", PixelWidth);
  printf("Maximal number of iterations = iterMax = %ld \n", iterMax);
  printf("iNumberOfUknknown  = %d \n", iNumberOfUknknown);
  printf("ratio of image  = %f \n", ratio);

  return 0;
}

/* -----------------------------------------  main   -------------------------------------------------------------*/
int main()
{
  setup(iPeriodParent, iPeriodChild);
  int n; 

 
   
    
 
  printf("components of Fatou set  : \n");
  ComputeFatouComponents(data, iterMax);
  SaveArray2PGMFile( data, 10000); 
  printf("done \n");

  ComputeBoundaries(data, julia);
  SaveArray2PGMFile( julia, 20000); 
  printf(" Julia set \n");
  CopyBoundaries(julia, data);
  SaveArray2PGMFile( data, 30000);



  
  
  for(n = 0; n<=25; ++n)
    { 
      ComputeZerosOfQnc(zero, n);
      ComputeBoundaries(zero,edge);
      CopyBoundaries(julia, edge); // julia is computed more precisely then edge 
      CopyBoundaries(edge, zero);
      SaveArray2PGMFile( zero, n); 
    }
  
 
 

  printf(" allways free memory  to avoid buffer overflow \n");
 
  free(data);
  free(edge);
  free(zero);
  free(julia);

  
  info();

  return 0;
}