iter  =GiveLastIteration(Zx, Zy );
 
   if (iter< iterMax ) 
    { color = iExterior; } // if point escapes = exterior 
    else // Filled Julia Set = bounded orbit = not escapes

 dX=Zx-_ZAx;
 dY=Zy-_ZAy;
 d=dX*dX+dY*dY; // d = square of distance from zn to alfa

 while ( d>_AR2) // check if z is inside target 
    {
      for(i=0;i<period ;++i) // iMax = period !!!!

  for (i = 0; i < period; ++i)
    {  if ((Angles[i]<angle) && (angle<Angles[i+1])) return  i;}
  if ((angle<Angles[0]) || (Angles[iMax-1]< angle)) return (iMax-1);

 angles[j] = (double)i/iMax  - Offset; // this turn offset is computed with Maxima CAS file

 int NumberOfPetal =GiveNumberOfPetal(Zx,Zy, Cx, Cy, period, AR2, creal(alfa), cimag(alfa));
      color = Colors[NumberOfPetal];

/*

  Adam Majewski
  fraktal.republika.pl

  c console progam using 
  * symmetry
  * openMP

  draw  parabolic Julia set 

  and saves it to pgm file

  period TurnOffset
  5	0.17298227404701
  6	0.04
  7	0.092



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

  File big_0.092000.pgm saved. 
  InternalAngle  = 0.142857 
  Cx  = 0.367375 
  Cy  = 0.147184 
  alfax  = 0.311745 
  alfay  = 0.390916 
  iHeight  = 501 
  PixelWidth  = 0.006000 
  TargetRadiusInPixels  = 7.477612 
  AR  = 0.044866 
  TurnOffset  = 0.092000 
  TurnOffset/InternalAngle  = 0.644000 
  distorsion of image  = 1.000000 

  real	47m4.555s
  user	84m48.390s
  sys	0m6.240s

  File 0.053055.pgm saved. 
  InternalAngle  = 0.142857 
  Cx  = 0.367375 
  Cy  = 0.147184 
  alfax  = 0.311745 
  alfay  = 0.390916 
  iHeight  = 2001 
  PixelWidth  = 0.001500 
  TargetRadiusInPixels  = 29.865672 
  AR  = 0.044799 
  TurnOffset  = 0.053055 
  TurnOffset/InternalAngle  = 0.371382 
  distorsion of image  = 1.000000 

  real	681m8.312s


*/



#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


/* --------------------------------- global variables and constans ------------------------------------------------------------ */
//unsigned int period = 7; // period of secondary component joined by root point
#define period 7

//  unsigned int denominator;  denominator = period;
double InternalAngle;
unsigned char Colors[period]; //={255,230,180, 160,140,120,100}; // NumberOfPetal of colors >= period
unsigned char iExterior = 245;

// 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 iyAbove ; // var, measured from 1 to (iyAboveAxisLength -1)
unsigned int iyAboveMin = 1 ; //
unsigned int iyAboveMax ; //
unsigned int iyAboveAxisLength ; //
unsigned int iyBelowAxisLength ; //
unsigned int iHeight = 2000; //  odd NumberOfPetal !!!!!! = (iyMax -iyMin + 1) = iyAboveAxisLength + iyBelowAxisLength +1
// The size of array has to be a positive constant integer 
unsigned int iSize ; // = iWidth*iHeight; 


// memmory 1D arrays 
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 */
const double ZxMin=-1.5;
const double ZxMax=1.5;
const double ZyMin=-1.5;
const double ZyMax=1.5;
double PixelWidth; // =(ZxMax-ZxMin)/iXmax;
double PixelHeight; // =(ZyMax-ZyMin)/iYmax;
 
double ratio ;


// angles 
double TwoPi=2*M_PI;
// angles measured in turns 
double TurnOffset; 
double Angles[period]; // array of angles = repelling directions in target set 


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

double complex alfa; // alfa fixed point alfa=f(alfa)

unsigned long int iterMax  = 100; //iHeight*100;
// target set for escaping points is a exterior of circle with center in origin
double ER = 2.0; // Escape Radius for bailout test 
double ER2;
// target set for points falling into alfa fixed point
// is a circle around alfa fixed point
// with radius = AR
double AR ; // radius of target set around alfa fixed point in world coordinate AR = PixelWidth*TargetWidth; 
double AR2; // =AR*AR;
double TargetRadiusInPixels; // radius of target set in pixels ; 







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

*/
double complex GiveC(double InternalAngleInTurns, double InternalRadius, unsigned int iPeriod)
{
  //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 ( iPeriod ) // 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 period  there are 2^(period-1) roots. 
    default: // higher periods : to do
      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 // general form  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;
}



// computes angles of repelling directions
// and saves it to array
int InitAngles( double angles[], double Offset)
{
  // repelling directions create period-arm symmetric star
  // star is slightly turned so this offset 
  // theory : Complex Dynamics by  Lennart Carleson,Theodore Williams Gamelin  page 40
  // TurnOffset = carg(a)/q
  // where a = coeff( z^(q+1))

  //fprintf(stderr,"Initialized of array by assignment  \n");
  int i,j;
  int iMax = period; // uses global var 

  for (i = 0; i < iMax; ++i)
    { if (Offset<0.00001) // no offset 
	j=i;
      else {
	// ascending order of angles 
	if (i==0) 
	  j=iMax-1 ;
	else j=i-1;
      }
    
      angles[j] = (double)i/iMax  - Offset;
      if (angles[j]<0) angles[j] = angles[j] + 1.0; //   argument in turns from 0 to 1.0
      // printf("i= %d angle= %f  \n",i, angles[i]);
    }

  for (i = 0; i < iMax; ++i) printf("i= %d angle= %f  \n",i, angles[i]);
   
  return 0;
}


// colors of components interior = shades of gray
int InitColors()
{
  int i;
  int iMax = period; // uses global var period and Colors
  unsigned int iStep;

  iStep=150/period;

  for (i = 1; i <= iMax; ++i)
    {Colors[i-1] = iExterior -i*iStep; 
      printf("i= %d color = %i  \n",i-1, Colors[i-1]);}
  return 0;
}


/* -------------------------------------------------- SETUP --------------------------------------- */
int setup()
{
  /* 2D array ranges */
  if (!(iHeight % 2)) iHeight+=1; // it sholud be even NumberOfPetal (variable % 2) or (variable & 1)
  iWidth = iHeight;
  iSize = iWidth*iHeight; // size = NumberOfPetal 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].

  TargetRadiusInPixels = iHeight/67.0; // = 15.0; // radius of target set in pixels ; Maybe increase to 20 = 1000/50
  InternalAngle = 1.0/((double) period); // angle for 3/7 is different then angle for 2/7
  c = GiveC(InternalAngle, 1.0, 1) ;
  alfa = GiveAlfaFixedPoint(c);
  TurnOffset=0.053054566595992; // 1.0/(4.0*period); // find by guess and check method
  //it is smaller the 1/period so here = 1/6  = 0.1666666...
  // aproximation = 1.0/(3.0*period)

  
  Cx=creal(c);
  Cy=cimag(c);

  /* 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;
  AR = PixelWidth*TargetRadiusInPixels; // !!!! important value  ( and precision and time of creation of the pgm image ) 
  AR2= AR*AR;
   
  /* create dynamic 1D arrays for colors ( shades of gray ) */
  data = malloc( iSize * sizeof(unsigned char) );
  edge = malloc( iSize * sizeof(unsigned char) );
  if (edge == NULL || edge == NULL )
    {
      fprintf(stderr," Could not allocate memory\n");
      return 1;
    }
  else fprintf(stderr," memory is OK \n");

  InitAngles( Angles, TurnOffset);
  InitColors();
  
 
  return 0;

}






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

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


// forward iteration of complex quadratic polynomial
// fc(z) = z*z +c
// initial point is a z0 = critical point 
// checks if points escapes to infinity = exterior of filled julia set
// uses global vars : ER2, Cx, Cy, iterMax
long int GiveLastIteration(double Zx, double Zy )
{ 
  long int iter;
  // z= Zx + Zy*i
  double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
 
  /* initial value of orbit = critical point Z= 0 */
                       
  Zx2=Zx*Zx;
  Zy2=Zy*Zy;
  // for iter from 0 to iterMax because of ++ after last loop
  for (iter=0; iter<iterMax && ((Zx2+Zy2)<ER2); ++iter)
    {
      Zy=2*Zx*Zy + Cy;
      Zx=Zx2-Zy2 + Cx;
      Zx2=Zx*Zx;
      Zy2=Zy*Zy;
    };
  
  return iter;
}




double GiveTurn(double complex z)
{
  double argument;

  argument = carg(z); //   argument in radians from -pi to pi
  if (argument<0) argument=argument + TwoPi; //   argument in radians from 0 to 2*pi
  return argument/TwoPi ; // argument in turns from 0.0 to 1.0
}



/* 
   Checks to which petal ( = part of target set)  point of interior falls 
   First bail-out test should be done ( = check if points escapes )
   used for coloring interior of Julia set
   and/or finding Julia set
*/
int GiveNumberOfPetal(double _Zx0, double _Zy0,double C_x, double C_y, int iMax, double _AR2, double _ZAx, double _ZAy )
{ 
  int i;
  double Zx, Zy; /* z = zx+zy*i */
  double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
  double d, dX, dY; /* distance from z to Alpha  */
  double angle;
 
  Zx=_Zx0; /* initial value of orbit  */
  Zy=_Zy0;
  
  Zx2=Zx*Zx;
  Zy2=Zy*Zy;
  dX=Zx-_ZAx;
  dY=Zy-_ZAy;
  d=dX*dX+dY*dY; // distance to alfa fixed point
 
  // iterate until point z will fall into taget set 
  while ( d>_AR2) // while z is outside target set 
    {
      for(i=0;i<period ;++i) // iMax = period !!!!
	{ // ieration z(n+1) = fc(zn) = zn*zn + c
	  Zy=2*Zx*Zy + C_y;
	  Zx=Zx2-Zy2 +C_x;
	  Zx2=Zx*Zx;
	  Zy2=Zy*Zy;
	}



      dX=Zx-_ZAx;
      dY=Zy-_ZAy;
      d=dX*dX+dY*dY;
    };

  
  angle =GiveTurn(dX +dY*I) ; 

  // divide interior into components 
  // find to which petal ( sector of target set ) angle belong
  // angles are measured in turns 
  for (i = 0; i < period; ++i)
    {  if ((Angles[i]<angle) && (angle<Angles[i+1])) return  i;}
  if ((angle<Angles[0]) || (Angles[iMax-1]< angle)) return (iMax-1);

  return 0;

}


 



unsigned char GiveColor(unsigned int ix, unsigned int iy)
{ 
  double Zx, Zy; //  Z= Zx+ZY*i;
  int iter;
  unsigned char color; // shades of gray =  from 0 to 255 

  // from screen to world coordinate 
  Zx = GiveZx(ix);
  Zy = GiveZy(iy);
  
  iter  =GiveLastIteration(Zx, Zy );
  if (iter< iterMax ) // uses global var iterMax 
    { color = iExterior; } // if point escapes = exterior 
  else // Filled Julia Set = bounded orbit = not escapes
    {   // divide interior into components according to which petal it falls
      int NumberOfPetal =GiveNumberOfPetal(Zx,Zy, Cx, Cy, period, AR2, creal(alfa), cimag(alfa));
      color = Colors[NumberOfPetal]; 
    } 


  return color;
}


/* -----------  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 i; /* index of 1D array */
  i = Give_i(ix,iy); /* compute index of 1D array from indices of 2D array */
  data[i] = iColor;

  return 0;
}


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


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


