/*

  Adam Majewski
  fraktal.republika.pl
https://gitlab.com/adammajewski/parabolic_chesboard_using_triangle_in_c

  c console progam using 
  * symmetry
  * openMP

  draw  julia sets

  gcc i.c -lm -Wall -fopenmp -march=native 
  time ./a.out
  time ./a.out > info.txt

  How to tune up parameter ? (  )
   
parabolic checkerboard
use target set ( 4 triangles) defined by  points : 
* critical point Zcr=0
* parabolic alfa fixed point Za 
* preimage of alfa 
* 2 precritical points :
* + f^{-3} (Zcr) 
* - f^{-3} (Zcr)

*/

#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> // OpenMP; needs also -fopenmp

/* --------------------------------- global variables and consts ------------------------------------------------------------ */
#define iPeriodChild 3 // iPeriodChild of secondary component joined by root point
int iPeriodParent = 1;
unsigned int denominator;
double InternalAngle;

// virtual 2D array and integer ( screen) coordinate
// Indexes of array starts from 0 not 1 
unsigned int ix, iy; // var
unsigned int ixMin = 0; // Indexes of array starts from 0 not 1
unsigned int ixMax ; //
unsigned int iWidth ; // horizontal dimension of array
unsigned int ixAxisOfSymmetry  ; // 
unsigned int iyMin = 0; // Indexes of array starts from 0 not 1
unsigned int iyMax ; //
unsigned int iyAxisOfSymmetry  ; // 
unsigned int iyAbove ; // var, measured from 1 to (iyAboveAxisLength -1)
unsigned int iyAboveMin = 1 ; //
unsigned int iyAboveMax ; //
unsigned int iyAboveAxisLength ; //
unsigned int iyBelowAxisLength ; //
unsigned int iHeight = 1000; //  odd number !!!!!! = (iyMax -iyMin + 1) = iyAboveAxisLength + iyBelowAxisLength +1
// The size of array has to be a positive constant integer 
unsigned int iSize ; // = iWidth*iHeight; 

// memmory 1D array 
unsigned char *data;
unsigned char *edge;

// unsigned int i; // var = index of 1D array
unsigned int iMin = 0; // Indexes of array starts from 0 not 1
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 */
const double ZxMin=-1.5;
const double ZxMax=1.5;
const double ZyMin=-1.5;
const double ZyMax=1.5;
double PixelWidth; // =(ZxMax-ZxMin)/iXmax;
double PixelWidth2; // =  PixelWidth*PixelWidth;
double PixelHeight; // =(ZyMax-ZyMin)/iYmax;

double ratio ;

// complex numbers of parametr plane 
double Cx; // c =Cx +Cy * i
double Cy;
double complex c; // 

double complex Za; // alfa fixed point alfa=f(alfa)
double Zax, Zay;

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

// points defining target sets 
double Zlx, Zly, Zrx, Zry; 

// critical point Zcr
double Zcrx = 0.0;
double Zcry=0.0;

unsigned char  iColorsOfInterior[iPeriodChild]={110, 160,210}; // number of colors >= iPeriodChild
static unsigned char iColorOfExterior = 245;
static unsigned char iColorOfUnknown = 100;
unsigned char iJulia = 0;
/* ------------------------------------------ functions -------------------------------------------------------------*/

/* find c in component of Mandelbrot set 
 
   uses code by Wolf Jung from program Mandel
   see function mndlbrot::bifurcate from mandelbrot.cpp
   http://www.mndynamics.com/indexp.html

*/
double complex GiveC(double InternalAngleInTurns, double InternalRadius, unsigned int Period)
{
  //0 <= InternalRay<= 1
  //0 <= InternalAngleInTurns <=1
  double t = InternalAngleInTurns *2*M_PI; // from turns to radians
  double R2 = InternalRadius * InternalRadius;
  //double Cx, Cy; /* C = Cx+Cy*i */
  switch ( Period ) // of component 
    {
    case 1: // main cardioid
      Cx = (cos(t)*InternalRadius)/2-(cos(2*t)*R2)/4; 
      Cy = (sin(t)*InternalRadius)/2-(sin(2*t)*R2)/4; 
      break;
    case 2: // only one component 
      Cx = InternalRadius * 0.25*cos(t) - 1.0;
      Cy = InternalRadius * 0.25*sin(t); 
      break;
      // for each iPeriodChild  there are 2^(iPeriodChild-1) roots. 
    default: // higher periods : to do, use newton method 
      Cx = 0.0;
      Cy = 0.0; 
      break; }

  return Cx + Cy*I;
}

