/*

  Adam Majewski
  adammaj1 aaattt o2 dot pl  // o like oxygen not 0 like zero 
  
  
  console program in c programing language 
  
  ==============================================
  
  
  ---------------------------------
  indent e.c 
  default is gnu style 
  -------------------



  c console progam 
  
		
  	gcc e.c -lm -Wall -march=native 
  	time ./a.out > b.txt


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


  time ./a.out

  time ./a.out >e.txt

  ----------------------
  
 real	0m19,809s
user	2m26,763s
sys	0m0,161s





gnuplot> set xlabel "distance"
gnuplot> set ylabel "distance"
gnuplot> set xlabel "turn"
gnuplot> set output "data.png"
gnuplot> plot "data.txt"
gnuplot> set title "distance from ray to the landing point after n iteration"
gnuplot> set output "data.png"
gnuplot> plot "data.txt"
gnuplot> unset key
gnuplot> set output "data.png"
gnuplot> plot "data.txt"



*/

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


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



// 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;	//  proportional to size of image file !!!!!!


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

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


double ZxMin ; //= -1.7; //-2.0;	//-0.05;
double ZxMax ; //=  1.7; // 2.0;	//0.75;
double ZyMin ;// //= -1.7; //-2.0;	//-0.1;
double ZyMax ; // =  1.7; //2.0;	//0.7;

double ImageWidth;

double plane_radius = 2.0 ;
double complex plane_center = 0.0;

static double PixelWidth;	// =(ZxMax-ZxMin)/ixMax;
static double PixelHeight;	// =(ZyMax-ZyMin)/iyMax;
static double ratio;


// complex numbers of parametr plane 
double complex c;		// parameter of function fc(z)=z^2 + c
double complex a; // alfa fixed point

int child_period = 11;
int t_numerator = 5;
int t_denominator; // = child period
int e_angle_numerator0 = 341;
int e_angle_denominator = 2047;


// for external rays


int bCounted = 0; // boolean
double MinDistanceToLandingPoint = 1.0; // initial value  
//int Period = 2;


unsigned long int iterMax_ray      = 2000000000000;	//it is also used for file name so it shlould be different values !!
unsigned long int iterMax_cr_orbit = 1000000000000;	//

/* colors = shades of gray from 0 to 255 */
unsigned char iColorOfExterior = 245;
unsigned char iColorOfInterior = 55;
unsigned char iColorOfInterior1 = 210;
unsigned char iColorOfInterior2 = 180;
unsigned char iColorOfBoundary = 0;
unsigned char iColorOfUnknown = 30;









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


int SetPlane(complex double center, double radius, double AspectRatio){


	ZxMin = creal(center) - radius*AspectRatio; //-2.0;	//-0.05;
  	ZxMax = creal(center) + radius*AspectRatio; // 2.0;	//0.75;
	ZyMin = cimag(center) - radius; //-2.0;	//-0.1;
	ZyMax = cimag(center) + radius; //2.0;	//0.7;
	ImageWidth = ZxMax - ZxMin;

	return 0;

}



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


/* 
   gives position ( index) in 1D virtual array  of 2D point Z 
   without bounds check !!
*/
int Give_i_from_d(complex double Z){ // double version of Give_k
  /* translate from world to screen coordinate */

  //  iY=(ZyMax-Zy)/PixelHeight; /*  */
  int ix=(creal(Z)-ZxMin)/PixelWidth;
  int iy=(ZyMax - cimag(Z))/PixelHeight; /* reverse Y  axis */		

	
  return Give_i(ix,iy);

}

 


void ColorPixel(int iColor, int i, unsigned char A[])
{
  A[i]   = iColor;
  
}




