/*

  https://web.archive.org/web/20161024194536/http://www.ijon.de/mathe/julia/some_julia_sets_3.html

  a: (0,15+0,15i)
  b: (-3+3i)
  f: 1/(a*z^5+ b*z^3 + z);
  dz = -(5*a*z^4+3*b*z^2+1)/(a*z^5+b*z^3+z)^2


  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 d.c -lm -Wall -march=native -fopenmp
  time ./a.out > b.txt


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


  time ./a.out

  time ./a.out >i.txt
  time ./a.out >e.txt
  
  
  
  
  
  
  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 ------------------------------------------------------------ */


//FunctionType
typedef enum  {Fatou = 0, IntLSM =1 , ExtLSM = 2, LSM = 3, DEM = 4, Unknown = 5 , BD = 6, MBD = 7 , SAC = 8, DLD = 9, ND = 10 , NP= 11, POT = 12 , Blend = 13, White= 14, Fatou_ab = 15, Fatou_abi = 16, LSM_m = 17 
		
} FunctionTypeT; 
// FunctionTypeT FunctionType;

// virtual 2D array and integer ( screen) coordinate
// Indexes of array starts from 0 not 1 
//unsigned int ix, iy; // var
static unsigned int ixMin = 0;	// Indexes of array starts from 0 not 1
static unsigned int ixMax;	//
static unsigned int iWidth;	// horizontal dimension of array

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

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

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

// unsigned int i; // var = index of 1D array
//static unsigned int iMin = 0; // Indexes of array starts from 0 not 1
static unsigned int iMax;	// = i2Dsize-1  = 
// The size of array has to be a positive constant integer 
// unsigned int i1Dsize ; // = i2Dsize  = (iMax -iMin + 1) =  ;  1D array with the same size as 2D array



// see SetPlane

double radius = 3.0; 
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; 


/*
  ER = pow(10,ERe);
  AR = pow(10,-ARe);
*/
//int ARe ;			// increase ARe until black ( unknown) points disapear 
//int ERe ;
double ER;
double ER2;			//= 1e60;
double AR; // bigger values do not works
double AR2;

double AR_max;
//double AR12;



int IterMax = 100000;
int IterMax_LSM = 1000;
int IterMax_DEM = 100000;

/* colors = shades of gray from 0 to 255 

   unsigned char colorArray[2][2]={{255,231},    {123,99}};
   color = 245;  exterior 
*/
unsigned char iColorOfExterior = 245;
unsigned char iColorOfInterior = 99;
unsigned char iColorOfInterior1 = 100;
unsigned char iColorOfInterior2 = 183;
unsigned char iColorOfBoundary = 0;
unsigned char iColorOfUnknown = 5;

// pixel counters
unsigned long long int uUnknown = 0;
unsigned long long int uInterior = 0;
unsigned long long int uExterior = 0;



// critical points




// critical points
complex double zcr1; //
complex double zcr2;// = -2.2351741790771484375e-08+9.4296410679817199707e-09*I;
complex double zcr3; //
complex double zcr4;// 

const complex double z_cr[4]= {0.389374737779177*I-0.9400337727919567,
  0.9400337727919567-0.389374737779177*I,
  -1.816005797447402*I-0.7522142306508819,
  1.816005797447402*I+0.7522142306508819};

const int period = 4;

// periodic points = attractors
//complex double z1 =  0.0 ; //fixed point (period  1) =  attracting cycle

/*
  attracting periodic points : 
  2 period 4 cycles found by marcm200
  https://fractalforums.org/fractal-mathematics-and-new-theories/28/rational-function/4279/msg29227#msg29227
*/
const complex double zp4a[4]= { -0.753562725059588878195881989086-1.81926135093768714945383635495*I, 0.110964246025498911030204851613-0.0459628956422664519676501981849*I , 0.0951541874285335154137754898329-0.0394141549494900420014253938916*I, -0.877971698675046430260238139454-2.11961118232104128722426139575*I};
const complex double zp4b[4]= {0.753562725059588878195881989086+1.81926135093768714945383635495*I, -0.0951541874285335154137754898329+0.0394141549494900420014253938916*I,  0.877971698675046430260238139454+2.11961118232104128722426139575*I, -0.110964246025498911030204851613+0.0459628956422664519676501981849*I};

