/*

  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 c.c -lm -Wall -march=native 
  time ./a.out


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


  time ./a.out

  time ./a.out >a.txt

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









  period 2 cycle :
  z = -0.500000000000000  +1.000000000000000 i
  z = -0.500000000000000  -1.000000000000000 i



*/

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

//#define iPeriodChild 0		// Period of secondary component joined by root point with the parent component
//int iPeriodParent = 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 = 4000;	//  
// The size of array has to be a positive constant integer 
static unsigned int iSize;	// = iWidth*iHeight; 

// memmory 1D array 
unsigned char *data;


// 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.15;	//-0.05;
static const double ZxMax = 1.15;	//0.75;
static const double ZyMin = -1.15;	//-0.1;
static const double ZyMax = 1.15;	//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 = 1000000;	//iHeight*100;

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


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

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

static unsigned char iColorOfExterior = 250;
//static unsigned char iColorOfInterior = 150;
unsigned char iColorOfUnknown = 50;




int NoOfUnknownPoints = 0;


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


}



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





int IsInTarget(complex double z){
  // here target set is a square with side r
  // with fixed  point on it's boundary

  double s = 0.25;
	
  double zx = creal(z);
  //double dx; // = fabs(creal(z) - r - 0.5);
  double dy; // = fabs(cimag(z));
	
  // fixed point is z=0.5
  // target set is any 
  if (zx > 0.50) return 0; //
  //dx = fabs(zx - r - 0.5);
  dy = fabs(cimag(z));
  if (	zx <= 0.5  &&
	zx > 0.0  &&
	dy< s) return 1;
	
  return 0;
	
	



}


double GiveDistance2ToAlfa(complex double z){

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

}


double GiveDistance2MaxLocal(complex double z){
	
  double d2 =  GiveDistance2ToAlfa(z); // d2 = d*d = d without sqrt
  double d2max = 0.0;
		
  while (d2>d2max) {
		
    d2max = d2; // save new local max
    z = csqrt(z - 0.25); // inverse iteration : complex quadratic polynomial
    d2 =  GiveDistance2ToAlfa(z);  
		
  }
	
	
  return d2max;
}


unsigned char GiveColorOfInterior(complex double z){
	
	
  double D2MaxLocal = GiveDistance2MaxLocal(z);
	 
  unsigned char iColor = 255.0*D2MaxLocal/D2MaxGlobal; 
	
	
  //printf(" z = %.16f ; %.16f \t D2MaxLocal = %.16f \tcolor = %.16f = %d \n ", creal(z), cimag(z), D2MaxLocal, 255.0*D2MaxLocal/D2MaxGlobal, iColor);
	
  return iColor;
}




unsigned char ComputeColorSepal(complex double z){



  int nMax = iterMax;
	
  int n;

  for (n=0; n < nMax; n++){

    if (cabs2(z)>ER2) return iColorOfExterior;
    if (IsInTarget(z)) return GiveColorOfInterior(z);
    z = z*z +c ; /* forward iteration : complex quadratic polynomial */
  
    
  }
  printf("unknown point ( probably periodic) : z = %.16f; %.16f\n", creal(z), cimag(z));
  NoOfUnknownPoints +=1;
  return iColorOfUnknown;
}













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




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

  return 0;
}


// save "A" array to pgm file 
int
SaveArray2PGMFile (unsigned char A[], double k)
{

  FILE *fp;
  const unsigned int MaxColorComponentValue = 255;	/* color component is coded from 0 to 255 ;  it is 8 bit color file */
  char name[30];		/* name of file */
  sprintf (name, "%.0f", k);	/*  */
  char *filename = strcat (name, ".pgm");
  char *comment = "# ";		/* comment should start with # */

  /* 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 A array to the file in one step */
  printf ("File %s saved. \n", filename);
  fclose (fp);

  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 = ( %f ; %f ) \n", creal(c), cimag(c));
  printf ("Maximal  Distance to alfa =  %.16f  \n", D2MaxGlobal);

  printf ("Image Width = %f \n", ZxMax - ZxMin);
  printf ("PixelWidth = %f \n", PixelWidth);
  // 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\n");




  c = 0.25;



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




  D2MaxGlobal = GiveDistance2ToAlfa(0.5*I); // manually chooosen


  // for numerical optimisation in iteration
  ER2 = ER * ER;




  /* create dynamic 1D arrays for colors ( shades of gray ) */
  data = malloc (iSize * sizeof (unsigned char));
  

  if (data == NULL )
    {
      fprintf (stderr, " Could not allocate memory");
      getchar ();
      return 1;
    }









  printf (" end of setup \n");

  return 0;

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




int end(){


  printf (" allways free memory  to avoid buffer overflow \n");

  free (data);
  	
  info ();
  return 0;

}

/* -----------------------------------------  main   -------------------------------------------------------------*/
int main () {
  setup ();
  // 
  DrawSepals(data);
  SaveArray2PGMFile (data, iHeight);
  //
  end();

  return 0;
}