File:Julia set DEM c = -0.7432918908524300+0.1312405523087980*I.png

/*


  

  Adam Majewski
  adammaj1 aaattt o2 dot pl  // o like oxygen not 0 like zero 
  
  --------------------------------------------------------------------
 
  Structure of a program or how to analyze the program 
  
  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 n.c -lm -Wall -march=native -fopenmp
  time ./a.out > j.txt


  gcc n.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 PROGRAM_VERSION 20230409

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


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


DynamicTypeT DynamicType; // it is set manually 


// see ComputeColor_FunctionType_DynamicType belwo



// see SetPlane

double radius = 1.2; 
complex double center = 0.0 ;
double  DisplayAspectRatio  = 1.25; // 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;

double R_BDM;
double R2_BDM;


complex double parabolic_trap_center;
double PixelWidth2;


int IterMax = pow(10,7);
int IterMax_LSM = pow(10,7);
int IterMax_BDM = pow(10,9);
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;
const int period_child = 11;
const int t_numerator = 5;
int t_denominator = period_child;
double t0; // for MBD https://www.youtube.com/watch?v=JttLtB0Gkdk&t=894s
double t2 = 0.435;

/* critical point */
complex double zcr = 0.0; //
complex double c = -0.743291890852430 +0.131240552308798*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.0;  // one of periodic pointes , alfa
complex double zcr_last = -0.3230313566972330 +0.1160747675433939*I ; //


/*
  for parabolic triangle trap ( computed in other program: last point of first and last ray ( ordered by t )
 
  Maximal number of iterations = iterMax_ray = 2000000000000 
  p = 0 	 t = 0.1665852467024914 	z = -0.4170188752901364 +0.1506018053351233*I 	 distance to landing point = 0.0634786084515321 = 70.4965212747848113* PixelWidth = 0.0352658935841845*ImageWidth
  p = 1 	 t = 0.3331704934049829 	z = -0.5388360314336247 +0.1504188918297714*I 	 distance to landing point = 0.0598567181450451 = 66.4742108733029085* PixelWidth = 0.0332537323028029*ImageWidth
  p = 2 	 t = 0.6663409868099658 	z = -0.4223414442666910 +0.1139242453319452*I 	 distance to landing point = 0.0634130274946685 = 70.4236899788012778* PixelWidth = 0.0352294597192603*ImageWidth
 
  Maximal number of iterations = iterMax_ray = 20000 
  p = 0 	 t = 0.1665852467024914 	z = -0.3350594976618946 +0.1652825185133397*I 	 distance to landing point = 0.1467326739625789 = 162.9547862506640570* PixelWidth = 0.0815181522014327*ImageWidth
  p = 1 	 t = 0.3331704934049829 	z = -0.6051132306923424 +0.1652675081250403*I 	 distance to landing point = 0.1277193817845747 = 141.8394689929804144* PixelWidth = 0.0709552121025415*ImageWidth
  p = 2 	 t = 0.6663409868099658 	z = -0.3512112430996852 +0.0760154138719138*I 	 distance to landing point = 0.1439685504118509 = 159.8850734851610866* PixelWidth = 0.0799825280065838*ImageWidth
  
  z2 = parabolic fixed point
  
  triangle trap for components are different ( bigger) then triangle trap for BDM !!!!
  triangle (z1,zp,z3) 
  z1
  zp		zcr_last
  z3
  
*/
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 ; //  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 -------------------------------------------------------------*/

// complex function
// 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(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;

}


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

/*
  z1
  zp		zcr_last
  z3
*/

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

  return PointInTriangle (creal(z), cimag(z), creal(z1), cimag(z1), creal(zp), cimag(zp), creal(z3), cimag(z3));

}

/*  divide circle around critical point into 4 equal segments , manual tuning 
    t0 = 5/11
    tn = n/4

    t0+1/4					0.18				0.175					0.1786329062544335
																		
    0.0					0.0				0.0							
																	
    1.0*2/4		1*3/4			0.43		0.93		0.425		0.925		0.4286329062544335	0.9286329062544335
												
    1*					0.7				0.675				0.6786329062544335
	
    z = -0.0559603266109820 -0.1163177312664140
    t1 = 0.1786329062544335 = t2 - 0.25
    t2 = 0.4286329062544335 = t3 - 0.25  ** 0.429 0.44 ( tu much )  0.435
    t3 = 0.6786329062544335 = c_turn(z)
    t4 = 0.9286329062544335 = t3 + 0.25   **0.929


    t2 - t1 = 0.2500000000000000
    t3 - t2 = 0.2500000000000000
    t4 - t3 = 0.2500000000000000
    t1+ (1.0 - t4) = 0.2500000000000001

*/
int IsPointIn_BDM_TriangleB(const double complex z){
	
  int b  = PointInTriangle (creal(z), cimag(z), creal(z1), cimag(z1), creal(zp), cimag(zp), creal(zcr_last), cimag(zcr_last));
  int bm = PointInTriangle (creal(z), cimag(z), creal(-z1), cimag(-z1), creal(-zp), cimag(-zp), creal(-zcr_last), cimag(-zcr_last));
	
  int c = 0;
  if (cabs2(z) < R2_BDM ){
		
    double turn = c_turn(z);
    //c =  ((turn > 0.75 + t0 && turn < 0 + t0) || (turn > 0.25 + t0  && turn <  0.5 + t0)); // divide circle into 4 equal segments rotated by t0
    c =  ((turn > 0.1786 && turn < t2) || (turn > 0.6786 && turn <  0.935 )); // divide circle into 4 equal segments rotated by t0      		}
  }		
	
	
  return (b || bm || c);

}



