/*

  Adam Majewski
  fraktal.republika.pl


  c console progam using 
  * symmetry
  * openMP

  draw  julia sets



  gcc f.c -lm -Wall -fopenmp -march=native 
  time ./a.out
  time ./a.out > info.txt


   

forward iteration with 2 tests : 
* for exterior =  escape to infinity ( bailout test ) with exter of circle with center 0.0 and radius 2.0
* attraction to parabolic fixed point alfa with target set =  wide triangle inside last component of immediate basin



wide triangle aproximates last component of filled Julia set near alfa fixed point
wide triangle is defined by 3 points : 
* parabolic alfa fixed point Za 
* 2 points on the boundary of Julia set. These points are found in procedure ComputeRays


narrow triangle is defined by 3  points : 
* parabolic alfa fixed point Za 
* point of critical orbit inside 0 componet
* 1 point on the boundary of last componet of the Julia set. This points is found in procedure ComputeRays
 

last componet = component  to the left of component including critical point zcr = 0.0
 this component is almost not changing when iPeriodChild is increasing , see 

https://commons.wikimedia.org/wiki/File:Parabolic_rays_landing_on_fixed_point.ogv




another option is 
polygon test !!!!
http://stackoverflow.com/questions/11716268/point-in-polygon-algorithm
http://apodeline.free.fr/FAQ/CGAFAQ/CGAFAQ-3.html
http://alienryderflex.com/polyspline/
https://www.ecse.rpi.edu/~wrf/Research/Short_Notes/pnpoly.html
http://geomalgorithms.com/a03-_inclusion.html


memcpy(dest, src, sizeof(int) * 10);
ULONG_LONG_MAX = 2^x 



iterMax  = 10000000; real	76m27.826s


last point of the ray for angle = 1 / 1073741823 = 0.0000000009313226 is = (0.3249187638558812 ;  0.1396868612816465 ) ; Distance to fixed = 112  pixels 
last point of the ray for angle = 2 / 1073741823 = 0.0000000018626452 is = (0.3467472200853008 ;  0.0930454493000004 ) ; Distance to fixed = 95  pixels 
last point of the ray for angle = 4 / 1073741823 = 0.0000000037252903 is = (0.3722636151216278 ;  0.0667981865256394 ) ; Distance to fixed = 82  pixels 
last point of the ray for angle = 8 / 1073741823 = 0.0000000074505806 is = (0.3948056374899506 ;  0.0520047535676791 ) ; Distance to fixed = 72  pixels 
last point of the ray for angle = 16 / 1073741823 = 0.0000000149011612 is = (0.4138544330353308 ;  0.0433352245497908 ) ; Distance to fixed = 64  pixels 
last point of the ray for angle = 32 / 1073741823 = 0.0000000298023224 is = (0.4300849860618463 ;  0.0381406343691283 ) ; Distance to fixed = 59  pixels 
last point of the ray for angle = 64 / 1073741823 = 0.0000000596046448 is = (0.4442058232194739 ;  0.0350791132213122 ) ; Distance to fixed = 55  pixels 
last point of the ray for angle = 128 / 1073741823 = 0.0000001192092897 is = (0.4567757051281010 ;  0.0334363775867163 ) ; Distance to fixed = 52  pixels 
last point of the ray for angle = 256 / 1073741823 = 0.0000002384185793 is = (0.4682134893147842 ;  0.0328175348022398 ) ; Distance to fixed = 49  pixels 
last point of the ray for angle = 512 / 1073741823 = 0.0000004768371586 is = (0.4788343167618370 ;  0.0330029099389567 ) ; Distance to fixed = 48  pixels 
last point of the ray for angle = 1024 / 1073741823 = 0.0000009536743173 is = (0.4888805465402919 ;  0.0338775367387197 ) ; Distance to fixed = 47  pixels 
last point of the ray for angle = 2048 / 1073741823 = 0.0000019073486346 is = (0.4985439369887921 ;  0.0353958224973381 ) ; Distance to fixed = 46  pixels 
last point of the ray for angle = 4096 / 1073741823 = 0.0000038146972692 is = (0.5079806286332533 ;  0.0375644305208307 ) ; Distance to fixed = 46  pixels 
last point of the ray for angle = 8192 / 1073741823 = 0.0000076293945384 is = (0.5173206684734973 ;  0.0404356910951789 ) ; Distance to fixed = 46  pixels 
last point of the ray for angle = 16384 / 1073741823 = 0.0000152587890767 is = (0.5266730648027788 ;  0.0441081224435497 ) ; Distance to fixed = 47  pixels 
last point of the ray for angle = 32768 / 1073741823 = 0.0000305175781534 is = (0.5361264266276797 ;  0.0487328049437455 ) ; Distance to fixed = 48  pixels 
last point of the ray for angle = 65536 / 1073741823 = 0.0000610351563068 is = (0.5457440949600917 ;  0.0545255739865563 ) ; Distance to fixed = 50  pixels 
last point of the ray for angle = 131072 / 1073741823 = 0.0001220703126137 is = (0.5555510148725334 ;  0.0617857048638559 ) ; Distance to fixed = 52  pixels 
last point of the ray for angle = 262144 / 1073741823 = 0.0002441406252274 is = (0.5655068927012731 ;  0.0709219068721949 ) ; Distance to fixed = 56  pixels 
last point of the ray for angle = 524288 / 1073741823 = 0.0004882812504547 is = (0.5754555647084196 ;  0.0824853391349234 ) ; Distance to fixed = 59  pixels 
last point of the ray for angle = 1048576 / 1073741823 = 0.0009765625009095 is = (0.5850327116561950 ;  0.0972049795890196 ) ; Distance to fixed = 64  pixels 
last point of the ray for angle = 2097152 / 1073741823 = 0.0019531250018190 is = (0.5935019015039471 ;  0.1160078703463312 ) ; Distance to fixed = 70  pixels 
last point of the ray for angle = 4194304 / 1073741823 = 0.0039062500036380 is = (0.5994741169304216 ;  0.1399734680216506 ) ; Distance to fixed = 77  pixels 
last point of the ray for angle = 8388608 / 1073741823 = 0.0078125000072760 is = (0.6004640809029653 ;  0.1700926269903841 ) ; Distance to fixed = 86  pixels 
last point of the ray for angle = 16777216 / 1073741823 = 0.0156250000145519 is = (0.5923130464369473 ;  0.2065407105412746 ) ; Distance to fixed = 97  pixels 
last point of the ray for angle = 33554432 / 1073741823 = 0.0312500000291038 is = (0.5688631155654750 ;  0.2469451995378435 ) ; Distance to fixed = 109  pixels 
last point of the ray for angle = 67108864 / 1073741823 = 0.0625000000582077 is = (0.5233107483736542 ;  0.2832277156244103 ) ; Distance to fixed = 122  pixels 
last point of the ray for angle = 134217728 / 1073741823 = 0.1250000001164153 is = (0.4543236362603006 ;  0.2987038999473556 ) ; Distance to fixed = 132  pixels 
last point of the ray for angle = 268435456 / 1073741823 = 0.2500000002328306 is = (0.3778733826407535 ;  0.2736881682116156 ) ; Distance to fixed = 135  pixels 
last point of the ray for angle = 536870912 / 1073741823 = 0.5000000004656613 is = (0.3285705161321115 ;  0.2091106321891394 ) ; Distance to fixed = 128  pixels 
Numerical approximation of parabolic Julia set for complex quadratic polynomial fc(z)= z^2 + c 
parameter c  = 0.2606874359562527 , 0.0022716846399296 
c is a root point between hyperbolic components of period 1 and 30  of Mandelbrot set 
combinatorial rotation number = internal angle  = 1 / 30 = 0.033333

image size in pixels = 2001 x 2001 
image size in world units : (ZxMin = -1.500000, ZxMax =  1.500000) ,  (ZyMin = -1.500000, ZyMax =  1.500000) 
ratio ( distortion) of image  = 1.000000 ; it should be 1.000 ...
PixelWidth = 0.001500 

parabolic alfa fixed point Za  = 0.4890738003669028 ; 0.1039558454088799 
critical point Zcr  = 0.0000000000000000 ; 0.0000000000000000 
precritical point Zl  = 0.3778733826407535 ; 0.2736881682116156 
precritical point Zr  = 0.3285705161321115 ; 0.2091106321891394 
Maximal number of iterations = iterMax = 100000000 
quality = 0.0000000000000000 = iNumberOfUknknown/ iNumberOfAllPixels  = 0 / 4004001 ; It should be 1.0  !!!! 
MaxDistance From Unknown 2 Fixed = 0.0000000000000000 [world units]= 0 [pixels] 

real	0m6.937s
user	0m6.940s
sys	0m0.000s


















*/