int  ColorPixel_d(complex double z, int iColor,  unsigned char A[]){
  int i = Give_i_from_d(z); // compute index of 1D array
  if ( i<0  || i>iSize) {fprintf(stderr, " bad i from color pixel\n");return 1;}
  ColorPixel(iColor, i, A);
  //printf("plot z = %f;%f ; i = %d\n",creal(z), cimag(z), i);
	
  return 0;

}






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

 /* 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(complex double point,unsigned char iColor, unsigned char A[] )
{

  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(ix, iy, iColor, A);
return 0;
}

/*
http://rosettacode.org/wiki/Bitmap/Bresenham%27s_line_algorithm
Instead of swaps in the initialisation use error calculation for both directions x and y simultaneously:
*/
void iDrawLine( int x0, int y0, int x1, int y1, unsigned char iColor, unsigned char A[]) 
{
  int x=x0; int y=y0;
  int dx = abs(x1-x0), sx = x0<x1 ? 1 : -1;
  int dy = abs(y1-y0), sy = y0<y1 ? 1 : -1; 
  int err = (dx>dy ? dx : -dy)/2, e2;

  for(;;){
    iDrawPoint(x, y, iColor, A);
    if (x==x1 && y==y1) break;
    e2 = err;
    if (e2 >-dx) { err -= dy; x += sx; }
    if (e2 < dy) { err += dx; y += sy; }
  }
}




int dDrawLineSegment(double Zx0, double Zy0, double Zx1, double Zy1, int color, unsigned char *array) 
{

 unsigned int ix0, iy0; // screen coordinate = indices of virtual 2D array 
 unsigned int ix1, iy1; // screen coordinate = indices of virtual 2D array

   // first step of clipping
   //if (  Zx0 < ZxMax &&  Zx0 > ZxMin && Zy0 > ZyMin && Zy0 <ZyMax 
    //  && Zx1 < ZxMax &&  Zx1 > ZxMin && Zy1 > ZyMin && Zy1 <ZyMax )
   ix0= (Zx0- ZxMin)/PixelWidth; 
   iy0 = (ZyMax - Zy0)/PixelHeight; // inverse Y axis 
   ix1= (Zx1- ZxMin)/PixelWidth; 
   iy1= (ZyMax - Zy1)/PixelHeight; // inverse Y axis 
   // second step of clipping
   if (ix0 >=ixMin && ix0<=ixMax && ix0 >=ixMin && ix0<=ixMax && iy0 >=iyMin && iy0<=iyMax 
      && iy1 >=iyMin && iy1<=iyMax )
   iDrawLine(ix0,iy0,ix1,iy1,color, array) ;

return 0;
}







// -------------------------------
// https://en.wikipedia.org/wiki/Cohen%E2%80%93Sutherland_algorithm

typedef int OutCode;

const int INSIDE = 0; // 0000
const int LEFT = 1;   // 0001
const int RIGHT = 2;  // 0010
const int BOTTOM = 4; // 0100
const int TOP = 8;    // 1000

// Compute the bit code for a point (x, y) using the clip rectangle
// bounded diagonally by (ZxMin, ZyMin), and (ZxMax, ZyMax)

// ASSUME THAT ZxMax, ZxMin, ZyMax and ZyMin are global constants.

OutCode ComputeOutCode(double x, double y)
{
	OutCode code;

	code = INSIDE;          // initialised as being inside of [[clip window]]

	if (x < ZxMin)           // to the left of clip window
		code |= LEFT;
	else if (x > ZxMax)      // to the right of clip window
		code |= RIGHT;
	if (y < ZyMin)           // below the clip window
		code |= BOTTOM;
	else if (y > ZyMax)      // above the clip window
		code |= TOP;

	return code;
}