/*
  z1
  zp		zcr_last
  z3
*/

int IsPointIn_BDM_TriangleA(const double complex z){

  int a  = PointInTriangle (creal(z), cimag(z), creal(zcr_last), cimag(zcr_last), creal(zp), cimag(zp), creal(z3), cimag(z3));
  int am = PointInTriangle (creal(z), cimag(z), creal(-zcr_last), cimag(-zcr_last), creal(-zp), cimag(-zp), creal(-z3), cimag(-z3));
	
  int c = 0;
  if (cabs2(z) < R2_BDM ){
		
    double turn = c_turn(z);
    //c =  ((turn > 0.5 + t0 && turn < 0.75 + t0) || (turn > 0 + t0  && turn <  0.25 + t0)); // divide circle into 4 equal segments rotated by t0
    c =  ((turn > 0.935  || turn < 0.1786) || (turn > t2  && turn <  0.6786)); // divide circle into 4 equal segments rotated by t0
  }
      				
	
	
  return (a || am || c);


}


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




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


// circular trap with zp as a center

int IsInsideTrap(int ix, int iy){


  complex double z = GiveZ(ix, iy);
  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 -----------------------------------------------------------*/
// ********************************************************************************************************************



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

  pixel counters work not good with OpenMP !!!

*/


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 zp= 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: attracting point zp is not equall to zero 
      if ( cabs2(zp-z) < 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;}
	
      // no Attraction basins 
      
		 
      			
     
      z = f(z);		//  iteration: z(n+1) = f(zn)
	

    }

  uUnknown += 1;
  return iColorOfUnknown;


}


unsigned char ComputeColor_FatouBasins_dendrite (complex double z)
{



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


		
      // infinity is superattracting here !!!!!	
      if ( cabs2(z) > ER2 ){ return iColorOfExterior;}
	
      // no Attraction basins 
     
		 
      			
     
      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)
{



  //double cabsz2;
  //double distance;


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


      //double cabs2z = cabs2(z);	

      // infinity is superattracting here , only one basin	
      if ( cabs2(z)      > ER2 ) 	{ return (13*i) % 255;} // 
      if ( cabs2(zp - z) < AR2 ) 	{ return 255- ((7*i) % 255);} //	
       
       
	 
      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)
    {


      //cabsz = cabs(z);	

      // infinity is superattracting here !!!!!	
      if ( cabs2(z) > ER2 )
	{ return (10*i) % 255;} // cabs2(zp-z) = cabs2(z) because zp = zcr = 0
       		
      //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;


}



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



  //double cabsz2;
  //double distance;


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


      //cabsz = cabs(z);	

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


}



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



  //double cabsz2;
  //double distance;


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


      //cabsz = cabs(z);	

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


}



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



unsigned char ComputeColor_BDM_superattracting (complex double z)
{


  double cabs2z;
  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  || cabs2z < AR2) // exterior
	{ 
	  if (cimag(z) > 0) // 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  ) // 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) // binary decomposition of target set
	    {  return 0;}
	  else {return 255; }
	}
		
		
      if (IsPointIn_BDM_TriangleB(z))
	{  return 100 ; } //50 + (i % period_child);}
      if (IsPointIn_BDM_TriangleA(z))
	{return 200; } // 255 - (i % period_child); }     		

	
      z = f(z);	

    }

  return iColorOfUnknown;


}


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

  double cabs2z;
  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) // binary decomposition of target set
	    {  return 0;}
	  else {return 255; }
	}
		
		     
	
      z = f(z);	

    }

  return iColorOfUnknown;


}




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

  double cabs2z;
  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) // binary decomposition of target set
	    {  return 0;}
	  else {return 255; }
	}
		
		     
	
      z = f(z);	

    }

  return iColorOfUnknown;


}