/*

  http://en.wikipedia.org/wiki/Periodic_points_of_complex_quadratic_mappings
  z^2 + c = z
  z^2 - z + c = 0
  ax^2 +bx + c =0 // ge3neral for  of quadratic equation
  so :
  a=1
  b =-1
  c = c
  so :

  The discriminant is the  d=b^2- 4ac 

  d=1-4c = dx+dy*i
  r(d)=sqrt(dx^2 + dy^2)
  sqrt(d) = sqrt((r+dx)/2)+-sqrt((r-dx)/2)*i = sx +- sy*i

  x1=(1+sqrt(d))/2 = beta = (1+sx+sy*i)/2

  x2=(1-sqrt(d))/2 = alfa = (1-sx -sy*i)/2

  alfa : attracting when c is in main cardioid of Mandelbrot set, then it is in interior of Filled-in Julia set, 
  it means belongs to Fatou set ( strictly to basin of attraction of finite fixed point )

*/
// uses global variables : 
//  ax, ay (output = alfa(c)) 
double complex GiveAlfaFixedPoint(double complex c)
{
  double dx, dy; //The discriminant is the  d=b^2- 4ac = dx+dy*i
  double r; // r(d)=sqrt(dx^2 + dy^2)
  double sx, sy; // s = sqrt(d) = sqrt((r+dx)/2)+-sqrt((r-dx)/2)*i = sx + sy*i
  double ax, ay;
 
  // d=1-4c = dx+dy*i
  dx = 1 - 4*creal(c);
  dy = -4 * cimag(c);
  // r(d)=sqrt(dx^2 + dy^2)
  r = sqrt(dx*dx + dy*dy);
  //sqrt(d) = s =sx +sy*i
  sx = sqrt((r+dx)/2);
  sy = sqrt((r-dx)/2);
  // alfa = ax +ay*i = (1-sqrt(d))/2 = (1-sx + sy*i)/2
  ax = 0.5 - sx/2.0;
  ay =  sy/2.0;
 

  return ax+ay*I;
}