#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; needs also -fopenmp
#include <limits.h> // INT_MAX

/* --------------------------------- global variables and consts ------------------------------------------------------------ */
#define iPeriodChild 30 // iPeriodChild of secondary component joined by root point
int iPeriodParent = 1;
// internal angle , 
unsigned int numerator = 1; // see ULONG_MAX and integer overflow 
unsigned int denominator;
double InternalAngle; // numerator/denominator

// virtual 2D array and integer ( screen) coordinate
// Indexes of array starts from 0 not 1 
//unsigned int ix, iy; // var
unsigned int ixMin = 0; // Indexes of array starts from 0 not 1
unsigned int ixMax ; //
unsigned int iWidth ; // horizontal dimension of array
unsigned int ixAxisOfSymmetry  ; // 
unsigned int iyMin = 0; // Indexes of array starts from 0 not 1
unsigned int iyMax ; //
unsigned int iyAxisOfSymmetry  ; // 
//
unsigned int iyAboveMin = 1 ; //
unsigned int iyAboveMax ; //
//unsigned int iyBelow ; // var, measured from 1 to (iyAboveAxisLength -1)
unsigned int iyBelowMin = 0 ; //
unsigned int iyBelowMax ; //
unsigned int iyAboveAxisLength ; //
unsigned int iyBelowAxisLength ; //
unsigned int iHeight = 1000; //  odd number !!!!!! = (iyMax -iyMin + 1) = iyAboveAxisLength + iyBelowAxisLength +1
// The size of array has to be a positive constant integer 
unsigned int iSize ; // = iWidth*iHeight; 


