/*

  Adam Majewski
  adammaj1 aaattt o2 dot pl  // o like oxygen not 0 like zero 
  
  
  console program in c programing language 
===============================================================
to do:
[[:File:Julia_IIM_6_basilica.png]]
Modified binary decomposition of dynamical plane


parabolic chesboard
[[:File:A_parabolic_checkerboard_for_rotation_number_1_over_2.png]]



https://arxiv.org/pdf/0802.3248.pdf





  
  ==============================================
  
  
  Structure of a program or how to analyze the program 
  
  
  ============== Image X ========================
  
  DrawImageOfX -> DrawPointOfX -> ComputeColorOfX 
  
  first 2 functions are identical for every X
  check only last function =  ComputeColorOfX
  which computes color of one pixel !
  
  

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

  
  ---------------------------------
  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 >a.txt

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


  

*/

#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

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


// see SetPlane

double radius = 1.4; 
complex double center = 0.0;
double  DisplayAspectRatio  = 1.0; // https://en.wikipedia.org/wiki/Aspect_ratio_(image)

// 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 = 10000;	//  
// The size of array has to be a positive constant integer 
static unsigned int iSize;	// = iWidth*iHeight; 

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

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


double ZxMin ;	//-0.05;
double ZxMax ;	//0.75;
double ZyMin ;	//-0.1;
double ZyMax ;	//0.7;
double PixelWidth;	// =(ZxMax-ZxMin)/ixMax;
double PixelHeight;	// =(ZyMax-ZyMin)/iyMax;
double ratio;


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




static unsigned long int iterMax = 1000000;	//iHeight*100;
unsigned long int iterMax_LSM = 255;


double ER = 200.0;		// EscapeRadius for bailout test 
double EscapeRadius=1000000; // = ER big !!!!
double ER_LSM ; // see GiveER_LSM  // 27.764 =  manually find value such that level curves of escape time cross critical point and it's  preimages
double ER_DLD ; // see GiveER_LSM  // 27.764 =  manually find value such that level curves of escape time cross critical point and it's  preimages


// SAC/J
double lnER; // ln(ER)
int i_skip = 2; // exclude (i_skip+1) elements from average
unsigned char s = 7; // stripe density

double BoundaryWidth = 3.0; // % of image width  
double distanceMax; //distanceMax = BoundaryWidth*PixelWidth;



//  ------------- DLD  ----------------------
const int N = 20;		// fixed number : maximal number of iterations
double p = 0.180; //0.01444322;		//
// DLD colors
//double me = 1.0;
double mi = 0.9;





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





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





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





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

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


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


}




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



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

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








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


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





// ***************************************************************************************************************************
// ************************** 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;
  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 (cabsz< PixelWidth) return iColorOfInterior; // falls into finite attractor = interior
  			
    dz = 2.0*z * dz; 
    z = z*z +c ; /* forward iteration : complex quadratic polynomial */ 
  }
  
  
  distance = 2.0 * cabsz* log(cabsz)/ cabs(dz);
  if (distance <distanceMax) return iColorOfBoundary; // distanceMax = BoundaryWidth*PixelWidth;
  // else
  
  return iColorOfExterior;

 
}



// plots raster point (ix,iy) 
int DrawPointOfDEMJ (unsigned char A[], 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 = ComputeColorOfDEMJ(z);
  A[i] = iColor ;		// interior
  
  return 0;
}




// fill array 
// uses global var :  ...
// scanning complex plane 
int DrawImagerOfDEMJ (unsigned char A[])
{
  unsigned int ix, iy;		// pixel coordinate 

  	//printf("compute image \n");
 	// for all pixels of image 
	#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax)
  	for (iy = iyMin; iy <= iyMax; ++iy){
    		printf (" %d from %d \r", iy, iyMax);	//info 
    		for (ix = ixMin; ix <= ixMax; ++ix)
      			DrawPointOfDEMJ(A, ix, iy);	//  
  }

  return 0;
}




// ***************************************************************************************************************************
// ************************** only boundary by  DEM/J*****************************************
// ****************************************************************************************************************************

