/*

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

  https://plus.google.com/116648956837292097606/posts/b6J6z2u8soL

  how to show sepals inside main box of parabolic chessboard ?
  - compute full orbit ( forward and backward of every point)
  - for each whole orbit ( not point) compute maximal distance from orbit to fixed point alfa
  - normalize distance ( dustance/ distance max ) so it will have value from 0 to 1.0
  _ use such normalized distance for coloring
  - then one can see orbits 

  -------------------------------
  cd existing_folder
  git init
  git remote add origin git@gitlab.com:adammajewski/SepalsOfCauliflower.git
  git add .
  git commit
  git push -u origin master
  ---------------------------------
  indent d.c 
  defaoult is gnu style 
  -------------------

  c console progam 

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

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

  time ./a.out

  time ./a.out >a.txt

----------------------- splint --------------------------
splint d.c
Splint 3.1.2 --- 03 May 2009

d.c:60:17: Cannot find include file omp.h on search path: /usr/include;/usr/incl
              ude
  Preprocessing error. (Use -preproc to inhibit warning)
d.c:132:106: Comment starts inside comment
  A comment open sequence  appears within a comment.  This usually means an
  earlier comment was not closed. (Use -nestcomment to inhibit warning)
Preprocessing error for file: /home/a/c/julia/parabolic/z2c/1over0/Sepals/d.c
*** Cannot continue.

gcc -C -E d.c -lm -Wall -march=native -fopenmp -o d.txt

gcc -C -E d.c -lm -Wall -march=native -fopenmp -o d.txt
d.c:47:28: warning within comment [-Wcomment]
   A comment open sequence  appears within a comment.  This usually means an
 
------------------------------ gcov -------------------------------

 gcc d.c -fprofile-arcs -ftest-coverage -lm -Wall -march=native -fopenmp 
 time ./a.out

 gcov d.c

gcov -a d.c

File 'd.c'
Lines executed:91.26% of 286
Creating 'd.c.gcov'

----------------------

setup
 end of setup 
components of Fatou set  : 
File 25.0.pgm saved. 
File 25.1.pgm saved. 
 find boundaries in A array using  Sobel filter
File 25.2.pgm saved. 
copy boundaries from edge array to data array 
File 25.3.pgm saved. 
File 25.5.pgm saved. 
File 25.6.pgm saved. 
 allways free memory  to avoid buffer overflow 
Numerical approximation of parabolic Julia set for fc(z)= z^2 + c 
iPeriodParent = 1 
iPeriodOfChild  = 0 
parameter c = ( 0.250000 ; 0.000000 ) 
Image Width = 3.000000 
PixelWidth = 0.006000 
Maximal number of iterations = iterMax = 10000 
ratio of image  = 1.000000 ; it should be 1.000 ...

real	0m0.963s
user	0m7.546s
sys	0m0.000s

======================================================================

   How about drawing orbits and coloring  the regions between different orbits?  or using an approximate Fatou
   coordinate?
  
   Actually,  you will need some approximation to avoid kinks in those places, 
   where the points of the same orbit have a relatively large distance from
   each other.
  
      Best regards,
  
      Wolf

*/

#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 0		// Period of secondary component joined by root point with the parent component
int iPeriodParent = 1;		// 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 = 2000;	//  
// The size of array has to be a positive constant integer 
static unsigned int iSize;	// = iWidth*iHeight; 

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

// 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 
   yrange = [-0.1,0.7],
   xrange = [-0.05,0.75],

*/

static const double ZxMin = -2.0;	//-0.05;
static const double ZxMax = 2.0;	//0.75;
static const double ZyMin = -2.0;	//-0.1;
static const double ZyMax = 2.0;	//0.7;
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;

/*
max distance of all arbits, computed manually

  double d = GiveAverageMaxDistanceOfFullOrbit(0.393075688878712,  +0.636009824757034 );
  printf("d = %.16f \n", d);
  d = GiveAverageMaxDistanceOfFullOrbit(0.0,  +0.5 );
  printf("d = %.16f \n", d);

*/
double DistanceMaxGlobal;	//= 0.0497920256372717 ;
//double DistanceMaxGlobal2  ;