const complex double zp4a0 = -0.753562725059588878195881989086-1.81926135093768714945383635495*I;
const complex double zp4b0 = 0.753562725059588878195881989086+1.81926135093768714945383635495*I;


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

/* 
   a: (0.15+0.15*%i);
   b: (-3+3*%i);
   f: 1/(a*z^5+ b*z^3 + z);
   dz = -(5*a*z^4+3*b*z^2+1)/(a*z^5+b*z^3+z)^2

*/
const complex double a =  0.15 + 0.15*I;
const complex double b = -3.0  + 3.0*I;

// complex function
complex double f(const complex double z0) {

  double complex z = z0;
  complex double z2 = z*z;
  complex double z3 = z2*z;
  complex double z5 = z3*z2;
  z = 1.0/(a*z5 + z3 + b*z);
  return  z;
}
	

// 
complex double dfz(const complex double z0) {


  // dz= -(5*a*z^4+3*z^2+b) / (a*z^5+z^3+b*z)^2
  double complex z = z0;
  complex double z2= z*z;
  complex double z4 = z2*z2;
  complex double numerator = -5.0*a*z4 - 3.0*z2 - b ;
  complex double denom = a*z4*z + z2*z + b*z;
  denom = denom*denom; // ^2
	
	
  return  numerator/denom;
	
}
		














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




}







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

double cabs2(complex double z){

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


}






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












/* 

is it possible to adjust AR so that level curves in interior have figure 8?

find such AR for internal LCM/J and LSM that level curves croses critical point and it's preimages
for attracting ( also weakly attracting = parabolic) dynamics

it may fail 
* if one iteration is bigger then smallest distance between periodic point zp and Julia set
* if critical point is attracted by another cycye ( then change periodic point zp)

Made with help of Claude Heiland-Allen

*/
double GiveTunedAR(const int iter_Max){

  fprintf(stdout, " GiveTunedAR\n");

  complex double z = zcr1; // initial point z0 = criical point 
  int iter;
  double r = 50 * PixelWidth; // initial value 
  double t;
  
  // iterate critical point
	for (iter=0; iter< iter_Max; ++iter ){
		// check attractor from first basin  
		t = cabs(zp4a0 - z);
		if ( t < r)
		{
      			r = t;
			break;
		}
		// check second basin
		t = cabs(zp4b0 - z);
		if (t < r)
		{
			r = t;
			break;
		}
  	
  	z = f(z); // forward iteration
  	
  
  }
  // check distance between zn = f^n(zcr) and periodic point zp
  
  fprintf(stdout, "  r = %f = %d * pixeWidth \n",  r, (int) (r/PixelWidth));
  
  
 // use it as a AR
 return r;
	
	
}



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


// ???????

int IsInterior(complex double z){

  if (
      cabs2(zp4a0-z) < AR2 || 
      cabs2(zp4b0-z) < AR2 
      )
    {return 1;}
  return 0;



}





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

*/
unsigned char ComputeColorOfFatou (complex double z)
{



	
	
  double r2;


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


		
	
      r2 =cabs2(z);
		
      if (r2 > ER2) // esaping = exterior
	{
	  uExterior += 1;
	  return iColorOfExterior;
	}
	
	
      if (IsInterior(z)) // 
	{
	  
	  return iColorOfInterior;
	}	
				
	
     
      z = f(z);		//  iteration: z(n+1) = f(zn)
	

    }

  
  return iColorOfUnknown;


}



/*
  3 basins 
  - exterior
  - interior ( consist of 2 attracting basins)
  - - basin 1 
  - - basin 2
  - unknown ( possibly empty set ) 

*/