unsigned char ComputeColorOfDEMJ_boundary(complex double z){
// https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Julia_set#DEM.2FJ


  
  int nMax = iterMax;
  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 (cabsz< PixelWidth) return iColorOfInterior; // falls into finite attractor = interior
  			
    dz = 2.0*z * dz; 
    z = z*z +c ; /* forward iteration : complex quadratic polynomial */ 
  }
  
  
  distance = 2.0 * cabsz* log(cabsz)/ cabs(dz);
  if (distance <distanceMax) return iColorOfBoundary; // distanceMax = BoundaryWidth*PixelWidth;
  // else
  
  return iColorOfExterior;

 
}



// plots raster point (ix,iy) 
int DrawPointOfDEMJ_boundary (unsigned char A[], 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 = ComputeColorOfDEMJ_boundary(z);
  if (iColor == iColorOfBoundary) 
  	{ A[i] = iColor ;} // draw only boundary without changing other parts 		
  
  return 0;
}




// fill array 
// uses global var :  ...
// scanning complex plane 
int DrawImagerOfDEMJ_boundary (unsigned char A[])
{
  unsigned int ix, iy;		// pixel coordinate 

  	//printf("compute image \n");
 	// for all pixels of image 
	#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax)
  	for (iy = iyMin; iy <= iyMax; ++iy){
    		printf (" %d from %d \r", iy, iyMax);	//info 
    		for (ix = ixMin; ix <= ixMax; ++ix)
      			DrawPointOfDEMJ_boundary(A, ix, iy);	//  
  }

  return 0;
}











// ***************************************************************************************************************************
// ************************** Unknown: boundary and slow dynamics *****************************************
// ****************************************************************************************************************************

unsigned char ComputeColorOfUnknown(complex double z){



  
  int nMax = 20; // very low value
  
  double cabsz;
	
  int n;

  for (n=0; n < nMax; n++){ //forward iteration
	cabsz = cabs(z);
    	if (cabsz > 10000000000*ER )  return iColorOfExterior; // big values
    	if (cabsz < (PixelWidth/100)) return iColorOfInterior; // falls into finite attractor = interior
  			
    
    z = z*z +c ; /* forward iteration : complex quadratic polynomial */ 
  }
  
  
  
  
  //printf("found \n");
  return iColorOfUnknown;

 
}



// plots raster point (ix,iy) 
int DrawPointOfUnknown (unsigned char A[], 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 = ComputeColorOfUnknown(z);
  A[i] = iColor ;		// interior
  
  return 0;
}




// fill array 
// uses global var :  ...
// scanning complex plane 
int DrawImagerOfUnknown (unsigned char A[])
{
  unsigned int ix, iy;		// pixel coordinate 

  	//printf("compute image \n");
 	// for all pixels of image 
	#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax)
  	for (iy = iyMin; iy <= iyMax; ++iy){
    		//printf (" %d from %d \r", iy, iyMax);	//info 
    		for (ix = ixMin; ix <= ixMax; ++ix)
      			DrawPointOfUnknown(A, ix, iy);	//  
  }

  return 0;
}





// ***************************************************************************************************************************
// ************************** LSM/J*****************************************
// ****************************************************************************************************************************

unsigned char ComputeColorOfLSM(complex double z){

 int nMax = iterMax_LSM;
  double cabsz;
  unsigned char iColor;
	
  int n;

  for (n=0; n < nMax; n++){ //forward iteration
	cabsz = cabs(z);
    	if (cabsz > ER_LSM) break; // esacping
    	//if (cabsz< PixelWidth) break; // fails into finite attractor = interior, but not for disconnected Julia sets, then critical point and its preimages  !!!!
  			
   
     	z = z*z +c ; /* forward iteration : complex quadratic polynomial */ 
  }
  
  // manually udjusted series of ordered colors ( shades of gray )
  iColor = 255 - 230.0*((double) n)/18.0; // nMax or lower values in denominator
  
  
  return iColor;


}



// plots raster point (ix,iy) 
int DrawPointOfLSM (unsigned char A[], 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 = ComputeColorOfLSM(z);
  A[i] = iColor ;		// interior
  
  return 0;
}




// fill array 
// uses global var :  ...
// scanning complex plane 
int DrawImagerOfLSM (unsigned char A[])
{
  unsigned int ix, iy;		// pixel coordinate 

  	//printf("compute image \n");
 	// for all pixels of image 
	#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax)
  	for (iy = iyMin; iy <= iyMax; ++iy){
    		printf (" %d from %d \r", iy, iyMax);	//info 
    		for (ix = ixMin; ix <= ixMax; ++ix)
      			DrawPointOfLSM(A, ix, iy);	//  
  }

  return 0;
}




