/*

https://www.facebook.com/photo.php?fbid=111908386634835&set=a.110008616824812&type=3&theater

center of period 85 componnet. Adress 1->5 -> 85
4K antialiased version. Max iterations 5000, and it was reached!

Paul Bourke I rendered a 30,000 pixel version

  Adam Majewski
  adammaj1 aaattt o2 dot pl  // o like oxygen not 0 like zero 
  
  
  console program in c programing language 
  
  ==============================================
  
  
  Structure of a program or how to analyze the program 
  
  
  ============== Image X ========================
  
  DrawImageOfX -> DrawPointOfX -> ComputeColorOfX -> IsInTheYTrap
  
  first 2 functions are identical for every X
  check only last function =  ComputeColorOfX
  which computes color of one pixel !
  
  

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

  -------------------------------
  cd existing_folder
  git init
  git remote add origin git@gitlab.com:adammajewski/SepalsOfCauliflower.git
  git add .
  git commit
  git push -u origin master
  ---------------------------------
  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

  ----------------------
  
  a@zelman:~/c/julia/demj/85$ gcc d.c -lm -Wall -march=native -fopenmp
a@zelman:~/c/julia/demj/85$ time ./a.out
setup start
 end of setup 
File 30001.0000000000.pgm saved . Comment = DEM/J 
 allways free memory (deallocate )  to avoid memory leaks 
Numerical approximation of Julia set for fc(z)= z^2 + c 
parameter c = ( 0.3555340000000000 ; -0.3372920000000000 ) 
Image Width = 2.400000 in world coordinate
PixelWidth = 0.000080 
distanceMax = 0.0000800026667556
Maximal number of iterations = iterMax = 50000 
ratio of image  = 1.000000 ; it should be 1.000 ...
gcc version: 7.3.0

real	2021m48,597s
user	16100m25,607s
sys	1m52,759s

  

*/

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



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


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


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


static unsigned long int iterMax = 50000;	//iHeight*100;

static double ER = 1e60;		// EscapeRadius for bailout test  , minimal =  2.0




double distanceMax;


//double D2MaxGlobal;	//= 0.0497920256372717 ;
//double DistanceMaxGlobal2  ;

/* colors = shades of gray from 0 to 255 */

static unsigned char iColorOfExterior = 250;
static unsigned char iColorOfInterior = 255;

unsigned char iColorOfBoundary = 0;


double BoundaryWidth = 1.0; //  distanceMax = BoundaryWidth*PixelWidth; // here boundary is changing with resolution, maybe % of ImageWidth would be better ? 



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





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