unsigned char ComputeColorOfFatou_ab (complex double z)
{



	
	
  double r2;


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


		
	
      r2 =cabs2(z);
		
      if (r2 > ER2) // esaping = exterior
	{
	  uExterior += 1;
	  return iColorOfExterior;
	}
	
	
      //Interior = 2 Attraction basins 
      if ( cabs2(zp4a0-z) < AR2 ){ return iColorOfInterior1;}
	 
      if (cabs2(zp4b0-z) < AR2 ) { return iColorOfInterior2;}			
	
     
      z = f(z);		//  iteration: z(n+1) = f(zn)
	

    }

  
  return iColorOfUnknown;


}


/*
  3 basins 
  - exterior
  - interior ( consist of 2 attracting basins) each basin is colored according to position in the cycle
  - - basin 1 
  - - basin 2
  - unknown ( possibly empty set ) 

*/

unsigned char ComputeColorOfFatou_abi (complex double z)
{



	
	
  double r2;


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


		
	
      r2 =cabs2(z);
		
      if (r2 > ER2) // esaping = exterior
	{
	  uExterior += 1;
	  return iColorOfExterior;
	}
	
	
      //Interior = 2 Attraction basins 
      if ( cabs2(zp4a0-z) < AR2 ){ return iColorOfInterior1 - (i % period)*20;}
	 
      if (cabs2(zp4b0-z) < AR2 ) { return iColorOfInterior2 - (i % period)*19;}			
	
     
      z = f(z);		//  iteration: z(n+1) = f(zn)
	

    }

  
  return iColorOfUnknown;


}






unsigned char ComputeColorOfLSM (complex double z)
{



	
	
  double r2;


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


		

      // complex iteration f(z)=z^3 + c
      r2 =cabs2(z);
		
      if (r2 > ER2) // esaping = exterior
	{
	  uExterior += 1;
	  return 255- (i % 255);
	}			
      if (IsInterior(z)) // 
	{
	  
	  return 255- ((2*i) % 255);
	}	
	
      z = f(z);	

    }

  return iColorOfUnknown;


}






unsigned char ComputeColorOfLSM_m (complex double z)
{



	
	
  double r2;


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


		

      // complex iteration f(z)=z^3 + c
      r2 =cabs2(z);
		
      if (r2 > ER2) // esaping = exterior
	{
	  uExterior += 1;
	  return 255- (i % 255);
	}			
	
      if (IsInterior(z)) // 
	{
	  
	  return  (i % 255);
	}	
	
      z = f(z);	

    }

  return iColorOfUnknown; //iColorOfUnknown;


}













// ***************************************************************************************************************************
// ************************** DEM/J*****************************************
// ****************************************************************************************************************************

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
    cabsz = cabs(z);
    if (cabsz > 1e60 || cabs(dz)> 1e60) { break; }// big values 
    if (IsInterior(z)) { return iColorOfInterior2;} // falls into finite attractor = interior
  			
    dz = dfz(z)*dz; 
    z = f(z) ; /* forward iteration : complex cubic polynomial */ 
  }
  
  
  distance = 2.0 * cabsz* log(cabsz)/ cabs(dz);
  if (distance <distanceMax) return iColorOfBoundary; // distanceMax = BoundaryWidth*PixelWidth;
  // else
  
  return iColorOfExterior;

 
}






