/*


  here are:
  * 1 critical point  z=0.0
  * 1 weakly attracting cycle period 4
  

 

  Adam Majewski
  adammaj1 aaattt o2 dot pl  // o like oxygen not 0 like zero 
  
  
  
  Structure of a program or how to analyze the program 
  
  
  ============== Image X ========================
  
  DrawImageOf -> DrawPointOf -> ComputeColorOf ( FunctionTypeT FunctionType , complex double z) -> ComputeColor
  
  
  check only last function  which computes color of one pixel for given Function Type
  
  

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

  
  ---------------------------------
  indent d.c 
  default is gnu style 
  -------------------



  c console progam 
  
  export  OMP_DISPLAY_ENV="TRUE"	
  gcc j.c -lm -Wall -march=native -fopenmp
  time ./a.out > j.txt


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


  time ./a.out

  time ./a.out >i.txt
  time ./a.out >e.txt
  
  make
  
  
  
  
  
  
  convert -limit memory 1000mb -limit disk 1gb dd30010000_20_3_0.90.pgm -resize 2000x2000 10.png

  
  
  
*/

#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
#include <limits.h>		// Maximum value for an unsigned long long int



// https://sourceforge.net/p/predef/wiki/Standards/

#if defined(__STDC__)
#define PREDEF_STANDARD_C_1989
#if defined(__STDC_VERSION__)
#if (__STDC_VERSION__ >= 199409L)
#define PREDEF_STANDARD_C_1994
#endif
#if (__STDC_VERSION__ >= 199901L)
#define PREDEF_STANDARD_C_1999
#endif
#endif
#endif




/* --------------------------------- global variables and consts ------------------------------------------------------------ */
#define VERSION 20230217

int NumberOfImages = 0;




// virtual 2D array and integer ( screen) coordinate
// Indexes of array starts from 0 not 1 
//unsigned int ix, iy; // var

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

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




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

// memmory 1D array for 8 bit color 
unsigned char *data;
unsigned char *edge;
unsigned char *edge2;
//unsigned char *edge2;






// unsigned int i; // var = index of 1D array
//static 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



//FunctionType = representing functions
typedef enum  {FatouBasins = 0, FatouComponents = 2,  LSM = 3, LSM_m = 4, Unknown = 5 , BDM = 6, MBD = 7 , MBD2 = 8, SAC = 9, DLD = 10, ND = 11, NP= 12, POT = 13 , Blend = 14, DEM = 15, BDM_LSM = 16
		
} FunctionTypeT; 


typedef enum  {superattracting = 100, attracting = 200, parabolic = 300,  repelling = 400
		
} DynamicTypeT; 

// Dynamic type of the interior, exterior is allways superattracting 
DynamicTypeT DynamicType = parabolic; // it is set manually , check also main and setup if there are other changes


// see ComputeColor_FunctionType_DynamicType below



// see SetPlane

double radius = 1.2; 
complex double center = 0.0 ;
double  DisplayAspectRatio  = 1.0; // https://en.wikipedia.org/wiki/Aspect_ratio_(image)
// dx = dy compare setup : iWidth = iHeight;
double ZxMin; //= -1.3;	//-0.05;
double ZxMax;// = 1.3;	//0.75;
double ZyMin;// = -1.3;	//-0.1;
double ZyMax;// = 1.3;	//0.7;
double PixelWidth;	// =(ZxMax-ZxMin)/ixMax;
double PixelHeight;	// =(ZyMax-ZyMin)/iyMax;

// dem
double BoundaryWidth ; //= 1.0*iWidth/2000.0  ; //  measured in pixels ( when iWidth = 2000) 
double distanceMax ; //= BoundaryWidth*PixelWidth;


double ratio; 


double ER;
double ER2;			//= 1e60;
double AR; // bigger values do not works
double AR2;
complex double parabolic_trap_center;
double PixelWidth2;


int IterMax = 10000000;
int IterMax_LSM = 10000000;
int IterMax_DEM = 10000000;

/* colors = shades of gray from 0 to 255 */
unsigned char iColorOfBasin1 = 170;
unsigned char iColorOfInterior = 150;
unsigned char iColorOfExterior = 225;

unsigned char iColorOfBoundary = 0;
unsigned char iColorOfUnknown = 5;

// pixel counters; pixel counters work not good with OpenMP !!!
unsigned long long int uUnknown = 0;
unsigned long long int uInterior = 0;
unsigned long long int uExterior = 0;




// const int period = 1;
const int period_parent = 1;
#define period_child 4 
int internal_angle_numerator = 1;
int internal_angle_denominator ; // = period_child in setup

/* critical point */
complex double zcr = 0.0; //
complex double c =  0.2500000000000000 +0.5000000000000000*I  ; //1/4
complex double zprecr = -0.2288905993372869 -0.0151096456992677*I;





// 1 superattracting cycle period 1
double delta; // delta is a distance between fixed points
/*
alfa < 1/2 <beta
cabs(beta - alfa) = delta

alfa = ( 1 - delta)/2
beta = ( 1 + delta)/2
delta = sqrt(1- 4c)


*/
complex double zp = 0.0000000000000000 +0.5000000000000000*I ;  // one of periodic pointes , alfa , z_fixed
complex double zcr_last; //
complex double zprep = -0.228890599337288  -0.015109645699267*I; //  -0.4227200833935814 -0.0286616828394151*I ; // z_prefixed
complex double zpp[period_child]; 



/*
 for parabolic triangle trap ( computed in other program: last point of first and last ray ( ordered by t )
   
  triangle trap for components are different ( bigger) then triangle trap for BDM !!!!
 
  circle around fixed point zp = trap for interior
  biggest triangle (zp, zprep, izpprep) = trap for components
  2 smaller triangles (zp, zprecr, zcr) and (zp, zcr, -zprecr) = traps for BD
  
   
  		-z_prep
  		-zprecr
  zp		zcr 
  		zprecr
  		zprep
 */
complex double z1 =       -0.3350594976618946 +0.1652825185133397*I; //-0.4170188752901364 +0.1506018053351233*I ;

complex double z3 =       -0.3512112430996852 +0.0760154138719138*I ; // -0.4223414442666910 +0.1139242453319452*I;
complex double z13;

// for MBD
static double TwoPi=2.0*M_PI; // texture
double t0 ; // manually tuned t for MBD
// see https://www.youtube.com/watch?v=JttLtB0Gkdk&t=894s
// 

// update with f function 
const char *f_description = "Numerical approximation of Julia set for f(z)= z^2 + c "; // without /n !!!
/* ------------------------------------------ functions -------------------------------------------------------------*/

// main complex function f 
// upadte f_description also
complex double f(const complex double z0) {

  double complex z = z0;
  z = z*z + c;
  return  z;
}
	


double c_arg(complex double z)
{
 double arg;
 arg = carg(z);
 if (arg<0.0) arg+= TwoPi ; 
 return arg; 
}

double c_turn(complex double z)
{
 double arg;
 arg = c_arg(z);
 return arg/TwoPi; 
}


double cabs2(complex double z){

  return creal(z)*creal(z)+cimag(z)*cimag(z);


}



int is_z_outside_viewport(complex double z){

  if (creal(z) >ZxMax ||
      creal(z) <ZxMin ||
      cimag(z) >ZyMax ||
      cimag(z) <ZyMin)
    {return 1; } // is outside = true
      
    
  return 0; // is inside = false



}


complex double GiveFixed(complex double c){
/* 

Equation defining fixed points : z^2-z+c = 0
	z*2+c = z
	z^2-z+c = 0

coefficients of standard form ax^2+ bx + c  
 	a = 1 , b = -1 , c = c
 
The discriminant d is 

	d=b^2- 4ac 
	d = 1 - 4c
	
 
 alfa =  (1-sqrt(d))/2 
*/

	complex double d = 1-4*c;
	complex double z = (1-csqrt(d))/2.0;
	return z;

}






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

}








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



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