// fast cabs
double cabs2(complex double z) {
  return (creal(z) * creal(z) + cimag(z) * cimag(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;
	
	


}




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

    if (cabs(z)> ER || 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;
  cabsz = cabs(z);
  distance = 2.0 * cabsz* log(cabsz)/ cabs(dz);
  if (distance <distanceMax) {//printf(" distance = %f \n", distance); 
  		return iColorOfBoundary;}
  
  
  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[])
{
  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;
}









// *******************************************************************************************
// ********************************** 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, "%.10f", k); /*  */
  char *filename =strncat(name,".pgm", 4);
  
  
  
  // save image to the pgm file 
  fp= fopen(filename,"wb"); // create new file,give it a name and open it in binary mode 
  fprintf(fp,"P5\n # %s\n %u %u\n %u\n", comment, iWidth, iHeight, MaxColorComponentValue);  // write header to the file
  fwrite(A,iSize,1,fp);  // write array with image data bytes to the file in one step 
  fclose(fp); 
  
  // info 
  printf("File %s saved ", filename);
  if (comment == NULL || strlen(comment) ==0)  
    printf("\n");
  else printf (". Comment = %s \n", comment); 

  return 0;
}




int info ()
{

  
  // 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);
  
  if ( distanceMax<0.0 || distanceMax > ER ) printf("bad distanceMax\n");
	printf("distanceMax = %.16f\n",  distanceMax); 
  
  // image corners in world coordinate
  // center and radius
  // center and zoom
  // GradientRepetition
  printf ("Maximal number of iterations = iterMax = %ld \n", iterMax);
  printf ("ratio of image  = %f ; it should be 1.000 ...\n", ratio);
  //
  printf("gcc version: %d.%d.%d\n",__GNUC__,__GNUC_MINOR__,__GNUC_PATCHLEVEL__); // https://stackoverflow.com/questions/20389193/how-do-i-check-my-gcc-c-compiler-version-for-my-eclipse
  // OpenMP version is diplayed in the console 
  return 0;
}


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

int setup ()
{

  printf ("setup start\n");
  c = 0.355534-I*0.337292; //
  
  
  
  
  
	
  /* 2D array ranges */
  
  iWidth = iHeight;
  iSize = iWidth * iHeight;	// size = number of points in array 
  // iy
  iyMax = iHeight - 1;		// Indexes of array starts from 0 not 1 so the highest elements of an array is = array_name[size-1].
  //ix

  ixMax = iWidth - 1;

  /* 1D array ranges */
  // i1Dsize = i2Dsize; // 1D array with the same size as 2D array
  iMax = iSize - 1;		// Indexes of array starts from 0 not 1 so the highest elements of an array is = array_name[size-1].

  /* Pixel sizes */
  PixelWidth = (ZxMax - ZxMin) / ixMax;	//  ixMax = (iWidth-1)  step between pixels in world coordinate 
  PixelHeight = (ZyMax - ZyMin) / iyMax;
  ratio = ((ZxMax - ZxMin) / (ZyMax - ZyMin)) / ((float) iWidth / (float) iHeight);	// it should be 1.000 ...
	
   
	
  
  //ER2 = ER * ER; // for numerical optimisation in iteration
  
  
   	
  /* 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;
  }

  
 	
  
  
  distanceMax = BoundaryWidth*PixelWidth; // here boundary is changing with resolution, maybe % of ImageWidth would be better ? 
  
  
  
  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);
  info ();
  return 0;

}

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

int main () {
  setup ();
  // ******************************** DEM/J **********************************************************
  DrawImagerOfDEMJ(data);
  SaveArray2PGMFile (data, iWidth+BoundaryWidth, "DEM/J");
  
  
  end();

  return 0;
}

/*

https://www.facebook.com/photo.php?fbid=111908386634835&set=a.110008616824812&type=3&theater

center of period 85 componnet. Adress 1->5 -> 85
4K antialiased version. Max iterations 5000, and it was reached!

Paul Bourke I rendered a 30,000 pixel version

  Adam Majewski
  adammaj1 aaattt o2 dot pl  // o like oxygen not 0 like zero 
  
  
  console program in c programing language 
  
  ==============================================
  
  
  Structure of a program or how to analyze the program 
  
  
  ============== Image X ========================
  
  DrawImageOfX -> DrawPointOfX -> ComputeColorOfX -> IsInTheYTrap
  
  first 2 functions are identical for every X
  check only last function =  ComputeColorOfX
  which computes color of one pixel !
  
  
Vec3f(color[2], color[1], color[0]);
   
  ==========================================

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

  ----------------------
 ls -l *

=====


*/

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



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

int iColorSize = 3 ; // RGB = 3*(unsigned char)


// memmory 1D arrays 
double *dData; // double for dem ( double = 8 bytes per pixel)
unsigned char * ucColor; // 


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


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


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


static unsigned long int iterMax = 50000;	//iHeight*100;

static double ER = 160;		// EscapeRadius for bailout test  , minimal =  2.0




//double distanceMax;
double distanceMaxImage = 0.0;
double BoundaryWidth = 0.0;

//double D2MaxGlobal;	//= 0.0497920256372717 ;
//double DistanceMaxGlobal2  ;


// sum of points = iSize;
int iExterior = 0;
int iInterior = 0;
int iBoundary = 0;
int iAll =0;
double PixelWidthRatio = 1.0; //  BoundaryWidth= PixelWidthRatio*PixelWidth; // here boundary is changing with resolution, maybe % of ImageWidth would be better ? 



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





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



// fast cabs
double cabs2(complex double z) {
  return (creal(z) * creal(z) + cimag(z) * cimag(z));
}



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

// uses global constants
double GiveZy (int iy) {
  //return (ZyMin + iy * PixelHeight);
  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;
}





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

double ComputeDistance(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 = 0.0;
  double cabsz;
	
  int n;

  // compute distance from point z to the boundary of Julia set
  for (n=0; n < nMax; n++){ //forward iteration

    if (cabs(z)> ER || cabs(dz)> 1e60) break; // big values
    
  			
    dz = 2.0*z * dz; 
    z = z*z +c ; /* forward iteration : complex quadratic polynomial */ 
  }
  
  if (n==nMax) 
  	{distance = -1.0;} // 
  	else {
  		cabsz = cabs(z);
  		distance = 2.0 * cabsz* log(cabsz)/ cabs(dz);
  		}
  	
  		
   
  return distance;
}



// compute and save  raster point (ix,iy) data
int ComputePointData (double A[], int ix, int iy)
{
  int i;			/* index of 1D array */
  double distance;
  complex double z;


  i = Give_i (ix, iy);		/* compute index of 1D array from indices of 2D array */
  z = GiveZ(ix,iy);
  distance = ComputeDistance(z);
  
  if (distance >distanceMaxImage ) distanceMaxImage = distance; // compute also max distance of the image
  A[i] = distance ;		// 
  
  return 0;
}




// fill array 
// uses global var :  ...
// scanning complex plane 
int FillDataArray (double A[])
{
  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)
      			ComputePointData(A, ix, iy);	//  
  }

  return 0;
}