// Cohen–Sutherland clipping algorithm clips a line from
// P0 = (x0, y0) to P1 = (x1, y1) against a rectangle with 
// diagonal from (xmin, ymin) to (xmax, ymax).
// CohenSutherlandLineClipAndDraw
void dDrawLine(double x0, double y0, double x1, double y1,unsigned char iColor, unsigned char A[])
{
	// compute outcodes for P0, P1, and whatever point lies outside the clip rectangle
	OutCode outcode0 = ComputeOutCode(x0, y0);
	OutCode outcode1 = ComputeOutCode(x1, y1);
	bool accept = false;

	while (true) {
		if (!(outcode0 | outcode1)) { // Bitwise OR is 0. Trivially accept and get out of loop
			accept = true;
			break;
		} else if (outcode0 & outcode1) { // Bitwise AND is not 0. (implies both end points are in the same region outside the window). Reject and get out of loop
			break;
		} else {
			// failed both tests, so calculate the line segment to clip
			// from an outside point to an intersection with clip edge
			double x, y;

			// At least one endpoint is outside the clip rectangle; pick it.
			OutCode outcodeOut = outcode0 ? outcode0 : outcode1;

			// Now find the intersection point;
			// use formulas y = y0 + slope * (x - x0), x = x0 + (1 / slope) * (y - y0)
			if (outcodeOut & TOP) {           // point is above the clip rectangle
				x = x0 + (x1 - x0) * (ZyMax - y0) / (y1 - y0);
				y = ZyMax;
			} else if (outcodeOut & BOTTOM) { // point is below the clip rectangle
				x = x0 + (x1 - x0) * (ZyMin - y0) / (y1 - y0);
				y = ZyMin;
			} else if (outcodeOut & RIGHT) {  // point is to the right of clip rectangle
				y = y0 + (y1 - y0) * (ZxMax - x0) / (x1 - x0);
				x = ZxMax;
			} else if (outcodeOut & LEFT) {   // point is to the left of clip rectangle
				y = y0 + (y1 - y0) * (ZxMin - x0) / (x1 - x0);
				x = ZxMin;
			}

			// Now we move outside point to intersection point to clip
			// and get ready for next pass.
			if (outcodeOut == outcode0) {
				x0 = x;
				y0 = y;
				outcode0 = ComputeOutCode(x0, y0);
			} else {
				x1 = x;
				y1 = y;
				outcode1 = ComputeOutCode(x1, y1);
			}
		}
	}
	if (accept) {
			
               // printf(	"x0 = %d, y0 = %d, x1 = %d, y1 =%d \n",x0, y0, x1, y1);			
		dDrawLineSegment(x0, y0, x1, y1, iColor, A);
	}
}

// -----------------------------------------------------------------------






int DrawClippedLineSegment(complex double z0, complex double z1, unsigned char iColor, unsigned char A[]){


	double x0 = creal(z0);
	double y0 = cimag(z0);
	double x1 = creal(z1);
	double y1 = cimag(z1);
	
	dDrawLine( x0, y0, x1, y1, iColor, A);
	return 0; 



}



































//------------------mapping between world and integer coordinate -----------------------------------------------------





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

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


complex double GiveZ( int ix, int iy){
  double Zx = GiveZx(ix);
  double Zy = GiveZy(iy);
	
  return Zx + Zy*I;
	
	


}




// ****************** DYNAMICS = trap tests ( target sets) ****************************



// bailout test
// z escapes when 
// abs(z)> ER or cabs2(z)> ER2 
// https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Julia_set#Boolean_Escape_time

//int Escapes(complex double z){
 // here target set (trap) is the exterior  circle with radsius = ER ( EscapeRadius) 
  // with ceter = origin z= 0
  // on the Riemann sphere it is a circle with point at infinity as a center  
   
  //if (cabs(z)>ER) return 1;
  //return 0;
//}





// compute alfa fixed point
// https://en.wikipedia.org/wiki/Periodic_points_of_complex_quadratic_mappings#Period-1_points_(fixed_points)
complex double GiveAlfa(complex double c)
{
	// d=1-4c 
	// alfa = (1-sqrt(d))/2 
	return (1.0-csqrt(1.0 - 4.0*c))/2.0 ;

}





 
 
// *************************************************************************************************
// ********************************* critical orbit ************************************************
// *********************************************************************************************

// fill array 
// uses global var :  ...