// memmory 1D array 
unsigned char *data0;
unsigned char *data;
unsigned char *edge;

// unsigned int i; // var = index of 1D array
unsigned int iMin = 0; // Indexes of array starts from 0 not 1
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


/* world ( double) coordinate = dynamic plane described by center and radius */
double ZxMin=-1.5;
double ZxMax=1.5;
double ZyMin=-1.5;
double ZyMax=1.5;
double PixelWidth; // =(ZxMax-ZxMin)/iXmax;
double PixelWidth2; // =  PixelWidth*PixelWidth;
double PixelHeight; // =(ZyMax-ZyMin)/iYmax;

double ratio ;

// complex numbers of parametr plane 
double Cx; // c =Cx +Cy * i
double Cy;
double complex c; // 

double complex Za; // alfa fixed point alfa=f(alfa)
double Zax, Zay;


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




// points defining target sets 
double Zlx, Zly, Zrx, Zry; 
complex double Zl;
complex double zz_rays[iPeriodChild];


// critical point Zcr
double Zcrx = 0.0;
double Zcry=0.0;



unsigned char  iColorsOfInterior[iPeriodChild]; //={110, 160,210, 223, 240}; // number of colors >= iPeriodChild
static unsigned char iColorOfExterior = 245;
static unsigned char iColorOfUnknown = 100;
unsigned char iJulia = 0;


// unfortunately , because of lazy= slow dynamic
// some points z need very long time to reach attractor 
static unsigned long int iterMax  = 100000000; // ~ iHeight*100;
unsigned int iNumberOfUnknown = 0;
double quality ;
double MaxDistanceFromU2F2 = 0.0; // Max Distance from unknown  to Fixed
double radius = 0.1; // radius of circle around fixed 
double radius2; 
/* ------------------------------------------ functions -------------------------------------------------------------*/

// colors of components interior = shades of gray
int InitColors(int iMax, unsigned char a[])
{
  int i;
  //int iMax = iPeriodChild; iPeriodChild and iColorsOfInterior
  unsigned int iStep;

  iStep=  150/iMax;

  for (i = 1; i <= iMax; ++i)
    {
      a[i-1] = iColorOfExterior -i*iStep; 
      // printf("i= %d color = %i  \n",i-1, iColors[i-1]); // debug
    }
  return 0;
}


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