/* 
 normalize values in  array D 
 normalize = change range to [0,1.0] 
 using transformation : d = d/dMax 
 to do it one muust know dMax = max (D)
 uses global var :  ...
 
*/ 
int NormalizeDataArray (double D[])
{
  int i=0;		// array index 

  	//printf("compute image \n");
 	// for all pixels of image 
	#pragma omp parallel for schedule(dynamic) private(i) shared(iAll, iBoundary, iExterior, iInterior,  D, iSize)
  	for (i = 0; i < iSize; ++i){
    		D[i] = D[i]/distanceMaxImage;
  }
  
  //printf(" i = %d\n", i);
	
  return 0;
}







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

	 FILE * fp;
  
  	char name [100]; /* name of file */
  	snprintf(name, sizeof name, "%.0f", k); /*  */
  	char *filename =strncat(name,".dat", 4);
  
  
  
  	// save image to the pgm file 
  	fp= fopen(filename,"wb"); // create new file,give it a name and open it in binary mode 
  	fwrite(A,iSize,1,fp);  // write array with image data bytes to the file in one step 
  	fclose(fp); 
  
  	// info 
  	printf("File %s saved ", filename);
  	if (comment == NULL || strlen(comment) ==0)  
    	printf("\n");
  	else printf (". Comment = %s \n", comment); 

  	return 0;

}






// *******************************************************************************************
// **********************************  pfm file ****************************
// *********************************************************************************************



/*

input :
* int i
* array D of double numbers ( distance)

output : array of rgb colors 

*/
int ComputePointColor(const int i, const double D[], unsigned char  C[] ){


	
	int iC = i*iColorSize; // compute index of F array
	// color channels from 0 to 255  
	unsigned char R;
	unsigned char G;
	unsigned char B;
	
	double Position = D[i]; // read distance form // distance/distanceMax;
	
	
	if (Position<0.0)// interior = white
		{
		  R = 255;
		  G = 255;
		  B = 255; 
		  iInterior +=1;
		} 
		else {
			if ( Position > BoundaryWidth) // exterior = color gradient
				{	Position = sqrt(Position); // Position*Position*Position;
					R = 0;  // red
					G = 255.0*Position; //x^3  = non linear gradient
					B = 0; // blue 
					iExterior +=1;
				} 
				else {	// boundary = black 
					 
					R = 0;
					G = 0;
					B = 0;
					iBoundary +=1;
					} // boundary
		} 
		
	//printf("R = %u \t G = %u \t B = %u \n", R, G, B); // info 
	// save color to the array
	C[iC] 	= R;
	C[iC+1]	= G;
	C[iC+2] = B;
	
	
	iAll +=1;
	
	
	
	
	
	
	
	return 0;

};







// fill array f using data from d array
// uses global var :  ...
int FillColorArray (double D[], unsigned char C[])
{
  int i=0;		// array index 

  	//printf("compute image \n");
 	// for all pixels of image 
	#pragma omp parallel for schedule(dynamic) private(i) shared(iAll, iBoundary, iExterior, iInterior,  D, C, iSize)
  	for (i = 0; i < iSize; ++i){
    		//printf (" %d from %d \n", i, iSize);	//info 
    		ComputePointColor(i, D, C);	//  
  }
  
  //printf(" i = %d\n", i);
	
  return 0;
}





int Z_IsInFirstQuadrant(const int ix, const int iy){

	double zx = GiveZx(ix);
	double zy = GiveZy(iy);
	int r = 0;
	
	if (zx >0.0 && zy > 0.0) { r = 1;}
	
	return r;
	
}




int InverseColor(const int i, unsigned char C[]){

	C[i]   = 255 - C[i];
	C[i+1] = 255 - C[i+1];
	C[i+2] = 255 - C[i+2];
	
	return 0;
}


// scanning complex plane 
int CheckArrayOrientation (unsigned char C[])
{
  int ix, iy;		// pixel coordinate 
  int i; 

  	//printf("compute image \n");
 	// for all pixels of image 
	//#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(F, ixMax , iyMax)
  	for (iy = iyMin; iy <= iyMax; ++iy){ 
    		for (ix = ixMin; ix <= ixMax; ++ix){
      			
      			if ( Z_IsInFirstQuadrant(ix, iy)) {
      			
      				i = Give_i(ix,iy);
      				InverseColor(i*iColorSize, C);
      			
      				}
      			
      			
      			}
      				//  
  }

  return 0;
}