// -------------------- parabolic triangle trap =---------------------------------------------------------- 

double IsTriangleCounterclockwise(double xa, double ya, double xb, double yb, double xc, double yc)
{return  ((xb*yc + xa*yb +ya*xc) - (ya*xb +yb*xc + xa*yc)); }

int DescribeTriangle(double xa, double ya, double xb, double yb, double xc, double yc)
{
 double t = IsTriangleCounterclockwise( xa,  ya, xb,  yb,  xc,  yc);
 
 if (t>0)  printf("this triangle is oriented counterclockwise,     determinent = %f \n", t);
 if (t<0)  printf("this triangle is oriented clockwise,            determinent = %f\n", t);
 if (t==0) printf("this triangle is degenerate: colinear or identical points, determinent = %f\n", t);

 return 0;
}


int Describe_Triangle(const double complex z1, const double complex z2, const double complex z3){

	return DescribeTriangle( creal(z1), cimag(z1), creal(z2), cimag(z2), creal(z3), cimag(z3));

}


//  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.
double side (double  x1, double y1, double x2,double y2,double x3, double y3)
{
    return (x1 - x3) * (y2 - y3) - (x2 - x3) * (y1 - y3);
}





// the triangle node numbering is counter-clockwise / clockwise  
int  PointInTriangle (double x, double y, double x1, double y1, double x2, double y2, double x3, double y3)
{
    int  b1, b2, b3;

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

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


int PointInTriangle_c(const complex double z, const complex double z1, const complex double z2, const complex double z3 ){

	return PointInTriangle (creal(z), cimag(z), creal(z1), cimag(z1), creal(z2), cimag(z2), creal(z3), cimag(z3));
}

// uses global var: z1, zp, z3 describing triangle(z1, zp, z3) 
int IsPointIn_Parabolic_Components_Triangle(const double complex z){

	return PointInTriangle_c(z, zp, zprecr, -zprecr);

}



int IsPointIn_BDM_Triangle_B(const double complex z){

	return PointInTriangle_c(z, zp, zcr, -zprecr);

}

int IsPointIn_BDM_Triangle_mB(const double complex z){

	return PointInTriangle_c(z, -zp, zcr, zprecr);

}


int IsPointIn_BDM_Triangle_A(const double complex z){

	return PointInTriangle_c(z, zp, zprecr, zcr);

}

int IsPointIn_BDM_Triangle_mA(const double complex z){

	return PointInTriangle_c(z, -zp, -zprecr, zcr);

}

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


int IsPointIn_BDM_Triangle_A_mA(const double complex z){

	//return (IsPointIn_BDM_Triangle_A(z) + IsPointIn_BDM_Triangle_mA(z));
	 int t = IsPointIn_BDM_Triangle_A(z); 
	 if (t) 
	 	{return t;}
	 	else { return IsPointIn_BDM_Triangle_mA(z);}
	
	
}

int IsPointIn_BDM_Triangle_B_mB(const double complex z){

	// return (IsPointIn_BDM_Triangle_B(z) + IsPointIn_BDM_Triangle_mB(z));
	int t = IsPointIn_BDM_Triangle_B(z); 
	 if (t) 
	 	{return t;}
	 	else { return IsPointIn_BDM_Triangle_mB(z);}
	
	
}

/*
 for parabolic triangle trap ( computed in other program: last point of first and last ray ( ordered by t )
   
  triangle trap for components are different ( bigger) then triangle trap for BDM !!!!
 
  circle around fixed point zp = trap for interior
  biggest triangle (zp, zprep, izpprep) = trap for components
  2 smaller triangles (zp, zprecr, zcr) and (zp, zcr, -zprecr) and it's preimages'= traps for BD
  
   
  reflection around origin: 
  		
  		-zprecr
  zp		zcr 		-zp
  		zprecr
  		
 */
















// circular trap with zp as a center

int IsInside_basin_trap_i(int ix, int iy){


  complex double z = GiveZ(ix, iy);
  if (  cabs2(z - zp) < AR2 )
    {return 1;}
  return 0;



}

int IsInside_basin_trap(complex double z){
if (  cabs2(z - zp) < AR2 )
    {return 1;}
  return 0;
}








/*


**************************************************  ComputeColor_FunctionType_DynamicType *********************************************

Make ComputeColor_FunctionType_DynamicType function for each combination of 2 enums
* FunctionTypeT
* DynamicTypeT


then update Compute8BitColor procedure : 


	 case FunctionType + DynamicType: {ComputeColor_FunctionType_DynamicType(z); break;}


run procedure inside MakeImages using DrawImage (array, FunctionType);	 
second enum DynamicType is updated manually inside main function

********************************************************************************************************************************************

*/






// ********************************************************************************************************************
/* ---------------------FatouBasins -----------------------------------------------------------*/
// ********************************************************************************************************************



/*
  
*/

unsigned char ComputeColor_FatouBasins_superattracting (complex double z)
{



	double cabs2z;
	
  int i;			// number of iteration
  for (i = 0; i < IterMax; ++i)
    {

	// 2 basins, both superattracting
	cabs2z = cabs2(z);
	// infinity is superattracting here !!!!!	
       if ( cabs2z > ER2 ){ return iColorOfExterior;}
	
      // second superAttraction basins. Here one of superattracting point is z= 0 
      if ( cabs2z < AR2 ){ uInterior += 1;return iColorOfInterior;}
		 
      			
     
      z = f(z);		//  iteration: z(n+1) = f(zn)
	

    }

  
  return iColorOfUnknown;


}


/*
  2 basins 
  
  - exterior
  - interior
  - unknown ( possibly empty set ) 

pixel counters work not good with OpenMP !!!

*/

unsigned char ComputeColor_FatouBasins_attracting (complex double z)
{



	
	
  int i;			// number of iteration
  for (i = 0; i < IterMax_LSM; ++i)
    {


		
      // infinity is superattracting here !!!!!	
       if ( cabs2(z) > ER2 ){ uExterior +=1;    return iColorOfExterior;}
	
      // 1 Attraction basins 
      if ( cabs2(zp-zpp[0]) < AR2 ){ uInterior += 1; return iColorOfInterior;}
		 
      			
     
      z = f(z);		//  iteration: z(n+1) = f(zn)
	

    }

  uUnknown += 1;
  return iColorOfUnknown;


}


// basin works bad for parabolic !!!!!!!
unsigned char ComputeColor_FatouBasins_parabolic (complex double z)
{



	
	
  int i;			// number of iteration
  for (i = 0; i < IterMax; ++i)
    {


		
      // infinity is superattracting here !!!!!	
       if ( cabs2(z) > ER2 ){ return iColorOfExterior;}
	
	// parabolic basins 
	if ( cabs2(zp-z) < AR2)
		{ uInterior += 1; return iColorOfInterior;}
		 
      z = f(z);		//  iteration: z(n+1) = f(zn)
	

    }

  uUnknown += 1;
  return iColorOfUnknown;


}


unsigned char ComputeColor_FatouBasins_repelling (complex double z)
{



	
	
  int i;			// number of iteration
  for (i = 0; i < IterMax; ++i)
    {


		
      // infinity is superattracting here !!!!!	
       if ( cabs2(z) > ER2 ){ return iColorOfExterior;}
	
      
      			
     
      z = f(z);		//  iteration: z(n+1) = f(zn)
	

    }

  uUnknown += 1;
  return iColorOfUnknown;


}




// ********************************************************************************************************************
/* ---------------------FatouComponents -----------------------------------------------------------*/
// ********************************************************************************************************************

unsigned char ComputeColor_FatouComponents_superattracting (complex double z)
{



	
	
  double cabs2z;


  int i;			// number of iteration
  for (i = 0; i < IterMax; ++i)
    {

	cabs2z = cabs2(z);
      // first superAttraction basin : infinity is superattracting for all polynomials
       if ( cabs2z > ER2 ){ return iColorOfExterior;}
	
  
      //second  superattraction basins 
      if ( cabs2z < AR2 ){ return iColorOfBasin1 - (i % period_child)*40;}
	 
		
     
      z = f(z);		//  iteration: z(n+1) = f(zn)
	

    }

  
  return iColorOfUnknown;


}



unsigned char ComputeColor_FatouComponents_attracting (complex double z)
{

  int i;			// number of iteration
  for (i = 0; i < IterMax; ++i)
    {


      // first basin is superAttraction: infinity is superattracting for all polynomials
       if ( cabs2(z) > ER2 ){ return iColorOfExterior;}
	
  
      //1 Attraction basins 
      if ( cabs2(zp-z) < AR2 ){ return iColorOfBasin1 - (i % period_child)*20;}
	 
		
     
      z = f(z);		//  iteration: z(n+1) = f(zn)
	

    }

  
  return iColorOfUnknown;


}

unsigned char ComputeColor_FatouComponents_parabolic (complex double z)
{

  int i;			// number of iteration
  for (i = 0; i < IterMax; ++i)
    {


      // superattracting basin ( infinity is superattracting here  )
       if ( cabs2(z) > ER2 ){ return iColorOfExterior;}
	
  
      // parabolic basins 
	if (IsPointIn_Parabolic_Components_Triangle(z) )
		{ return iColorOfBasin1 - (i % period_child)*11;}
	 
		
     
      z = f(z);		//  iteration: z(n+1) = f(zn)
	

    }

  
  return iColorOfUnknown;


}



unsigned char ComputeColor_FatouComponents_repelling (complex double z)
{


  int i;			// number of iteration
  for (i = 0; i < IterMax; ++i)
    {


		
      // there is onlt one component = basin of infinity, infinity is superattracting 
       if ( cabs2(z) > ER2 ){ return iColorOfExterior;}
	
      z = f(z);		//  iteration: z(n+1) = f(zn)
	

    }

  uUnknown += 1;
  return iColorOfUnknown;


}






// ********************************************************************************************************************
/* ---------------------Level Set Method = LSM -----------------------------------------------------------*/
// ********************************************************************************************************************


/*
2 basins
 exterior is basin of infinity
 interior is superattracting 
 julia set is connected 
*/
unsigned char ComputeColor_LSM_superattracting(complex double z)
{

  int i;			// number of iteration
  for (i = 0; i < IterMax_LSM; ++i)
    {


	double cabs2z = cabs2(z);	

      // infinity is superattracting here ,
      // if ( cabs2z > ER2 ) 	{ return (13*i) % 255;} // exterior
      // if ( cabs2z < AR2 ) 	{ return 255- ((7*i) % 255);} // interior	
       
        if ( cabs2z > ER2 || ( cabs2z < AR2 )) // extarior and interior
       		{ return (10*i) % 255;} // cabs2(zp-z) = cabs2(z) because zp = zcr = 0
	 
      z = f(z);	

    }

  return iColorOfUnknown;


}


/*
2 basins
 exterior is basin of infinity
 interior is attracting 
 julia set is connected 
*/

unsigned char ComputeColor_LSM_attracting(complex double z)
{

  int i;			// number of iteration
  for (i = 0; i < IterMax_LSM; ++i)
    {
      // infinity is superattracting here 
       if ( cabs2(z)      > ER2 ) 	{ return (10*i) % 255;} // exterior
       if ( cabs2(zp - z) < AR2 ) 	{ return 255- ((7*i) % 255);} // interior	
       
       
	 
      z = f(z);

    }

  return iColorOfUnknown;


}


/*
no interior = julia set is disconnected = only one basin here 
*/
unsigned char ComputeColor_LSM_repelling(complex double z)
{


  int i;			// number of iteration
  for (i = 0; i < IterMax_LSM; ++i)
    {

	
      // infinity is superattracting here , only one basin	
       if ( cabs2(z) > ER2 )
       		{ return (10*i) % 255;} // cabs2(zp-z) = cabs2(z) because zp = zcr = 0
       		
      z = f(z);	

    }

  return iColorOfUnknown;


}

/*
for child_period 1 and 2 it can work
for child period > 2 it not works, because target set is a triangle fragment of circle = not works good for parabolic case

z_n < center < z_p

here  AR = (z_p - z_n)/2

It is parabolic case: compute AR and change trap center in the local setup procedure


*/
unsigned char ComputeColor_LSM_parabolic(complex double z)
{

  int i;			// number of iteration
  for (i = 0; i < IterMax_LSM; ++i)
    {


	

      // infinity is superattracting here !!!!!	
       if ( cabs2(z) > ER2 )
       		{ return (10*i) % 255;} // 
       		
	//distance = cabs(z - zp);
	if (IsPointIn_Parabolic_Components_Triangle(z) ) // if z is inside target set ( orbit trap) = interior of cirlce with radius AR
		{ return (10*i) % 255;} //
     
	
      z = f(z);	

    }

  return iColorOfUnknown;


}







// ********************************************************************************************************************
/* ---------------------Binary Decomposition Method = BDM -----------------------------------------------------------*/
// ********************************************************************************************************************



unsigned char ComputeColor_BDM_superattracting (complex double z)
{



	
	
	double cabs2z;
	//double distance2;
	double turn;
	int i;			// number of iteration
	
	
	for (i = 0; i < IterMax_LSM; ++i)
	{
		cabs2z = cabs2(z); // numerical speed up : cabs2(zp-z) = cabs2(z) because zp = zcr = 0	

		//  if z is inside target set ( orbit trap) = exterior of circle with radius ER 
		if ( cabs2z > ER2  ) // exterior
       		{ 
       			if (cimag(z) > 0) // standard binary decomposition of target set
       				{  return 0;}
       				else {return 255; }
		}
		
		if ( cabs2z  < AR2 ) // if z is inside target set ( orbit trap) = interior of cirlce with radius AR
      		{
      			turn = c_turn(z);
      			if (turn < t0 || turn > t0+0.5) // modified binary decomposition of target set
      				{  return 0;}
       				else {return 255; }
     
      		
      		}
		
		
		
		
		
		z = f(z);	

	}

  return iColorOfUnknown;


}







unsigned char ComputeColor_BDM_attracting (complex double z)
{



	int i;			// number of iteration
	
	
	for (i = 0; i < IterMax_LSM; ++i)
	{
		//cabs2z = ; // numerical speed up : cabs2(zp-z) = cabs2(z) because zp = zcr = 0	

		//  if z is inside target set ( orbit trap) = exterior of circle with radius ER 
		if ( cabs2(z) > ER2  ) // exterior
       		{ 
       			if (cimag(z) > 0) // binary decomposition of target set
       				{  return 0;}
       				else {return 255; }
		}
		
		if ( cabs2(zp - z) < AR2  ) // interior
       		{ 
       			if (cimag(z) > 0) // binary decomposition of target set // ? modificate ?
       				{  return 0;}
       				else {return 255; }
		}     
	
		z = f(z);	

	}

  return iColorOfUnknown;


}






/*
no interior = julia set is disconnected = only one basin here 
*/
unsigned char ComputeColor_BDM_repelling (complex double z)
{

	int i;			// number of iteration
	
	
	for (i = 0; i < IterMax_LSM; ++i)
	{
		// numerical speed up : cabs2(zp-z) = cabs2(z) because zp = zcr = 0	
		//  if z is inside target set ( orbit trap) = exterior of circle with radius ER 
		if ( cabs2(z) > ER2  ) // exterior
       		{ 
       			if (cimag(z) > 0) // binary decomposition of target set
       				{  return 0;}
       				else {return 255; }
		}
		
		     
	
		z = f(z);	

	}

  return iColorOfUnknown;


}

unsigned char ComputeColor_BDM_parabolic (complex double z)
{

	int i;			// number of iteration
	
	
	for (i = 0; i < IterMax_LSM; ++i)
	{
		
		//  if z is inside target set ( orbit trap) = exterior of circle with radius ER 
		if ( cabs2(z) > ER2  ) // exterior
       		{ 
       			if (cimag(z) > 0) // BDM = binary decomposition (of target set) method
       				{  return 0;}
       				else {return 255; }
		}
		
		// interior: modified BDM
		if (IsPointIn_BDM_Triangle_B_mB(z)) { return 100; }//50 + (i % period_child);}
      		if (IsPointIn_BDM_Triangle_A_mA(z)) { return 200; }// - (i % period_child); }     		
      		
     
	
		z = f(z);	

	}

  return iColorOfUnknown;


}



// BDM for interior + LSM for exterior

unsigned char ComputeColor_BDM_LSM_parabolic(complex double z)
{

	int i;			// number of iteration
	
	
	for (i = 0; i < IterMax_LSM; ++i)
	{
		
		//  exterior : LSM
		if ( cabs2(z) > ER2 )
       			{ /*
       				int d = 8*i;
       				if (d>255){d = 255;}
       				return d;
       			*/
       			//return ((10*i) % 255);
       			return 250; // all pixels wth the same color	
       			} 
       			
		// interior: BDM 
      		if (IsPointIn_BDM_Triangle_B_mB(z)) { return 100; }//50 + (i % period_child);}
      		if (IsPointIn_BDM_Triangle_A_mA(z)) { return 200; }// - (i % period_child); }     		
      		
     
	
		z = f(z);	

	}

  return iColorOfUnknown;


}








/* 
==================================================================================================
============================= Draw functions ===============================================================
=====================================================================================================
	
	
	typedef enum  {FatouBasins = 0, FatouComponents = 2,  LSM = 3, LSM_m = 4, Unknown = 5 , BDM = 6, MBD = 7 , MBD2 = 8, SAC = 9, DLD = 10, ND = 11, NP= 12, POT = 13 , Blend = 14, DEM = 15} FunctionTypeT; 
	typedef enum  {superattracting = 100, attracting = 200, parabolic = 300,  repelling = 400} DynamicTypeT; 

	
*/ 
unsigned char Compute8BitColor(FunctionTypeT FunctionType, complex double z){

	unsigned char iColor;
	
	// not exery pair gives good result
	switch(DynamicType+FunctionType){
	
		// case FunctionType + DynamicType: {ComputeColor_FunctionType_DynamicType(z); break;}
  		case  FatouBasins + superattracting: 	{iColor = ComputeColor_FatouBasins_superattracting(z); 	break;}
		case  FatouBasins + attracting: 	{iColor = ComputeColor_FatouBasins_attracting(z); 	break;}
		case  FatouBasins + parabolic: 		{iColor = ComputeColor_FatouBasins_parabolic(z); 	break;}
		case  FatouBasins + repelling: 		{iColor = ComputeColor_FatouBasins_repelling(z); 	break;}
		
		case FatouComponents + superattracting: {iColor = ComputeColor_FatouComponents_superattracting(z); break;}
		case FatouComponents + attracting:      {iColor = ComputeColor_FatouComponents_attracting(z); break;}
		case FatouComponents + parabolic:       {iColor = ComputeColor_FatouComponents_parabolic(z); break;}
		case FatouComponents + repelling:       {iColor = ComputeColor_FatouComponents_repelling(z); break;}
  
		case LSM + superattracting: 	{iColor = ComputeColor_LSM_superattracting(z); 	break;}
		case LSM + attracting: 		{iColor = ComputeColor_LSM_attracting(z); 	break;}
  		case LSM + parabolic:  		{iColor = ComputeColor_LSM_parabolic(z); 	break;}
		case LSM + repelling:  		{iColor = ComputeColor_LSM_repelling(z);	break;}
  		
  		case BDM + superattracting: 	{iColor = ComputeColor_BDM_superattracting(z); 	break;}
		case BDM + attracting: 		{iColor = ComputeColor_BDM_attracting(z); 	break;}
  		case BDM + parabolic:  		{iColor = ComputeColor_BDM_parabolic(z); 	break;}
  		case BDM + repelling:  		{iColor = ComputeColor_BDM_repelling(z); 	break;}
  		
  		
  		case BDM_LSM + parabolic: 	{iColor = ComputeColor_BDM_LSM_parabolic(z); 	break;}
  
  	default: break;
	
	
  }
	
  return iColor;



}


// plots raster point (ix,iy) 
int DrawPoint ( unsigned char A[], FunctionTypeT FunctionType, int ix, int iy)
{
  int i;			/* index of 1D array */
  unsigned char iColor;
  complex double z;


  i = Give_i (ix, iy);		/* compute index of 1D array from indices of 2D array */
  if(i<0 && i> iMax)
    { return 1;}
  
  z = GiveZ(ix,iy);
  iColor = Compute8BitColor(FunctionType, z);
  A[i] = iColor ;		// 
  		
  	  
  return 0;
}




int DrawImage ( unsigned char A[], FunctionTypeT FunctionType)
{
  unsigned int ix, iy;		// pixel coordinate 

  fprintf (stderr, "compute image %d \t c = %.16f%+.16f*i \n",  FunctionType, creal(c), cimag(c));
  // for all pixels of image 
#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax, uUnknown, uInterior, uExterior)
  for (iy = iyMin; iy <= iyMax; ++iy)
    {
      fprintf (stderr, " %d from %d \r", iy, iyMax);	//info 
      for (ix = ixMin; ix <= ixMax; ++ix)
	DrawPoint(A, FunctionType, ix, iy);	//  
    }
  fprintf (stderr, "\n");	//info 
  return 0;
}







int PlotPoint(const complex double z, unsigned char A[]){

	
  unsigned int ix = (creal(z)-ZxMin)/PixelWidth;
  unsigned int iy = (ZyMax - cimag(z))/PixelHeight;
  unsigned int i = Give_i(ix,iy); /* index of _data array */
	
  	
	
  if(i>-1 && i< iMax)
    {A[i]= 0; // 255-A[i];
    }
	
	
  return 0;
	
}






int IsInsideCircle (int x, int y, int xcenter, int ycenter, int r){

	
  double dx = x- xcenter;
  double dy = y - ycenter;
  double d = sqrt(dx*dx+dy*dy);
  if (d<r) {    return 1;}
  return 0;
	  

} 

// Big point = disk 
int PlotBigPoint(const complex double z, double p_size, unsigned char A[]){

	
  unsigned int ix_seed = (creal(z)-ZxMin)/PixelWidth;
  unsigned int iy_seed = (ZyMax - cimag(z))/PixelHeight;
  unsigned int i;
  //unsigned char temp;
	
	
  if (  is_z_outside_viewport(z)) 
    {fprintf (stdout,"PlotBigPoint :  z= %.16f %+.16f*I is outside\n", creal(z), cimag(z)); return 1;} // do not plot	
	
  /* mark seed point by big pixel */
  int iSide =p_size*iWidth/2000.0 ; /* half of width or height of big pixel */
  int iY;
  int iX;
  for(iY=iy_seed-iSide;iY<=iy_seed+iSide;++iY){ 
    for(iX=ix_seed-iSide;iX<=ix_seed+iSide;++iX){ 
      if (IsInsideCircle(iX, iY, ix_seed, iy_seed, iSide)) {
	i= Give_i(iX,iY); /* index of _data array */
	//if(i>-1 && i< iMax)
	//temp = A[i];
	//if ( temp > 110 && temp < 160)
	//	{ A[i] = 0; }
	//	else {A[i]= 255 - temp  ; }
	A[i] = 0; 
      }
      // else {printf(" bad point \n");}
	
    }}
	
	
  return 0;
	
}






int DrawForwardOrbit(const complex double z0, const unsigned long long int i_Max, double p_size, unsigned char A[]){
 

  
  unsigned long long int i; /* nr of point of critical orbit */
  complex double z = z0;
  //complex double zt = f(f(z)); 
  fprintf(stdout, "draw forward orbit \n");
  fprintf(stderr, "draw forward orbit \n");
  PlotBigPoint(z, p_size, A);
  int iMax;
  
  switch ( DynamicType){
  	case superattracting: iMax = period_child; break;
  	default: iMax = i_Max*period_child ;
  	}
  
  
  
  
  
  /* forward orbit of critical point  */
  for (i=1;i< iMax; ++i)
    {
      z  = f(z);
      //if (cabs2(z - zp) < AR2) {break;} // 
      if (cabs2(z - zp) < PixelWidth2) 
      	{	fprintf (stdout,"last point of the orbit z= %.16f %+.16f*I \n After i = %llu iterations forward orbit of critical point reaches trap: circle with radius = PixelWidth around fixed point \n", creal(z), cimag(z), i );
      		break;} // 
	//if (cabs2(z-zt) > PixelWidth2)		
		PlotBigPoint(z, p_size/2 , A);
	//zt = z;
    }
    
    
  zcr_last = z;
  fprintf (stdout,"first point of the orbit z0= %.16f %+.16f*I \n", creal(z0), cimag(z0));
  // printf (stdout,"last point of the orbit z= %.16f %+.16f*I \n After i = %llu iterations forward orbit of critical point reaches trap: circle with radius = AR around fixed point \n", creal(z), cimag(z), i );
  fprintf (stdout,"last point of the orbit z= %.16f %+.16f*I after i = %llu iterations. \n", creal(z), cimag(z), i);
  fprintf (stdout,"distance between last point of the orbit and fixed point = %.16f  = %.16f * ImageWidth = %.1f * PixelWidth = \n ", cabs(z - zp), cabs(z-zp)/(ZxMax - ZxMin), cabs(z - zp)/PixelWidth); 
  fprintf (stdout,"\n \n "); 
  return 0;
 
}



// ***********************************************************************************************
// ********************** draw line segment ***************************************
// ***************************************************************************************************




// 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 */
	if (ix >=ixMin && ix<=ixMax && iy >=iyMin && iy<=iyMax )
		{A[Give_i(ix,iy)] = iColor;}
		else {fprintf (stderr,"iDrawPoint :   (%d; %d) is outside\n", ix,iy); }

  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 complex Z0, double complex Z1, int color, unsigned char *array) 
{

  double Zx0 = creal(Z0);
  double Zy0 = cimag(Z0);
  double Zx1 = creal(Z1);
  double Zy1 = cimag(Z1);
  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;
}