int DrawImage_CriticalOrbit (unsigned char A[])
{
	fprintf(stderr, "draw critical orbit \n");
  int i = 0; // iteration = number of the point
  int iMax = child_period*iterMax_cr_orbit;
  complex double z = 0.0;
  
   
	
  
  	for (i = 0; i < iMax; ++i){
  		fprintf(stderr, "\t%d / %d\r", i, iMax); // info
  		ColorPixel_d(z, 255, A );	//  draw point and check if point is outside image 
  		z = z*z+c; // forward iteration : complex quadratic polynomial
  }
  	
  return 0;
}

 
 
 
// fill array 
// uses global var :  ...
int DrawImage_CriticalOrbit_Approx (unsigned char A[])
{
	fprintf(stderr, "draw critical orbit \n");
  int i = 0; // iteration = number of the point
  int iMax = child_period*iterMax_cr_orbit;
  complex double z = 0.0;
  
   
	
  
  	for (i = 0; i < iMax; ++i){
  		fprintf(stderr, "\t%d / %d\r", i, iMax); // info
  		ColorPixel_d(z, 255, A );	//  draw point and check if point is outside image 
  		z = z*z+c; // forward iteration : complex quadratic polynomial
  }
  	// aproximation and info
  	double t = ((double)e_angle_numerator0)/e_angle_denominator; // first external angle in turns	
	for (i = 0; i < child_period; ++i){
		printf (" i = %d \t t = %.16f \tz = %.16f %+.16f*I \t distance to fixed point = %.16f = %.16f*PixelWidth = %.16f*ImageWidth\n ", i, t, creal(z), cimag(z), cabs(a - z), cabs(a-z)/PixelWidth, cabs(a-z)/ImageWidth );
		// approximate when landing point is known
      		DrawClippedLineSegment(z, a, 155, A) ;
      		// next angle 
  		t *= 2.0; // t = 2*t angle doubling map
      		if (t > 1.0) {t--;}  // t = t modulo 1 
      		// next z
		z = z*z+c; // forward iteration : complex quadratic polynomial
		}

  return 0;
}

 

 
  
// *************************************************************************************************
// ********************************* external ray ************************************************
// *********************************************************************************************



 
// check if period p of n/d is proper under period doubling map   
  // p times (2*n) mod d  should give n  
 
int TestPeriod(int n, int d, int period){
	int p;
	int pMax = period;
	int nn = n; // termporary value
	
	if (n<=0 || d<=0 || period <= 0) {
		fprintf(stderr, " input should be positive \n");
		return 1;
	
	}
	
	for(p=0; p<pMax; p++) 
		nn *=2;
	nn = nn % d; // remainder 
	if (nn-n ==0) 
		return 1; // true
	return 0; // false



}
 
 


// dot product 

double dot(complex double a, complex double b){

	// ax*bx + ay*by
	return creal(a)*creal(b) + cimag(a)*cimag(b);


}


int InverseIterationAndDrawSegment (int per, complex double zz[per], unsigned char A[]){
	
	int p;
	int pMax = per;
	complex double z;
	complex double zn; // next vsalue : zn = f^(-1)(z)
	
	
	
	// compute preimages
	for(p=0; p<pMax; ++p){
		z = zz[(p+1) % per]; //read point from the ray p+1 
		zn = csqrt(z-c); // inverse iteration = preimage;  fc(z) maps  z_{l,j} to z_{l-1,j+1}} so f^{-1} map to z_{l+1, j-1}
		
		// choose correct preimage using dot product 
		if (dot(zz[p],zn) < 0 ) zn = - zn; 
		
		// draw segment of ray p
		DrawClippedLineSegment(zz[p],zn, 255, A);
		zz[p] = zn; //save 
	}
	
	
	
	
		
	
  	return 0; 


}