int setup()
{

  
  

  denominator = iPeriodChild;
  InternalAngle = 1.0/((double) denominator);

  c = GiveC(InternalAngle, 1.0, 1) ;
  Cx=creal(c);
  Cy=cimag(c);
  Za = GiveAlfaFixedPoint(c);
  Zax = creal(Za);
  Zay = cimag(Za);

  //

  

 

  /* virtual 2D array ranges */
  if (!(iHeight % 2)) iHeight+=1; // it sholud be even number (variable % 2) or (variable & 1)
  iWidth = 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].
  iyAboveAxisLength = (iHeight -1)/2;
  iyAboveMax = iyAboveAxisLength ; 
  iyBelowAxisLength = iyAboveAxisLength; // the same 
  iyAxisOfSymmetry = iyMin + iyBelowAxisLength ; 
  // 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 
  ER2 = ER * ER;
  
  PixelWidth2 =  PixelWidth*PixelWidth;
 
  /* 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");

  
  // points defining target sets 
  //aaa = -0.229955135116281  -0.141357981605006 i
  Zlx = -0.229955135116281 ; // aaa  
  Zly = -0.141357981605006 ;
  Zrx = -Zlx; // aab = 0.229955135116281  +0.141357981605006 i
  Zry = -Zly;

  
 
  return 0;

}

// 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

// ============ http://stackoverflow.com/questions/2049582/how-to-determine-a-point-in-a-2d-triangle
// In general, the simplest (and quite optimal) algorithm is checking on which side of the half-plane created by the edges the point is.
// Kornel Kisielewicz
double sign (double  x1, double y1, double x2,double y2,double x3, double y3)
{
    return (x1 - x3) * (y2 - y3) - (x2 - x3) * (y1 - y3);
}

// Kornel Kisielewicz
int  PointInTriangle (double x, double y, double x1, double y1, double x2, double y2, double x3, double y3)
{
    double  b1, b2, b3;

    b1 = sign(x, y, x1, y1, x2, y2) < 0.0;
    b2 = sign(x, y, x2, y2, x3, y3) < 0.0;
    b3 = sign(x, y, x3, y3, x1, y1) < 0.0;

    return ((b1 == b2) && (b2 == b3));
}

int IfInsideTarget(double Zx, double Zy)
{
 
return PointInTriangle(Zx, Zy, Zax, Zay, Zrx, Zry,  Zlx, Zly)       ;
  
}

unsigned char GiveColor(unsigned int ix, unsigned int iy)
{ 
  // check behavour of z under fc(z)=z^2+c
  // using 2 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
  // 2. interior : triangle 

  double Zx2, Zy2;
  int i=0;
  //int j; // iteration = fc(z)
  
  double Zx, Zy;

  int t=0; // test , boolean value ; 0 = false 
  
  
  
  
  // from screen to world coordinate 
  Zx = GiveZx(ix);
  Zy = GiveZy(iy);
  
  t = IfInsideTarget( Zx, Zy);
  if (t) return iColorsOfInterior[2];
    
    

  // if not inside target set around attractor ( alfa fixed point )
  while (!t )
    { // then iterate 
     

    
      for(i=0;i<iPeriodChild ;++i) // iMax = period !!!!
	{  
	  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; 
	  //
	  t = IfInsideTarget( Zx, Zy);
	  if (t) return iColorsOfInterior[ i ];
     
   
	}
      
      
      
    }

  
  
  return iColorOfUnknown; // never
}

 

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

/* 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; }
//  ix = i % iWidth;
//  iy = (i- ix) / iWidth;
//  i  = Give_i(ix, iy);

// plots raster point (ix,iy) 
int PlotPoint(unsigned int ix, unsigned int iy, unsigned char iColor, unsigned char a[])
{
  unsigned i; /* index of 1D array */
  i = Give_i(ix,iy); /* compute index of 1D array from indices of 2D array */
  a[i] = iColor;

  return 0;
}

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

  // for all pixels of image 
  for(iy = iyMin; iy<=iyMax; ++iy) 
    { printf(" %d z %d\r", iy, iyMax); //info 
      for(ix= ixMin; ix<=ixMax; ++ix) PlotPoint(ix, iy, GiveColor(ix, iy) , a); //  
    } 
   
  return 0;
}

// fill array using symmetry of image 
// uses global var :  ...
int FillArraySymmetric(unsigned char a[] )
{
   
  unsigned char Color; // gray from 0 to 255 

  printf("axis of symmetry \n"); 
  iy = iyAxisOfSymmetry; 
#pragma omp parallel for schedule(dynamic) private(ix,Color) shared(ixMin,ixMax, iyAxisOfSymmetry)
  for(ix=ixMin;ix<=ixMax;++ix) {//printf(" %d from %d\n", ix, ixMax); //info  
    PlotPoint(ix, iy, GiveColor(ix, iy), a);
  }

  /*
    The use of ‘shared(variable, variable2) specifies that these variables should be shared among all the threads.
    The use of ‘private(variable, variable2)’ specifies that these variables should have a seperate instance in each thread.
  */

#pragma omp parallel for schedule(dynamic) private(iyAbove,ix,iy,Color) shared(iyAboveMin, iyAboveMax,ixMin,ixMax, iyAxisOfSymmetry)

  // above and below axis 
  for(iyAbove = iyAboveMin; iyAbove<=iyAboveMax; ++iyAbove) 
    {printf(" %d from %d\r", iyAbove, iyAboveMax); //info 
      for(ix=ixMin; ix<=ixMax; ++ix) 

	{ // above axis compute color and save it to the array
	  iy = iyAxisOfSymmetry + iyAbove;
	  Color = GiveColor(ix, iy);
	  PlotPoint(ix, iy, Color , a); 
	  // below the axis only copy Color the same as above without computing it 
	  PlotPoint(ixMax-ix, iyAxisOfSymmetry - iyAbove , Color , a); 
	} 
    }  
  return 0;
}

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");   
#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= 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]=255;} /* 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;
}

