/*
 
errors : why graphic file has 999 pixel not 1000 ? 

  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 


  forward 
  backward


===========================
Q: 
Let's take 
t =  0.3332223278292314 

find bifurcation point from period 1

so


c = -0.124396035791842  +0.649518736914556 i    period = 10000


then go to dynamic plane
increase iteration number to see better Julia set
then draw better critical orbit


then go back to critical point z=0

Then use
AA BAA AAA........

to stay on critical orbit ( and to draw rare points of critical orbits ) 

Is it possiible to find other sequences of inverse iterations ( A, B)
to draw rare points of such critical orbits ?
------------ Answer -----------------------
No,  there is only one backwards orbit contained in the boundary of the
Siegel disk.  Ideally it should be AAAAAAAAAAAAAAAAAAA... , but since
the commands  a  and  b  work only approximately, it might be different.

  Best regards,

  Wolf

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



  










*/

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

int n; 
#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;
unsigned char *data3;


/* */
double ZxMin = -0.75;
double ZxMax = 0.35;
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 iPeriodChild=3;

// unsigned int m;
// check rang of type for big numbers : iMax, i, NrOfCrOrbitPoints
unsigned long long int NrOfCrOrbitPoints ;
 
double InnerRadius= 0.05997604;
 
Real radius = 1.0;
Real tt(int n)
{
  // t(n) = [0;3,  10^n, golden_ratio]
  Real phi = (sqrt(Real(5)) + 1) / 2;
  Real tenn = pow(Real(10), n);
  return 0 +1/( 3  +1/( tenn +1/( phi )));
}

/* fc(z) = z*z + c */
/*   */
Cplx C; 
Real Cx; /* C = Cx + Cy*i   */
Real Cy ;
 
Cplx alfa; 
Real Ax;
Real Ay;

Cplx z3a, z3b, z3c;
 
