/*




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



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

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



int IterMax = 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 iColorOfInterior1 = 99;
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
complex double zc1a = -0.47068779553447764874 + 0.0026098248687148155323*I; //period 1 attract
complex double zc1p = 0.24496023578261891251 + 0.0050657319837705428595*I; // period 1  attracting from parabolic

// periodic points = attractors
complex double zp1a =-0.33036439123272171026-0.035156244692189524137*I ; //period 1 attract
complex double zp1p = 0.024367377494104072722 +0.043208283893844831591*I ; // period 1 attracting from parabolic



/*

  F(z,C):=z(1+z)(1+z)(C-(2C+1)z+(3+4C)zz/4)


  Then for any C, F(z,C) has the following properties:

  1. z=0 is a fixed point with multiplier C.
  2. z=-1 is a critical point and F(-1,C)=0.
  3. z=1 is another critical point, and F(1,C)=-1, thus F(F(1,C),C)=0




  C*z^5+(3*z^5)/4+z^4/2-2*C*z^3-(5*z^3)/4-z^2+C*z

  C*z^5
  +(3*z^5)/4
  +z^4/2
  -2*C*z^3
  -(5*z^3)/4
  -z^2
  +C*z

  ==============
  (%o12) C*z^5+(3*z^5)/4+z^4/2-2*C*z^3-(5*z^3)/4-z^2+C*z

  (%i13) coeff(f,z,5);
  (%o13) C+3/4 = 0.058*%i+1.77


  (%i14) coeff(f,z,4);
  (%o14) 1/2

  (%i15) coeff(f,z,3);
  (%o15) (-2*C)-5/4 = (-0.116*%i)-3.29


  (%i16) coeff(f,z,2);
  (%o16) -1

  (%i17) coeff(f,z,1);
  (%o17) C =  1.02 + 0.05*I;


  (%i18) coeff(f,z,0);
  (%o18) 0
 
  =============================

  coefficients read from input file kawahira_sc_c3.txt
  degree 5 coefficient = ( +1.7700000000000000 +0.0580000000000000*i) 
  degree 4 coefficient = ( +0.5000000000000000 +0.0000000000000000*i) 
  degree 3 coefficient = ( -3.2900000000000000 -0.1160000000000000*i) 
  degree 2 coefficient = ( -1.0000000000000000 +0.0000000000000000*i) 
  degree 1 coefficient = ( +1.0200000000000000 +0.0500000000000000*i) 
  degree 0 coefficient = ( +0.0000000000000000 +0.0000000000000000*i) 

  Input polynomial p(z)=(1.7700000000000000178+0.058000000000000002942i)*z^5+(0.5+0i)*z^4+(-3.2900000000000000355-0.11600000000000000588i)*z^3+(-1+0i)*z^2+(1.0200000000000000178+0.050000000000000002776i)*z^1

  3 critical points found

  cp#0: -0.47068779553447764874,0.0026098248687148155323 . It's critical orbit is bounded and enters cycle #0 length=1 and it's stability = |multiplier|=0.66225 =attractive 
  internal angle = 0.97054734997537162045
  cycle = {
  -0.33036439123272171026,-0.035156244692189524137 ; }

  cp#1: 0.24496023578261891251,0.0050657319837705428595 . It's critical orbit is bounded and enters cycle #1 length=1 and it's stability = |multiplier|=0.98594 =attractive 
  internal angle = 0.99082387912800862217
  cycle = {
  0.024367377494104072722,0.043208283893844831591 ; }

  cp#2: -1.0000287546513304537,-0.00068489313584920526092 . It's critical orbit is bounded  and enters cycle #1



*/



// C=1.02+0.05i
complex double C = 1.02 + 0.05*I;


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




}



double cabs2(complex double z){

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


}




// =====================
int IsPointInsideTrap1(complex double  z){

	
	 
	
  if ( cabs2(z - zp1a) < AR2) {return 1;} // circle with prabolic point zp on it's boundary
  return 0; // outside



}



// =====================
int IsPointInsideTrap2(complex double  z){

	
  if (cabs2(z - zp1p) <AR2) {return 1;} // circle around periodic point
	
  return 0; // outside



}


complex double F(complex double z, complex double C){


  return (z*(1+z)*(1+z)*(C-(2*C+1)*z+(3+4*C)*z*z/4));


};






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