// ***********************************************************************************************
// ********************** mark attractors, traps ***************************************
// ***************************************************************************************************

int DrawAttractors(const complex double zpp[], int kMax, double p_size, unsigned char A[]){

	// join points by line to create closed curve
  	for (int p = 0; p < period_child-1; ++p){
  		dDrawLineSegment(zpp[p], zpp[p+1],0,A);
  		PlotBigPoint(zpp[p], p_size, A);}
  	
  	dDrawLineSegment(zpp[period_child-1], zpp[0],0,A);
  	PlotBigPoint(zpp[period_child-1], p_size, A);
  	return 0;

}




int Mark_Parabolic_Components_Triangle_Trap(unsigned char A[]){

  unsigned int ix, iy;		// pixel coordinate 
  unsigned int i;
  unsigned char temp;


  fprintf (stderr, "Mark_Parabolic_Components_Triangle_Trap\n");
  // for all pixels of image 
#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax)
  for (iy = iyMin; iy <= iyMax; ++iy)
    {
      fprintf (stderr, " %d from %d \r", iy, iyMax);	//info 
      for (ix = ixMin; ix <= ixMax; ++ix){
	if ( IsPointIn_Parabolic_Components_Triangle(GiveZ(ix,iy))){
	
	
	  	i= Give_i(ix,iy); /* index of _data array */
	  	//A[i]= 255-A[i]; // inverse color
	  	temp = A[i];
		if ( temp > 110 && temp < 160)
			{ A[i] = 0; }
			else {A[i]= 255 - temp  ; }
	  	
	  	
	}}}
  return 0;
}