/* colors */
const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
const int iExterior = 150; /* exterior of Julia set */
const int iJulia = 0; /* border , boundary*/
const int iJuliaInv = 50; /* border , boundary*/
const int iInterior = 200;
const int iBackground = 255; 
 
 
/* ----------------------- 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;
}


/*

  http://en.wikipedia.org/wiki/Periodic_points_of_complex_quadratic_mappings
  z^2 + c = z
  z^2 - z + c = 0
  ax^2 +bx + c =0 // ge3neral for  of quadratic equation
  so :
  a=1
  b =-1
  c = c
  so :

  The discriminant is the  d=b^2- 4ac 

  d=1-4c = dx+dy*i
  r(d)=sqrt(dx^2 + dy^2)
  sqrt(d) = sqrt((r+dx)/2)+-sqrt((r-dx)/2)*i = sx +- sy*i

  x1=(1+sqrt(d))/2 = beta = (1+sx+sy*i)/2

  x2=(1-sqrt(d))/2 = alfa = (1-sx -sy*i)/2

  alfa : attracting when c is in main cardioid of Mandelbrot set, then it is in interior of Filled-in Julia set, 
  it means belongs to Fatou set ( strictly to basin of attraction of finite fixed point )

*/
// uses global variables : 
//  ax, ay (output = alfa(c)) 
Cplx GiveAlfaFixedPoint(Cplx c)
{
  Real dx, dy; //The discriminant is the  d=b^2- 4ac = dx+dy*i
  Real r; // r(d)=sqrt(dx^2 + dy^2)
  Real sx, sy; // s = sqrt(d) = sqrt((r+dx)/2)+-sqrt((r-dx)/2)*i = sx + sy*i
  Real ax, ay;
 
  // d=1-4c = dx+dy*i
  dx = 1 - 4*creal(c);
  dy = -4 * cimag(c);
  // r(d)=sqrt(dx^2 + dy^2)
  r = sqrt(dx*dx + dy*dy);
  //sqrt(d) = s =sx +sy*i
  sx = sqrt((r+dx)/2);
  sy = sqrt((r-dx)/2);
  // alfa = ax +ay*i = (1-sqrt(d))/2 = (1-sx + sy*i)/2
  ax = 0.5 - sx/2.0;
  ay =  sy/2.0;
 

  return ax+ay*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( Real Zx, Real Zy, unsigned char data[], int iColor)
{
  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", get_double(Zx), get_double(Zy)); return -1; } 
#endif
  
  iX = (int)get_double((Zx-ZxMin)*invPixelWidth);
  iY = (int)get_double((ZyMax-Zy)*invPixelHeight); // reverse Y axis
  index = iX + iY*iXmax;//f(iX,iY);
  
  
  data[index] = iColor;  /* 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;
 
}


/* mndyncxmics::root from mndyncxmo.cpp  by Wolf Jung (C) 2007-2014. */

// input = x,y
// output = u+v*I = sqrt(x+y*i) 
Cplx
GiveRoot (Cplx z)
{
  Real x = creal (z);
  Real y = cimag (z);
  Real u, v;

  v = sqrt (x * x + y * y);

  if (x > 0.0)
    {
      u = sqrt (0.5 * (v + x));
      v = 0.5 * y / u;
      return u + v * I;
    }
  if (x < 0.0)
    {
      v = sqrt(Real(0.5 * (v - x)));
      if (y < 0.0)
	v = -v;
      u = 0.5 * y / v;
      return u + v * I;
    }
  if (y >= 0.0)
    {
      u = sqrt(Real(0.5 * y));
      v = u;
      return u + v * I;
    }

  u = sqrt (Real (-0.5 * y));
  v = -u;
  return u + v * I;
}

// from mndlbrot.cpp  by Wolf Jung (C) 2007-2014. part of Madel 5.12 
// input : c, z , mode
// c = cx+cy*i where cx and cy are global variables defined in mndynamo.h
// z = x+y*i
// 
// output : z = x+y*i
Cplx InverseIteration(Cplx z, Cplx c)
{
  Real x = creal (z);
  Real y = cimag (z);
  Real cx = creal (c);
  Real cy = cimag (c);

  // f^{-1}(z) = inverse with principal value
  if (cx * cx + cy * cy < 1e-20)
    {
      z = GiveRoot (x - cx + (y - cy) * I);	// 2-nd inverse function = key b 
      //if (mode & 1) { x = -x; y = -y; } // 1-st inverse function = key a   
      return -z;
    }

  //f^{-1}(z) =  inverse with argument adjusted
  Real u, v;
  Cplx uv;
  Real w = cx * cx + cy * cy;

  uv = GiveRoot (-cx / w - (cy / w) * I);
  u = creal (uv);
  v = cimag (uv);
  //
  z = GiveRoot (w - cx * x - cy * y + (cy * x - cx * y) * I);
  x = creal (z);
  y = cimag (z);
  //
  w = u * x - v * y;
  y = u * y + v * x;
  x = w;

  //if (mode & 1) // mode = -1
  //  { x = -x; y = -y; } // 1-st inverse function = key a

  return x + y * I;		// 2-nd inverse function = key b 
}


 int DrawInverseOrbit(unsigned long long int iMax,  unsigned char A[], unsigned char iColor ){

    Cplx Z;

/* critical point z = 0 */
 
 
 Z = -  InverseIteration(Cplx(0.0,0.0), C); // A
 DrawPoint(get_double(creal(Z)),get_double(cimag(Z)),A, iColor);
 
 Z = -  InverseIteration(Z, C); // A
 DrawPoint(get_double(creal(Z)),get_double(cimag(Z)),A, iColor);

 Z =   InverseIteration(Z, C); // B
 DrawPoint(get_double(creal(Z)),get_double(cimag(Z)),A, iColor);
 for (unsigned long long int i = 0; i < iMax; ++i){
 Z = -  InverseIteration(Z, C); // A
 DrawPoint(get_double(creal(Z)),get_double(cimag(Z)),A, iColor);
 }
 


return 0; 

}


// forward iteration
// Cx and Cy = global
Cplx Forward(Cplx z){
  Real Zx,Zy, tmp;
  Zx = creal(z);
  Zy = cimag(z);

	  /* f(z) = z*z+c */
	  tmp = mul_pwr2(Zx*Zy, 2) + Cy;
	  Zx = sqr(Zx) - sqr(Zy) + Cx;
	  Zy = tmp;
   return Cplx(Zx,Zy);

}

int DrawForwardOrbit(Cplx Z, unsigned long long int iMax,  unsigned char A[], unsigned char iColor  )
{
  //unsigned long long int iMax10000 = iMax / 10000;
  unsigned long long int i; /* nr of point of critical orbit */
  Real Zx,Zy, tmp;
  int IsGood;  
 
  /* Z0 */
  Zx = creal(Z);
  Zy = cimag(Z);
  DrawPoint(Zx,Zy,A,iColor);
  
   
      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;}
      
	 // IsGood =
          DrawPoint(Zx,Zy,A,iColor); /* draws critical orbit */
#ifdef BOUNDS_CHECKS
	//  if (unlikely(IsGood<0)) { printf ("i = %llu \n", i); return 1;}
#endif
	}
    
    
   
  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 iMaxInverse, unsigned long long int iMaxForward, unsigned char A[] )
{
  
  unsigned long long int i; /* nr of point of critical orbit */
  Cplx Z0;
  Cplx Z; 
 
   
  /* critical point z = 0 */
  Z0 = Cplx(0.0,0.0);
  //printf("   point z = %f , %f   \n", get_double(creal(Z)), get_double(cimag(Z)));
  printf("%dC",n); // symbolic sequence starting from n , then  crritical ,  

  DrawForwardOrbit(Z0, iMaxForward, A, iJulia);



  // first
  Z =  InverseIteration(Z, C); 
  if (cimag(Z) >0.0) // manual adjust 
                  {printf("B"); }
             else {Z = -Z; printf("A"); } 
  //printf("   point z = %f , %f   \n", get_double(creal(Z)), get_double(cimag(Z)));
  DrawPoint( get_double(creal(Z)), get_double(cimag(Z)),A,iJuliaInv);


  // second
  //DrawForwardOrbit(Z, 2, A);
  Z =  InverseIteration(Z, C); 
       if ( get_double(cimag(Z)) >0.0) 
                  {printf("B"); }
             else {Z = -Z; printf("A"); } 
  //printf("   point z = %f , %f   \n", get_double(creal(Z)), get_double(cimag(Z)));
 // DrawForwardOrbit(Z, 2, A);
  DrawPoint( get_double(creal(Z)), get_double(cimag(Z)),A,iJuliaInv);




  // third 
  Z =  InverseIteration(Z, C); 
       if ( get_double(cimag(Z)) >0.0) 
                  {printf("B"); }
             else {Z = -Z; printf("A"); } 

  //printf("   point z = %f , %f   \n", get_double(creal(Z)), get_double(cimag(Z)));



   for (i=1;i<=iMaxInverse ;i++){
     //DrawForwardOrbit(Z,2, A);
       DrawPoint( get_double(creal(Z)), get_double(cimag(Z)),A,iJuliaInv);
       Z =  InverseIteration(Z, C); 
       if ( get_double(cimag(Z)) >0.0) 
                  {printf("B"); }
             else {Z = -Z; printf("A"); } 
      // if ((get_double(creal(Z)) < ZxMin || ZxMax  < get_double(creal(Z)) || get_double(cimag(Z)) < ZyMin || ZyMax < get_double(cimag(Z)))) {  printf("   point z = %f , %f out of bounds  \n", get_double(creal(Z)), get_double(cimag(Z))); return -1; } 
    }
   printf("\n"); // symbolic sequence 
  return 0;
 
}