// Check Orientation of image : first quadrant in 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;
}

int IfInsideTarget2a(double Zx, double Zy)
{
 
return PointInTriangle(Zx, Zy, Zax, Zay,  Zrx, Zry, Zcrx, Zcry)       ;
  
}

int IfInsideTarget2b(double Zx,double Zy)
{
 
return PointInTriangle(Zx, Zy, Zax, Zay, Zcrx, Zcry, Zlx, Zly)       ;
  
}

int IsInside(double Zx, double Zy)
{
  // inside : check 4 triangles 
 // check one triangle and its preimage  
 if  (PointInTriangle(Zx, Zy, Zax, Zay,  Zrx, Zry, Zcrx, Zcry))   return 1;      
 if  (PointInTriangle(Zx, Zy, -Zax, -Zay,  Zlx, Zly, Zcrx, Zcry)) return 1;
 // check one triangle and its preimage
 if  (PointInTriangle(Zx, Zy, Zax, Zay,  Zlx, Zly, Zcrx, Zcry))   return 2;      
 if  (PointInTriangle(Zx, Zy, -Zax, -Zay,  Zrx, Zry, Zcrx, Zcry)) return 2;
 // if it is not in the target then not inside 
 return 0; 

  
 
  
}

// 
unsigned char GiveColor2(unsigned int ix, unsigned int iy)
{ 
  

  double Zx2, Zy2;
  long int i=0;
    
  double Zx, Zy, Zx0, Zy0;
  int flag=0;
  
  
  
  
  
  // from screen to world coordinate 
  Zx0 = GiveZx(ix);
  Zy0 = GiveZy(iy);
  Zx= Zx0;
  Zy = Zy0;
  
  
  flag = IsInside( Zx, Zy);
  if (flag ) 
   { if (flag==1) return iColorsOfInterior[2];
     if (flag==2) return iColorsOfInterior[2]+20;}
   
    
    

  // if not inside target set around attractor ( alfa fixed point )
  while (1 )
    { // then iterate 
     

    
      for(i=0;i<iPeriodChild ;++i) // iMax = period !!!!
	{  
	  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; 
	  //
	  flag = IsInside( Zx, Zy);
          if (flag ) 
            { if (flag==1) return iColorsOfInterior[i];
              if (flag==2) return iColorsOfInterior[i]+20;}
       }
      
       
              
      
    }

  
  
  return iColorOfUnknown; // it should never happen
}

int MakeInternalTilingS(unsigned char a[] )
{
     
  unsigned char Color; // gray from 0 to 255 

  printf("axis of symmetry \n"); 
  iy = iyAxisOfSymmetry; 
#pragma omp parallel for schedule(dynamic) private(ix,Color) shared(ixMin,ixMax, iyAxisOfSymmetry)
  for(ix=ixMin;ix<=ixMax;++ix) {//printf(" %d from %d\n", ix, ixMax); //info  
    PlotPoint(ix, iy, GiveColor2(ix, iy),a);
  }

  /*
    The use of ‘shared(variable, variable2) specifies that these variables should be shared among all the threads.
    The use of ‘private(variable, variable2)’ specifies that these variables should have a seperate instance in each thread.
  */

#pragma omp parallel for schedule(dynamic) private(iyAbove,ix,iy,Color) shared(iyAboveMin, iyAboveMax,ixMin,ixMax, iyAxisOfSymmetry)

  // above and below axis 
  for(iyAbove = iyAboveMin; iyAbove<=iyAboveMax; ++iyAbove) 
    {printf(" %d from %d\r", iyAbove, iyAboveMax); //info 
      for(ix=ixMin; ix<=ixMax; ++ix) 

	{ // above axis compute color and save it to the array
	  iy = iyAxisOfSymmetry + iyAbove;
	  Color = GiveColor2(ix, iy);
	  PlotPoint(ix, iy, Color, a ); 
	  // below the axis only copy Color the same as above without computing it 
	  PlotPoint(ixMax-ix, iyAxisOfSymmetry - iyAbove , Color, a ); 
	} 
    }  
  return 0;
}

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

  // for all pixels of image 
  for(iy = iyMin; iy<=iyMax; ++iy) 
    { printf(" %d z %d\r", iy, iyMax); //info 
      for(ix= ixMin; ix<=ixMax; ++ix) 
           PlotPoint(ix, iy, GiveColor2(ix, iy),a ); //  
    } 
   
  return 0;
}