int Mark_Parabolic_BDM_Traps(unsigned char A[]){

  unsigned int ix, iy;		// pixel coordinate 
  unsigned int i;
  //unsigned char temp;
complex double z;

  fprintf (stderr, "Mark_Parabolic_Components_Triangle_Trap\n");
  // for all pixels of image 
#pragma omp parallel for schedule(dynamic) private(ix,iy,i, z) shared(A, ixMax , iyMax, )
  for (iy = iyMin; iy <= iyMax; ++iy)
    {
      fprintf (stderr, " %d from %d \r", iy, iyMax);	//info 
      for (ix = ixMin; ix <= ixMax; ++ix){
      	i= Give_i(ix,iy); /* index of _data array */
      	z = GiveZ(ix,iy);
      	
	if ( IsPointIn_BDM_Triangle_B_mB(z)) { A[i]= 150;} // 
	if ( IsPointIn_BDM_Triangle_A_mA(z)) { A[i]= 250;} // 
		
	  	
	  	
	}}
  return 0;
}





int Mark_Basin_Trap(unsigned char A[]){

  unsigned int ix, iy;		// pixel coordinate 
  unsigned int i;
  unsigned char temp;


  fprintf (stderr, "Mark_Basin_Trap\n");
  // for all pixels of image 
#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax)
  for (iy = iyMin; iy <= iyMax; ++iy)
    {
      fprintf (stderr, " %d from %d \r", iy, iyMax);	//info 
      for (ix = ixMin; ix <= ixMax; ++ix){
	if ( IsInside_basin_trap(zpp[0])){
	
	
	  	i= Give_i(ix,iy); /* index of _data array */
	  	//A[i]= 255-A[i]; // inverse color
	  	temp = A[i];
		if ( temp > 110 && temp < 160)
			{ A[i] = 0; }
			else {A[i]= 255 - temp  ; }
	  	
	  	
	}}}
  return 0;
}