// ***************************************************************************************************************************
// ************************** binary decomposition BD/J*****************************************
// ****************************************************************************************************************************

unsigned char ComputeColorOfBD(complex double z){

 int nMax = iterMax_LSM;
  double cabsz;
  unsigned char iColor;
	
  int n;

  for (n=0; n < nMax; n++){ //forward iteration
	cabsz = cabs(z);
    	if (cabsz > ER_LSM) break; // esacping
    	//if (cabsz< PixelWidth) break; // fails into finite attractor = interior but not for disconnected Julia sets, then critical point and its preimages  !!!!
  			
   
     	z = z*z +c ; /* forward iteration : complex quadratic polynomial */ 
  }
  
  if (cimag(z)>0.0) 
  	iColor = 255; 
  	else iColor = 0;
  
  
  return iColor;


}



// plots raster point (ix,iy) 
int DrawPointOfBD (unsigned char A[], 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 = ComputeColorOfBD(z);
  A[i] = iColor ;		// interior
  
  return 0;
}




// fill array 
// uses global var :  ...
// scanning complex plane 
int DrawImagerOfBD (unsigned char A[])
{
  unsigned int ix, iy;		// pixel coordinate 

  	//printf("compute image \n");
 	// for all pixels of image 
	#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax)
  	for (iy = iyMin; iy <= iyMax; ++iy){
    		printf (" %d from %d \r", iy, iyMax);	//info 
    		for (ix = ixMin; ix <= ixMax; ++ix)
      			DrawPointOfBD(A, ix, iy);	//  
  }

  return 0;
}





// ***************************************************************************************************************************
// ************************** modified binary decomposition BD/J*****************************************
// ****************************************************************************************************************************

unsigned char ComputeColorOfMBD(complex double z){
// const number of iterations
 int nMax = 7;
  //double cabsz;
  unsigned char iColor;
	
  int n;

  for (n=0; n < nMax; n++){ //forward iteration
	//cabsz = cabs(z);
    	//if (cabsz > ER) break; // esacping
    	//if (cabsz< PixelWidth) break; // falls into finite attractor = interior
  			
   
     	z = z*z +c ; /* forward iteration : complex quadratic polynomial */ 
  }
  
  //if (cabs(z) > 2.0)
  	{ // exterior
  		if (creal(z)>0.0) 
  			iColor = 255; 
  			else iColor = 0;
  	}
  //	else iColor = iColorOfInterior;
  	
  return iColor;


}



// plots raster point (ix,iy) 
int DrawPointOfMBD (unsigned char A[], 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 = ComputeColorOfMBD(z);
  A[i] = iColor ;		// interior
  
  return 0;
}




// fill array 
// uses global var :  ...
// scanning complex plane 
int DrawImagerOMfBD (unsigned char A[])
{
  unsigned int ix, iy;		// pixel coordinate 

  	//printf("compute image \n");
 	// for all pixels of image 
	#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax)
  	for (iy = iyMin; iy <= iyMax; ++iy){
    		printf (" %d from %d \r", iy, iyMax);	//info 
    		for (ix = ixMin; ix <= ixMax; ++ix)
      			DrawPointOfMBD(A, ix, iy);	//  
  }

  return 0;
}




// ***********************************************************************************************
//*************************************** SAC/J **************************************************
// *****************************************************************************************
// https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/stripeAC
// SAC = Stripe Average Coloring

//

// the addend function
// input : complex number z
// output : double number t 
double Give_t(double complex z){

  return 0.5+0.5*sin(s*carg(z));

}