// ************************** DEM/J*****************************************
/* 
it can be used for 
* whole image thru Compute8BitColor function
* only drawing boundary thru 
*/
unsigned char ComputeColorOfDEMJ(complex double z){
// https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Julia_set#DEM.2FJ


  int nMax = IterMax_DEM;
  complex double dz = 1.0; //  is first derivative with respect to z.
  double distance;
  double cabsz;
	
  int n;

  for (n=0; n < nMax; n++){ //forward iteration

    if (cabs2(z)> ER2 || cabs(dz)> 1e60) break; // big values
    
  			
    dz = 2.0*z * dz; 
    z = z*z +c ; /* forward iteration : complex quadratic polynomial */ 
  }
  
  if (n==nMax) return iColorOfInterior; // not escaping
  // escaping and boundary
  cabsz = cabs(z);
  distance = 2.0 * cabsz* log(cabsz)/ cabs(dz);
  double g = tanh(distance /PixelWidth);
  
  
  //if (distance <distanceMax) {//printf(" distance = %f \n", distance); 
  //		return iColorOfBoundary;}
  
  
  return 255*g; // iColorOfExterior;
}




/* 
   ==================================================================================================
   ============================= 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;
	
	
  switch(DynamicType+FunctionType){
	
    // case FunctionType + DynamicType: {ComputeColor_FunctionType_DynamicType(z); break;}
    //    FatouBasins
  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  FatouBasins + dendrite: 		{iColor = ComputeColor_FatouBasins_dendrite(z); 	break;}
		
    //   FatouComponents
  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;}
		
    //   LSM	
  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 LSM + dendrite:  		{iColor = ComputeColor_LSM_dendrite(z);		break;}
		
    //   BDM
  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 + dendrite:  		{iColor = ComputeColor_BDM_dendrite(z);	break;}
  
	
	case DEM+ dendrite:		{iColor = ComputeColorOfDEMJ(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;
}



// this function changes color of all image pixels 
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(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);
  
  
  /* forward orbit of critical point  */
  for (i=1;i<i_Max ; ++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 traps ***************************************
// ***************************************************************************************************


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_Triangle_Traps(unsigned char A[]){

  unsigned int ix, iy;		// pixel coordinate 
  unsigned int i;
  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){
	z = GiveZ(ix,iy);
	i= Give_i(ix,iy); /* index of _data array */
	// mark traps
	if (IsPointIn_BDM_TriangleB(z)){ A[i] = 200;}
	if (IsPointIn_BDM_TriangleA(z)){ A[i] = 100;}
	  	
      }}
  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 !!!!
int  MarkImmediateBasin( unsigned char A[]){
  fprintf (stderr, "mark immediate basin of attracting cycle \n");
	
  //printf("  \n");
  unsigned char iColor = 100;
  //for (int i=0;i<period ; ++i){
  FillContour(zp, 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;
}


// *******************************************************************************************
// ********************************** Boundary = DEM/J         ****************************
// *********************************************************************************************


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

  fprintf (stderr, "compute image DEM/J c = %.16f%+.16f*i \n",   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);	//  
	{
		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 = ComputeColorOfDEMJ(z);
  		if (iColor < iColorOfInterior)
  			{ A[i] = iColor ; }	
	
	
	
	}
	
    }
  fprintf (stderr, "\n");	//info 
  return 0;
}








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

int SaveArray2PGMFile (unsigned char A[], complex double c, char * n, unsigned int iHeight, 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*I_%d_%.3f",n, creal(c),cimag(c), iHeight, t2);	/* 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;
}








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


  // display info messages
  fprintf (stdout, "Program version = %d  \n", PROGRAM_VERSION ); // https://stackoverflow.com/questions/12509038/can-i-print-defines-given-their-values-in-c
  fprintf (stdout, "%s  \n", f_description );
  fprintf (stdout, "c =  %.16f %+.16f*i  \n", creal (c), cimag (c));
  	
  double ImageWidth = (ZxMax - ZxMin);
  	
  	
  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: 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, "\tbasin type is attracting \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, "\tbasin type is superattracting \n");
    fprintf (stdout, "\there superattracting periodic point ie equall to critical point : zcr =  %.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, "\tbasin type  is parabolic \n");
    fprintf (stdout, "\t parabolic trap for Julia set components is a triangle: ( z1,z2,z3) = part of the circle around parabolic fixed point\n");
    fprintf (stdout, "\t z2 and z3 are computed in other program \n");
    fprintf (stdout, "\tz1 =  %.16f %+.16f*I\n", creal (z1), cimag(z1));// 
    fprintf (stdout, "\tz2 is a parabolic fixed point zp = %.16f %+.16f*I\n", creal (zp), cimag(zp));// 
    fprintf (stdout, "\tz3 =  %.16f %+.16f*I\n", creal (z3), cimag(z3));// 
    fprintf (stdout, "\t for parabolic trap for Julia set components see procedure: IsPointIn_Parabolic_Components_Triangle\n");
    break;
	case dendrite:
	case siegel:
	default: break;
  }		
			
			
  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, "t0 = %.16f = %d / %d\n", t0, t_numerator, t_denominator);
  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 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