// ***********************************************************************************************
// ********************** mark immediate basin of attracting cycle***************************************
// ***************************************************************************************************


int FillContour(complex double seed,  unsigned char color, unsigned char _data[])
{ 
  /* 
     fills contour with black border ( color = iColorOfBoundary)  using seed point inside contour 
     and horizontal lines 
     it starts from seed point, saves max right( iXmaxLocal) and max left ( iXminLocal) interior points of horizontal line,
     in new line ( iY+1 or iY-1) it computes new interior point  : iXmidLocal=iXminLocal + (iXmaxLocal-iXminLocal)/2;
     result is stored in _data array : 1D array of 1-bit colors ( shades of gray)
     it does not check if index of _data array is good  so memory error is possible 
     
     it need array with components boundaries mrked by iColorOfBoundary
     
  */
	double dXseed = creal(seed);
	double dYseed = cimag(seed);
	// from 
  	int iXseed = (int)((dXseed - ZxMin)/PixelWidth);
  	int iYseed = (int)((ZyMax - dYseed )/PixelHeight); // reversed Y axis
  	
  	
  	
  	
  	int iX; /* seed integer coordinate */
    	int iY = iYseed;
    	/* most interior point of line iY */
    	int iXmidLocal=iXseed;
    	/* min and max of interior points of horizontal line iY */
    	int iXminLocal;
    	int iXmaxLocal; 
  	int i ; /* index of _data array */;


	//fprintf (stderr, "FillContour seed = %.16f %+.16f = %d %+d\n",creal(seed), cimag(seed), iXseed,iYseed);
  
  	/* ---------  move up --------------- */ 
  do{
    iX=iXmidLocal;
    i =Give_i(iX,iY); /* index of _data array */;
  
    /* move to right */
    while (_data[i] != iColorOfBoundary) 
      { _data[i]=color;
	iX+=1; 
	i=Give_i(iX,iY);  
      }
    iXmaxLocal=iX-1;

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

    iY+=1; /* move up */
    iXmidLocal=iXminLocal + (iXmaxLocal-iXminLocal)/2; /* new iX inside contour */
    i=Give_i(iXmidLocal,iY); /* index of _data array */;
    if ( _data[i] == iColorOfBoundary)  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 _data array */;
  
    /* move to right */
    while (_data[i] != iColorOfBoundary) /*  */
      { _data[i]=color;
	iX+=1;
	i=Give_i(iX,iY);  
      }
    iXmaxLocal=iX-1;

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

	//fprintf (stderr, "FillContour done \n");
  return 0;
}


// fill countours of componnets of immediate basin of attraction 
// with color 
// needs zpp and period global var
// it needs componnets boundaris in A array !!!!
// for attracting and superattracting basins

int  MarkImmediateBasin( unsigned char A[]){
	fprintf (stderr, "mark immediate basin of attracting cycle \n");
	
	//printf("  \n");
	unsigned char iColor = 100;
	for (int p=0; p<period_child ; ++p){
  		  FillContour(zpp[p], iColor , A);
  		 
  		}
 	return 0;
 	}






// ***********************************************************************************************
// ********************** edge detection usung Sobel filter ***************************************
// ***************************************************************************************************

// 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
  memset(D, iColorOfBasin1, iSize*sizeof(*D)); // for heap-allocated arrays, where N is the number of elements = FillArrayWithColor(D , iColorOfBasin1);
 
  // 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)
  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;
}