// fill array using symmetry of image 
// uses global var :  ...
int FillArraySymmetric(unsigned char data[] )
{
   
  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));
  }


  /*
    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 ); 
	  // below the axis only copy Color the same as above without computing it 
	  PlotPoint(ixMax-ix, iyAxisOfSymmetry - iyAbove , Color ); 
	} 
    }  
  return 0;
}

int AddBoundaries(unsigned char data[])
{

  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; 
 
  


  printf(" find boundaries in data 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= data[Give_i(iX-1,iY+1)] + 2*data[Give_i(iX,iY+1)] + data[Give_i(iX-1,iY+1)] - data[Give_i(iX-1,iY-1)] - 2*data[Give_i(iX-1,iY)] - data[Give_i(iX+1,iY-1)];
      Gh= data[Give_i(iX+1,iY+1)] + 2*data[Give_i(iX+1,iY)] + data[Give_i(iX-1,iY-1)] - data[Give_i(iX+1,iY-1)] - 2*data[Give_i(iX-1,iY)] - data[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) {edge[i]=255;} /* background */
      else {edge[i]=0;}  /* boundary */
    }
  }

  // copy boundaries from edge array to data array 
  //for(iY=1;iY<iyMax-1;++iY){ 
  //   for(iX=1;iX<ixMax-1;++iX){i= Give_i(iX,iY); if (edge[i]==0) data[i]=0;}}



  return 0;
}