/*
The distance between (x1, y1) and (x2, y2) is given by:
https://en.wikipedia.org/wiki/Distance

*/

Real GiveD(Real x1, Real y1, Real x2, Real y2 ){

return sqrt((x1-x2)*(x1-x2) +(y1-y2)*(y1-y2));


}

int MarkSets(unsigned long long int iterMax,  unsigned char A[] )
{
  int ix, iy; // pixel coordinate 
  Real Zx, Zy; //  Z= Zx+ZY*i;
  Real Zx2, Zy2;
  Real d; 
  
  int i; /* index of 1D array */
  unsigned int iter;
  for(iy=0;iy<=iYmax;++iy) 
    {
      //Zy =ZyMax - iy*PixelHeight;
      for(ix=0;ix<=iXmax;++ix) 
	{
          i = f(ix, iy); /* compute index of 1D array from indices of 2D array */
          if ( A[i] != iJulia) {

	  // from screen to world coordinate
          Zy =ZyMax - iy*PixelHeight; 
	  Zx = ZxMin + ix*PixelWidth ;
	  
          Zx2=Zx*Zx;
          Zy2=Zy*Zy;

          for (iter=1;iter<=iterMax ;iter++)
	{
	  /* f(z) = z*z+c */
	  Zy=2*Zx*Zy + Cy;
          Zx=Zx2-Zy2 +Cx;
          Zx2=Zx*Zx;
          Zy2=Zy*Zy;
       
        if (Zx2+Zy2>4) {A[i]=iExterior; break;}
        d = GiveD(Zx, Zy, Ax, Ay);
        if ( d< InnerRadius) {A[i]=iInterior; break;}
        

       } //iter
        
      } // if ( A[i] != iJulia)
	  

	} // ix
    } //iy
   
    
   
  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]=iBackground;
  return 0;
}