/* 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]={255,231,180, 160}; //
static unsigned char iColorOfExterior = 245;
static unsigned char iColorOfInterior = 200;
unsigned char iColorOfUnknown = 133;
//int iNumberOfUknknown = 0;

//static double TwoPi=2*M_PI;

// number of point used to compute average maximal distance of full orbit
// n point befor + n point after + point with maximal distance = 2*n+1
// see function GiveAverageMaxDistanceOfFullOrbit
#define NrOfPoints 1000
// #define NrOfPointsUsedForComputation (2*NrOfPoints + 1)

// number of point used to compute average maximal distance of full orbit
// n point befor + n point after + point with maximal distance = 2*n+1
// see function GiveAverageMaxDistanceOfFullOrbit
 // double MaxDistances[NrOfPointsUsedForComputation]; // one must have more space to find NrOfPoints because we do not know when we find max distance

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

unsigned int
Give_is (unsigned int ix, unsigned int iy)
{				// double Zy = - GiveZy(iy); // symmetric
  //unsigned int iys =  (int)((ZyMax - Zy )/PixelHeight); 

  int dy = iy - (iHeight / 2);
  iy = (iHeight / 2) - dy;

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

/* mndyncxmics::root from mndyncxmo.cpp  by Wolf Jung (C) 2007-2014. */

// input = x,y
// output = u+v*I = sqrt(x+y*i) 
complex double
GiveRoot (complex double z)
{
  double x = creal (z);
  double y = cimag (z);
  double u, v;

  v = sqrt (x * x + y * y);

  if (x > 0.0)
    {
      u = sqrt (0.5 * (v + x));
      v = 0.5 * y / u;
      return u + v * I;
    }
  if (x < 0.0)
    {
      v = sqrt (0.5 * (v - x));
      if (y < 0.0)
	v = -v;
      u = 0.5 * y / v;
      return u + v * I;
    }
  if (y >= 0.0)
    {
      u = sqrt (0.5 * y);
      v = u;
      return u + v * I;
    }

  u = sqrt (-0.5 * y);
  v = -u;
  return u + v * I;
}

// from mndlbrot.cpp  by Wolf Jung (C) 2007-2014. part of Madel 5.12 
// input : c, z , mode
// c = cx+cy*i where cx and cy are global variables defined in mndynamo.h
// z = x+y*i
// 
// output : z = x+y*i
complex double
InverseIteration (complex double z, complex double c)
{
  double x = creal (z);
  double y = cimag (z);
  double cx = creal (c);
  double cy = cimag (c);

  // f^{-1}(z) = inverse with principal value
  if (cx * cx + cy * cy < 1e-20)
    {
      z = GiveRoot (x - cx + (y - cy) * I);	// 2-nd inverse function = key b 
      //if (mode & 1) { x = -x; y = -y; } // 1-st inverse function = key a   
      return -z;
    }

  //f^{-1}(z) =  inverse with argument adjusted
  double u, v;
  complex double uv;
  double w = cx * cx + cy * cy;

  uv = GiveRoot (-cx / w - (cy / w) * I);
  u = creal (uv);
  v = cimag (uv);
  //
  z = GiveRoot (w - cx * x - cy * y + (cy * x - cx * y) * I);
  x = creal (z);
  y = cimag (z);
  //
  w = u * x - v * y;
  y = u * y + v * x;
  x = w;

  //if (mode & 1) // mode = -1
  //  { x = -x; y = -y; } // 1-st inverse function = key a

  return x + y * I;		// 2-nd inverse function = key b 
}

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;

    }

  return iColorOfInterior;	//
}

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

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

int
ComputeBoundariesIn (unsigned char A[])	// compute in A, but save to edge
{

  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 edge array ( global var )

  // clear all pixels
  for (i = 1; i < iSize; ++i)
    edge[i] = iColorOfExterior;

  printf (" find boundaries in A 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 =
	    A[Give_i (iX - 1, iY + 1)] + 2 * A[Give_i (iX, iY + 1)] +
	    A[Give_i (iX - 1, iY + 1)] - A[Give_i (iX - 1, iY - 1)] -
	    2 * A[Give_i (iX - 1, iY)] - A[Give_i (iX + 1, iY - 1)];
	  Gh =
	    A[Give_i (iX + 1, iY + 1)] + 2 * A[Give_i (iX + 1, iY)] +
	    A[Give_i (iX - 1, iY - 1)] - A[Give_i (iX + 1, iY - 1)] -
	    2 * A[Give_i (iX - 1, iY)] - A[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 */
	}
    }

  return 0;
}