int DrawCriticalOrbit(unsigned char A[], unsigned int IterMax)
{
 
  unsigned int ix, iy; // pixel coordinate 
  // initial point z0 = critical point
  double Zx=0.0; 
  double Zy=0.0; //  Z= Zx+ZY*i;
  double Zx2=0.0;
  double Zy2=0.0;
  unsigned int i; /* index of 1D array */
  unsigned int j;

  // 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 1; // if escaping stop iteration and return error code
       
      // 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;
}

// fill array 
// uses global var :  ...
// scanning complex plane 
int FillAndMarkTarget(unsigned char a[] )
{
  unsigned int ix, iy; // pixel coordinate 
  double Zx, Zy; //  Z= Zx+ZY*i;
  unsigned i; /* index of 1D array */
  int flag;

  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 ( IfInsideTarget2a(Zx,Zy)) data[i]=data[i]-43;
          //if ( IfInsideTarget2b(Zx,Zy)) data[i]=data[i]+43;   // 

          flag = IsInside( Zx, Zy);
          if (flag ) 
           { if (flag==1)  a[i]=a[i] - 43;
              if (flag==2)  a[i]=a[i]+43;}
          
	}
    }
   
  return 0;
}

// save data 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 [50]; /* 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 image data bytes 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 complex quadratic polynomial fc(z)= z^2 + c \n");
  printf("parameter c  = %.16f , %.16f \n", Cx, Cy);
  printf("c is a root point between hyperbolic components of period %d and %d  of Mandelbrot set \n", iPeriodParent,  iPeriodChild);
  printf("parabolic alfa fixed point Za  = %.16f ; %.16f \n", Zax, Zay);
  printf("image size in pixels = %d x %d \n", iWidth, iHeight);
  printf("image size in world units : (ZxMin = %f, ZxMax =  %f) ,  (ZyMin = %f, ZyMax =  %f) \n", ZxMin , ZxMax ,  ZyMin , ZyMax);
  printf("ratio ( distortion) of image  = %f ; it should be 1.000 ...\n", ratio);
  printf("PixelWidth = %f \n", PixelWidth);
  printf("critical point Zcr  = %.16f ; %.16f \n", Zcrx, Zcry);
  printf("precritical point Zl  = %.16f ; %.16f \n", Zlx, Zly);
  printf("precritical point Zr  = %.16f ; %.16f \n", Zrx, Zry);
  return 0;
}

/* -----------------------------------------  main   -------------------------------------------------------------*/
int main()
{
  

  setup();

  // here are procedures for creating image file
  
  // compute colors of pixels = image
  //FillArray( data ); // no symmetry
  FillArraySymmetric(data); 
  AddBoundaries(data);
  SaveArray2PGMFile(data , 0); // save array (image) to pgm file 

  FillAndMarkTarget(data);
  SaveArray2PGMFile(data , 1); // save array (image) to pgm file

  DrawCriticalOrbit( data, 1000000);
  SaveArray2PGMFile(data , 2); // save array (image) to pgm file 

  SaveArray2PGMFile(edge , 3); // save array (image) to pgm file 

  MakeInternalTilingS(data);
  SaveArray2PGMFile(data , 4); // save array (image) to pgm file 

 // DrawCriticalOrbit( data, 1000000);
//  SaveArray2PGMFile(data , 5); // save array (image) to pgm file 

  AddBoundaries(data);
  SaveArray2PGMFile(data , 6); // save array (image) to pgm file 

  FillAndMarkTarget(data);
  SaveArray2PGMFile(data , 7); // save array (image) to pgm file

  CheckOrientation(data);
  SaveArray2PGMFile(data , 8); // save array (image) to pgm file

  free(data);
  free(edge);
  info();
  return 0;
}