*/
double complex GiveC(double InternalAngleInTurns, double InternalRadius, unsigned int Period)
{
  //0 <= InternalRay<= 1
  //0 <= InternalAngleInTurns <=1
  double t = InternalAngleInTurns *2*M_PI; // from turns to radians
  double R2 = InternalRadius * InternalRadius;
  //double 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)) 
double complex GiveAlfaFixedPoint(double complex c)
{
  double dx, dy; //The discriminant is the  d=b^2- 4ac = dx+dy*i
  double r; // r(d)=sqrt(dx^2 + dy^2)
  double sx, sy; // s = sqrt(d) = sqrt((r+dx)/2)+-sqrt((r-dx)/2)*i = sx + sy*i
  double 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;
}





/* 
   principal square  root of complex number 
   http://en.wikipedia.org/wiki/Square_root

   z1= I;
   z2 = root(z1);
   printf("zx  = %f \n", creal(z2));
   printf("zy  = %f \n", cimag(z2));
*/
double complex root(double x, double y)
{ 
  
  double u;
  double v;
  double r = sqrt(x*x + y*y); 
  
  v = sqrt(0.5*(r - x));
  if (y < 0) v = -v; 
  u = sqrt(0.5*(r + x));
  return u + v*I;
}









// This function only works for periodic  angles.
// You must know the iPeriodChild n before calling this function.
// draws all "iPeriodChild" external rays 
// http://commons.wikimedia.org/wiki/File:Backward_Iteration.svg
// based on the code by Wolf Jung from program Mandel
// http://www.mndynamics.com/

int ComputeRays( //unsigned char A[],
			    int n, //iPeriodChild of ray's angle under doubling map
			    int iterMax,
                            complex double C, 
                            complex double Alfa,
                            double PixelWidth,
                            complex double zz[iPeriodChild] // output array

			    )
{
  

   
  //
  double Cx = creal(C); 
  double Cy = cimag(C);
  // 
  double dAlfaX = creal(Alfa);
  double dAlfaY = cimag(Alfa);


  double xNew; // new point of the ray
  double yNew;
  
  const double R = 10000; // very big radius = near infinity
  int j; // number of ray 
  int iter; // index of backward iteration

  double t,t0; // external angle in turns 
  double num, den; // t = num / den
  
  double complex zPrev;
  double u,v; // zPrev = u+v*I
 
  int iDistance ; // dDistance/PixelWidth = distance to fixed in pixels

  
  /* dynamic 1D arrays for coordinates ( x, y) of points with the same R on preperiodic and periodic rays  */
  double *RayXs, *RayYs;
  int iLength = n+2; // length of arrays ?? why +2



  //  creates arrays :  RayXs and RayYs  and checks if it was done
  RayXs = malloc( iLength * sizeof(double) );
  RayYs = malloc( iLength * sizeof(double) );
  if (RayXs == NULL || RayYs==NULL)
    {
      fprintf(stderr,"Could not allocate memory");
      getchar(); 
      return 1; // error
    }
  




   // external angle of the first ray 
   num = 1.0;
   den = pow(2.0,n) -1.0;
   t0 = num/den; // http://fraktal.republika.pl/mset_external_ray_m.html
   t=t0;
   // printf(" angle t = %.0f / %.0f = %f in turns \n", num, den, t0);


  //  starting points on preperiodic and periodic rays 
  //  with angles t, 2t, 4t...  and the same radius R
  for (j = 0; j < n; j++)
    { // z= R*exp(2*Pi*t)
      RayXs[j] = R*cos((2*M_PI)*t); 
      RayYs[j] = R*sin((2*M_PI)*t);
      //
      // printf(" %d angle t = = %.0f / %.0f =  %.16f in turns \n", j, num , den,  t); 
      //
      num *= 2.0;
      t *= 2.0; // t = 2*t
      if (t > 1.0) t--; // t = t modulo 1
      
    }
  //zNext = RayXs[0] + RayYs[0] *I;

  // printf("RayXs[0]  = %f \n", RayXs[0]);
  // printf("RayYs[0]  = %f \n", RayYs[0]);

  // z[k] is n-periodic. So it can be defined here explicitly as well.
  RayXs[n] = RayXs[0]; 
  RayYs[n] = RayYs[0];
  

  //   backward iteration of each point z
  for (iter = -10; iter <= iterMax; iter++)
    { 
     	
      for (j = 0; j < n; j++) // n +preperiod
	{ // u+v*i = sqrt(z-c)   backward iteration in fc plane 
	  zPrev = root(RayXs[j+1] - Cx , RayYs[j+1] - Cy ); // , u, v
	  u=creal(zPrev);
	  v=cimag(zPrev);
                
	  // choose one of 2 roots: u+v*i or -u-v*i
	  if (u*RayXs[j] + v*RayYs[j] > 0) 
	    { xNew = u; yNew = v; } // u+v*i
	  else { xNew = -u; yNew = -v; } // -u-v*i
                

	  // draw part of the ray = line from zPrev to zNew
	 // dDrawLine(A, RayXs[j], RayYs[j], xNew, yNew, j, 255);
                
	  
	  //  
	  RayXs[j] = xNew; RayYs[j] = yNew;
                
	       
                  

                
	} // for j ...

          //RayYs[n+k] cannot be constructed as a preimage of RayYs[n+k+1]
      RayXs[n] = RayXs[0]; 
      RayYs[n] = RayYs[0];
          
      // convert to pixel coordinates 
      //  if z  is in window then draw a line from (I,K) to (u,v) = part of ray 
   
      // printf("for iter = %d cabs(z) = %f \n", iter, cabs(RayXs[0] + RayYs[0]*I));
     
    }

  // aproximate end of ray by straight line to its landing point here = alfa fixed point
 // for (j = 0; j < n + 1; j++)
  //  dDrawLine(A, RayXs[j],RayYs[j], dAlfaX, dAlfaY,j, 255 ); 
 
  

  // this check can be done only from inside this function
  t=t0;
  num = 1.0;
  for (j = 0; j < n ; j++)
    {

      zz[j] = RayXs[j] +  RayYs[j] * I; // save to the output array
      // aproximate end of ray by straight line to its landing point here = alfa fixed point
      //dDrawLine(RayXs[j],RayYs[j], creal(alfa), cimag(alfa), 0, data); 
      iDistance = (int) round(sqrt((RayXs[j]-dAlfaX)*(RayXs[j]-dAlfaX) +  (RayYs[j]-dAlfaY)*(RayYs[j]-dAlfaY))/PixelWidth);
      printf("last point of the ray for angle = %.0f / %.0f = %.16f is = (%.16f ;  %.16f ) ; Distance to fixed = %d  pixels \n",num, den, t, RayXs[j], RayYs[j],  iDistance);
      num *= 2.0;
      t *= 2.0; // t = 2*t
      if (t > 1) t--; // t = t modulo 1
    } // end of the check 

 

  // free memmory
  free(RayXs);
  free(RayYs);

  
  

  

  return  0; //  
}




// ============ http://stackoverflow.com/questions/2049582/how-to-determine-a-point-in-a-2d-triangle
// In general, the simplest (and quite optimal) algorithm is checking on which side of the half-plane created by the edges the point is.
double side (double  x1, double y1, double x2,double y2,double x3, double y3)
{
    return (x1 - x3) * (y2 - y3) - (x2 - x3) * (y1 - y3);
}





// the triangle node numbering is counter-clockwise / clockwise  
int  PointInTriangle (double x, double y, double x1, double y1, double x2, double y2, double x3, double y3)
{
    int  b1, b2, b3;

    b1 = side(x, y, x1, y1, x2, y2) < 0.0;
    b2 = side(x, y, x2, y2, x3, y3) < 0.0;
    b3 = side(x, y, x3, y3, x1, y1) < 0.0;

    return ((b1 == b2) && (b2 == b3));
}


int InsideNarrowTriangle(double Zx, double Zy)
{
 
return PointInTriangle(Zx, Zy, Zax, Zay, Zrx, Zry,  Zlx, Zly)       ;
  
}




int InsideNarrowTriangle2a(double Zx, double Zy)
{
 
return PointInTriangle(Zx, Zy, Zax, Zay,  Zrx, Zry, Zcrx, Zcry)       ;
  
}

int InsideNarrowTriangle2b(double Zx,double Zy)
{
 
return PointInTriangle(Zx, Zy, Zax, Zay, Zcrx, Zcry, Zlx, Zly)       ;
  
}


int IsInsideWideTriangle(double Zx, double Zy)
{
 
return PointInTriangle(Zx, Zy, Zax, Zay, Zrx, Zry,  Zlx, Zly)       ;
  
}








int setup()
{

  
  

  denominator = iPeriodChild;
  InternalAngle = numerator/((double) denominator);

  c = GiveC(InternalAngle, 1.0, iPeriodParent) ;
  Cx=creal(c);
  Cy=cimag(c);
  Za = GiveAlfaFixedPoint(c);
  Zax = creal(Za);
  Zay = cimag(Za);



  //

  

 

   /* virtual 2D array ranges */
  if (!(iHeight % 2)) iHeight+=1; // it sholud be even number (variable % 2) or (variable & 1)
  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].
  iyAboveAxisLength = (iHeight -1)/2;
  iyAboveMax = iyAboveAxisLength ; 
  iyBelowAxisLength = iyAboveAxisLength; // the same 
  iyAxisOfSymmetry = iyMin + iyBelowAxisLength ; 
  // 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 ...
  
  

  // for numerical optimisation 
  ER2 = ER * ER;
  radius2 = radius*radius;
  PixelWidth2 =  PixelWidth*PixelWidth;
 
  /* create dynamic 1D arrays for colors ( shades of gray ) */
  data0 = malloc( iSize * sizeof(unsigned char) );
  data = malloc( iSize * sizeof(unsigned char) );
  edge = malloc( iSize * sizeof(unsigned char) );
  if (data0 == NULL || data == NULL || edge == NULL)
    {
      fprintf(stderr," Could not allocate memory\n");
      return 1;
    }
  else fprintf(stderr," memory is OK \n");

  
  
  

  

  // fill array iColorsOfInterior with iPeriodChild colors ( shades of gray )
  InitColors(iPeriodChild, iColorsOfInterior);



   // points defining target sets 
  ComputeRays( iPeriodChild, iterMax, c, Za, PixelWidth, zz_rays ) ;
  Zlx = creal(zz_rays[iPeriodChild-2]);   
  Zly = cimag(zz_rays[iPeriodChild-2]);
  Zrx = creal(zz_rays[iPeriodChild-1]); 
  Zry = cimag(zz_rays[iPeriodChild-1]);

  fprintf(stderr," setup done \n");
 
  return 0;

}



