/*

  Adam Majewski
  adammaj1 aaattt o2 dot pl  // o like oxygen not 0 like zero 
  

  https://plus.google.com/116648956837292097606/posts/b6J6z2u8soL






  how to show sepals inside main box of parabolic chessboard ?
  - compute full orbit ( forward and backward of every point)
  - for each whole orbit ( not point) compute maximal distance from orbit to fixed point alfa
  - normalize distance ( dustance/ distance max ) so it will have value from 0 to 1.0
  _ use such normalized distance for coloring
  - then one can see orbits 
  ==========================================

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

  gcc b.c -lm -Wall -march=native 
  time ./a.out


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


  time ./a.out

  time ./a.out >a.txt

  ----------------------

*/

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

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

// https://mrob.com/pub/muency/child.html
int ChildPeriod = 2; 		// Period of secondary component joined by root point with the parent component
int ParentPeriod = 1;		// main cardioid of Mandelbrot set    



// 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 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.6;	//-0.1;
static const double ZyMax =  1.6;	//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
double complex a; // alfa fixed point


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

static double ER = 2.0;		// Escape Radius for bailout test 
static double ER2;

double radius; //=  1.0-cabs(1.0-csqrt(1.0-4.0*c)) ; //0.1; // half of distance between critical point and fixed point

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

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

static unsigned char iColorOfExterior = 250;
static unsigned char iColorOfInterior = 60;
unsigned char ColorStep; // (240- iColorOfInterior)/ChildPeriod
unsigned char iColorOfUnknown = 50;




int NoOfUnknownPoints = 0;


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





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

/* 
   c functions using complex type numbers
   computes c from  component  of Mandelbrot set */
complex double Give_c( int Period,  int p, int q , double InternalRadius )
{
  
  complex double w;  // point of reference plane  where image of the component is a unit disk 
 // alfa = ax +ay*i = (1-sqrt(d))/2 ; // result 
  double t; // InternalAngleInTurns
  
  t  = (double) p/q; 
  t = t * M_PI * 2.0; // from turns to radians
  
  w = InternalRadius*cexp(I*t); // map to the unit disk 
  
  switch ( Period ) // of component 
    {
    case 1: // main cardioid = only one period 1 component
      c = w/2 - w*w/4; // https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Mandelbrot_set/boundary#Solving_system_of_equation_for_period_1
      break;
    case 2: // only one period 2 component 
      c = (w-4)/4 ; // https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Mandelbrot_set/boundary#Solving_system_of_equation_for_period_2
      break;
      // period > 2 
    default: 
      printf("higher periods : to do, use newton method \n");
      printf("for each q = Period of the Child component  there are 2^(q-1) roots \n");
      c = 10000.0; // bad value 
       
      break; }
  return c;
}


// compute alfa fixed point
// https://en.wikipedia.org/wiki/Periodic_points_of_complex_quadratic_mappings#Period-1_points_(fixed_points)
complex double GiveAlfa(complex double c)
{
	// d=1-4c 
	// alfa = (1-sqrt(d))/2 
	return (1.0-csqrt(1.0 - 4.0*c))/2.0 ;

}



// angle in turns
// https://en.wikipedia.org/wiki/Turn_(geometry)
double GiveTurn( double complex z){
double t;

  t =  carg(z);
  t /= 2*M_PI; // now in turns
  if (t<0.0) t += 1.0; // map from (-1/2,1/2] to [0, 1) 
  return (t);
}

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


}



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


// bailout test
// escapes = abs(z)> ER
int Escapes(complex double z){

  if (cabs2(z)>ER2) return 1;
  return 0;
}

int IsInTarget(complex double z){
  // here target set is a circle inside immediate basin component containing  critical point 
  // with fixed  point on it's boundary
  // attracting petal
  complex double center = a+radius;
  

  if (cabs(z-center) <= radius) return 1;
  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, ER2)
  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;
}





// ***************************************************************************************************************************
// ************************** Interior : components of Immediate Basin of Attraction *****************************************
// ****************************************************************************************************************************

