/*

Adam Majewski
fraktal.republika.pl

c console progam using 
* symmetry
* openMP





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







*/



#include <stdio.h>
#include <stdlib.h> // malloc
#include <string.h> // strcat
#include <math.h> // M_PI; needs -lm also 
#include <stdbool.h>
#include <complex.h>
#include <omp.h> // OpenMP; needs also -fopenmp


/* --------------------------------- global variables and consts ------------------------------------------------------------ */
unsigned int period = 2; // period of secondary component joined by root point

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

// 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=-2.5;
  const double ZxMax=2.5;
  const double ZyMin=-2.5;
  const double ZyMax=2.5;
  double PixelWidth; // =(ZxMax-ZxMin)/ixMax;
  
  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 alfa; // alfa fixed point alfa=f(alfa)
double complex beta; // beta fixed point alfa=f(alfa)

unsigned long int iterMax  = 100; //iHeight*100;

double ER = 2.0; // Escape Radius for bailout test 
double ER2;
double AR,AR2; // AR2 = AR*AR where AR is a radius around attractor

/* colors = shades of gray from 0 to 255 */
unsigned char iExterior=245;
unsigned char iInteriorRightUp=231;
unsigned char iInteriorRightDown=99;
unsigned char iInteriorLeftUp=123;
unsigned char iInteriorLeftDown=255;
// {{255,231},{123,99}};
unsigned char iBoundary=0;

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

double complex GiveBetaFixedPoint(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);
 // beta = 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;
}


// distance2 = distance*distance
double GiveDistance2Between(double complex z1, double z2x, double z2y )
{double dx,dy;
 
 dx = creal(z1) - z2x;
 dy = cimag(z1) - z2y;
 return (dx*dx+dy*dy);
 
} 



int setup()
{

  
  unsigned int denominator;
  double InternalAngle;
  

  denominator = period;
  InternalAngle = 1.0/((double) denominator);

  c = GiveC(InternalAngle, 0, period) ; // internal radius= o gives center of component 
  Cx=creal(c);
  Cy=cimag(c);
  alfa = GiveAlfaFixedPoint(c);
  beta = GiveBetaFixedPoint(c);

  /* 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 in iteration
  ER2 = ER * ER;
  AR =  PixelHeight; // radius of the target set around fixed attractor 
  AR2 = AR*AR; 
  
    
  
/* create dynamic 1D arrays for colors ( shades of gray ) */
  
  data = malloc( iSize * sizeof(unsigned char) );
   if (data == NULL )
    {
      fprintf(stderr," Could not allocate memory");
      getchar(); 
      return 1;
    }
   
  
 
  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


// forward iteration of complex quadratic polynomial
// fc(z) = z*z +c
// z0 = initial point 
// uses global var 
/*
 Main loop : forward iteration of initial point 


 */
unsigned char dGiveColor(double Zx, double Zy, double Cx, double Cy ,  int iter_max)
 { 
 int i;
 double  x = Zx, /* Z = x+y*i */
         y = Zy, 
         x2,
         y2 ; 

    x2 = x*x;
    y2 = y*y;
    if (x2+y2 > ER2) return iExterior; // escapes = exterior
 
 for (i = 1; i <= iter_max; i++)
  { 
    /* z = z*z + c = x+y*i */
    
    y = 2*x*y + Cy; 
    x = x2 - y2 + Cx; 
    x2 = x*x;
    y2 = y*y;
    if (GiveDistance2Between(alfa,x,y)<AR2) break; // interior
    if (x2+y2 > ER2) return iExterior; // escapes = exterior
  } // for   
    
  // if not escapes then z is in a filled Julia set 
  // interior color :  tiling
 if (x>creal(alfa))
      {if (y>cimag(alfa))   return iInteriorRightUp;
                       else return iInteriorRightDown;}
 else /* x<=creal(alfa) */ 
    if (y>cimag(alfa))  return iInteriorLeftUp;
 // last case 
 return iInteriorLeftDown; 
 }


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

  // from screen to world coordinate 
  Zx = GiveZx(ix);
  Zy = GiveZy(iy);
  
  color = dGiveColor(Zx, Zy,  Cx,  Cy ,iterMax);


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




// plots raster point (ix,iy) 
int iDrawPoint(unsigned int ix, unsigned int iy, unsigned char iColor, unsigned char A[])
{ 

 /* i =  Give_i(ix,iy) compute index of 1D array from indices of 2D array */
 A[Give_i(ix,iy)] = iColor;

return 0;
}


// draws point to memmory array data
// uses complex type so #include <complex.h> and -lm 
int dDrawPoint(complex double point,unsigned char iColor, unsigned char A[] )
{

  unsigned int ix, iy; // screen coordinate = indices of virtual 2D array
  //unsigned int i; // index of 1D array
  
  ix = (creal(point)- ZxMin)/PixelWidth; 
  iy = (ZyMax - cimag(point))/PixelHeight; // inverse Y axis 
  iDrawPoint(ix, iy, iColor, A);
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;

  for(;;){
    iDrawPoint(x, y, iColor, A);
    if (x==x1 && y==y1) break;
    e2 = err;
    if (e2 >-dx) { err -= dy; x += sx; }
    if (e2 < dy) { err += dx; y += sy; }
  }
}




int dDrawLineSegment(double Zx0, double Zy0, double Zx1, double Zy1, int color, unsigned char *array) 
{

 unsigned int ix0, iy0; // screen coordinate = indices of virtual 2D array 
 unsigned 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 >=ixMin && ix0<=ixMax && ix0 >=ixMin && ix0<=ixMax && iy0 >=iyMin && iy0<=iyMax 
      && iy1 >=iyMin && iy1<=iyMax )
   iDrawLine(ix0,iy0,ix1,iy1,color, array) ;

return 0;
}