// f(z)=1+z−3z2−3.75z3+1.5z4+2.25z5
unsigned char
ComputeColor_Fatou (complex double z, int IterMax)
{



	
	
  double r2;


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


		

      z = F(z,C);		// complex iteration f(z)=z^6+A*z+c
      r2 =cabs2(z);
		
      if (r2 > ER2) // esaping = exterior
	{
	  uExterior += 1;
	  return iColorOfExterior;
	}			
	
      // solid color for each Fatou components
	
      if ( IsPointInsideTrap1(z)) {
	uInterior +=1;
	return iColorOfInterior1;
      } // 50 + (i % 114); }
	
      if (IsPointInsideTrap2(z)){
	uInterior +=1;
	return iColorOfInterior2;}

	

    }

  uUnknown += 1;
  return iColorOfUnknown;


}





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


  i = Give_i (ix, iy);		/* compute index of 1D array from indices of 2D array */
  z = GiveZ (ix, iy);
  iColor = ComputeColor_Fatou (z, IterMax);
  A[i] = iColor;		// interior

  return 0;
}




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

  fprintf (stdout, "compute Fatou image \n");
  // 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)
	DrawFatouPoint (A, ix, iy, IterMax);	//  
    }

  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(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 =1.0*iWidth/4000.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]= 255-A[i];}}}
	
	
  return 0;
	
}


// fill array 
// uses global var :  ...
// scanning complex plane 
int MarkAttractors (unsigned char A[])
{
  
	
	
	
  fprintf (stderr, "mark attractors \n");
  
  PlotBigPoint(zp1a, A); // period 114  cycle
  PlotBigPoint(zp1p, A);	// period 2 attracting cycle
    		 
      	

  return 0;
}


// =====================
int IsPointInsideTraps(unsigned int ix, unsigned int iy){

	
  complex double  z = GiveZ (ix, iy);
	
  if ( IsPointInsideTrap1(z)) {return 1;} // circle with prabolic point on it's boundary
	
  if (IsPointInsideTrap2(z)) {return 1;}
	
  return 0; // outside



}





int MarkTraps(unsigned char A[]){

  unsigned int ix, iy;		// pixel coordinate 
  unsigned int i;


  fprintf (stderr, "Mark traps \n");
  // 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){
	if (IsPointInsideTraps(ix, iy)) {
	  i= Give_i(ix,iy); /* index of _data array */
	  A[i]= 255-A[i]; // inverse color
	}}}
  return 0;
}






int PlotPoint(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]= 255-A[i]; // Mark point with inveres color
	
	
  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 */
  
 
  PlotBigPoint(z, A);
  
  /* forward orbit of critical point  */
  for (i=1;i<iMax ; ++i)
    {
      z  = F(z,C);
      if (cabs2(z - zp1p) > 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[], int a, int b,  int c, 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, "%d_%d_%d", a, b, c );	/*  */
  char *filename = strcat (name, ".pgm");
  char long_comment[200];
  sprintf (long_comment, "fc(z)=z*(1+z)*(1+z)*(C-(2C+1)*z+(3+4*C)*z*z/4)  %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,C) =  z*(1+z)*(1+z)*(C-(2C+1)*z+(3+4*C)*z*z/4) \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 = %.16f \n", PixelWidth);
  printf ("AR = %.16f = %f *PixelWidth\n", AR, AR / PixelWidth);


  printf("pixel counters\n");
  printf ("uUnknown = %llu\n", uUnknown);
  printf ("uExterior = %llu\n", uExterior);
  printf ("uInterior = %llu\n", uInterior);
  printf ("Sum of pixels  = %llu\n", uInterior+uExterior + uUnknown);
  printf ("all pixels of the array = iSize = %llu\n", iSize);


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

}


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

  ER = 2.0; // 
  ER2 = ER*ER;
  AR = 18.5*PixelWidth*iWidth/2000.0 ; // adjust first number 
  AR2 = AR * AR;
  //AR12 = AR/2.0;
  
  
  	
  
  



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


  DrawFatouImage (data, IterMax);	// first find Fatou
  SaveArray2PGMFile (data, iWidth, IterMax, 0, "Fatou, name = iWidth_IterMax_n");
  
  ComputeBoundaries(data,edge);
  SaveArray2PGMFile (edge, iWidth, IterMax, 1, "Boundaries of Fatou; name = iWidth_IterMax_n"); 
  
  CopyBoundaries(edge,data);
  SaveArray2PGMFile (data, iWidth, IterMax, 2, "Fatou with boundaries; name = iWidth_IterMax_n"); 
  
  //MarkAttractors(data);
  MarkTraps(data);
  SaveArray2PGMFile (data, iWidth, IterMax, 4, "Fatou with boundaries and traps; name = iWidth_IterMax_n"); 
  	
  DrawForwardOrbit(zc1p, 10000, data);
  SaveArray2PGMFile (data, iWidth, IterMax, 5, "Fatou with boundaries and traps, critical orbit; name = iWidth_IterMax_n");  
  	

  end ();

  return 0;
}