// uses global var :  ...
/*
"In the dynamic plane, external rays can be drawn by backwards iteration. It is most effective for a periodic or preperiodic angle.

You must keep track of points on the finite collection of rays with angles phi ,2phi ,4phi ...

Say z_{l,j} corresponds to a radius  R^{1/(2l)}  and the angle  2j*phi .

Then  fc(z) maps  z_{l,j} to z_{l-1,j+1}}

This point, which was constructed before, has two preimages under {\displaystyle f_{c}(z)} {\displaystyle f_{c}(z)} .

The one that is closer to {\displaystyle z_{l-1,j}} {\displaystyle z_{l-1,j}} is the correct one. 
This criterion was proved by Thierry Bousch. The ray will look better when you introduce intermediate points." Wolf Jung


this procedure works only for periodic angles 

*/ 
int DrawExternalDynamicRaysBI (int n, int m, int period, int iMax, unsigned char A[])
{ // In the dynamic plane, external rays can be drawn by backwards iteration. This procedure is only for the periodic angle
  // https://commons.wikimedia.org/wiki/File:Backward_Iteration.svg
  fprintf(stderr, "draw external ray \n");
  
  int i = 0; // iteration = number of points
  //int iMax = 10000000;
  double r = 10000.0; // very big radius = near infinity where  z=w so one can swith to dynamical plane ( Boettcher conjugation )
  complex double zz[period]; // zz is an array of z  
  
  double t;
  int p;
  int pMax = period; // number of rays to draw 
  int result ;
  
  
  // check if period is proper 
  // p times (2*n) mod m  should give n 
  if (TestPeriod(n, m, period)!=1 ) {
  	fprintf(stderr, "bad input to the TestPeriod from external ray : error \n");
  	return 1;
  }
  fprintf(stderr, "good input to the TestPeriod from external ray \n");
  
  
  // to do 
  // find landing point  : periodic point  
  // compute how many times ray rotates around it's landing point 
  // compute distance to tha landing point of last point of the ray ( along the ray ? = infinity  so to the circle with fixed radius around landing point) 
  
  
  
  
  
  // initial points  on rays
  t = (double)n/m; // first external angle in turns	
  for(p=0; p<pMax; ++p){
  	// initial point
  	zz[p] =  r*cexp(2.0*I * M_PI * t ); // Euler's formula
  	// next angle 
  	t *= 2.0; // t = 2*t angle doubling map
      	if (t > 1.0) t--;  // t = t modulo 1 
  
  
  }
  
  
	for (i = 0; i < iMax; ++i){
    		// inverse iteration of  complex quadratic polynomial:  z = csqrt(z-c)  with proper choose of preimage for one point on every ray 
    		fprintf(stderr, "\t%d / %d\r", i, iMax); // info
    		result = InverseIterationAndDrawSegment(period, zz, A);
    		if (result != 0) {fprintf (stderr, " lost precision\n"); return 1;} 
    		}
  			
	fprintf(stderr, "end of drawing external ray \n");
  return 0;
}