unsigned char ComputeColorOfImmediateBasin(complex double z){



  int nMax = iterMax;
	
  int n;

  for (n=0; n < nMax; n++){ //forward iteration

    if (Escapes(z)) return iColorOfExterior;
    if (IsInTarget(z)) return iColorOfInterior + (n % ChildPeriod)* ColorStep;	// immediate basin of attraction and it's preimages 
  			
   
    z = z*z +c ; /* forward iteration : complex quadratic polynomial */ 
  }
  printf("unknown point : z = %.16f; %.16f\n", creal(z), cimag(z));
  NoOfUnknownPoints +=1;
  return iColorOfUnknown;
}



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




// fill array 
// uses global var :  ...
// scanning complex plane 
int DrawImagerOfImmediateBasin (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, iColorOfUnknown)
  	for (iy = iyMin; iy <= iyMax; ++iy){
    		//printf (" %d from %d \r", iy, iyMax);	//info 
    		for (ix = ixMin; ix <= ixMax; ++ix)
      			DrawPointOfImmediateBasin(A, ix, iy);	//  
  }

  return 0;
}





// *******************************************************************************************************
//******************** Interior:  Level Sets of Attraction time ******************************************
// *******************************************************************************************************
unsigned char ComputeColorOfInteriorLevelSets(complex double z){

  int nMax = iterMax;
	
  int n;
  int p;
  int pMax = ChildPeriod; // 

  for (n=0; n < nMax; n++){ //forward iteration

    if (Escapes(z)) return iColorOfExterior;
    for (p=0; p < pMax; p++){ //forward iteration
      	if (IsInTarget(z)) return iColorOfInterior + (n % ChildPeriod)* ColorStep;	// immediate basin of attraction and it's preimages 
      	z = z*z +c ; /* forward iteration : complex quadratic polynomial */ 
    }
  }
  //
  printf("unknown point : z = %.16f; %.16f\n", creal(z), cimag(z));
  NoOfUnknownPoints +=1;
  return iColorOfUnknown;
}


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




// fill array 
// uses global var :  ...
// scanning complex plane 
int DrawImageOfInteriorLevelSets (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, iColorOfUnknown)
  for (iy = iyMin; iy <= iyMax; ++iy){
    //printf (" %d from %d \r", iy, iyMax);	//info 
    for (ix = ixMin; ix <= ixMax; ++ix)
      DrawPointOfInteriorLevelSets(A, ix, iy);	//  
  }

  return 0;
}



//******************** Interior:  chessboard  ******************************************
// *******************************************************************************************************
unsigned char ComputeColorOfChessboard (complex double z){

  int nMax = iterMax;
	
  int n;
  int p;
  int pMax = ChildPeriod; // 

  for (n=0; n < nMax; n++){ //forward iteration

    if (Escapes(z)) return iColorOfExterior;
    for (p=0; p < pMax; p++){ //forward iteration
      	if (IsInTarget(z)) 
      		{if (cimag(z)>0.0) 
      			return 20 ;	//  above critical orbit
      			else return 243;	//  below critical orbit 
      			}
      	z = z*z +c ; /* forward iteration : complex quadratic polynomial */ 
    }
  }
  //
  printf("unknown point : z = %.16f; %.16f\n", creal(z), cimag(z));
  NoOfUnknownPoints +=1;
  return iColorOfUnknown;
}


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




// fill array 
// uses global var :  ...
// scanning complex plane 
int DrawImageOfChessboard  (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, iColorOfUnknown)
  for (iy = iyMin; iy <= iyMax; ++iy){
    	//printf (" %d from %d \r", iy, iyMax);	//info 
    	for (ix = ixMin; ix <= ixMax; ++ix)
      		DrawPointOfChessboard(A, ix, iy);	//  
  }

  return 0;
}