// -------------------------------
// https://en.wikipedia.org/wiki/Cohen%E2%80%93Sutherland_algorithm

typedef int OutCode;

const int INSIDE = 0; // 0000
const int LEFT = 1;   // 0001
const int RIGHT = 2;  // 0010
const int BOTTOM = 4; // 0100
const int TOP = 8;    // 1000

// Compute the bit code for a point (x, y) using the clip rectangle
// bounded diagonally by (ZxMin, ZyMin), and (ZxMax, ZyMax)

// ASSUME THAT ZxMax, ZxMin, ZyMax and ZyMin are global constants.

OutCode ComputeOutCode(double x, double y)
{
	OutCode code;

	code = INSIDE;          // initialised as being inside of [[clip window]]

	if (x < ZxMin)           // to the left of clip window
		code |= LEFT;
	else if (x > ZxMax)      // to the right of clip window
		code |= RIGHT;
	if (y < ZyMin)           // below the clip window
		code |= BOTTOM;
	else if (y > ZyMax)      // above the clip window
		code |= TOP;

	return code;
}

// Cohen–Sutherland clipping algorithm clips a line from
// P0 = (x0, y0) to P1 = (x1, y1) against a rectangle with 
// diagonal from (xmin, ymin) to (xmax, ymax).
// CohenSutherlandLineClipAndDraw
void dDrawLine(double x0, double y0, double x1, double y1,unsigned char iColor, unsigned char A[])
{
	// compute outcodes for P0, P1, and whatever point lies outside the clip rectangle
	OutCode outcode0 = ComputeOutCode(x0, y0);
	OutCode outcode1 = ComputeOutCode(x1, y1);
	bool accept = false;

	while (true) {
		if (!(outcode0 | outcode1)) { // Bitwise OR is 0. Trivially accept and get out of loop
			accept = true;
			break;
		} else if (outcode0 & outcode1) { // Bitwise AND is not 0. (implies both end points are in the same region outside the window). Reject and get out of loop
			break;
		} else {
			// failed both tests, so calculate the line segment to clip
			// from an outside point to an intersection with clip edge
			double x, y;

			// At least one endpoint is outside the clip rectangle; pick it.
			OutCode outcodeOut = outcode0 ? outcode0 : outcode1;

			// Now find the intersection point;
			// use formulas y = y0 + slope * (x - x0), x = x0 + (1 / slope) * (y - y0)
			if (outcodeOut & TOP) {           // point is above the clip rectangle
				x = x0 + (x1 - x0) * (ZyMax - y0) / (y1 - y0);
				y = ZyMax;
			} else if (outcodeOut & BOTTOM) { // point is below the clip rectangle
				x = x0 + (x1 - x0) * (ZyMin - y0) / (y1 - y0);
				y = ZyMin;
			} else if (outcodeOut & RIGHT) {  // point is to the right of clip rectangle
				y = y0 + (y1 - y0) * (ZxMax - x0) / (x1 - x0);
				x = ZxMax;
			} else if (outcodeOut & LEFT) {   // point is to the left of clip rectangle
				y = y0 + (y1 - y0) * (ZxMin - x0) / (x1 - x0);
				x = ZxMin;
			}

			// Now we move outside point to intersection point to clip
			// and get ready for next pass.
			if (outcodeOut == outcode0) {
				x0 = x;
				y0 = y;
				outcode0 = ComputeOutCode(x0, y0);
			} else {
				x1 = x;
				y1 = y;
				outcode1 = ComputeOutCode(x1, y1);
			}
		}
	}
	if (accept) {
			
               // printf(	"x0 = %d, y0 = %d, x1 = %d, y1 =%d \n",x0, y0, x1, y1);			
		dDrawLineSegment(x0, y0, x1, y1, iColor, A);
	}
}