int ClearArray3( unsigned char A[] )
{
  int index; /* index of 1D array  */

  int k ; //
  for(index=0;index<iLength-1;++index) {

    k = index * 3;
    //printf("index = %d k = %d \n", index,k);
    A[k]   = iBackground;
    A[k+1] = iBackground;
    A[k+2] = iBackground;
   }
  return 0;
}


// https://commons.wikimedia.org/wiki/File:Eight_point_neighborhood.svg
int GiveNeighbours(unsigned char S[], int i){

int neighbours =0;
int iX, iY;
div_t result; // ang. result = iloraz 
  
 // i = iX + iY*iXmax;
 // compute iX and iY from i
 
 result = div(i, iXmax);
 iY = result.quot;
 iX = result.rem;

 if ( iX<1 || iY <1 || iX>iXmax-1 || iY> iYmax-1 ) {printf("out of bounds in GiveNeighbours\n"); return -1;}
 
 // iY-1
  if (S[f(iX-1,iY-1)]==iJulia) neighbours+=1;
  if (S[f(iX,iY-1)]==iJulia) neighbours+=1;
  if (S[f(iX+1,iY-1)]==iJulia) neighbours+=1;
 // iY
  if (S[f(iX-1,iY)]==iJulia) neighbours+=1;
  if (S[f(iX+1,iY)]==iJulia) neighbours+=1;
 // iY+1
  if (S[f(iX-1,iY+1)]==iJulia) neighbours+=1;
  if (S[f(iX,iY+1)]==iJulia) neighbours+=1;
  if (S[f(iX+1,iY+1)]==iJulia) neighbours+=1;




return neighbours;


}



/*
 S = source = 1 bit color
 D = destination = 3 bit color
*/