/* ==================================================================================================
   ============================= Draw functions ===============================================================
   =====================================================================================================
*/ 
unsigned char ComputeColor(FunctionTypeT FunctionType, complex double z){

  unsigned char iColor;
	
	
	
  switch(FunctionType){
  
  case Fatou :{iColor = ComputeColorOfFatou(z); break;}
  	
  case Fatou_ab :{iColor = ComputeColorOfFatou_ab(z); break;}
  	
  case Fatou_abi :{iColor = ComputeColorOfFatou_abi(z); break;}
  
 
    // case IntLSM :{iColor = ComputeColorOfIntLSM(z); break;}
	
    // case ExtLSM :{iColor = ComputeColorOfExtLSM(z); break;}
  
  case LSM :{iColor = ComputeColorOfLSM(z); break;}
  
  	
  case LSM_m :{iColor = ComputeColorOfLSM_m(z); break;}
		
  case DEM : {iColor = ComputeColorOfDEMJ(z); break;}
	
    /*	
  	case Unknown : {iColor = ComputeColorOfUnknown(z); break;}
		
  	case BD : {iColor = ComputeColorOfBD(z); break;}
		
  	case MBD : {iColor = ComputeColorOfMBD(z); break;}
		
  	case SAC : {iColor = ComputeColorOfSAC(z); break;}
  
  	case DLD : {iColor = ComputeColorOfDLD(z); break;}
		
  	case ND : {iColor = ComputeColorOfND(z); break;}
		
  	case NP : {iColor = ComputeColorOfNP(z); break;}
		
  	case POT : {iColor = ComputeColorOfPOT(z); break;}
		
  	case Blend : {iColor = ComputeColorOfBlend(z); break;}
    */	
  	
  case White : {iColor = 255;}
  	
  	
	
  default: {}
	
	
  }
	
  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 */
  
  z = GiveZ(ix,iy);
  

  iColor = ComputeColor(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 \n", FunctionType);
  // 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 */
	
	
  A[i]= 0; //255-A[i]; // Mark point with inveres color
	
	
  return 0;
	
}






int IsInside (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;
	  

} 


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

	
  unsigned int ix_seed = (creal(z)-ZxMin)/PixelWidth;
  unsigned int iy_seed = (ZyMax - cimag(z))/PixelHeight;
  unsigned int i;
	
	
  /* mark seed point by big pixel */
  int iSide =4.0*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 (IsInside(iX, iY, ix_seed, iy_seed, iSide)) {
	i= Give_i(iX,iY); /* index of _data array */
	A[i]= 0; //255-A[i];
      }
      else {printf(" bad point \n");}
	
    }}
	
	
  return 0;
	
}



int PlotAllPoints(const complex double zz[], int kMax, unsigned char A[]){

  int k;
	
	
  printf("kMax = %d \n",kMax);
	

  for (k = 0; k < kMax; ++k)
    {
      //fprintf(stderr, "z = %+f %+f \n", creal(zz[k]),cimag(zz[k]));
      PlotBigPoint(zz[k], A);}
  return 0;





}




