/*
 
  c++ console program
  It can be compiled and run under Linux, windows, Mac 
  It needs gcc
 
  --------------------------------------
  draws critical orbit for f(z)=z*z+c 
  
 
 
 
  ------------------------------------------
  one can change :
  
  - n 
  - iSide ( width of image = iXmax = (iSide) 
  - NrOfCrOrbitPoints = ; // check rang of type for big numbers : iMax, i, NrOfCrOrbitPoints


 %lld and %llu for print long long int
  -----------------------------------------
  1.pgm file code is  based on the code of Claudio Rocchini
  http://en.wikipedia.org/wiki/Image:Color_complex_plot.jpg
  create 8 bit color graphic file ,  portable gray map file = pgm 
  see http://en.wikipedia.org/wiki/Portable_pixmap
  to see the file use external application ( graphic viewer)
  I think that creating graphic can't be simpler
  ---------------------------
  2. first it creates data array which is used to store color values of pixels,
  fills tha array with data and after that writes the data from array to pgm file.
  It alows free ( non sequential) acces to "pixels"
  -------------------------------------------
  Here are 4 items : 
  1. complex plane Z with points z = zx+zy*i ( dynamic plane and world coordinate )
  2. virtual 2D array indexed by iX and iYmax ( integer or screen coordinate )
  3. memory 1D array data  indexed by i =f(iX,iY)
  4. pgm file 
 
 
 
  Adam Majewski   fraktal.republika.pl 
 
 
  to compile : 
  gcc d.c -lm -Wall -march=native
  to run ( Linux console) :
  time ./a.out



  convert -version
 
  convert data0.pgm -convolve "-1,-1,-1,-1,8,-1,-1,-1,-1"  -threshold 5% -negate data0.png
  convert data0.pgm -convolve "0,1,0,1,1,1,0,1,0" -threshold 5% -negate data0c.png
convert data0.pgm -convolve "-0.125,0,0.125,  -0.25,0,0.25,  -0.125,0,0.125" -threshold 5% -negate data0g.png

  convert data0.pgm -edge 3 -negate data0e.png
  convert data0.pgm -edge 3 -negate data0e.png
http://www.imagemagick.org/Usage/transform/#vision
"As you can see, the edge is added only to areas with a color gradient that is more than 50% white! I don't know if this is a bug or intentional, but it means that the edge in the above is located almost completely in the white parts of the original mask image. This fact can be extremely important when making use of the results of the "-edge" operator. For example if you are edge detecting an image containing an black outline, the "-edge" operator will 'twin' the black lines, producing a weird result." 

  convert data0.pgm -negate -edge 3 data0n.png
  convert data0n.png -edge 3 -negate data0e.png


http://unix.stackexchange.com/questions/299218/imagemagick-how-to-thicken-lines

convert 4.pgm -negate 4n.pgm
convert 4n.pgm -morphology Dilate Octagon 4nf.pgm
convert 4n.pgm -morphology Thicken '3x1+2+0:1,0,0' 4nfb.pgm
convert 4n.pgm -morphology Thicken ConvexHull 4nfc.pgm
convert 4nf.pgm -negate 4thick.pgm


--------------------
check if curve is closed : 
if x(1) == x(end) && y(1) == y(end)
  % It's closed
else
  % It's not closed
end


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

http://www.scholarpedia.org/article/File:InfoldingSiegelDisk.gif
1,000 × 1,000 pixels, file size: 91 KB, MIME type: image/gif, looped, 9 frames, 12s)
The Siegel disks have been translated in the plane so that the critical point is always at the same point on the screen (near the top). 


for n :
 0,1   I have used  1.0e5 = pow(10,5) points 
 n= 2                    1.0e6 = pow(10,6)
 n = 3                   1.0e7 = pow(10,7)
 n = 4                   1.0e8 = pow(10,8) // good result
 n= 5                    1.0e9 = pow(10,9) // not good 
 
For n=5 I have to try pow(10,12).
 
Do you use such high number of iterations or other method ?

I think in this particular case, I iterated a lot. However, in some other pictures, I used an acceleration method, an approximation of a high iterate of f.



--------------------------------------------------------
n="" "0" "" ; t= "" "0.2956859994078892" "" ; c = "" "0.6153124581224951*%i+0.06835556662164869" "

"n="" "1" "" ; t= "" "0.2875617458610296" "" ; c = "" "0.599810068302661*%i+0.1057522049785167" "

"n="" "2" "" ; t= "" "0.285916253540726" "" ; c = "" "0.5963646240901801*%i+0.1130872227062027" "

"n="" "3" "" ; t= "" "0.2857346725405881" "" ; c = "" "0.5959783359361234*%i+0.1138915132131216" "

"n="" "4" "" ; t= "" "0.2857163263170416" "" ; c = "" "0.5959392401496606*%i+0.1139727184036316" "

"n="" "5" "" ; t= "" "0.2857144897937824" "" ; c = "" "0.5959353258460209*%i+0.1139808467588366" "

"n="" "6" "" ; t= "" "0.2857143061224276" "" ; c = "" "0.5959349343683512*%i+0.1139816596727942" "

"n="" "7" "" ; t= "" "0.2857142877551018" "" ; c = "" "0.5959348952201112*%i+0.1139817409649742" "

"n="" "8" "" ; t= "" "0.2857142859183673" "" ; c = "" "0.5959348913052823*%i+0.1139817490942002" "

"n="" "9" "" ; t= "" "0.2857142857346939" "" ; c = "" "0.5959348909137995*%i+0.1139817499071228" "
"n="" "10" "" ; t= "" "0.2857142857163265" "" ; c = "" "0.5959348908746511*%i+0.1139817499884153" "






-------------------------------------
cd existing_folder
git init
git remote add origin git@gitlab.com:adammajewski/InfoldingSiegelDisk_in_c_2over7.git
git add d.c
git commit -m ""
git push -u origin master
----------------------------------------------------





*/