int CopyArray( unsigned char S[], unsigned char D[] )
{
  int i; /* index of 1D array  */
  //int Neighbours;
  int k ; //
  for(i=0;i<iLength-1;++i) {

     k = i*3;
    if (S[i]==iJulia){

   // Neighbours = GiveNeighbours(S,i);
    
   
    //printf("index = %d k = %d \n", index,k);

    //if (Neighbours<2) // disconnected 
       {D[k]=255;
        D[k+1]=0;
        D[k+2]=0;}
    //  else { // connected
    //    D[k]=0;
    //    D[k+1]=0;
   //     D[k+2]=255;}
    } // julia


     if (S[i]==iJuliaInv){

    //Neighbours = GiveNeighbours(S,i);
    
   
    //printf("index = %d k = %d \n", index,k);

    //if (Neighbours<2) // disconnected 
       {D[k]=0;
        D[k+1]=0;
        D[k+2]=255;}
    //  else { // connected
     //   D[k]=0;
     //   D[k+1]=255;
      //  D[k+2]=255;}
    } // julia



    
    if (S[i]==iInterior){
         D[k]=iInterior;
        D[k+1]=iInterior;
        D[k+2]=iInterior;}// interior
    


    if (S[i]==iExterior){
         D[k]=iExterior;
        D[k+1]=iExterior;
        D[k+2]=iExterior;}// interior
     //
} //
   
  return 0;
}

/*
 S = source = 1 bit color
 D = destination = 3 bit color
*/

int CopyTestArray( unsigned char S[], unsigned char D[] )
{
  int i; /* index of 1D array  */
  int Neighbours;
  int k ; //
  for(i=0;i<iLength-1;++i) {

     k = i*3;
    if (S[i]==iJulia){

    Neighbours = GiveNeighbours(S,i);
    
   
    //printf("index = %d k = %d \n", index,k);

    if (Neighbours<2) // disconnected 
       {D[k]=255;
        D[k+1]=0;
        D[k+2]=0;}
      else { // connected
        D[k]=0;
        D[k+1]=0;
        D[k+2]=255;}
    } // julia


     if (S[i]==iJuliaInv){

    Neighbours = GiveNeighbours(S,i);
    
   
    //printf("index = %d k = %d \n", index,k);

    if (Neighbours<2) // disconnected 
       {D[k]=255;
        D[k+1]=0;
        D[k+2]=255;}
      else { // connected
        D[k]=0;
        D[k+1]=255;
        D[k+2]=255;}
    } // julia



    
    if (S[i]==iInterior){
         D[k]=iInterior;
        D[k+1]=iInterior;
        D[k+2]=iInterior;}// interior
    


    if (S[i]==iExterior){
         D[k]=iExterior;
        D[k+1]=iExterior;
        D[k+2]=iExterior;}// interior
     //
} //
   
  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;
}

/*
http://rosettacode.org/wiki/Bitmap/Bresenham%27s_line_algorithm
Instead of swaps in the initialisation use error calculation for both directions x and y simultaneously:
*/
void iDrawLine( int x0, int y0, int x1, int y1, unsigned char iColor, unsigned char A[]) 
{
  int x=x0; int y=y0;
  int dx = abs(x1-x0), sx = x0<x1 ? 1 : -1;
  int dy = abs(y1-y0), sy = y0<y1 ? 1 : -1; 
  int err = (dx>dy ? dx : -dy)/2, e2;
  int index; 
  e2 = err;
  for(;;){
    //iDrawPoint(A, x, y,color);

   index = x + y*iXmax;//f(iX,iY);
   A[index] = iColor;  /* draw */




   if (x==x1 && y==y1) break;
   // if (x<0 || x>iXmax || y<0 || y > iYmax ) break; // cliping 
    
    if (e2 >-dx) { err -= dy; x += sx; }
    if (e2 < dy) { err += dx; y += sy; }
    e2 = err;
    
  }
}
/*

paritition a dynamic plane into k=PeriodChild sectors
with center = alf afixed point 
when point z from Siegel disc if iterated 
it moves from one sector to other
*/