// uses global var :  ...
/*
"In the dynamic plane, external rays can be drawn by backwards iteration. It is most effective for a periodic or preperiodic angle.

You must keep track of points on the finite collection of rays with angles phi ,2phi ,4phi ...

Say z_{l,j} corresponds to a radius  R^{1/(2l)}  and the angle  2j*phi .

Then  fc(z) maps  z_{l,j} to z_{l-1,j+1}}

This point, which was constructed before, has two preimages under {\displaystyle f_{c}(z)} {\displaystyle f_{c}(z)} .

The one that is closer to {\displaystyle z_{l-1,j}} {\displaystyle z_{l-1,j}} is the correct one. 
This criterion was proved by Thierry Bousch. The ray will look better when you introduce intermediate points." Wolf Jung


this procedure works only for periodic angles 

*/ 
int DrawExternalDynamicRaysBI_approx (int n, int m, int period, int iMax, unsigned char A[])
{ // In the dynamic plane, external rays can be drawn by backwards iteration. This procedure is only for the periodic angle
  // https://commons.wikimedia.org/wiki/File:Backward_Iteration.svg
  fprintf(stderr, "draw external ray \n");
  
  int i = 0; // iteration = number of points
  //int iMax = 10000000;
  double r = 10000.0; // very big radius = near infinity where  z=w so one can swith to dynamical plane ( Boettcher conjugation )
  complex double zz[period]; // zz is an array of z  
  
  double t;
  int p;
  int pMax = period; // number of rays to draw 
  int result ;
  
  
  // check if period is proper 
  // p times (2*n) mod m  should give n 
  if (TestPeriod(n, m, period)!=1 ) {
  	fprintf(stderr, "bad input to the TestPeriod from external ray : error \n");
  	return 1;
  }
  fprintf(stderr, "good input to the TestPeriod from external ray \n");
  // to do 
  // find landing point  : periodic point  
  // compute how many times ray rotates around it's landing point 
  // compute distance to tha landing point of last point of the ray ( along the ray ? = infinity  so to the circle with fixed radius around landing point) 
  
  
  
  
  
  // initial points  on rays
  t = (double)n/m; // first external angle in turns	
  for(p=0; p<pMax; ++p){
  	// initial point
  	zz[p] =  r*cexp(2.0*I * M_PI * t ); // Euler's formula
  	// next angle 
  	t *= 2.0; // t = 2*t angle doubling map
      	if (t > 1.0) t--;  // t = t modulo 1 
  
  
  }
  
  
	for (i = 0; i < iMax; ++i){
    		// inverse iteration of  complex quadratic polynomial:  z = csqrt(z-c)  with proper choose of preimage for one point on every ray 
    		fprintf(stderr, "\t%d / %d\r", i, iMax); // info
    		result = InverseIterationAndDrawSegment(period, zz, A);
    		if (result != 0) {fprintf (stderr, " lost precision\n"); return 1;} 
    		}
  			

	// print and approximate if one knows landing point !
	printf("last point of the ray\n");
  	t = (double)n/m; // first external angle in turns	
	for(p=0; p<pMax; ++p){
		printf (" p = %d \t t = %.16f \tz = %.16f %+.16f*I \t distance to landing point = %.16f = %.16f* PixelWidth = %.16f*ImageWidth\n ", p, t, creal(zz[p]), cimag(zz[p]), cabs(a - zz[p]), cabs(a- zz[p])/PixelWidth, cabs(a- zz[p])/ImageWidth  );
		// next angle 
  		t *= 2.0; // t = 2*t angle doubling map
      		if (t > 1.0) {t--;}  // t = t modulo 1 
      		// approximation when landing point is known
    		DrawClippedLineSegment(zz[p], a, 155, A) ;
		}
		
	  fprintf(stderr, "end of drawing external ray \n");
  return 0;
}

 
 
 
 
 






// *******************************************************************************************
// ********************************** save A array to pgm file ****************************
// *********************************************************************************************

int SaveArray2PGMFile( unsigned char A[], double k, char* comment )
{
  
  FILE * fp;
  const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
  char name [200]; /* name of file */
  snprintf(name, sizeof name, "z%.0f", k); /*  */
  char *filename =strcat(name,".pgm");
  
  
  
  // 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 array with image data bytes to the file in one step 
  fclose(fp); 
  
  // info 
  fprintf(stderr, "File %s saved ", filename);
  if (comment == NULL || strlen(comment) ==0)  
    fprintf(stderr,"\n");
  else fprintf (stderr, ". Comment = %s \n", comment); 

  return 0;
}