int Save_PPM( unsigned char Cl[], double k, char* comment )
{
  
  FILE * fp;
  
  char name [100]; /* name of file */
  snprintf(name, sizeof name, "%.0f", k); /*  */
  char *filename =strncat(name,".ppm", 4);
  
  
  
  // save image to the pgm file 
  fp= fopen(filename,"wb"); // create new file,give it a name and open it in binary mode 
  if (fp ) {
  	fprintf(fp,"P6\n%d %d\n255\n",  iWidth, iHeight);  // write header to the file
  	fwrite(Cl, iSizeC, 1,  fp);  // write array with image data bytes to the file in one step 
  	fclose(fp); 
  
  	// info 
  	printf("File %s saved ", filename);
  	if (comment == NULL || strlen(comment) ==0)  
    		printf("\n");
  		else printf ("Comment = %s \n", comment); 
	}
	else printf( "ERROR saving `%s'\n", filename);
  return 0;
}




int info ()
{

  	int iSum =  iExterior+iBoundary+iInterior;
  	
	// 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);
  
  	//if ( distanceMax<0.0 || distanceMax > ER ) printf("bad distanceMax\n");
	//printf("distanceMax = %.16f\n",  distanceMax); 
	
  	printf("distanceMaxImage = %.16f\n",  distanceMaxImage); 
  	printf("BoundaryWidth = %f = %f *distanceMaxImage = %f * PixelWidth\n" , BoundaryWidth , BoundaryWidth/((double) distanceMaxImage), PixelWidthRatio);	
	
		
  	printf("iSize = number of points in the  dData array = %d\n", iSize );
  	//
  	printf("iSizeC = number of points in the ucColor array = %d\n", iSizeC);
  	printf(" iExterior = %d \t iBoundary = %d \t iInterior = %d \t iSum = %d\n", iExterior, iBoundary,iInterior, iSum);
  	printf("iAll = %d \n", iAll);
  	if (iSum != iSize) { printf("Error : iSum != iSize \n");}
  
  	
  	printf("estimated ppm file size = %d B\n ", iSizeC*((int) sizeof(unsigned char)) +16); // header ( 16 bytes ) +  series of three 4-byte IEEE 754 single precision floating point numbers for each pixel
  
  // image corners in world coordinate
  // center and radius
  // center and zoom
  // GradientRepetition
  	printf ("Maximal number of iterations = iterMax = %ld \n", iterMax);
  	printf ("ratio of image  = %f ; it should be 1.000 ...\n", ratio);
  
  
  //
  	printf("gcc version: %d.%d.%d\n",__GNUC__,__GNUC_MINOR__,__GNUC_PATCHLEVEL__); // https://stackoverflow.com/questions/20389193/how-do-i-check-my-gcc-c-compiler-version-for-my-eclipse
  	// OpenMP version is diplayed in the console 
  	return 0;
}


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

int setup ()
{

  printf ("setup start\n");
  c = 0.355534-I*0.337292; //
  
  
  
  
  
	
  /* 2D array ranges */
  
  iWidth = iHeight;// square image
  iSize = iWidth * iHeight;	// size = number of points in the dData array 
  iSizeC = iSize * iColorSize ; // size of the fRGB 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].

  /* Pixel sizes */
  PixelWidth = (ZxMax - ZxMin) / ixMax;	//  ixMax = (iWidth-1)  step between pixels in world coordinate 
  PixelHeight = (ZyMax - ZyMin) / iyMax;
  ratio = ((ZxMax - ZxMin) / (ZyMax - ZyMin)) / ((float) iWidth / (float) iHeight);	// it should be 1.000 ...
	
   
	
  
  //ER2 = ER * ER; // for numerical optimisation in iteration
  
  
   	
  /* create dynamic 1D arrays for data  ( dData) and  colors ( ucColor ) */
  dData = malloc (iSize * sizeof (double));
  ucColor = malloc (iSizeC * sizeof (unsigned char));
  	
  if (dData == NULL || ucColor == NULL){
    printf (" Could not allocate memory");
    return 1;
  }

  
 	
  
  
  BoundaryWidth =  PixelWidthRatio*PixelWidth; // here boundary is changing with resolution, maybe % of ImageWidth would be better ? 
  
  
  
  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 (dData);
  free(ucColor);
  info ();
  return 0;

}

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

int main () {
  setup ();
  
  //  DEM/J 
  FillDataArray(dData); // compute distasnce of each pixel and also find max distance of the image
  NormalizeDataArray(dData); // convert distance  normalized distance = distance/distanceMaxImage
  SaveDataFile(dData, iHeight+PixelWidthRatio, "") ; 
  
  // 
  FillColorArray(dData, ucColor);
  Save_PPM(ucColor, iHeight+PixelWidthRatio, "");
  //
  // CheckArrayOrientation(fRGB);
  // Save_fRGB_Array_2_PFM(fRGB, iHeight, "");
  
  
  
  end();

  return 0;
}