// dot product 
double dot(complex double a, complex double b){

	// ax*bx + ay*by
	return creal(a)*creal(b) + cimag(a)*cimag(b);
}


int InverseIterationAndDrawSegment (int per, complex double zz[per], unsigned char A[]){
	
	int p;
	int pMax = per;
	complex double z;
	
	
	// compute preimages
	for(p=0; p<pMax; ++p){
		z = zz[(p+1) % per]; //read point from the ray p+1 
		z = csqrt(z-c); // inverse iteration = preimage;  fc(z) maps  z_{l,j} to z_{l-1,j+1}} so f^{-1} map to z_{l+1, j-1}
		
		// choose correct preimage using dot product 
		if (dot(zz[p],z) < 0 ) z = - z; 
		// draw segment of ray p
		DrawClippedLineSegment(zz[p],z, 255, A);
		zz[p] = z; //save 
	}
	
  	return 0; 

}




int DrawExternalDynamicRaysBI (int n, int m, int period, int iMax, unsigned char A[])
{ 
  	int i = 0; // iteration = number of points
  	double r = 10000.0; // very big radius = near infinity where  z=w so one can swith to dynamical plane ( Boettcher conjugation )
  	complex double zz[period]; // zz is an array of z  
  
  	double t;
  	int p;
  	int pMax = period; // number of rays to draw 
  
   
  	// initial points  on pMax rays
  	t = (double)n/m; // first external angle in turns	
  	
  	for(p=0; p<pMax; ++p){
  		
  		zz[p] =  r*cexp(2.0*I * M_PI * t ); // Euler's formula gives initial point on the p ray
  		t *= 2.0; // t = 2*t  = angle doubling map ; gives next angle 
      	  	}
  
	for (i = 0; i < iMax; ++i)
    		 
    		InverseIterationAndDrawSegment(period, zz, A); // inverse iteration of  complex quadratic polynomial:  with proper choose of preimage 
      	return 0;
}