int dDrawLine(double Zx0, double Zy0, double Zx1, double Zy1, unsigned char iColor, unsigned char A[]) 
{

  int ix0, iy0; // screen coordinate = indices of virtual 2D array 
  int ix1, iy1; // screen coordinate = indices of virtual 2D array

   // first step of clipping
   if (  Zx0 < ZxMax &&  Zx0 > ZxMin && Zy0 > ZyMin && Zy0 <ZyMax 
      && Zx1 < ZxMax &&  Zx1 > ZxMin && Zy1 > ZyMin && Zy1 <ZyMax )
   ix0= (Zx0- ZxMin)/PixelWidth; 
   iy0 = (ZyMax - Zy0)/PixelHeight; // inverse Y axis 
   ix1= (Zx1- ZxMin)/PixelWidth; 
   iy1= (ZyMax - Zy1)/PixelHeight; // inverse Y axis 
   // second step of clipping
   if (ix0 >=0 && ix0<=iXmax && ix1 >=0 && ix1<=iXmax && iy0 >=0 && iy0<=iYmax && iy1 >=0 && iy1<=iYmax )

   //
   iDrawLine(ix0,iy0,ix1,iy1, iColor, A) ;

return 0;
}

/* 
from mndlbrot.cpp  part of Mandel 5.12 by Wolf Jung (C) 2007-2014.  
 the GNU General Public License as published by the Free Software Foundation; either
   version 3, or (at your option) any later version.

int find(int sg, unsigned int preper, unsigned int per, double cx, double cy, double &x, double &y) 
{  double a = cx, b = cy, fx, fy, px, py, w; 
   unsigned int i, j;
   for (i = 0; i < 30; i++)
   {  if (sg > 0) { a = x; b = y; }
      if (!preper)
      {  if (sg > 0) { fx = 0; fy = 0; px = 0; py = 0; }
         else { fx = -x; fy = -y; px = -1; py = 0; }
      }
      else
      {  fx = x; fy = y; px = 1.0; py = 0;
         for (j = 1; j < preper; j++)
         {  if (px*px + py*py > 1e100) return 1;
            w = 2*(fx*px - fy*py); py = 2*(fx*py + fy*px);
            px = w; if (sg > 0) px++;
            w = fx*fx - fy*fy + a; fy = 2*fx*fy + b; fx = w;
         }
      }
      double Fx = fx, Fy = fy, Px = px, Py = py;
      for (j = 0; j < per; j++)
      {  if (px*px + py*py > 1e100) return 2;
         w = 2*(fx*px - fy*py); py = 2*(fx*py + fy*px);
         px = w; if (sg > 0) px++;
         w = fx*fx - fy*fy + a; fy = 2*fx*fy + b; fx = w;
      }
      fx += Fx; fy += Fy; px += Px; py += Py;
      w = px*px + py*py; if (w < 1e-100) return -1;
      x -= (fx*px + fy*py)/w; y += (fx*py - fy*px)/w;
   }
   return 0;
}



*/













Cplx  find( int preper, int period, Real a, Real b) 
{  
   Real x=0.0;
   Real y=0.0;
   Real  fx=0.0;
   Real  fy=0.0;
   Real px=0.0;
   Real py=0.0;
   Real w=0.0; 
   int i, j;

   Real Px = 0.0;
   Real Py = 0.0;
   Real Fx=0.0;
   Real Fy=0.0;


    for (i = 0; i < 30; i++)
   { 
      if (!preper)
      { { fx = -x; fy = -y; px = -1.0; py = 0.0; }
      }
      else
      {  fx = x; fy = y; px = 1.0; py = 0.0;
         for (j = 1; j < preper; j++)
         {  if (px*px + py*py > 1e100) return Cplx(0.0,0.0);
            w = 2*(fx*px - fy*py); py = 2*(fx*py + fy*px);
            px = w; 
            w = fx*fx - fy*fy + a; fy = 2*fx*fy + b; fx = w;
         }
      }
      Real Fx = fx, Fy = fy, Px = px, Py = py;
      for (j = 0; j < period; j++)
      {  if (px*px + py*py > 1e100) Cplx(1.0,0.0);
         w = 2*(fx*px - fy*py); py = 2*(fx*py + fy*px);
         px = w;  
         w = fx*fx - fy*fy + a; fy = 2*fx*fy + b; fx = w;
      }
      fx += Fx; fy += Fy; px += Px; py += Py;
      w = px*px + py*py; if (w < 1e-100) return Cplx(-1,0.0);
      x -= (fx*px + fy*py)/w; y += (fx*py - fy*px)/w;
   }
   return Cplx(x,y);
}