// -----------------------------------------------------------------------






// 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) iDrawPoint(ix, iy, GiveColor(ix, iy), data ); //  
   } 
   
 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  
                              iDrawPoint(ix, iy, GiveColor(ix, iy), data);
}


/*
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.
*/
printf("symmetric parts :\n"); 
#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);
    iDrawPoint(ix, iy, Color, data ); 
    // below the axis only copy Color the same as above without computing it 
    iDrawPoint(ixMax-ix, iyAxisOfSymmetry - iyAbove , Color, data ); 
   } 
}  
printf("\ndone\n"); 
 return 0;
}








int AddEdges(unsigned char data[] )
{
  // memmory 1D array 
  unsigned char *edge;

  /* sobel filter */
  unsigned char G, Gh, Gv; 
  unsigned int i; /* index of 1D array  */
  printf("find boundaries in data array using  Sobel filter : ");  

 /* create dynamic 1D arrays for colors ( shades of gray ) */
    edge = malloc( iSize * sizeof(unsigned char) );
  if (edge==NULL)
    {
      fprintf(stderr," Could not allocate memory for the edge array.Stop the program. \n");
      return 1;
    }
    
 
  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 */
    }
  }
 
    //printf(" copy boundaries from edge to data array \n");
    for(iy=1;iy<iyMax-1;++iy){ 
     for(ix=1;ix<ixMax-1;++ix)
      {i= Give_i(ix,iy); /* compute index of 1D array from indices of 2D array */
    if (edge[i]==0) data[i]=0;}}


   free(edge);
   printf(" done\n"); 

 return 0;
}


// Check Orientation of image : mark first quadrant 
// it should be in the 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;
}






/* 
 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 complex z)
{ 
  double x = creal(z);
  double y = cimag(z);
  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;
}

double complex preimage(double complex z1, double complex z2,  double complex c)
{ 
  double complex zPrev;

  zPrev = root(creal(z1) - creal(c) + ( cimag(z1) - cimag(c))*I);
  // choose one of 2 roots 
  if (creal(zPrev)*creal(z2) + cimag(zPrev)*cimag(z2) > 0) 
        return zPrev ;     //  u+v*i
       else return -zPrev; // -u-v*i
}



// This function only works for periodic or preperiodic angles.
// You must determine the period n and the preperiod k before calling this function.
// draws all "period" external rays 

double complex DrawRay(double t0, // external angle in turns 
		       int n, //period of ray's angle under doubling map
		       int k, // preperiod
                       int iterMax,
                       double complex c
                          )
{
  double xNew; // new point of the ray
  double yNew;
  double xend ; // re of the endpoint of the ray
  double yend; // im of the endpoint of the ray
  const double R = 4000; // very big radius = near infinity
  int j; // number of ray 
  int iter; // index of backward iteration
  double t=t0;
  
  
  double complex zPrev;
  double u,v; // zPrev = u+v*I
 // double complex zNext;


printf(" preperiod = %d \n" , k);


   /* dynamic 1D arrays for coordinates ( x, y) of points with the same R on preperiodic and periodic rays  */
  double *RayXs, *RayYs;
  int iLength = n+k+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
    }
  


 //  starting points on preperiodic and periodic rays 
 //  with angles t, 2t, 4t...  and the same radius R
  for (j = 0; j < n + k; j++)
  { // z= R*exp(2*Pi*t)
    RayXs[j] = R*cos((2*M_PI)*t); 
    RayYs[j] = R*sin((2*M_PI)*t);
    
    t *= 2; // t = 2*t
    if (t > 1) 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+k] = RayXs[k]; 
  RayYs[n+k] = RayYs[k];
  

  //   backward iteration of each point z
  for (iter = -10; iter <= iterMax; iter++)
    { 
     	
        for (j = 0; j < n+k; j++) // period +preperiod
         { // u+v*i = sqrt(z-c)   backward iteration in fc plane 
		zPrev = root(RayXs[j+1] - creal(c)+(RayYs[j+1] - cimag(c))*I ); // , 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
                //
                dDrawLine(xNew, yNew, RayXs[j], RayYs[j], 0, data);  
                RayXs[j] = xNew; RayYs[j] = yNew;
                
         } // for j ...

          //RayYs[n+k] cannot be constructed as a preimage of RayYs[n+k+1]
          RayXs[n+k] = RayXs[k]; 
          RayYs[n+k] = RayYs[k];
          
          // 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));
     
  }

 
 // last point of a ray 0
 xend = RayXs[0];
 yend = RayYs[0];
  printf("landing point of ray for angle = %f is = %f ; %f \n",t0, RayXs[0], RayYs[0]);
 


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