#define BOUNDS_CHECKS
// uncomment next line for tiny speedup but less safety (possible invalid memory writes)
//#undef BOUNDS_CHECKS

#ifdef __cplusplus

#include <complex>

#ifdef USE_QD_REAL
#define QD_INLINE
#include <qd/qd_real.h>
typedef qd_real Real;
#define get_double(z) ((z).x[0])
#define pi (Real("3.141592653589793238462643383279502884197169399375105820974944592307816406286208"))
#else
#ifdef USE_DD_REAL
#define QD_INLINE
#include <qd/dd_real.h>
typedef dd_real Real;
#define get_double(z) ((z).x[0])
#define pi (Real("3.141592653589793238462643383279502884197169399375105820974944592307816406286208"))
#else
typedef double Real;
#define get_double(z) (z)
#define pi (3.141592653589793238462643383279502884197169399375105820974944592307816406286208)
static inline Real sqr(Real x) { return x * x; }
static inline Real mul_pwr2(Real x, double y) { return x * y; }
#endif
#endif

typedef std::complex<Real> Cplx;

#define creal(c) (std::real(c))
#define cimag(c) (std::imag(c))
#define I (Cplx(0,1))

#else

#include <complex.h>

typedef double Real;
typedef double _Complex Cplx;

#endif


#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>


#define NMAX 13
#define L  (NMAX +1)

#define unlikely(x) (__builtin_expect(!!(x),0))
 
/* iXmax/iYmax =  */
const int iSide = 1000;
int iXmax ; /* width of image in pixels = (15*iSide); */
int iYmax ;
int iLength ;

/* dynamic 1D arrays for colors ( shades of gray )  */
  
  unsigned char *data;



/* */
 double ZxMin = -0.6;
 double ZxMax = 0.4;
double ZyMin  = -0.15;
 double ZyMax = 0.85;
/* (ZxMax-ZxMin)/(ZyMax-ZyMin)= iXmax/iYmax  */
 
 
double PixelWidth ;
double PixelHeight ;
double invPixelWidth ;
double invPixelHeight ;
 
 unsigned int period=1;
 // unsigned int m;
 // check rang of type for big numbers : iMax, i, NrOfCrOrbitPoints
  unsigned long long int NrOfCrOrbitPoints ;
 
 
 
 Real radius = 1.0;
 Real tt(int n)
{
  // t(n) = [0;3, 2 , 10^n, golden_ratio]
  Real phi = (sqrt(Real(5)) + 1) / 2;
  Real tenn = pow(Real(10), n);
  return 0 +1/( 3 +1/( 2 +1/( tenn +1/( phi ))));
}

/* fc(z) = z*z + c */
/*   */
Cplx C; 
Real Cx; /* C = Cx + Cy*i   */
Real Cy ;
 
 
 
/* colors */
const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
const int iExterior = 245; /* exterior of Julia set */
const int iJulia = 0; /* border , boundary*/
const int iInterior = 230;
 
 
 
/* ----------------------- functions ---------------------------------------- */