// from screen to world coordinate ; linear mapping
// uses global cons
double GiveZx(unsigned int ix)
{ return (ZxMin + ix*PixelWidth );}

// uses globaal cons
double GiveZy(unsigned int iy)
{ return (ZyMax - iy*PixelHeight);} // reverse y axis




unsigned char GiveColor(unsigned int ix, unsigned int iy)
{ 
  // check behavour of z under fc(z)=z^2+c
  // using 2 target set:
  // 1. exterior or circle (center at origin and radius ER ) 
  // as a target set containing infinity = for escaping points ( bailout test)
  // for points of exterior of julia set
  // 2. interior : triangle 

  double Zx0, Zy0;
  double Zx2, Zy2;
  int i=0;
  int j=0; // iteration = fc(z)
  
  double Zx, Zy;
  //double distance2 ;

  // from screen to world coordinate 
  Zx0 = GiveZx(ix);
  Zy0 = GiveZy(iy);

  //
  Zx = Zx0;
  Zy = Zy0;
  
  
 
  if (IsInsideWideTriangle( Zx, Zy)) 
       return iColorsOfInterior[iPeriodChild-1];
    
    

  while ( j<iterMax)
    { // then iterate and check behaviour
      for(i=0;i<iPeriodChild ;++i) // iMax = period !!!!
	{  
	  Zx2 = Zx*Zx; 
	  Zy2 = Zy*Zy;
       
	  // bailout test : if escaping ( = exterior ) then stop iteration
	  if (Zx2 + Zy2 > ER2) return iColorOfExterior; 
       
	  // if not escaping then iterate 
	  // new z : Z(n+1) = Zn * Zn  + C
	  Zy = 2*Zx*Zy + Cy; 
	  Zx = Zx2 - Zy2 + Cx; 
	  //  attracting = interior 
	  if (IsInsideWideTriangle( Zx, Zy)) return iColorsOfInterior[i];
          j+=1;
   	} // for(i
      }// while


  // unknown = not escaping and not attracting 
  //distance2 = (Zx0-Zax)*(Zx0-Zax) + (Zy0-Zay)*(Zy0-Zay);
  //if ( distance2 < radius2 ) 
  //{ if (distance2 > MaxDistanceFromU2F2) MaxDistanceFromU2F2 = distance2;
   //printf("a unknown ix = %d, iy = %d, Z  = %.16f ; %.16f ; Distance2Fixed= %f \n",ix, iy,  Zx0, Zy0, sqrt(distance2)); 
 // }
  iNumberOfUnknown +=1;
  //printf("b unknown ix = %d, iy = %d, Z  = %.16f ; %.16f ; Distance2Fixed= %f \n",ix,iy, Zx0, Zy0, sqrt(distance2));
  return iColorOfUnknown; // 
}


 




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