// FillAllArrayWithColor
//memset (data, 255, sizeof (unsigned char ) * iSize);








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

int SaveArray2PGMFile (unsigned char A[], complex double c, char * n, 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[100];		/* name of file */
  snprintf (name, sizeof name, "%s_%.16f%+.16f",n, creal(c),cimag(c));	/* radius and iHeght are global variables */
  char *filename = strcat (name, ".pgm");
  char long_comment[200]; // to long comment can cause: "*** stack smashing detected ***: terminated"
  sprintf (long_comment, "%s %s", f_description , comment); // f_description is global var





  // save image array 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", long_comment, iWidth, iHeight, MaxColorComponentValue);	// write header to the file
  size_t rSize = fwrite (A, sizeof(A[0]), iSize, fp);	// write whole array with image data bytes to the file in one step 
  fclose (fp);

  // info 
  if ( rSize == iSize) 
    {
      printf ("File %s saved ", filename);
      if (long_comment == NULL || strlen (long_comment) == 0)
	printf ("\n");
      else { printf (". Comment = %s \n", long_comment); }
    }
  else {printf("wrote %zu elements out of %llu requested\n", rSize,  iSize);}
  
  NumberOfImages +=1; // count images using global variable	

  return 0;
}



// *******************************************************************************************
// ********************************** info ****************************
// *********************************************************************************************





int PrintCInfo ()
{

  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 displayed in the console : export  OMP_DISPLAY_ENV="TRUE"

  printf ("__STDC__ = %d\n", __STDC__);
  printf ("__STDC_VERSION__ = %ld\n", __STDC_VERSION__);
  printf ("c dialect = ");
  switch (__STDC_VERSION__)
    {				// the format YYYYMM 
    case 199409L:
      printf ("C94\n");
      break;
    case 199901L:
      printf ("C99\n");
      break;
    case 201112L:
      printf ("C11\n");
      break;
    case 201710L:
      printf ("C18\n");
      break;
      //default : /* Optional */

    }

  return 0;
}


int
PrintProgramInfo ()
{


 	fprintf(stdout,"Program info\n");
  	fprintf (stdout, "%s  \n", f_description );
  	fprintf (stdout, "c =  %.16f %+.16f*i  \n", creal (c), cimag (c));
  	
  	double ImageWidth = (ZxMax - ZxMin);
  	
  	fprintf (stdout, "Dynamic type of the Julia set interior. Julia set exterior is allways superattracting for polynomials. \n")
  	fprintf (stdout, "DynamicType value is setup manually; One can do it also numerically ( from multiplier of fixed point alfa or from some other properities)\n");
  		switch ( DynamicType){
		case repelling: 
				fprintf (stdout, "\tThere is only one Fatou basin: superattracting basin of infinity \n");
				fprintf (stdout, "\tthere is no interior = Julia set is disconnected \n");
				fprintf (stdout, "\tcritical point z=0 is repelling = attracted to infinity \n");
				break;
		case attracting: 
	  			fprintf (stdout, "\tinterior basin type is attracting = critical point is attracted by periodic point\n");
	  			fprintf (stdout, "\tzcr_last =  %.16f \talfa fixed point zp = %.16f\n", creal (zcr_last), creal(zp));// 
	  			fprintf (stdout, "\tdelta =  %.16f is the distance between fixed points\n", delta);// 
	  			fprintf (stdout, "\tAtracting Radius AR is set manually  = %.16f = %f * PixelWidth = %f * ImageWidth \n", AR, AR / PixelWidth, AR /ImageWidth );
	  			break;
  			
		case superattracting: 
				fprintf (stdout, "\tinterior basin type is superattracting  =  critical point is als superattracting periodic point\n");
	  			fprintf (stdout, "\tzcr =  %.16f  = zp = %.16f\n", creal (zcr), creal(zp));// 
	  			fprintf (stdout, "\tAtracting Radius AR is set manually  = %.16f = %f *PixelWidth = %f *ImageWidth \n", AR, AR / PixelWidth, AR /ImageWidth);
	  			break;
		
		case parabolic:
				fprintf (stdout, "\tinterior basin type  is parabolic : all  periodic points are equal to fixed point\n"); 
				fprintf (stdout, "\t internal angle = %d/%d\n", internal_angle_numerator, internal_angle_denominator);
				fprintf (stdout, "\tc is a root point between period %d and period %d componnets\n", internal_angle_numerator, internal_angle_denominator);//
				
				
				fprintf (stdout, "\t trap for componnets is a triangle: ( zp,zprecr,-zprecr) =  curvilinear sector ( triangular part) of the circle around parabolic fixed point\n");
				fprintf (stdout, "\t zprecr is computed in other program \n");
				fprintf (stdout, "\t-zprecrz =  %.16f %+.16f*I\n", -creal (zprecr), -cimag(zprecr));// 
				fprintf (stdout, "\tzp is a parabolic fixed point zp = %.16f %+.16f*I\n", creal (zp), cimag(zp));// 
				fprintf (stdout, "\tzprecr =  %.16f %+.16f*I\n", creal (zprecr), cimag(zprecr));// 

				fprintf (stdout, "\t traps for BDM are the triangles: ( zp, zprecr, zcr), ( zp, zcr, - zprecr) =  curvilinear sectors ( triangular parts) of the circle around parabolic fixed point\n");
				fprintf (stdout, "\t zprecr is computed in other program \n");
				fprintf (stdout, "\t-zprecrz =  %.16f %+.16f*I\n", -creal (zprecr), -cimag(zprecr));// 
				fprintf (stdout, "\tzp is a parabolic fixed point zp = %.16f %+.16f*I\n", creal (zp), cimag(zp));// 
				fprintf (stdout, "\tzprecr =  %.16f %+.16f*I\n", creal (zprecr), cimag(zprecr));// 
				
				fprintf (stdout," ASCI graphic: target sets for BDM ( 4 triangles, )\n");
				fprintf (stdout,"\t\t -zprecr\n");
				fprintf (stdout,"zp\t\t		zcr\t\t 		-zp\n");
				fprintf (stdout,"\t\t zprecr\n");
				
				fprintf (stdout," ASCI graphic: target set for componnets ( triangle)\n");
				fprintf (stdout,"\t\t -zprecr\n");
				fprintf (stdout,"zp\n");
				fprintf (stdout,"\t\t zprecr\n");
				
				

				break;
		default: 
	
	
	
	
	}		
			
			
			
			
  	fprintf (stdout, "Image Width = %f in world coordinate\n", ImageWidth);
  	fprintf (stdout, "PixelWidth = %.16f \n", PixelWidth);
  	


  	fprintf (stdout, "plane description \n");
  	fprintf (stdout, "\tcenter z =  %.16f %+.16f*i  and radius = %.16f \n", creal (center), cimag (center), radius);
  	// center and radius
  	// center and zoom
  	// GradientRepetition
  	fprintf (stdout, "Maximal number of iterations = iterMax = %d \n", IterMax);
  	fprintf (stdout, "Maximal number of iterations = iterMax_LSM = %d \n", IterMax_LSM);
  	fprintf (stdout, "ratio of image  = %f ; it should be 1.000 ...\n", ratio);
  
   	fprintf (stdout, "\tAttracting Radius = AR = %.16f = %f *PixelWidth = %f * ImageWidth \n", AR, AR / PixelWidth, AR /ImageWidth);
	fprintf (stdout, "\tEscaping Radius = ER = %.16f = %f *PixelWidth = %f * ImageWidth \n", ER, ER / PixelWidth, ER /ImageWidth);
   	fprintf(stdout, " periodic point ");
   	//for (int i=0;i<period ; ++i){
  		fprintf(stdout, "z =  %.16f %+.16f*i  \n", creal (zp), cimag (zp));
  		//}
  	fprintf (stdout, "Unknown pixels = %llu = %.16f * iSize \n", uUnknown, ((double) uUnknown)/iSize  );
	fprintf (stdout, "Exterior pixels = %llu = %.16f * iSize \n", uExterior, ((double) uExterior)/iSize  );
	fprintf (stdout, "Interior pixels = %llu = %.16f * iSize \n", uInterior, ((double) uInterior)/iSize  );
	//printf("Number of images = %d \n", NumberOfImages);	


  return 0;
}



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

  ZxMin = creal(center) - radius*a_ratio;	
  ZxMax = creal(center) + radius*a_ratio;	//0.75;
  ZyMin = cimag(center) - radius;	// inv
  ZyMax = cimag(center) + radius;	//0.7;
  return 0;

}