/* find c in component of Mandelbrot set 
 
   uses code by Wolf Jung from program Mandel
   see function mndlbrot::bifurcate from mandelbrot.cpp
   http://www.mndynamics.com/indexp.html

*/
Cplx GiveC(Real InternalAngleInTurns, Real InternalRadius, unsigned int Period)
{
  //0 <= InternalRay<= 1
  //0 <= InternalAngleInTurns <=1
  Real t = InternalAngleInTurns *2*pi; // from turns to radians
  Real R2 = InternalRadius * InternalRadius;
  //Real Cx, Cy; /* C = Cx+Cy*i */
  switch ( Period ) // of component 
    {
    case 1: // main cardioid
      Cx = (cos(t)*InternalRadius)/2-(cos(2*t)*R2)/4; 
      Cy = (sin(t)*InternalRadius)/2-(sin(2*t)*R2)/4; 
      break;
    case 2: // only one component 
      Cx = InternalRadius * 0.25*cos(t) - 1.0;
      Cy = InternalRadius * 0.25*sin(t); 
      break;
      // for each iPeriodChild  there are 2^(iPeriodChild-1) roots. 
    default: // higher periods : to do, use newton method 
      Cx = 0.0;
      Cy = 0.0; 
      break; }

  return Cx + Cy*I;
}




 
 
inline unsigned int f(unsigned int iX, unsigned int iY)
/* 
   gives position of point (iX,iY) in 1D array  ; uses also global variables 
   it does not check if index is good  so memory error is possible 
*/
{return (iX + iY*iXmax );}
 
 
 
inline int DrawPoint( double Zx, double Zy, unsigned char data[])
{
  unsigned int iX,iY; /* indices of 2D virtual array (image) = integer coordinate */
  unsigned int index; /* index of 1D array  */ 
 
#ifdef BOUNDS_CHECKS
  if (unlikely(Zx < ZxMin || ZxMax  < Zx || Zy < ZyMin || ZyMax < Zy)) {  printf("   point z = %f , %f out of bounds  \n", Zx, Zy); return -1; } 
#endif
  
  iX = (int)((Zx-ZxMin)*invPixelWidth);
  iY = (int)((ZyMax-Zy)*invPixelHeight); // reverse Y axis
  index = iX + iY*iXmax;//f(iX,iY);
  
  
  data[index] = iJulia;  /* draw */
  return 0;
}
 


int TestCriticalOrbit(unsigned long long int iMax )
{
  unsigned long long int i; /* nr of point of critical orbit */
  Real Zx,Zy, tmp;
   
 
  /* critical point z = 0 */
  Zx = 0.0;
  Zy = 0.0;
  //
  ZxMin = 0.0;
  ZxMax = 0.0;
  ZyMin = 0.0;
  ZyMax = 0.0;
  
  /* forward orbit of critical point  */
  for (i=1;i<=iMax ;i++)
    {
      /* f(z) = z*z+c */
      tmp = mul_pwr2(Zx*Zy, 2) + Cy;
      Zx = sqr(Zx) - sqr(Zy) + Cx;
      Zy = tmp;
       
 
      
     // if (Zx2+Zy2>4) { printf("   point z = %f , %f escapes \n",Zx, Zy); break;}
     if (Zx>ZxMax) ZxMax=get_double(Zx);
     if (Zx<ZxMin) ZxMin=get_double(Zx);
     if (Zy>ZyMax) ZyMax=get_double(Zy); 
     if (Zy<ZyMin) ZyMin=get_double(Zy);         

    }

   printf(" ZxMin = %.16f  ZxMax = %.16f \n", ZxMin, ZxMax);       
   printf(" ZyMin = %.16f  ZyMax = %.16f \n", ZyMin, ZyMax);  
  return 0;
 
}



 
 
/*

 check rang of type for big numbers : iMax, i, NrOfCrOrbitPoints
x := x^2 - y^2 + cx
y := 2 x y + cy
*/


int DrawCriticalOrbit(unsigned long long int iMax,  unsigned char A[] )
{
  unsigned long long int iMax100 = iMax / 100;
  unsigned long long int i; /* nr of point of critical orbit */
  Real Zx,Zy, tmp;
  int IsGood;  
 
  /* critical point z = 0 */
  Zx = 0.0;
  Zy = 0.0;
  DrawPoint(get_double(Zx),get_double(Zy),A);
  
  /* forward orbit of critical point  */
  for (int j = 0; j < 100; ++j)
  {
    printf("%4d \r", j);
    fflush(stdout);
    for (i=1;i<=iMax100 ;i++)
    {
      /* f(z) = z*z+c */
      tmp = mul_pwr2(Zx*Zy, 2) + Cy;
      Zx = sqr(Zx) - sqr(Zy) + Cx;
      Zy = tmp;
       
 
      
     // if (Zx2+Zy2>4) { printf("   point z = %f , %f escapes \n",Zx, Zy); break;}
      
     IsGood = DrawPoint(get_double(Zx),get_double(Zy),A); /* draws critical orbit */
#ifdef BOUNDS_CHECKS
     if (unlikely(IsGood<0)) { printf ("i = %llu \n", i); return 1;}
#endif
    }
  }

    
   
  return 0;
 
}