int DrawStar(unsigned char A[] ){

  dDrawLine(get_double(Ax), get_double(Ay),get_double(creal(z3a)), get_double(cimag(z3a)), 0, A);
  dDrawLine(get_double(Ax), get_double(Ay),get_double(creal(z3b)), get_double(cimag(z3b)), 0, A);
  dDrawLine(get_double(Ax), get_double(Ay),get_double(creal(z3c)), get_double(cimag(z3c)), 0, A);

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




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

  FILE * fp;
  char name [20]; /* name of file */
  sprintf(name,"%u",n); /*  */
  char *filename =strcat(name,".ppm");
  //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,"P6\n %u %u\n %u\n",iXmax,iYmax,MaxColorComponentValue);  /*write header to the file*/
  fwrite(A,iLength*3,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 ;

  //i = 10000000000;
  // i = fabs(get_double(7 / (2 - 7 * tt(n))));

  unsigned long long int i ; 
  switch( n )
  { case  0 : i =      100000; break;
    case  1 : i =      100000; break;  
    case  2 : i =      1000000; break;
    case  3 : i =      10000000; break;
    case  4 : i =      100000000; break;
    case  5 : i =      1000000000; break;
    case  6 : i =      10000000000; break;
    case  7 : i =      10000000000; break;
    case  8 : i =      100000000000; break;
    case  9 : i =      100000000000; break;
    case 10 : i =      1000000; break;
    case 11 : i =      1000000; break;
    default: i = 100000; break;
    }
  return i;


  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;
 
 
  // 1 byte color 
  data = (unsigned char *) malloc( iLength * sizeof(unsigned char) );
  if (data == NULL )
    {
      fprintf(stderr," Could not allocate memory");
      return 1;
    }
  
          
  

  // 3 byte color
  data3 = (unsigned char *) malloc( 3*iLength * sizeof(unsigned char) );
  if (data3 == NULL )
    {
      fprintf(stderr," Could not allocate memory");
      return 1;
    }
  ClearArray( data );
  ClearArray3( data3 );

   
  //
  C = GiveC(tt(n), radius, period);
  Cx = creal(C); 
  Cy = cimag(C); 
  //
  alfa = GiveAlfaFixedPoint(C);
  Ax = creal (alfa);
  Ay = cimag (alfa);
  // 
  NrOfCrOrbitPoints = GiveNrOfCrOrbitPoints(n); //pow(10,5+n);
  //

  z3a = find(0,iPeriodChild,Cx, Cy); 
  z3b = Forward(z3a);
  z3c = Forward(z3b);

  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("   alfa = %.16f , %.16f  \n",get_double(Ax), get_double(Ay));
  printf(" periodic orbit : \n");
  printf("   z3a = %.16f , %.16f  \n",get_double(creal(z3a)), get_double(cimag(z3a)));
  printf("   z3a = %.16f , %.16f  \n",get_double(creal(z3b)), get_double(cimag(z3b)));
  printf("   z3a = %.16f , %.16f  \n",get_double(creal(z3c)), get_double(cimag(z3c)));

  printf("   NrOfCrOrbitPoints = %.llu  \n",NrOfCrOrbitPoints);
  printf("   Inner Radius = %.16f \n", InnerRadius);
  //
  printf("   iSide = %d  \n", iSide);
  printf("   iXmax = %d  \n", iXmax);
  printf("   iLength = %d  \n", iLength);


  return 0;
}
 






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

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

    


   // SaveArray2ppm(data3,100+n);
    //CheckOrientation(data);
    //SaveArray2pgm(data, n+1000);

    //DrawStar(data);
    //SaveArray2pgm(data, 100+n);


     info(n); 

  }
  

  //
  free(data);
  free(data3);
 
 
 
 
 
  return 0;
}