int
CopyBoundariesTo (unsigned char A[])	// copy boundaries from edge to A
{

  unsigned int iX, iY;		/* indices of 2D virtual array (image) = integer coordinate */
  unsigned int i;		/* index of 1D array  */

  printf ("copy boundaries from edge array to data array \n");
  for (iY = 1; iY < iyMax - 1; ++iY)
    for (iX = 1; iX < ixMax - 1; ++iX)
      {
	i = Give_i (iX, iY);
	if (edge[i] == 0)
	  A[i] = 0;
      }

  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)
{

  double Zx2, Zy2;
  int i = 0;			// number of the iteration = fc(z)

  unsigned char iColor;
  double Zx, Zy;

  // from screen to world coordinate 
  Zx = GiveZx (ix);
  Zy = GiveZy (iy);

  while (i < iMax)
    {				// then iterate 

      Zx2 = Zx * Zx;
      Zy2 = Zy * Zy;

      // new z : Z(n+1) = Zn * Zn  + C
      Zy = 2 * Zx * Zy + Cy;
      Zx = Zx2 - Zy2 + Cx;
      //
      i += 1;
    }

  /*

     From program Mandel by Wolf Jung http://www.mndynamics.com/indexp.html
     On parameter plane for complex quadratic map algorithm " 9 is showing the location of centers / p.p. The period is set with q. (Use 9 also to check the displayed period, which might be inaccurate.)"
   */
  // check position near fixed point after n iterations 
  if (Zx > 0 && 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 (Zx > 0 && 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);
  if (data[i] == iColorOfInterior)
    A[i] = iColor + 15;		// interior
  // exterior
  else
    {
      if (n < 10)
	A[i] = iColor;		// exterior , only for low n 
      else
	A[i] = iColorOfExterior;
    }

  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(A, 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;
}

/*  gives sign of number */
double
sign (double d)
{
  if (d < 0)
    {
      return -1.0;
    }
  else
    {
      return 1.0;
    };
};

/*
  full orbit = forward and backward orbit 
  maximal distance of full orbit
  distance = distance between point z ( of orbit) and parabolic fixed point alfa 

*/
double
GiveMaxDistanceOfFullOrbit (double Zx0, double Zy0)
{

  double Distance2MaxLocal = 0.0;
  double temp;

  double Zx = Zx0;
  double Zy = Zy0;
  double Zx2, Zy2;

  double NewZx, NewZy;

  Zx2 = Zx * Zx;
  Zy2 = Zy * Zy;
  //
  temp = (Zx - 0.5) * (Zx - 0.5) + Zy * Zy;

  //forward orbit in case when one knows that orbit 
  // will allways fall into alfa !!!
  while (temp > Distance2MaxLocal)
    {

      Distance2MaxLocal = temp;
      // then iterate 

      Zx2 = Zx * Zx;
      Zy2 = Zy * Zy;

      // new z : Z(n+1) = Zn * Zn  + C
      Zy = 2 * Zx * Zy + Cy;
      Zx = Zx2 - Zy2 + Cx;
      //
      temp = (Zx - 0.5) * (Zx - 0.5) + Zy * Zy;
    }

  // go back to initial point
  Zx = Zx0;
  Zy = Zy0;
  // and make backward iterations

  /* Zn*Zn=Z(n+1)-c */
  Zx = Zx - Cx;
  Zy = Zy - Cy;
  /* sqrt of complex number algorithm from Peitgen, Jurgens, Saupe: Fractals for the classroom */
  if (Zx > 0.0)
    {
      NewZx = sqrt ((Zx + sqrt (Zx * Zx + Zy * Zy)) / 2);
      NewZy = Zy / (2 * NewZx);
    }
  else				/* ZX <= 0 */
    {
      if (Zx < 0.0)
	{
	  NewZy = sign (Zy) * sqrt ((-Zx + sqrt (Zx * Zx + Zy * Zy)) / 2);
	  NewZx = Zy / (2 * NewZy);
	}
      else			/* Zx=0 */
	{
	  NewZx = sqrt (fabs (Zy) / 2);
	  if (NewZx > 0)
	    NewZy = Zy / (2 * NewZx);
	  else
	    NewZy = 0.0;
	}
    };

  temp = (NewZx - 0.5) * (NewZx - 0.5) + NewZy * NewZy;
  Zx = NewZx;
  Zy = NewZy;

  //
  while (temp > Distance2MaxLocal)
    {
      Distance2MaxLocal = temp;

      // and make backward iterations

      /* Zn*Zn=Z(n+1)-c */
      Zx = Zx - Cx;
      Zy = Zy - Cy;
      /* sqrt of complex number algorithm from Peitgen, Jurgens, Saupe: Fractals for the classroom */
      if (Zx > 0.0)
	{
	  NewZx = sqrt ((Zx + sqrt (Zx * Zx + Zy * Zy)) / 2);
	  NewZy = Zy / (2 * NewZx);
	}
      else			/* ZX <= 0 */
	{
	  if (Zx < 0.0)
	    {
	      NewZy = sign (Zy) * sqrt ((-Zx + sqrt (Zx * Zx + Zy * Zy)) / 2);
	      NewZx = Zy / (2 * NewZy);
	    }
	  else			/* Zx=0 */
	    {
	      NewZx = sqrt (fabs (Zy) / 2);
	      if (NewZx > 0)
		NewZy = Zy / (2 * NewZx);
	      else
		NewZy = 0.0;
	    }
	};

      temp = (NewZx - 0.5) * (NewZx - 0.5) + NewZy * NewZy;
      Zx = NewZx;
      Zy = NewZy;

    }

  return sqrt (Distance2MaxLocal);
}

/*
  full orbit = forward and backward orbit 
  maximal distance of full orbit = max distance / number of tested points 
  distance = distance between point z ( of orbit) and parabolic fixed point alfa
  average distance = (sum of max distances )/ NrOfPoints
   number of point used to compute average maximal distance of full orbit
   n_points_before + point_with_maximal_distance  + n_points_after  = 2*n+1  points used to compute average distance 

*/
double
GiveAverageMaxDistanceOfFullOrbit (double Zx0, double Zy0)
{

  double Distance2MaxLocal = 0.0;
  double D2;

  double Zx = Zx0;
  double Zy = Zy0;
  double Zx2, Zy2;
  // Z with max distance from fixed point
  double ZxMaxLocal;
  double ZyMaxLocal;

  //double NewZx, NewZy;
  complex double z;

  double SumOfDistances;	// used for computation of max aver distance

  double NumberOfPointsUsed = 2.0 * NrOfPoints + 1.0;
  int i;

  //
  D2 = (Zx - 0.5) * (Zx - 0.5) + Zy * Zy;	// fixed point z = 0.5

  //start with forward orbit 
  // orbit will allways fall into alfa !!!
  while (D2 > Distance2MaxLocal)
    {

      Distance2MaxLocal = D2;
      ZxMaxLocal = Zx;
      ZyMaxLocal = Zy;

      // then iterate 
      Zx2 = Zx * Zx;
      Zy2 = Zy * Zy;

      // new z : Z(n+1) = Zn * Zn  + C
      Zy = 2 * Zx * Zy + Cy;
      Zx = Zx2 - Zy2 + Cx;
      //
      D2 = (Zx - 0.5) * (Zx - 0.5) + Zy * Zy;

    }

  // go back to initial point
  // and check backward orbit ( maybe z0 is  after max distance)

  /* Zn*Zn=Z(n+1)-c */
  z = InverseIteration (Zx0 + Zy0 * I, c);
  Zx = creal (z);
  Zy = cimag (z);

  D2 = (Zx - 0.5) * (Zx - 0.5) + Zy * Zy;

  //
  while (D2 > Distance2MaxLocal)
    {
      Distance2MaxLocal = D2;
      ZxMaxLocal = Zx;
      ZyMaxLocal = Zy;

      // and make backward iterations

      /* Zn*Zn=Z(n+1)-c */
      z = InverseIteration (Zx + Zy * I, c);
      Zx = creal (z);
      Zy = cimag (z);

      D2 = (Zx - 0.5) * (Zx - 0.5) + Zy * Zy;

    }

  // z with max distance is found
  // compute average distance 

  SumOfDistances = sqrt (Distance2MaxLocal);

  // forward iteration of ZMaxLocal
  Zx = ZxMaxLocal;
  Zy = ZyMaxLocal;

  for (i = 0; i < NrOfPoints; i++)
    {

      // then iterate 
      Zx2 = Zx * Zx;
      Zy2 = Zy * Zy;

      // new z : Z(n+1) = Zn * Zn  + C
      Zy = 2 * Zx * Zy + Cy;
      Zx = Zx2 - Zy2 + Cx;
      //
      D2 = (Zx - 0.5) * (Zx - 0.5) + Zy * Zy;
      SumOfDistances += sqrt (D2);
    }

  // backward iteration of ZMaxLocal

  Zx = ZxMaxLocal;
  Zy = ZyMaxLocal;

  for (i = 0; i < NrOfPoints; i++)
    {
      /* Zn*Zn=Z(n+1)-c */
      /* Zn*Zn=Z(n+1)-c */
      z = InverseIteration (Zx + Zy * I, c);
      Zx = creal (z);
      Zy = cimag (z);
      D2 = (Zx - 0.5) * (Zx - 0.5) + Zy * Zy;

      SumOfDistances += sqrt (D2);
    }

  return (SumOfDistances / NumberOfPointsUsed);
}

unsigned char
GiveColorOfDistance (unsigned int ix, unsigned int iy)
{

  double Zx, Zy;
  double Distance;
  double fraction;
  unsigned char color;		// shade of gray

  // from screen to world coordinate 
  Zx = GiveZx (ix);
  Zy = GiveZy (iy);
  //
  //Distance = GiveMaxDistanceOfFullOrbit(Zx,Zy);
  Distance = GiveAverageMaxDistanceOfFullOrbit (Zx, Zy);
  Distance = Distance / DistanceMaxGlobal;	// normalisatio
  Distance = 55.0 * Distance;	// repetition of gradients 
  fraction = Distance - (int) Distance;
  color = 255 * fraction;

  if ((int) Distance % 2)
    color = 255 - color;	// smooth connection of 2 gradients 

  return color;

}

int
FillContour (double dXseed, double dYseed, unsigned char iNeededColor,
	     unsigned char A[])
{
  /* 
     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 A array : 1D array of 1-bit colors ( shades of gray)
     it does not check if index of A array is good  so memory error is possible 
   */

  int iXseed = (int) ((dXseed - ZxMin) / PixelWidth);
  int iYseed = (int) ((ZyMax - dYseed) / PixelHeight);	// (ZyMax - cimag(point))/PixelHeight; // inverse Y axis 

  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 A array */ ;
  int is;			// symmetric i 
  unsigned char iColor;

  /* ---------  move up --------------- */
  do
    {
      iX = iXmidLocal;
      i = Give_i (iX, iY); /* index of A array */ ;
      is = Give_is (iX, iY); /* index of A array */ ;
      // printf(" i =%d is = %d\n", i, is);
      /* move to right */
      while (A[i] == iNeededColor)
	{
	  iColor = GiveColorOfDistance (iX, iY);
	  A[i] = iColor;
	  A[is] = iColor;	// error , some lines are not drawe
	  //printf("found %d,%d \n", iX,iY);
	  iX += 1;
	  i = Give_i (iX, iY);
	  is = Give_is (iX, iY); /* index of A array */ ;
	  // printf(" i =%d is = %d\n", i, is);
	}
      iXmaxLocal = iX - 1;

      /* move to left */
      iX = iXmidLocal - 1;
      i = Give_i (iX, iY);
      is = Give_is (iX, iY); /* index of A array */ ;
      while (A[i] == iNeededColor)
	{
	  iColor = GiveColorOfDistance (iX, iY);
	  A[i] = iColor;
	  A[is] = iColor;
	  //printf("found %d,%d \n", iX,iY);
	  iX -= 1;
	  i = Give_i (iX, iY);
	  is = Give_is (iX, iY); /* index of A array */ ;
	  // printf(" i =%d is = %d\n", i, is);
	}
      iXminLocal = iX + 1;

      iY += 1;			/* move up */
      iXmidLocal = iXminLocal + (iXmaxLocal - iXminLocal) / 2;	/* new iX inside contour */
      i = Give_i (iXmidLocal, iY); /* index of A array */ ;
      is = Give_is (iXmidLocal, iY); /* index of A array */ ;
      if (A[i] == 0)
	break;			/*  it should not cross the border */

    }
  while (iY < iyMax);

  /* ------  move down ----------------- */
  iXmidLocal = iXseed;
  iY = iYseed - 1;

  do
    {
      iX = iXmidLocal;
      i = Give_i (iX, iY); /* index of A array */ ;
      is = Give_is (iX, iY); /* index of A array */ ;
      /* move to right */
      while (A[i] == iNeededColor)	/*  */
	{
	  iColor = GiveColorOfDistance (iX, iY);
	  A[i] = iColor;
	  A[is] = iColor;
	  //printf("found %d,%d \n", iX,iY);
	  iX += 1;
	  i = Give_i (iX, iY);
	  is = Give_is (iX, iY); /* index of A array */ ;
	}
      iXmaxLocal = iX - 1;

      /* move to left */
      iX = iXmidLocal - 1;
      i = Give_i (iX, iY);
      is = Give_is (iX, iY); /* index of A array */ ;
      while (A[i] == iNeededColor)	/*  */
	{
	  iColor = GiveColorOfDistance (iX, iY);
	  A[i] = iColor;
	  A[is] = iColor;
	  //printf("found %d,%d \n", iX,iY);
	  iX -= 1;
	  i = Give_i (iX, iY);
	  is = Give_is (iX, iY); /* index of A array */ ;
	}
      iXminLocal = iX + 1;

      iY -= 1;			/* move down */
      iXmidLocal = iXminLocal + (iXmaxLocal - iXminLocal) / 2;	/* new iX inside contour */
      i = Give_i (iXmidLocal, iY); /* index of A array */ ;
      is = Give_is (iXmidLocal, iY); /* index of A array */ ;
      if (A[i] == 0)
	break;			/*  it should not cross the border */
    }
  while (0 < iY);

  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, "%.1f", 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 ("DistanceMaxGlobal =  %.16f  \n", DistanceMaxGlobal);

  printf ("Image Width = %f \n", ZxMax - ZxMin);
  printf ("PixelWidth = %f \n", PixelWidth);
  printf ("Maximal number of iterations = iterMax = %ld \n", iterMax);
  printf ("ratio of image  = %f ; it should be 1.000 ...\n", ratio);
  return 0;
}

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

int
setup (int ParentPeriod, int ChildPeriod)
{

  printf ("setup\n");

  Cx = 0.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 = 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 ...

  DistanceMaxGlobal = GiveAverageMaxDistanceOfFullOrbit (0.0, +0.5);

  // for numerical optimisation in iteration
  ER2 = ER * ER;

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

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

  printf (" end of setup \n");

  return 0;

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

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

  printf ("components of Fatou set  : \n");
  ComputeFatouComponents (data, iterMax);
  SaveArray2PGMFile (data, NrOfPoints + 0.0);
  //
  ComputeZerosOfQnc (zero, 20);
  SaveArray2PGMFile (zero, NrOfPoints + 0.1);
  //  
  ComputeBoundariesIn (zero);	// from zero to edge
  SaveArray2PGMFile (edge, NrOfPoints + 0.2);
  // 
  CopyBoundariesTo (zero);	// copy from edge to zero
  SaveArray2PGMFile (zero, NrOfPoints + 0.3);
  //

  //
  FillContour (0.4, 0.4, 165, zero);
  // FillContour(0.4, -0.4, 215, zero); fill contour use symmetry
  SaveArray2PGMFile (zero, NrOfPoints + 0.6);

  printf (" allways free memory  to avoid buffer overflow \n");

  free (data);
  free (edge);
  free (zero);

  info ();

  return 0;
}