// Check Orientation of z-plane image : mark first quadrant of complex plane 
// it should be in the upper right position
// uses global var :  ...
int CheckZPlaneOrientation(unsigned char A[] )
{
 
  double Zx, Zy; //  Z= Zx+ZY*i;
  unsigned i; /* index of 1D array */
  unsigned int ix, iy;		// pixel coordinate 
	
  fprintf(stderr, "compute image CheckOrientation\n");
  // for all pixels of image 
#pragma omp parallel for schedule(dynamic) private(ix,iy, i, Zx, Zy) shared(A, ixMax , iyMax) 
  for (iy = iyMin; iy <= iyMax; ++iy){
    //fprintf (stderr, " %d from %d \r", iy, iyMax);	//info 
    for (ix = ixMin; ix <= ixMax; ++ix){
      // from screen to world coordinate 
      Zy = GiveZy(iy);
      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;
}







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

void usage(){

	fprintf (stdout, "Program is ready to use. To get different setting change importatnt input values, computed in other program\n" ); // 
	fprintf (stdout, "\tDynamicType\n" );
	fprintf (stdout, "\tDynamicType value is setup manually; One can do it also numerically ( from multiplier of fixed point alfa or from some other properities)\n");
	fprintf (stdout, "\tc : parameter of function ( map) f \n" );
	fprintf (stdout, "\tzp : fixed point of function ( map) f , z : zp = f(zp)  \n" );

	switch ( DynamicType){
		case repelling: 
				
			break;
		case attracting: 
			
	  		break;
  			
		case superattracting:  
			
	  		break;
		
		case parabolic:
			fprintf (stdout, "z points describing triangle traps for components, BDM: precritical point and prefixed point\n" );
			fprintf (stdout, "period_parent, period_child, internal_angle_numerator, internal_angle_denominator\n");

				break;
		default: 
	}


}



int local_setup(){
	
	
	complex double z;
	
	switch ( DynamicType){
		case repelling: // no  interior = no attracting fixed point = only escaping points
				
			break;
		case attracting: 
			delta = csqrt(1.0 - 4.0*c);  // delta is a distance between alfa and beta fixed points
			AR =  (delta /3.5) * iWidth / 2000;
			fprintf(stderr," error : compute periodic point: array zpp !!!\n");
			return 1;
	  		break;
  			
		case superattracting: // cabs(zp - zcr_last ) < PixelWidth 
			AR = 30.0* PixelWidth * iWidth / 2000 ; //
			// compute periodic point: array zpp
			z = zcr;
			for (int p = 0; p < period_child; ++p){
				zpp[p] = z;
				z = f(z);} 
			
			 
	  		break;
		
		case parabolic:
				/* Target set for the components for child period:  
					* 1 and 2: target set can be circle with parabolic fixed point on it's boundary ( target set is inside component )
					* > 2: target set is a triangle fragment of circle with center at parabolic fixed point

				
				 parabolic trap is a triangle: ( z1,z2,z3) = part of the circle around parabolic fixed point 
				 z2 and z3 are computed in other program
				 */
				Describe_Triangle(zp, zprecr, -zprecr); // test if triangle is oriented counterclockwise
				// all periodic point ( for period parent and child) are the same : zp = zpp[0] = zpp[1] ...
				AR = cabs(zp-zcr)/10.0;
				break;
		default: 
	}
	
	
	
	
	
	AR2 = AR*AR;
	
	
	return 0;
}




int general_setup()
{

  fprintf (stderr, "setup start\n");


	zp = GiveFixed(c);
	center = 0.0; //zp; 
	radius = 1.5; //0.24;

	z13 = (z1+z3)/2.0;
	 
  /* 2D array ranges */

  iWidth = iHeight* DisplayAspectRatio ;
  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].

  
  SetPlane( center, radius,  DisplayAspectRatio );	
  /* Pixel sizes */
  PixelWidth = (ZxMax - ZxMin) / ixMax;	//  ixMax = (iWidth-1)  step between pixels in world coordinate 
  PixelHeight = (ZyMax - ZyMin) / iyMax;
  ratio = ((ZxMax - ZxMin) / (ZyMax - ZyMin)) / ((double) iWidth / (double) iHeight);	// it should be 1.000 ...
  PixelWidth2 = PixelWidth*PixelWidth;
    
  
  // LSM  
  // escape radius ( of circle around infinity 
  ER = 200.0; // it can be 2.0: but then there is no level curves near Julia set;  but the small detailes will be visible 
  ER2 = ER*ER;
  
  
  
  //AR = 30* PixelWidth ; //GiveTunedAR(IterMax- 300); 
  //AR2 = AR * AR;
  
  // for MBD
  internal_angle_denominator = period_child;
  t0 = ((double) internal_angle_numerator)/ internal_angle_denominator; // Is it iternal angle from internal adress  ???
  
  
  // DEM
 // BoundaryWidth = 0.5*iWidth/2000.0  ; //  measured in pixels ( when iWidth = 2000) 
  //distanceMax = BoundaryWidth*PixelWidth;



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




 


  fprintf (stderr, " end of setup \n");

  return 0;

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




int end(void)
{


	fprintf(stderr, " allways free memory (deallocate )  to avoid memory leaks \n");	// https://en.wikipedia.org/wiki/C_dynamic_memory_allocation
	fprintf(stderr, " warning: too long comment in SaveArray2PGMFile can cause: *** stack smashing detected ***: terminated\n");
	free (data);
	free(edge);
	free(edge2);
	
	PrintCInfo ();
	return 0;

}



int MakeImages( ){

// warning: to long comment in SaveArray2PGMFile can cause: "*** stack smashing detected ***: terminated"

	
		

	/*
	DrawImage (data, FatouBasins);	 
	SaveArray2PGMFile (data, c, "FatouBasins" , "FatouBasins ");
	
	ComputeBoundaries(data,edge);
	SaveArray2PGMFile (edge, c, "FatouBasins_LCM" , "FatouBasins_LCM ");

	CopyBoundaries(edge, data);
	SaveArray2PGMFile (data, c,  "FatouBasins_LSCM" , "FatouBasins_LSCM");
	
	//DrawForwardOrbit(zcr, IterMax_LSM, 20, data);
	//DrawAttractors(zpp, period_child, 20, data);
	//SaveArray2PGMFile (data, c, "FatouBasins_attr" , "boundaries of Level set method ( LSM) =  Level Curve Method (LCM) and attractors ");
	
	//Mark_Parabolic_Components_Triangle_Trap(data);
	//SaveArray2PGMFile (data,  c, "FatouBasins_trap" , "FatouBasins  trap");

*/	
	

	


	DrawImage (data, FatouComponents);	 
	SaveArray2PGMFile (data, c, "FatouComponents" , "FatouComponents ");
    
    
	ComputeBoundaries(data,edge);
	SaveArray2PGMFile (edge, c, "FatouComponents_LCM" , "FatouComponents_LCM ");
    
	CopyBoundaries(edge, data);
	SaveArray2PGMFile (data, c,  "FatouComponents_LSCM" , "FatouComponents_LSCM");
	
	// MarkImmediateBasin(edge); // for attracting and superattracting, 
	// SaveArray2PGMFile (edge, c,  "FatouComponents_immmediate" , "FatouComponents_LCM + immediate basin");
	Mark_Parabolic_Components_Triangle_Trap(data);
	SaveArray2PGMFile (data,  c, "LSM_trap" , "LSM  trap");
	
	DrawForwardOrbit(zcr, 1000, 15, data);
	//DrawAttractors(zpp, period_child, 20, edge); // for attracting and superattracting 
	SaveArray2PGMFile (data, c, "FatouComponents_LCM_cr" , "FatouComponents_LSCM, critical orbit");
	
	
	


	
	//Mark_Parabolic_Components_Triangle_Trap(data);
	//SaveArray2PGMFile (data, c, "FatouComponents_LSCM_CRO_trap" , "FatouComponents_LSCM, trap and critical orbit ");
	
    
   /* 
  Mark_Parabolic_Components_Triangle_Traps(data);
  SaveArray2PGMFile (data,  "FatouComponents_LSCM_trap" , "FatouComponents_LSCM_trap");
    
    
 
 // DrawAttractors(zp3, period, 10, data);
  //SaveArray2PGMFile (data,  "FatouComponents_LSCM_zp" , "FatouComponents_LSCM_zp");_
     
    DrawImage (data, FatouBasins);	 
    ComputeBoundaries(data,edge);
   MarkImmediateBasin(edge);
    SaveArray2PGMFile (edge,  "FatouBasins_LCM_immediate" , "FatouBasins_LCM_immediate");
    
    //DrawAttractors(zp3, period, 10, edge);
  //SaveArray2PGMFile (edge,  "FatouBasins_LCM_immediate_zp" , "FatouBasins_LCM_immediate_zp");
   
    */


/*
	DrawImage (data, LSM);	
	SaveArray2PGMFile (data, c, "LSM" , " Level sets of integer escape time of ET  ");

	ComputeBoundaries(data,edge2);
	SaveArray2PGMFile (edge2, c, "LSM_LC" , "boundaries of integer Escape Time =  Level Curve ( LC) of integer Escape Time");
	
	CopyBoundaries(edge2,  data);
	SaveArray2PGMFile (data, c, "LSM_LCM" , " Level sets and it's baundaries ( Level Curves = LC) of integer escape time of ET  ");
	
*/	
	
	
	
	
		
	
	DrawImage (data, BDM);	
	SaveArray2PGMFile (data, c,  "BDM" , "BDM = Binary Decomposition Method for both exterior and interior = Level Sets of BDM");
	
	
	
	ComputeBoundaries(data,edge);
	SaveArray2PGMFile (edge, c, "BDM_LC", "boundaries of Binary Decomposition Method (LC of  BD) = LCBD");
	
	CopyBoundaries( edge, data);
	SaveArray2PGMFile (data, c,  "BDM_LCM" , " Binary Decomposition Method for both exterior and interior ");
	
	
	
	
	//Mark_Parabolic_Components_Triangle_Trap(data);
	//SaveArray2PGMFile (data, c,  "BDM_LCM_cro_tr" , " Binary Decomposition Method for both exterior and interior and critical orbit, com trap");
	
	Mark_Parabolic_BDM_Traps(data);
	SaveArray2PGMFile (data, c,  "BDM_LCM_traps" , "BDM and traps" );
	
	DrawForwardOrbit(zcr, 1000, 10, data);
	SaveArray2PGMFile (data, c,  "BDM_LCM_cro_traps" , " Binary Decomposition Method for both exterior and interior and critical orbit");
	
	//CopyBoundaries(edge2,  edge);
	//SaveArray2PGMFile (edge, c, "BDM_LC2_only ", "boundaries of Binary Decomposition Method and Level sets of Integer Escape Time");

	//CopyBoundaries( edge, data);
	//SaveArray2PGMFile (data, c,  "BDM_LC2" , " Binary Decomposition Method for both exterior and interior ");
	
	

	

    
  	//CopyBoundaries(edge2, data);
 	//SaveArray2PGMFile (data,  "LSCM" , "LSCM");
    
  
    	DrawImage (data, BDM_LSM);	
	SaveArray2PGMFile (data, c,  "BDM_LSM" , "BDM for  exterior and solid color for interior");
	
	
	
	ComputeBoundaries(data,edge);
	SaveArray2PGMFile (edge, c, "BDM_LSM_LC", "boundaries of BDM for  exterior and solid color for interior");
	
	CopyBoundaries( edge, data);
	SaveArray2PGMFile (data, c,  "BDM_LSM_LCM" , " boundaries and LSM of BDM for  exterior and solid color for interior ");
	
	Mark_Parabolic_BDM_Traps(data);
	SaveArray2PGMFile (data, c,  "BDM_LSM_LCM_traps" , "BDM (LS and LC) and traps for interior, solid for exterior" );
	
	DrawForwardOrbit(zcr, 1000, 10, data);
	SaveArray2PGMFile (data, c,  "BDM_LSM_LCM_cro" , " Binary Decomposition Method for exterior and solid color for interior and critical orbit");
  
  
  /*
    
  CopyBoundaries(edge, data);
  SaveArray2PGMFile (data,  "MBD_LSCM" , "MBD_LSCM");
  
  CopyBoundaries(edge2, edge);
  SaveArray2PGMFile (edge,  "MBD_LSM_LCM" , "MBD_LSM_LCM ");
    
    CopyBoundaries(edge, data);
  SaveArray2PGMFile (data,  "MBD_LSM_LSCM" , "MBD_LSM_LSCM");
  
  
	DrawImage (data, FL);	
	SaveArray2PGMFile (data,  "FL" , "FL");
*/


	return 0;



}


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

int main(void)
{

	//int i;
	//int i_max = 10;
	
	//double cx_min = 0.0;
	//double cx_max = 0.250001;
	//double dcx = (cx_max - cx_min)/ i_max;
	//DynamicType = superattracting;  // setup DynamicType value manually  ; Once can do it also numerically ( from multiplier or from some properities)
	general_setup();

	
	//for (i = 0; i < i_max; ++i){
		local_setup();
		
		MakeImages();
		PrintProgramInfo();
	//	cx -= dcx;
		
	//	}
	end();
	return 0;
}