// Check Orientation of image : first quadrant in upper right position
// uses global var :  ...
int CheckOrientation(unsigned char data[] )
{
  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) data[i]=255-data[i];   // check the orientation of Z-plane by marking first quadrant */

	}
    }
   
  return 0;
}



// save data array to pgm file 
int SaveArray2PGMFile( unsigned char data[], double t)
{
  
  FILE * fp;
  const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
  char name [10]; /* name of file */
  sprintf(name,"%f", t); /*  */
  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(data,iSize,1,fp);  /*write image data bytes to the file in one step */
  printf("File %s saved. \n", filename);
  fclose(fp);

  return 0;
}


int info()
{
  // diplay info messages
  printf("InternalAngle  = %f \n", InternalAngle);
  printf("Cx  = %f \n", Cx); 
  printf("Cy  = %f \n", Cy);
  // 
  printf("alfax  = %f \n", creal(alfa));
  printf("alfay  = %f \n", cimag(alfa));
  printf("iHeight  = %d \n", iHeight);
  printf("PixelWidth  = %f \n", PixelWidth);
  printf("TargetRadiusInPixels  = %f \n", TargetRadiusInPixels);
  printf("AR  = %f \n", AR);
  printf("TurnOffset  = %f \n", TurnOffset);
  printf("TurnOffset/InternalAngle  = %f \n", TurnOffset/InternalAngle);
  printf("distorsion of image  = %f \n", ratio);
  return 0;
}




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

  
  setup(); 

  // here are procedures for creating image file
  
  //FillArray( data ); // no symmetry
  FillArraySymmetric(data); 
  AddBoundaries(data);
  //  CheckOrientation( data );
  SaveArray2PGMFile(edge , TurnOffset); // save edge array (only boundaries) to pgm file 

  //
  free(data);
  free(edge);
  
  //
  info();

  return 0;
}