/*

  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 ========================
  
  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 >i.txt
  time ./a.out >e.txt
  
  =======================
  # gnuplot "i.plt"
set terminal svg enhanced background rgb 'white'
set xlabel "re(z)"
set ylabel "DLD"
set title "Relation between z and DLD in the interior of Julia set for c = -1"
set output "interior.svg"
plot "i.txt" with lines

  ----------------------
  
*/

#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 ; //= 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 *bin;

// 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.6;	//-0.05;
static const double ZxMax =  1.6;	//0.75;
static const double ZyMin = -1.0;	//-0.1;
static const double ZyMax =  1.0;	//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 = -0.75;		// parameter of function fc(z)=z^2 + c
double complex zf = -0.5;

double ER = 1e60;
double AR = 1e-10; //pek = p*10^k = p*pow(10.0, k) 



const int N = 1000; // fixed number : maximal number of iterations
double p  = 0.015625;  // 1/64 // 0.25; 
double m = 2.0; // density of curves


int iInterior = 0;
int iExterior = 0;


/* colors = shades of gray from 0 to 255 */
unsigned char iColorOfExterior = 250;
unsigned char iColorOfInterior = 200;
unsigned char iColorOfBoundary = 0;





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


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







// ***************************************************************************************************************************
// ************************** 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
	
	
	//if (cabs(z) < AR || cabs(z +1)< AR) return 5.0; // for z= 0.0 d = inf
	
	
	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 ) {
			iExterior +=1;
			d = -1.0; // (double)i/iMax; // escape time
			break; // exterior : big values produces NAN error in ppnorm computing 
			}
		// it not works ????
		if (cabs(z -zf) < AR ) 
			{ // interior
			  iInterior +=1;
				//d = -d;
				break; 
				
			}
			
		
	}
	 
	
	
	
	if (d<0.0) {// exterior = escape time
			//d = -d;
		
		
		
		}
		else { // interior = DLD
			d =  d/((double)i); // averaging 
			d = m* d;
			d = d - (int)d; // fractional part}
		}
		
	
	
	return d; 
	



}





unsigned char ComputeColorOfDLD(complex double z){

 	
  	//double cabsz;
  	int iColor;
  	double d;
  
	
  	d = lagrangian(z,c, N,p);
  	
  	if ( d<0.0) 
  		{iColor = iColorOfExterior;}
   		else {iColor = (int)(d*255.0)  % 255;} // interior
  
  
  return (unsigned char) iColor;


}



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




// fill array 
// uses global var :  ...
// scanning complex plane 
int DrawImagerOfDLD (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)
      			DrawPointOfDLD(A, ix, iy);	//  
  }

  return 0;
}





// test how values changes to tune color 
int test_interior(){

// choose 2 points such that color is changing the most
	complex double z = zf;
	complex double z2 = 0.3*I;
	int iMax = 20;
	complex double dz = (zf- z2)/iMax;
	printf("dz = %.16f ; %.16f\n", creal(dz), cimag(dz));
	int i;
	
	printf("z = %.16f ; %.16f\n", creal(z), cimag(z));
	printf("# z d\n"); // gnuplot
	for (i=0; i<iMax; ++i){
	
		double d = lagrangian(z, c, N, p);
		//int iColor = ComputeColorOfDLD(z);
		 
		// printf(" z = %.16f d = %.16f color = %d \n",creal(z), d, iColor);
		printf("%d %.16f %.16f %.16f\t %d\n",i ,  creal(z), cimag(z),d, (int)(d*255.0)  % 255); // gnuplot 
		z = z -  dz;
		}
		
	//		
	double d0 = lagrangian(zf, c, N, p);
	double db = lagrangian(z2, c, N, p);	
	double dd = d0 - db;
	printf("d0 - db  = %.16f - %.16f = %.16f\n",d0, db, dd);
	
		
	return 0;


}
 
 
 

// test how values changes to tune color 
int test_exterior(){

	complex double z;
	complex double z0 = zf;
	complex double z1 = 3.0;
	complex double dz = cabs(z1 - z0)/20;
	
	
	z = z0;
	printf("# z d color\n"); // gnuplot
	while (creal(z) < creal(z1)){
	
		double d = lagrangian(z, c, N, p);
		//int iColor = ComputeColorOfDLD(z);
		 
		// printf(" z = %.16f d = %.16f color = %d \n",creal(z), d, iColor);
		printf(" %.16f\t %.16f \t%d\n",creal(z), d, (int)(d*255.0)  % 255); // gnuplot 
		z += dz;
		}
		
	//		
	double d0 = lagrangian(z0, c, N, p);
	double d1 = lagrangian(z1, c, N, p);	
	double dd = d0 - d1;
	printf("d0 - d1  = %.16f - %.16f = %.16f\n",d0, d1, dd);
		
	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 =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 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("iInterior = %d \n", iInterior);
  printf("iExterior = %d \n", iExterior);
 
  // image corners in world coordinate
  // center and radius
  // center and zoom
  // GradientRepetition
  printf ("Maximal number of iterations = iterMax = %d \n", N);
  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");
   
  
  
  
  
	
  /* 2D array ranges */
  
  iHeight = 5*2000; 
  iWidth = 5*3200;
  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)) / ((double) iWidth / (double) iHeight);	// it should be 1.000 ...
	
   
	
  
   	
  /* create dynamic 1D arrays for colors ( shades of gray ) */
  data = malloc (iSize * sizeof (unsigned char));
  edge = malloc (iSize * sizeof (unsigned char));
  bin = malloc (iSize * sizeof (unsigned char));    
  	
  if (data == NULL  || edge == NULL || bin == NULL){
    fprintf (stderr, " Could not allocate memory");
    return 1;
  }

  
 	
  
  
  
  
  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(bin);
  
  
  PrintInfoAboutProgam();
  return 0;

}

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

int main () {
  setup ();
  
  
  
  DrawImagerOfDLD(data);
  SaveArray2PGMFile (data, iWidth+m+p, "DLD/J");

	ComputeBoundaries(data, edge);
	SaveArray2PGMFile (edge, iWidth+100+m+p, "boundaries of DLD/J");

  test_exterior();
  test_interior(); 
  
  end();

  return 0;
}