/* 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; }
//  ix = i % iWidth;
//  iy = (i- ix) / iWidth;
//  i  = Give_i(ix, iy);




// plots raster point (ix,iy) 
int PlotPoint(unsigned int ix, unsigned int iy, unsigned char iColor, unsigned char a[])
{
  unsigned i; /* index of 1D array */
  i = Give_i(ix,iy); /* compute index of 1D array from indices of 2D array */
  a[i] = iColor;

  return 0;
}


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


  // for all pixels of image 
  for(iy = iyMin; iy<=iyMax; ++iy) 
    { printf(" %d z %d\r", iy, iyMax); //info 
      for(ix= ixMin; ix<=ixMax; ++ix) PlotPoint(ix, iy, GiveColor(ix, iy) , a); //  
    } 
   

  
  return 0;
}



// fill array using symmetry of image 
// uses global var :  ...
int FillArraySymmetric(unsigned char a[] )
{
  unsigned int ix, iy; // var
  unsigned int iyAbove ; // var, measured from 1 to (iyAboveAxisLength -1) 
  unsigned char Color; // gray from 0 to 255 

  printf("axis of symmetry \n"); 
  iy = iyAxisOfSymmetry; 
#pragma omp parallel for schedule(dynamic) private(ix,Color) shared(ixMin,ixMax, iyAxisOfSymmetry)
  for(ix=ixMin;ix<=ixMax;++ix) {//printf(" %d from %d\n", ix, ixMax); //info  
    PlotPoint(ix, iy, GiveColor(ix, iy), a);
  }


  /*
    The use of ‘shared(variable, variable2) specifies that these variables should be shared among all the threads.
    The use of ‘private(variable, variable2)’ specifies that these variables should have a seperate instance in each thread.
  */

#pragma omp parallel for schedule(dynamic) private(iyAbove,ix,iy,Color) shared(iyAboveMin, iyAboveMax,ixMin,ixMax, iyAxisOfSymmetry)

  // above and below axis 
  for(iyAbove = iyAboveMin; iyAbove<=iyAboveMax; ++iyAbove) 
    {printf(" %d from %d\r", iyAbove, iyAboveMax); //info 
      for(ix=ixMin; ix<=ixMax; ++ix) 

	{ // above axis compute color and save it to the array
	  iy = iyAxisOfSymmetry + iyAbove;
	  Color = GiveColor(ix, iy);
	  PlotPoint(ix, iy, Color , a); 
	  // below the axis only copy Color the same as above without computing it 
	  PlotPoint(ixMax-ix, iyAxisOfSymmetry - iyAbove , Color , a); 
	} 
    }  
  return 0;
}













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