/*
  input :
  - complex number
  - intege
  output = average
 
*/
double Give_Arg(double complex z , int iMax)
{
  int i=0; // iteration 
   
   
  //double complex Z= 0.0; // initial value for iteration Z0
  double A = 0.0; // A(n)
  double prevA = 0.0; // A(n-1)
  double R; // =radius = cabs(Z)
  double d; // smooth iteration count
  double complex dz = 1.0; // first derivative with respect to z
  double de; // Distance Estimation from DEM/J  
   
    
  // iteration = computing the orbit
  for(i=0;i<iMax;i++)
    { 
    
      dz = 2.0 * z * dz ; 
      z = z*z + c; // https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/qpolynomials
      
      if (i>i_skip) A += Give_t(z); // 
      
      R = cabs(z);
      // if(R > EscapeRadius) break; // exterior of M set
  	if (R > 1e60 || cabs(dz)> 1e60) break; // prevent NAN 	 	
      prevA = A; // save value for interpolation
        
    } // for(i=0
   
   
  if (i == iMax) 
    A = -1.0; // interior 
  else { // exterior
    de = 2 * R * log(R) / cabs(dz);
    if (de < distanceMax) A = FP_ZERO; //  boundary
    else {
      // computing interpolated average
      A /= (i - i_skip) ; // A(n)
      prevA /= (i - i_skip - 1) ; // A(n-1) 
      // smooth iteration count
      d = i + 1 + log(lnER/log(R))/M_LN2;
      d = d - (int)d; // only fractional part = interpolation coefficient
      // linear interpolation
      A = d*A + (1.0-d)*prevA;
     }   
  }
    
  return A;  
}
 
 
 
 
 
unsigned char ComputeColorOfSAC(complex double z){

  unsigned char iColor;
  double arg;
  
   
   
  	arg = Give_Arg( z, 2500); //   N in wiki 
	
   	// color is proportional to arg 
	if (arg < 0.0)
           
		iColor = 0;  // interior                        
    
		else //  
			{if (arg == FP_ZERO) 
     				iColor = 255; // boundary     
        			else iColor = (unsigned char) (255 - 255*arg );// exterior
      			}
      
    
  return iColor;


}



// plots raster point (ix,iy) 
int DrawPointOfSAC (unsigned char A[], 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 = ComputeColorOfSAC(z);
  A[i] = iColor ;		//   
  return 0;
}




// fill array 
// uses global var :  ...
// scanning complex plane 
int DrawImagerOMfSAC (unsigned char A[])
{
  unsigned int ix, iy;		// pixel coordinate 

  	//printf("compute image \n");
 	// for all pixels of image 
	#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax)
  	for (iy = iyMin; iy <= iyMax; ++iy){
    		printf ("SAC/J :  %d from %d \r", iy, iyMax);	//info 
    		for (ix = ixMin; ix <= ixMax; ++ix)
      			DrawPointOfSAC(A, ix, iy);	//  
  }

  return 0;
}



 

// ***************************************************************************************************************************
// ************************** DLD/J*****************************************
// ****************************************************************************************************************************



/* partial pnorm 
   input: z , zn = f(z), p
   output ppn
   
   
*/
double
ppnorm (complex double z, complex double zn, double p)
{

  double s[2][3];		// array for 2 points on the Riemann sphere
  int j;
  double d;			// denominator 
  double x;
  double y;

  double ds;
  double ppn = 0.0;

  // map from complex plane to riemann sphere
  // z
  x = creal (z);
  y = cimag (z);
  d = x * x + y * y + 1.0;

  s[0][0] = (2.0 * x) / d;
  s[0][1] = (2.0 * y) / d;
  s[0][2] = (d - 2.0) / d;	// (x^2 + y^2 - 1)/d

  // zn
  x = creal (zn);
  y = cimag (zn);
  d = x * x + y * y + 1.0;
  s[1][0] = (2.0 * x) / d;
  s[1][1] = (2.0 * y) / d;
  s[1][2] = (d - 2.0) / d;	// (x^2 + y^2 - 1)/d

  // sum 
  for (j = 0; j < 3; ++j)
    {
      ds = fabs (s[1][j] - s[0][j]);
      //  normal:  neither zero, subnormal, infinite, nor NaN
      //if (fpclassify (ds) !=FP_INFINITE)
      //if (isnormal(ds)) 
      // it is solved by if (cabs(z) > 1e60 ) break; procedure in parent function 
      ppn += pow (ds, p);	// |ds|^p
      //      else {ppn = 10000.0; printf("ds = infty\t");} // 

    }


  return ppn;







}