// **************************************************************************************

int local_setup(){
	
	
	
	
  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;
			
    break;
  			
  case superattracting: // cabs(zp - zcr_last ) < PixelWidth 
    AR = 30.0* PixelWidth * iWidth / 2000 ; // 
    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(z1, zp, z3); // test if triangle is oriented counterclockwise
				
    AR = cabs(zp-z1);
    R_BDM = cabs(-0.0559603266109820 -0.1163177312664140*I); // it is a point from BDM trap // 20.0* PixelWidth * iWidth / 2000 ; // 

				
    break;
  default: break;
  }
	
	
	
	
	
  AR2 = AR*AR;
  R2_BDM = R_BDM* R_BDM;
	
  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 = 2000.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;
  
  
  t0 = ((double) t_numerator)/t_denominator; //1.0 / period; // Is it iternal angle from internal adress  ???
  
	// DEM
	BoundaryWidth = 1.0*iWidth/2000.0;
	distanceMax = BoundaryWidth*PixelWidth; // here boundary is changing with resolution, maybe % of ImageWidth would be better ? 
  
  
  
  
  /* 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, 11*IterMax_LSM, 2, data);
    SaveArray2PGMFile (data, c, "FatouBasins_CRO" , "boundaries of Level set method ( LSM) =  Level Curve Method (LCM) and critical orbit ");
	
    //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");
	
    DrawForwardOrbit(zcr, 11*IterMax_LSM, 4, data);
    SaveArray2PGMFile (data, c, "FatouComponents_LSCM_CRO" , "FatouComponents_LSCM and 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, "LSET" , iHeight, " Level sets of integer escape time of ET  ");

     ComputeBoundaries(data,edge2);
     SaveArray2PGMFile (edge2, c, "LCET" , iHeight, "boundaries of integer Escape Time =  Level Curve ( LC) of integer Escape Time");
	
     CopyBoundaries(edge2,  data);
     SaveArray2PGMFile (data, c, "LSLCET" , iHeight, " Level sets and it's baundaries ( Level Curves = LC) of integer escape time of ET  ");
	


  DrawImage (data, BDM);	
  SaveArray2PGMFile (data, c,  "BDM" , iHeight, "BDM = Binary Decomposition Method for both exterior and interior = Level Sets of BDM");
	
  ComputeBoundaries(data,edge);
  SaveArray2PGMFile (edge, c, "BDM_LC", iHeight, "boundaries of Binary Decomposition Method (LC of  BD) = LCBD");
	
  CopyBoundaries( edge, data);
  SaveArray2PGMFile (data, c,  "BDM_LC_LS" , iHeight, " Binary Decomposition Method for both exterior and interior ");
	
	*/ 
  //CopyBoundaries(edge2,  data); //  if LSM was done !!!!
  //SaveArray2PGMFile (data, c, "BDM_LSLCET" , " Level sets and it's baundaries ( Level Curves = LC) of integer escape time of ET  ");
	
 // Mark_parabolic_BDM_Triangle_Traps(data);
  //SaveArray2PGMFile (data, c, "BDM_LC_LS_traps", iHeight, "BDM and traps");
	
	
  //CopyBoundaries(edge2,  edge);
  //SaveArray2PGMFile (edge, c, "LCBDET ", "boundaries of Binary Decomposition Method and Level sets of Integer Escape Time");

  //Mark_Parabolic_Components_Triangle_Trap(edge2);
  //SaveArray2PGMFile (edge2,  c, "LSM_CRO_trap" , "LSM, critical orbit and trap");

  //CopyBoundaries(edge2, data);
  //SaveArray2PGMFile (data,  "LSCM" , "LSCM");
  
  
  DrawImage (data, DEM);	
  SaveArray2PGMFile (data, c,  "DEM" , iHeight, "DEM/J = Distance Estimation Method for both exterior and interior of Julia set");
  
  DrawImage (data, BDM);
  DrawBoundary_using_DEMJ ( data );
  SaveArray2PGMFile (data, c,  "DEM_BDM" , iHeight, "BDM and DEM/J = Distance Estimation Method for both exterior and interior of Julia set");  
  
  
  return 0;
}


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

int main(void)
{

  DynamicType = dendrite;  // setup DynamicType value manually  ; One can do it also numerically ( from multiplier or from some properities)
  general_setup();

	
  local_setup();
		
  MakeImages();
  PrintProgramInfo();
			
  end();
  return 0;
}