int DrawForwardOrbit(complex double z, unsigned long long int iMax,  unsigned char A[] )
{
  
  unsigned long long int i; /* nr of point of critical orbit */
  printf("draw forward orbit \n");
 
  PlotBigPoint(z, A);
  
  /* forward orbit of critical point  */
  for (i=1;i<iMax ; ++i)
    {
      z  = f(z);
      //if (cabs2(z - z2a) > 2.0) {return 1;} // escaping
      PlotBigPoint(z, 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, iColorOfExterior, iSize*sizeof(*D)); // for heap-allocated arrays, where N is the number of elements = FillArrayWithColor(D , iColorOfExterior);
 
  // 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;
}

















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

int SaveArray2PGMFile (unsigned char A[],  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", n );	/*  */
  char *filename = strcat (name, ".pgm");
  char long_comment[200];
  sprintf (long_comment, "Julia set f(z)=1/((0,15+0,15i)z5+z3+(-3+3i)z) on [-3;3]x[-3;3]. Location by Michael Becker %s", comment);





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

  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
  printf ("Numerical approximation of Julia set for F(z) =   z*c) \n");
  //printf ("parameter C = ( %.16f ; %.16f ) \n", creal (C), cimag (C));
  

  printf ("Image Width = %f in world coordinate\n", ZxMax - ZxMin);
  printf ("PixelWidth = %.16f \n", PixelWidth);
  //printf ("AR = %.16f = %f *PixelWidth = %f %% of ImageWidth \n", AR, AR / PixelWidth, AR / ZxMax - ZxMin);



  // image corners in world coordinate
  // center and radius
  // center and zoom
  // GradientRepetition
  printf ("Maximal number of iterations = iterMax = %d \n", IterMax);
  printf ("ratio of image  = %f ; it should be 1.000 ...\n", ratio);
  //




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

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






  /* 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 ...
  
  // critical points
  zcr1 = z_cr[0];
  zcr2 = z_cr[1];
  zcr3 = z_cr[2];
  zcr4 = z_cr[3];
  
  
  // LSM  
  // escape radius ( of circle around infinity 
  ER = 200.0; // 
  ER2 = ER*ER;
  
  
  // attracting radius of ciurcel arounf finite attractor
  //AR_max = 5*PixelWidth*iWidth/2000.0 ; // adjust first number 
  // GiveTunedAR(const int i_Max, const complex double zcr, const double c, const double zp){
  //AR = 50*PixelWidth; // 0.03; // 10*0.0006 = 0.006
  AR = GiveTunedAR(10000); 
  
  AR2 = AR * AR;
  //AR12 = AR/2.0;
  
  
  
  
  
  // 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));
  if (data == NULL || edge == NULL)
    {
      fprintf (stderr, " Could not allocate memory");
      return 1;
    }





 


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

  return 0;

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




int end ()
{


  fprintf (stderr, " allways free memory (deallocate )  to avoid memory leaks \n");	// https://en.wikipedia.org/wiki/C_dynamic_memory_allocation
  free (data);
  free(edge);


  PrintProgramInfo ();
  PrintCInfo ();
  return 0;

}

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

int main ()
{
  setup ();
  
  /*
    DrawImage (data, Fatou);	
    SaveArray2PGMFile (data,  "Fatou" , "Fatou ");

    ComputeBoundaries(data,edge);
    SaveArray2PGMFile (edge,  "LCM_Fatou" , "LCM of Fatou ");  
  
  
    DrawImage (data, Fatou_ab);	 
    SaveArray2PGMFile (data,  "Fatou_ab" , "Fatou_ab ");
  
    DrawImage (data, Fatou_abi);	 
    SaveArray2PGMFile (data,  "Fatou_abi" , "Fatou_abi ");
  */
  
    DrawImage (data, LSM);	// first 
    //SaveArray2PGMFile (data,  "LSM" , "LSM");
  
    ComputeBoundaries(data,edge);
    SaveArray2PGMFile (edge,  "LCM" , "LCM ");
    
    DrawForwardOrbit(zcr1, 200,  edge);
    SaveArray2PGMFile (edge,  "LCM_cr" , "LCM and critical orbit");
    
  /*
    CopyBoundaries(edge, data);
    SaveArray2PGMFile (data,  "LSCM" , "LSCM");
  
  
 
 

 
    PlotAllPoints(z_cr, 4, data);
    SaveArray2PGMFile (data,  "Fatou_cr" , "Fatou_cr");
  
  
 
    DrawImage (data, Fatou);	 
    PlotBigPoint(zp4a0, data);
    PlotBigPoint(zp4b0, data);
    SaveArray2PGMFile (data,  "Fatou_zp4" , "Fatou_zp4");
  
  
  
  
  
    DrawImage (data, LSM_m);	// first 
    SaveArray2PGMFile (data,  "LSM_m2" , "LSM_m ");
  
    ComputeBoundaries(data,edge);
    SaveArray2PGMFile (edge,  "LCM_m2_LS" , "LCM_m ");
  
    CopyBoundaries(edge, data);
    SaveArray2PGMFile (data,  "LSCM_mm2" , "LSCM mm");
  
  
  DrawImage (data, DEM);	// first 
  SaveArray2PGMFile (data,  "DEM" , "DEM ");
  */

  
  end ();

  return 0;
}