// DLD = Discret Lagrangian Descriptior
double
lagrangian (complex double z0, complex double c, int iMax, double p)
{

  int i;			// number of iteration
  double d = 0.0;		// DLD = sum
  double ppn;			// partial pnorm
  complex double z = z0;
  complex double zn;		// next z

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




      zn = z * z + c;		// complex iteration
      ppn = ppnorm (z, zn, p);
      d += ppn;			// sum
      //
      z = zn;

      //if (! isnormal(d)) { return 0.0; } // not works
      if (cabs (z) > ER_DLD ) //1e6)
	break;			// exterior : big values produces artifacts on the image  



    }





  //if (d<0.0) {// interior
  // d(z1a) - d(z21) = -0.0804163521959989        
  //      d = - d;
  //      d = (db - d) /dd ; // normalize, see test_interior
  //d = d*d;
  //if (d>1.0) {printf("d int > 1.0\n");
  ///     }
  //      else {

  d = d / ((double) i);		// averaging not summation
  //d = d*me;} // exterior

  return d;




}





unsigned char
ComputeColor_DLD (complex double z, int FatouType)
{


  //double cabsz;
  int iColor;
  double d;

  if (FatouType == 1)
    {				// interior
      d = lagrangian (z, c, N, p);
      // modify gradient position

      //{d = d - (int)d;} // only fractional part
      d = d * d * mi;
      //if ( d< 1.0 ) d = 0.0;

    }				//  
  else
    {
      d = lagrangian (z, c, 10 * N, p);
    }

  iColor = (int) (d * 255) % 255;	// nMax or lower walues in denominator



  return (unsigned char) iColor;


}



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


  i = Give_i (ix, iy);		/* compute index of 1D array from indices of 2D array */
  z = GiveZ (ix, iy);
  iColor = A[i];		// read color = read the information about Fatou component type ( interior/exterior)
  if (iColor == iColorOfInterior)
    {
      FatouType = 1;
    }				// tru = interior
  else
    {
      FatouType = 0;
    }

  iColor = ComputeColor_DLD (z, FatouType);	// compute new color 
  A[i] = iColor;		// save new colr to the array         

  return 0;
}