// ************************************************************************************
//******************** Interior:  both Levels sets and chesboard ******************************************
// *******************************************************************************************************
unsigned char ComputeColorOfBoth (complex double z){

  int nMax = iterMax;
	
  int n;
  int p;
  int pMax = ChildPeriod; // 
  double angle; // in turns 
  unsigned char color;

  for (n=0; n < nMax; n++){ //forward iteration

    if (Escapes(z)) return iColorOfExterior;
    for (p=0; p < pMax; p++){ //
      	if (IsInTarget(z)) 
      		{ 
      		  
      		  angle = GiveTurn(z - a); // now in (0,1) range
      		  
      		  // !!!!!!
      		  //if (angle> 7.0/8.0 ) angle = (angle)*8.0;
      		  //if (angle<1.0/8.0)  
      		  angle = angle*4.0; // repeated gradient
      		  //printf("angle = %.16f\n", angle);
      		  color = angle* 255; // now in (0,255) range
      		  return color;
      			}
      	z = z*z +c ; /* forward iteration : complex quadratic polynomial */ 
    }
  }
  //
  printf("unknown point : z = %.16f; %.16f\n", creal(z), cimag(z));
  NoOfUnknownPoints +=1;
  return iColorOfUnknown;
}


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




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

  // radius /=10.0;	
  //printf("compute image \n");
  // for all pixels of image 
#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax, iColorOfUnknown)
  for (iy = iyMin; iy <= iyMax; ++iy){
    	//printf (" %d from %d \r", iy, iyMax);	//info 
    	for (ix = ixMin; ix <= ixMax; ++ix)
      		DrawPointOfBoth(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 =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 parabolic 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 ("is a root point between period %d and %d components \n", ChildPeriod, ParentPeriod);
  printf ("alfa fixed point z = ( %.16f ; %.16f ) \n", creal(a), cimag(a));
  printf ("Image Width = %f in world coordinate\n", ZxMax - ZxMin);
  printf ("PixelWidth = %f \n", PixelWidth);
  printf("radius of attracting circular petal = %.16f\n", radius);
  // 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("NoOfUnknownPoints  = %d NoOfAllPoints = %d so ratio unknown/all = %f \n", NoOfUnknownPoints, iSize, (double) NoOfUnknownPoints/ iSize);
  return 0;
}


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

int setup ()
{

  printf ("setup start\n");
  c = Give_c(ParentPeriod, 1, ChildPeriod, 1.0);
  a = GiveAlfa(c); // -0.5; //  alfa fixed point
  radius =  cabs(c)/4.8 ; // choose such value that level sets cross at z=0
	
  /* 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;
  }

  
	
  ColorStep =  (243 - iColorOfInterior)/(ChildPeriod-1);
  if (ColorStep <1) {printf("error from setup : ColorStep < 0 ; It should be greater\n"); return 1; } // check
  
  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 ();
  // ******************************** components O f immediate basin **********************************************************
  DrawImagerOfImmediateBasin(data);
  SaveArray2PGMFile (data, iHeight, "Interior : components of immediate basin of attraction (IBA) and it's preimages");
  //
  ComputeBoundaries(data,edge);
  SaveArray2PGMFile (edge, iHeight+1.0, "only boundary of components");
  //
  CopyBoundaries(edge,data);
  SaveArray2PGMFile (data, iHeight+2.0, "components with boundaries");
  	
   // ***************** attraction time ***************************************************
  DrawImageOfInteriorLevelSets (data);
  SaveArray2PGMFile (data, iHeight+3.0+radius, "Interior: level sets of attraction time to the parabolic fixed point");
  //
  ComputeBoundaries(data,edge);
  SaveArray2PGMFile (edge, iHeight+4.0+radius, "only boundaries of level sets");
  //
  CopyBoundaries(edge,data);
  SaveArray2PGMFile (data, iHeight+5.0+radius, "level sets with boundaries");
  	
  // ***************** parabolic chessboard ***************************************************
  DrawImageOfChessboard(data);
  SaveArray2PGMFile (data, iHeight+6.0+radius, "Interior : parabolic chessboard");
  //
   ComputeBoundaries(data,edge);
  SaveArray2PGMFile (edge, iHeight+7.0+radius, "only boundaries of parabolic chessboard");
  //
  CopyBoundaries(edge,data);
  SaveArray2PGMFile (data, iHeight+8.0+radius, "parabolic chessboard with boundaries");
  
   // *****************  Interior:  both Levels sets and chesboard ***************************************************
  DrawImageOfBoth(data);
  SaveArray2PGMFile (data, iHeight+9.0+radius, "Interior:  both Levels sets and chesboard ");
  
  
  
  
  end();

  return 0;
}