int PrintInfoAboutProgam()
{

  
  // display info messages
  printf ("\n\nNumerical approximation of Julia set for fc(z)= z^2 + c \n");
  //printf ("iPeriodParent = %d \n", iPeriodParent);
  //printf ("iPeriodOfChild  = %d \n", iPeriodChild);
  printf ("parameter c = ( %.16f ; %.16f ) \n", creal(c), cimag(c));
  printf ("fixed point alfa z = a = ( %.16f ; %.16f ) \n", creal(a), cimag(a));
  
  
  printf ("iSize = %d\n", iSize);
  printf ("Image Width = %.16f in world coordinate\n", ZxMax - ZxMin);
  // https://en.wikipedia.org/wiki/Aspect_ratio_(image)
  printf ("Image aspect ratio = %.16f\n", ratio);
  printf ("PixelWidth = %.16f \n", PixelWidth);
  
  printf ("\nexternal ray \n");
  printf ("Maximal number of iterations = iterMax_ray = %ld \n", iterMax_ray);
  printf ("Maximal number of iterations = iterMax_cr_orbit = %ld \n", iterMax_cr_orbit);
  printf ("MinDistanceToLandingPoint = %.16f \n", MinDistanceToLandingPoint);
 
  
  
  printf ("\n plane : \n"); 
  printf (" radius = %.16f \n", plane_radius);
  printf (" center = z = %.16f %+.16f*I )\n ", creal(plane_center), cimag(plane_center) );
  // image corners in world coordinate
  // center and radius
  // center and zoom
  // GradientRepetition
  
  printf ("ratio of image  = %f ; it should be 1.000 ...\n", ratio);
  //
  printf("gcc version: %d.%d.%d\n",__GNUC__,__GNUC_MINOR__,__GNUC_PATCHLEVEL__); // https://stackoverflow.com/questions/20389193/how-do-i-check-my-gcc-c-compiler-version-for-my-eclipse
  // OpenMP version is diplayed in the console 
  return 0;
}






// *****************************************************************************
//;;;;;;;;;;;;;;;;;;;;;;  setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
// **************************************************************************************

int setup ()
{

  fprintf (stderr, "setup start ");
  c = -0.690059870015044  +0.276026482784614*I; //    t = 5/11
  a = GiveAlfa(c);
  t_denominator = child_period;
  plane_center = a;
  plane_radius = 0.9;
  
	
  /* 2D array ranges */
  
  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].
  
  ratio = (double) iWidth / iHeight;
  
  SetPlane(plane_center , plane_radius, ratio); // 


  /* Pixel sizes */
  PixelWidth = (ZxMax - ZxMin) / ixMax;	//  ixMax = (iWidth-1)  step between pixels in world coordinate 
  PixelHeight = (ZyMax - ZyMin) / iyMax;
  
  
  //((ZxMax - ZxMin) / (ZyMax - ZyMin)) / ((double) iWidth / (double) iHeight);	// it should be 1.000 ...
	
   
	
  
  //ER2 = ER * ER; // for numerical optimisation in iteration
  //lnER = log(EscapeRadius); // ln(ER) 
  
   	
  /* create dynamic 1D arrays for colors ( shades of gray ) */
  data = malloc (iSize * sizeof (unsigned char));
  edge = malloc (iSize * sizeof (unsigned char));
  edge2 = malloc (iSize * sizeof (unsigned char));
  	
  if (data == NULL || edge == NULL || edge2 == NULL){
    fprintf (stderr, " Could not allocate memory");
    return 1;
  }

  
 	
  
  //BoundaryWidth = 6.0*iWidth/2000.0  ; //  measured in pixels ( when iWidth = 2000) ; such function is stable when iWidth is changing
  //distanceMax = BoundaryWidth*PixelWidth; // distance to the boundary from exterior
  
  //InternalSiegelDiscRadius = GiveInternalSiegelDiscRadius(c,a);
  //ExternalSiegelDiscRadius = GiveExternalSiegelDiscRadius(c,a);
  
  
  fprintf (stderr, "and end\n");
	
  return 0;

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




int end(){


  // printf (" allways free memory (deallocate )  to avoid memory leaks \n"); // https://en.wikipedia.org/wiki/C_dynamic_memory_allocation
  free (data);
  free(edge);
  free(edge2);
  PrintInfoAboutProgam();
  return 0;

}

// ********************************************************************************************************************
/* -----------------------------------------  main   -------------------------------------------------------------*/
// ********************************************************************************************************************

int main () {
  setup ();
  
  
  
  DrawImage_CriticalOrbit_Approx(data);
  SaveArray2PGMFile (data, child_period*iterMax_cr_orbit, "critical orbit");
  
  DrawExternalDynamicRaysBI_approx(e_angle_numerator0,e_angle_denominator,child_period, child_period*iterMax_ray, data);
  SaveArray2PGMFile (data, child_period*iterMax_ray, "ray");
  
  
  
  
  end();

  return 0;
}