// fill array 
// uses global var :  ...
// scanning complex plane 
int
DrawDLDImage (unsigned char A[])
{
  unsigned int ix, iy;		// pixel coordinate 

  printf ("compute DLD image \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)
    {
      printf (" %d from %d \r", iy, iyMax);	//info 
      for (ix = ixMin; ix <= ixMax; ++ix)
	DrawDLDPoint (A, ix, iy);	//  
    }

  return 0;
}



//=========================================

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

int
SaveArray2PGMFile (unsigned char A[], double k, char *comment)
{

  FILE *fp;
  const unsigned int MaxColorComponentValue = 255;	/* color component is coded from 0 to 255 ;  it is 8 bit color file */
  char name[100];		/* name of file */
  snprintf (name, sizeof name, "%.3f", k );	/*  */
  char *filename = strcat (name, ".pgm");
  char long_comment[200];
  sprintf (long_comment, "fc(z)=z^2+ c l c = (%f %+f );  %s", creal(c), cimag(c),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 %u requested\n", rSize,  iSize);}

  return 0;
}
















int PrintInfoAboutProgam()
{

  
  // display info messages
  printf ("Numerical approximation of Julia set for fc(z)= z^2 + c \n");
  //printf ("iPeriodParent = %d \n", iPeriodParent);
  //printf ("iPeriodOfChild  = %d \n", iPeriodChild);
  printf ("parameter c = ( %.16f ; %.16f ) \n", creal(c), cimag(c));
  
  printf ("Image Width = %f in world coordinate\n", ZxMax - ZxMin);
  printf ("PixelWidth = %f \n", PixelWidth);
  
  printf("for DEM/J \n");
  if ( distanceMax<0.0 || distanceMax > ER ) printf("bad distanceMax\n");
	printf("Max distance from exterior to the boundary =  distanceMax = %.16f = %f pixels\n",  distanceMax, BoundaryWidth); 
  
  // image corners in world coordinate
  // center and radius
  // center and zoom
  // GradientRepetition
  printf ("Maximal number of iterations = iterMax = %ld \n", iterMax);
  
  printf ("For LSM/J \n");
  printf ("Maximal number of iterations = iterMax_LSM = %ld \n", iterMax_LSM);
  printf ("Escape Radius = ER_LSM = %f \n", ER_LSM);
  
  
  printf ("ratio of image  = %f ; it should be 1.000 ...\n", ratio);
  //
  printf("gcc version: %d.%d.%d\n",__GNUC__,__GNUC_MINOR__,__GNUC_PATCHLEVEL__); // https://stackoverflow.com/questions/20389193/how-do-i-check-my-gcc-c-compiler-version-for-my-eclipse
  // OpenMP version is diplayed in the console 
  return 0;
}





int PrintInfoAboutPoint(complex double z){

	//unsigned int ix, iy;		// pixel coordinate
	// to do 
	
	double arg;
	unsigned char iColor;
	
	arg = Give_Arg( z, 2500); //   N in wiki
	iColor = ComputeColorOfSAC(z);
	
	printf ("parameter z = ( %.16f ; %.16f ) \n", creal(z), cimag(z));
	printf ("SAC/J : arg = %.16f ; iColor = %d  \n", arg, iColor);
	
	

	return z; 

}

// find such ER for LSM/J that level curves croses critical point and it's preimages
double GiveER(int i_Max){

	complex double z= 0.0; // criical point
	int i;
	 ; // critical point escapes very fast here. Higher valus gives infinity
	for (i=0; i< i_Max; ++i ){
		z=z*z +c; 
	 
	 }
	 
	 return cabs(z);
	
	
}





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;

}





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

int setup ()
{

  printf ("setup start\n");
  c = 0.35; // parabolic parameter   
  
  
  
  
	
  /* 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 ...
	
   
	
  
  //ER2 = ER * ER; // for numerical optimisation in iteration
  lnER = log(EscapeRadius); // ln(ER) 
  ER_LSM = GiveER(10); // find such ER for LSM/J that level curves croses critical point and it's preimages
  ER_DLD = GiveER(15);
  
  
   	
  /* create dynamic 1D arrays for colors ( shades of gray ) */
  data = malloc (iSize * sizeof (unsigned char));
  edge = malloc (iSize * sizeof (unsigned char));
  edge2 = malloc (iSize * sizeof (unsigned char));
  	
  if (data == NULL || edge == NULL || edge2 == NULL){
    fprintf (stderr, " Could not allocate memory");
    return 1;
  }

  
 	
  
  BoundaryWidth = 6.0*iWidth/2000.0  ; //  measured in pixels ( when iWidth = 2000) 
  distanceMax = BoundaryWidth*PixelWidth;
  
  
  
  printf (" end of setup \n");
	
  return 0;

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




int end(){


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

}

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

int main () {
  setup ();
  
  
   
  
  DrawImagerOfDEMJ(data);
  SaveArray2PGMFile (data, iWidth+0.1, "boundary using DEM/J");
  
  
  
  DrawImagerOfBD(data);
  SaveArray2PGMFile (data, iWidth+0.2, "BD/J");
  
  ComputeBoundaries(data, edge);
  SaveArray2PGMFile (edge, iWidth+0.3, "boundaries of BD/J");
  
  DrawImagerOMfBD(data);
  SaveArray2PGMFile (data, iWidth+0.4, "MBD/J");
  
  ComputeBoundaries(data, edge2);
  SaveArray2PGMFile (edge2, iWidth+0.5, "boundaries of MBD/J");
  
  
  
 
  DrawImagerOfLSM(data);
  SaveArray2PGMFile (data, iWidth+0.6, "LSM/J");
  
  
  ComputeBoundaries(data, edge);
  SaveArray2PGMFile (edge, iWidth+0.7, "boundaries of LSM/J");
  
  CopyBoundaries(edge, data);
  SaveArray2PGMFile (data, iWidth+0.8, "LSM + boundaries of LSM/J");
  
  CopyBoundaries(edge, edge2);
  SaveArray2PGMFile (edge2, iWidth+0.9, "boundaries of LSM/J and MBD");
  
  
  DrawImagerOfUnknown(data);
  SaveArray2PGMFile (data, iWidth+1.0, "Unknown : boundary and slow dynamics");
  
  
  DrawImagerOMfSAC(data);
  SaveArray2PGMFile (data, iWidth+1.1, "SAC/J + DEM/J");
  /*
  
  p = 0.17;
  int n;
  for (n=0; n< 100; ++n )
  {
  p += 0.001;
  */
  DrawDLDImage(data);
  DrawImagerOfDEMJ_boundary(data);
  SaveArray2PGMFile (data, iWidth+1.2, "DLD/J + boundary by DEM");
  
 //} 
  
  //PrintInfoAboutPoint(ZxMin+ZyMax*I);
  
  end();

  return 0;
}