// Check Orientation of image : first quadrant in upper right position
// uses global var :  ...
int CheckOrientation(unsigned char a[] )
{
  unsigned int ix, iy; // pixel coordinate 
  double Zx, Zy; //  Z= Zx+ZY*i;
  unsigned i; /* index of 1D array */
  for(iy=iyMin;iy<=iyMax;++iy) 
    {
      Zy = GiveZy(iy);
      for(ix=ixMin;ix<=ixMax;++ix) 
	{

	  // from screen to world coordinate 
	  Zx = GiveZx(ix);
	  i = Give_i(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;
}





// 
unsigned char GiveColor2(unsigned int ix, unsigned int iy)
{ 
  

  double Zx2, Zy2;
  int i=0;
  long int j=0;  
  double Zx, Zy, Zx0, Zy0;
  int flag=0;
  
  
  
  
  
  // from screen to world coordinate 
  Zx0 = GiveZx(ix);
  Zy0 = GiveZy(iy);
  Zx= Zx0;
  Zy = Zy0;
  
  
  flag = IsInsideWideTriangle( Zx, Zy);
  if (flag ) 
   { if (flag==1) return iColorsOfInterior[iPeriodChild-1];
     if (flag==2) return iColorsOfInterior[iPeriodChild-1]+11;}
   
    
    

  // if not inside target set around attractor ( alfa fixed point )
   while (!flag && j<iterMax)
    { // then iterate 
     

    
      for(i=0;i<iPeriodChild ;++i) // iMax = period !!!!
	{  
	  Zx2 = Zx*Zx; 
	  Zy2 = Zy*Zy;
       
	  // bailout test 
	  if (Zx2 + Zy2 > ER2) return iColorOfExterior; // if escaping stop iteration
       
	  // if not escaping or not attracting then iterate = check behaviour
	  // new z : Z(n+1) = Zn * Zn  + C
	  Zy = 2*Zx*Zy + Cy; 
	  Zx = Zx2 - Zy2 + Cx; 
	  //
	  flag = IsInsideWideTriangle( Zx, Zy);
          if (flag ) 
            { if (flag==1) return iColorsOfInterior[i];
              if (flag==2) return iColorsOfInterior[i]+11;}
        j+=1;
       }
      
       
              
      
    }

  
  
  return iColorOfUnknown; // it should never happen
}







int MakeInternalTilingS(unsigned char a[] )
{
  unsigned int ix, iy; // pixel coordinate 

  unsigned char Color; // gray from 0 to 255 

  printf("axis of symmetry \n"); 
  iy = iyAxisOfSymmetry; 
#pragma omp parallel for schedule(dynamic) private(ix,Color) shared(ixMin,ixMax, iyAxisOfSymmetry)
  for(ix=ixMin;ix<=ixMax;++ix) {//printf(" %d from %d\n", ix, ixMax); //info  
    PlotPoint(ix, iy, GiveColor2(ix, iy),a);
  }


  /*
    The use of ‘shared(variable, variable2) specifies that these variables should be shared among all the threads.
    The use of ‘private(variable, variable2)’ specifies that these variables should have a seperate instance in each thread.
  */

#pragma omp parallel for schedule(dynamic) private(iy,ix,Color) shared(iyAboveMin, iyAboveMax,ixMin,ixMax, iyAxisOfSymmetry)

  // above and below axis 
  for(iy = iyAboveMin; iy<=iyAboveMax; ++iy) 
    {printf(" %d from %d\r", iy, iyAboveMax); //info 
      for(ix=ixMin; ix<=ixMax; ++ix) 

	{ // above axis compute color and save it to the array
	  iy = iyAxisOfSymmetry + iy;
	  Color = GiveColor2(ix, iy);
	  PlotPoint(ix, iy, Color, a ); 
	  // below the axis only copy Color the same as above without computing it 
	  PlotPoint(ixMax-ix, iyAxisOfSymmetry - iy , Color, a ); 
	} 
    }  
  return 0;
}


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


  // for all pixels of image 
  for(iy = iyMin; iy<=iyMax; ++iy) 
    { printf(" %d z %d\r", iy, iyMax); //info 
      for(ix= ixMin; ix<=ixMax; ++ix) 
           PlotPoint(ix, iy, GiveColor2(ix, iy),a ); //  
    } 
   
  return 0;
}


int DrawCriticalOrbit(unsigned char A[], unsigned int IterMax)
{
 
  unsigned int ix, iy; // pixel coordinate 
  // initial point z0 = critical point
  double Zx=0.0; 
  double Zy=0.0; //  Z= Zx+ZY*i;
  double Zx2=0.0;
  double Zy2=0.0;
  unsigned int i; /* index of 1D array */
  unsigned int j;


  // draw critical point  
  //compute integer coordinate  
  if ( Zx<=ZxMax && Zx>=ZxMin && Zy>=ZyMin && Zy<=ZyMax ){
  ix = (int)((Zx-ZxMin)/PixelWidth);
  iy = (int)((ZyMax-Zy)/PixelHeight); // reverse y axis
  i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */
  if (i>=0 && i<iSize) A[i]=255-A[i]; }

  // iterate
  for (j = 1; j <= IterMax; j++) //larg number of iteration s
    {  Zx2 = Zx*Zx; 
      Zy2 = Zy*Zy;
       
      // bailout test 
      if (Zx2 + Zy2 > ER2) return 1; // if escaping stop iteration and return error code
       
      // if not escaping iterate
      // Z(n+1) = Zn * Zn  + C
      Zy = 2*Zx*Zy + Cy; 
      Zx = Zx2 - Zy2 + Cx;
      //compute integer coordinate  
      if ( Zx<=ZxMax && Zx>=ZxMin && Zy>=ZyMin && Zy<=ZyMax )
      {ix = (int)round((Zx-ZxMin)/PixelWidth);
      iy = (int)round((ZyMax-Zy)/PixelHeight); // reverse y axis
      i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */
      if (i>=0 && i<iSize) A[i]=255-A[i];   // mark the critical orbit
      }
    }
  return 0;
}







// mark circle with center = (Zax, Zay ) and radius = sqrt(radius2)
int MarkCircle(unsigned char a[] )
{
  unsigned int ix, iy; // pixel coordinate 
  double Zx, Zy; //  Z= Zx+ZY*i;
  unsigned i; /* index of 1D array */
  double distance2;


  for(iy=iyMin;iy<=iyMax;++iy) 
    {
      Zy = GiveZy(iy);
      for(ix=ixMin;ix<=ixMax;++ix) 
	{

	  // from screen to world coordinate 
	  Zx = GiveZx(ix);
	  i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */
	  // 
          distance2 = (Zx-Zax)*(Zx-Zax) + (Zy-Zay)*(Zy-Zay);


          if (distance2 < radius2 )
          { if (a[i] == iColorOfUnknown)
              //printf("mark unknown ix= %d, iy = %d ; Z  = %.16f ; %.16f ; Distance2Fixed= %f \n",ix,iy, Zx, Zy, sqrt(distance2)); 
           a[i]= 255-a[i]; // mark circle 
          }
                        
	}
    }
   
  return 0;
}


int MarkWideTriangle(unsigned char a[] )
{
  unsigned int ix, iy; // pixel coordinate 
  double Zx, Zy; //  Z= Zx+ZY*i;
  unsigned i; /* index of 1D array */
  


  for(iy=iyMin;iy<=iyMax;++iy) 
    {
      Zy = GiveZy(iy);
      for(ix=ixMin;ix<=ixMax;++ix) 
	{

	  // from screen to world coordinate 
	  Zx = GiveZx(ix);
	  
	   if (IsInsideWideTriangle( Zx, Zy))
         
           {   i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */
               a[i]= 255 - a[i] ;
              }
          
	}
    }
   
  return 0;
}






// save data 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, "%.0f", 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 image data bytes to the file in one step */
  
  //
  printf("File %s saved. ", filename);
  if (comment == NULL)  printf ("empty comment \n");
                   else printf (" comment = %s \n", comment); 
  fclose(fp);

  return 0;
}











int info()
{
  // diplay info messages
  printf("Numerical approximation of parabolic Julia set for complex quadratic polynomial fc(z)= z^2 + c \n");
  printf("parameter c  = %.16f , %.16f \n", Cx, Cy);
  printf("c is a root point between hyperbolic components of period %d and %d  of Mandelbrot set \n", iPeriodParent,  iPeriodChild);
  printf("combinatorial rotation number = internal angle  = %d / %d = %f\n",iPeriodParent,  iPeriodChild, InternalAngle);
  printf("\n");
  printf("image size in pixels = %d x %d \n", iWidth, iHeight);
  printf("image size in world units : (ZxMin = %f, ZxMax =  %f) ,  (ZyMin = %f, ZyMax =  %f) \n", ZxMin , ZxMax ,  ZyMin , ZyMax);
  printf("ratio ( distortion) of image  = %f ; it should be 1.000 ...\n", ratio);
  printf("PixelWidth = %f \n", PixelWidth);
  printf("\n");
  printf("parabolic alfa fixed point Za  = %.16f ; %.16f \n", Zax, Zay);
  printf("critical point Zcr  = %.16f ; %.16f \n", Zcrx, Zcry);
  printf("precritical point Zl  = %.16f ; %.16f \n", Zlx, Zly);
  printf("precritical point Zr  = %.16f ; %.16f \n", Zrx, Zry);
  printf("Maximal number of iterations = iterMax = %ld \n", iterMax);
  printf("quality = %.16f = iNumberOfUknknown/ iNumberOfAllPixels  = %d / %d ; It should be 1.0  !!!! \n",quality,  iNumberOfUnknown, iSize);
  printf("MaxDistance From Unknown 2 Fixed = %.16f [world units]= %d [pixels] \n", sqrt(MaxDistanceFromU2F2), (int)round(sqrt(MaxDistanceFromU2F2)/PixelWidth));
  return 0;
}




   


/* -----------------------------------------  main   -------------------------------------------------------------*/
int main()
{
  

  // bounds checking
  if (iPeriodChild > log2(ULONG_MAX) ) 
     {printf("possible integer overflow, change type of numerator and denominator\n"); return 1;}


  setup();


  // here are procedures for creating image files
  
  //FillArray(data0); 
  FillArraySymmetric(data0); // some problems 
  SaveArray2PGMFile(data0 , iterMax +0, "components of filled Julia set - first aproximation"); // save array (image) to pgm file 

  ComputeBoundaries(data0, edge);
  SaveArray2PGMFile(edge , iterMax +1, "only boundaries"); // save array (image) to pgm file 
  

  CopyBoundaries(edge,data0);
  SaveArray2PGMFile(data0 , iterMax +2, "with boundaries "); // save array (image) to pgm file 
 //
  MarkWideTriangle(data0);
  SaveArray2PGMFile(data0 , iterMax +3, "marked wide triangle"); // save array (image) to pgm file 
  
 

  

  
  CheckOrientation(data0);
  SaveArray2PGMFile(data0 , iterMax+9, "marked first quadrant should be up and right"); // save array (image) to pgm file




  free(data);
  free(edge);
  free(data0);
  info();
  
  return 0;
}