/* 
close the curve by filing gaps by streight lines
curve is the simple closed 2d plane curve 

*/
int CloseTheCurve( unsigned char A[] ){
(void) A;
return 0;
}


int ClearArray( unsigned char A[] )
{
  int index; /* index of 1D array  */
  for(index=0;index<iLength-1;++index) 
                A[index]=iExterior;
  return 0;
}


// Check Orientation of image : first quadrant in upper right position
// uses global var :  ...
int CheckOrientation(unsigned char A[] )
{
  int ix, iy; // pixel coordinate 
  Real Zx, Zy; //  Z= Zx+ZY*i;
  int i; /* index of 1D array */
  for(iy=0;iy<=iYmax;++iy) 
    {
      Zy =ZyMax - iy*PixelHeight;
      for(ix=0;ix<=iXmax;++ix) 
	{

	  // from screen to world coordinate 
	  Zx = ZxMin + ix*PixelWidth ;
	  i = f(ix, iy); /* compute index of 1D array from indices of 2D array */
	  if (Zx>0 && Zy>0) A[i] = 255-A[i];   // check the orientation of Z-plane by marking first quadrant */

	}
    }
   
  return 0;
}





int SaveArray2pgm(unsigned char A[], unsigned int n)
{

FILE * fp;
  char name [20]; /* name of file */
  sprintf(name,"%u",n); /*  */
  char *filename =strcat(name,".pgm");
  const char *comment="# C= ";/* 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,iXmax,iYmax,MaxColorComponentValue);  /*write header to the file*/
  fwrite(A,iLength,1,fp);  /*write image data bytes to the file in one step */
  printf("File %s saved. \n", filename);
  fclose(fp);
  return 0;
}
 



unsigned long long int GiveNrOfCrOrbitPoints( int n){

  unsigned long long int i = fabs(get_double(7 / (2 - 7 * tt(n))));
  return i;
}

int setup(int n){

  

/* unsigned int iX,iY;  indices of 2D virtual array (image) = integer coordinate */
  iXmax = iSide-1; /* height of image in pixels */
  iYmax = iSide-1;
  iLength = (iSide*iSide);
 
 //
  PixelWidth =  ((ZxMax-ZxMin)/iXmax);
  PixelHeight = ((ZyMax-ZyMin)/iYmax);
  invPixelWidth = 1 / PixelWidth;
  invPixelHeight = 1 / PixelHeight;
 
 
  //
  data = (unsigned char *) malloc( iLength * sizeof(unsigned char) );
  if (data == NULL )
    {
      fprintf(stderr," Could not allocate memory");
      return 1;
    }
  
          
  ClearArray( data );
  
  //
  C = GiveC(tt(n), radius, period);
  Cx = creal(C); 
  Cy = cimag(C); 
  // 
  NrOfCrOrbitPoints = GiveNrOfCrOrbitPoints(n); //pow(10,5+n);
 

return 0;
}


 int info(int n){

  printf("   period = %d  \n", period);
  printf("   n = %d  \n", n);
  printf("   t = %.16f \n", get_double(tt(n)));
  printf("   c = %.16f , %.16f  \n",get_double(Cx), get_double(Cy));
  printf("   NrOfCrOrbitPoints = %.llu  \n",NrOfCrOrbitPoints);
  return 0;
}
 






/* --------------------------------------------------------------------------------------------------------- */
 
int main(){
  
int n;   
 
  for (n=10; n<=12; ++n){ 

   setup(n);
   //TestCriticalOrbit(NrOfCrOrbitPoints);
  /* ------------------ draw ----------------------- */
  printf(" for n = %d draw %llu points of critical orbit to array \n",n, NrOfCrOrbitPoints);       
  DrawCriticalOrbit(NrOfCrOrbitPoints, data);
  printf(" save  data array to the pgm file \n");
 SaveArray2pgm(data, n);

  //CheckOrientation(data);
  //SaveArray2pgm(data, n+1000);
   
 // info(n); 

 }
  
  free(data);
  
 
 
 
 
 
  return 0;
}