return  xend + yend*I; // return last point or ray for angle t 
}










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

  /* save image to the pgm file  */      
  fp= fopen(filename,"wb"); /*create new file,give it a name and open it in binary mode  */
  fprintf(fp,"P5\n %s\n %u %u\n %u\n",comment,iWidth,iHeight,MaxColorComponentValue);  /*write header to the file*/
  fwrite(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("Cx  = %f \n", Cx); 
  printf("Cy  = %f \n", Cy);
  // 
 
  printf("alfax  = %f \n", creal(alfa));
  printf("alfay  = %f \n", cimag(alfa));
  printf("betax  = %f \n", creal(beta));
  printf("betay  = %f \n", cimag(beta));
  printf("target set around fixed attractor has radius AR  = %f  = %f pixels wide \n", AR, AR/PixelWidth);
  printf("ratio of image  = %f ; it should be 1.000 ...\n", ratio);
return 0;
}

double GiveAngleInTurns(int numerator, int denominator)
{ 
printf("draw ray for angle = %d / %d ;  " , numerator,  denominator);
return ((double)(numerator % denominator))/((double)denominator);}


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

  int preperiod=0; // preperiod
  // external angle of dynamic ray in turns 
  int n=1; // numerator of angle
  // denominator of angle 
  int d=( (int)pow(2.0,period) -1) * (int)pow(2.0,preperiod); // http://fraktal.republika.pl/mset_external_ray_m.html
  
  setup();

  // here are procedures for creating image file
  // compute colors of pixels = image
  //FillArray( data ); // no symmetry
  FillArraySymmetric(data); 
  AddEdges(data);
  //  CheckOrientation( data );


  // external rays for fixed point and its preimages
  // preperiod k=0 ; period 2 = two rays
  DrawRay(GiveAngleInTurns(n,d) , period, preperiod, 1000,c);

  for (preperiod=1;preperiod<5; ++preperiod)
{

   // second new angle = (angle+1)/2
   DrawRay(GiveAngleInTurns(n+d,d*2) , period, preperiod, 1000,c);
   DrawRay(GiveAngleInTurns((d-n)+d,d*2) , period, preperiod, 1000,c);
   // first  new angle = angle/2
   d*=2;
   DrawRay(GiveAngleInTurns(n,d) , period, preperiod, 1000,c);
   DrawRay(GiveAngleInTurns(d-n,d) ,period , preperiod, 1000,c);
   
 } 
  
  //preperiod 5 ; width = 2/24
 DrawRay(GiveAngleInTurns(11,24) ,period , 5, 1000,c);
 DrawRay(GiveAngleInTurns(13,24) ,period , 5, 1000,c);
 //
 DrawRay(GiveAngleInTurns(7,24) ,period , 5, 1000,c);
 DrawRay(GiveAngleInTurns(5,24) ,period , 5, 1000,c);
 //
 DrawRay(GiveAngleInTurns(17,24) ,period , 5, 1000,c);
 DrawRay(GiveAngleInTurns(19,24) ,period , 5, 1000,c);

 // preperiod 6 ; width = 2/48
 DrawRay(GiveAngleInTurns(5,48) ,period , 6, 1000,c);
 DrawRay(GiveAngleInTurns(7,48) ,period , 6, 1000,c);

 //
 DrawRay(GiveAngleInTurns(43,48) ,period , 6, 1000,c);
 DrawRay(GiveAngleInTurns(41,48) ,period , 6, 1000,c);
  //
 DrawRay(GiveAngleInTurns(29,48) ,period , 6, 1000,c);
 DrawRay(GiveAngleInTurns(31,48) ,period , 6, 1000,c);
 //
 DrawRay(GiveAngleInTurns(17,48) ,period , 6, 1000,c);
 DrawRay(GiveAngleInTurns(19,48) ,period , 6, 1000,c);
 //
 DrawRay(GiveAngleInTurns(23,48) ,period , 6, 1000,c);

 //preperiod 7 ; width = 2/96
 DrawRay(GiveAngleInTurns(5,96) ,period , 7, 1000,c);
 DrawRay(GiveAngleInTurns(7,96) ,period , 7, 1000,c);
 //
 DrawRay(GiveAngleInTurns(53,96) ,period , 7, 1000,c);
 DrawRay(GiveAngleInTurns(55,96) ,period , 7, 1000,c);
 //
 DrawRay(GiveAngleInTurns(43,96) ,period , 7, 1000,c);
 DrawRay(GiveAngleInTurns(41,96) ,period , 7, 1000,c);
 //
 DrawRay(GiveAngleInTurns(91,96) ,period , 7, 1000,c);
 DrawRay(GiveAngleInTurns(89,96) ,period , 7, 1000,c);
  
 // 
 preperiod = 14;
  SaveArray2PGMFile( data,preperiod-1); // save array (image) to pgm file 
  free(data);
  info();

  return 0;
}