Jump to content

File:Julia set p(z)= z^3+(1.0149042485835864102+0.10183008497976470119i)*z; (zoom).png

Page contents not supported in other languages.
This is a file from the Wikimedia Commons
From Wikibooks, open books for an open world

Original file (2,000 × 2,000 pixels, file size: 1.65 MB, MIME type: image/png)

Summary

Description
English: Julia set p(z)= z^3+(1.0149042485835864102+0.10183008497976470119i)*z. Zoom and critical orbit. Location by Michael Becker[1] Image made with LCM ( level curves for interior and exterior). Level curves cross at critical point and its preimages. One can see spiral from attracting fixed point to repelling fixed point ( z= 0) which is a place with high density of level curves. Points of critical orbit ( including crirital point and attractor) are on the level curves like notes on the musical staff
Deutsch: f(z)=z3+dz mit d=1,02*e0,1i, dargestellt auf [-0.5;0.5]x[-0.5;0.5]. The Julia set (boundary of filled-n Julia set) itself is not drawn: we see it as the locus of points where the level curves (= the boundaries of level sets) are especially close to each other = a place with high density of level curves.
Date
Source Own work
Author Adam majewski
Other versions

Licensing

I, the copyright holder of this work, hereby publish it under the following license:
w:en:Creative Commons
attribution share alike
This file is licensed under the Creative Commons Attribution-Share Alike 4.0 International license.
You are free:
  • to share – to copy, distribute and transmit the work
  • to remix – to adapt the work
Under the following conditions:
  • attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
  • share alike – If you remix, transform, or build upon the material, you must distribute your contributions under the same or compatible license as the original.


c source code

/*

  Adam Majewski
  adammaj1 aaattt o2 dot pl  // o like oxygen not 0 like zero 
  
  
  console program in c programing language 
  ===============================================================

coefficients read from input file ijon_b011.txt
	degree 3 coefficient = ( +1.0000000000000000 +0.0000000000000000*i) 
	degree 1 coefficient = ( +1.0149042485835864 +0.1018300849797647*i) 

Input polynomial p(z)=(1+0i)*z^3+(1.0149042485835864102+0.10183008497976470119i)*z^1

derivative dp/dz = (3+0i)*z^2+(1.0149042485835864102+0.10183008497976470119i)

2 critical points found

	cp#0: 0.029142613176165576422,-0.58236647233272664792 . It's critical orbit is bounded and enters cycle #0 length=1 and it's stability = |multiplier|=0.99134 =parabolic 
	internal angle = 0.96706874577169499307
cycle = {
0.20977461555364212975,-0.24271307734496930242 ; }

	cp#1: -0.029142613176165569483,0.58236647233272664792 . It's critical orbit is bounded and enters cycle #1 length=1 and it's stability = |multiplier|=0.99134 =parabolic 
	internal angle = 0.96706874577169499307
cycle = {
-0.20977461555364212975,0.24271307734496930242 ; }



  
  ==============================================
  
  
  Structure of a program or how to analyze the program 
  
  
  ============== Image of (FunctionType, PlaneInversion) ========================
  
  DrawImageOf(X,y) -> DrawPoint(X,y) -> ComputeColorOf(X,y)
  
  so check  only compute color  function 
  
  
  

   
  ==========================================

  
  ---------------------------------
  indent d.c 
  default is gnu style 
  -------------------



  c console progam 
  
  export  OMP_DISPLAY_ENV="TRUE"	
  gcc d.c -lm -Wall -march=native -fopenmp
  time ./a.out > b.txt


  gcc d.c -lm -Wall -march=native -fopenmp


  time ./a.out

  time ./a.out >a.txt
  
  
  ./g.sh
  
  make

  ----------------------
  
  real	0m19,809s
  user	2m26,763s
  sys	0m0,161s


  

*/

#include <stdio.h>
#include <stdlib.h>		// malloc
#include <string.h>		// strcat
#include <math.h>		// M_PI; needs -lm also
#include <complex.h> 		// complex numbers : https://stackoverflow.com/questions/6418807/how-to-work-with-complex-numbers-in-c
#include <omp.h>		// OpenMP




/* --------------------------------- global variables and consts ------------------------------------------------------------ */


//  
double complex c;		// parameter of function fc(z)=z^3 + c*z
// attractors
double complex za0 =  0.20977461555364212975 -0.24271307734496930242 *I; // first attracting point of period 1 = cycle #0
double complex za1 = -0.20977461555364212975 +0.24271307734496930242 *I; // second attracting point of period 1 = cycle #1
// for bd
double za0im;
double za1im;


double complex zcr0 =  0.029142613176165576422-0.58236647233272664792 *I; // first critical = cycle #0
double complex zcr1 = -0.029142613176165569483+0.58236647233272664792 *I; // second critical point = cycle #1

















int NumberOfImages = 0;


//FunctionType
typedef enum  {Fatou, IntLSM, ExtLSM , LSM, DEM, Unknown, BD, MBD , SAC, DLD, ND, NP, POT, Blend
		
} FunctionTypeT; 
// FunctionTypeT FunctionType;

int PlaneInversion = 0; // boolean 1 = w = 1/z  plane; 0 = z plane 


// virtual 2D array and integer ( screen) coordinate
// Indexes of array starts from 0 not 1 
//unsigned int ix, iy; // var
static unsigned int ixMin = 0;	// Indexes of array starts from 0 not 1
static unsigned int ixMax;	//
static unsigned int iWidth;	// horizontal dimension of array

static unsigned int iyMin = 0;	// Indexes of array starts from 0 not 1
static unsigned int iyMax;	//

static unsigned int iHeight = 20000;	//  
// The size of array has to be a positive constant integer 
static unsigned int iSize;	// = iWidth*iHeight; 

// ----------memmory 1D arrays ==================
// unsigned char = for 1 byte ( 8 bit) colors 
unsigned char *data;
unsigned char *edge;
unsigned char *edge2;
 

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





// see SetZPlane

double radius = 0.5; 
complex double center = 0.0;
double  DisplayAspectRatio  = 1.0; // https://en.wikipedia.org/wiki/Aspect_ratio_(image)



// z plane = dynamic plane
double ZxMin ;	//-0.05;
double ZxMax ;	//0.75;
double ZyMin ;	//-0.1;
double ZyMax ;	//0.7;


double PixelWidth;	// =(ZxMax-ZxMin)/ixMax;
double PixelHeight;	// =(ZyMax-ZyMin)/iyMax;
double ratio;





// w plane = 1/z plane
double WxMin = -2000;	//-0.05;
double WxMax = 2000;	//0.75;
double WyMin = -2000;	//-0.1;
double WyMax = 2000;	//0.7;
double wPixelWidth;	// =(WxMax-WxMin)/ixMax;

double wPixelHeight;	// =(WyMax-WyMin)/iyMax;

static unsigned long int iterMax = 1000000;	//iHeight*100;
const long int iterMax_LSM = 1000;
const int iterMax_DLD = 200; // N in wiki = fixed number : maximal number of iterations
const int iterMax_pot = 400; // potential 
int iterMax_DEM = 1000;

double ER = 200.0;		// EscapeRadius for bailout test 
double EscapeRadius=1000000; // = ER big !!!!
double ER_LSM ; // see GiveER_LSM  // 27.764 =  manually find value such that level curves of escape time cross critical point and it's  preimages
double ER_DLD ; // see GiveER_LSM  // 27.764 =  manually find value such that level curves of escape time cross critical point and it's  preimages
double ER_NP = 100.0; 
double ER_pot = 100000.0;  // sqrt(1e24)
double ER_DEM = 1000.0;  // sqrt(1e24)

double loger; // = log(ER_LSM); // for texture
static double TwoPi=2.0*M_PI; // texture
double MaxFinalRadius;



double AR; // Attracting Radius = radius of circle around attractor ( trap)
double AR_max;

// SAC/J
double lnER; // ln(ER)
int i_skip = 2; // exclude (i_skip+1) elements from average
unsigned char s = 7; // stripe density

double BoundaryWidth = 3.0; // % of image width  
double distanceMax; //distanceMax = BoundaryWidth*PixelWidth;



//  ------------- DLD  ----------------------

double p = 0.180; //0.01444322;		//
// DLD colors
//double me = 1.0;
double mi = 0.9;

// potential
double MaxImagePotential;



/* colors = shades of gray from 0 to 255 */
unsigned char iColorOfExterior = 250;
unsigned char iColorOfInterior = 200;
unsigned char iColorOfInterior1 = 210;
unsigned char iColorOfInterior2 = 180;
unsigned char iColorOfBoundary = 0;
unsigned char iColorOfUnknown = 30;





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

/**
 * Find maximum between two numbers.
 https://codeforwin.org/2016/02/c-program-to-find-maximum-and-minimum-using-functions.html
*/
double max(double n1, double n2)
{
  return (n1 > n2 ) ? n1 : n2;
}



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

double min(double n1, double n2)
{
  return (n1 < n2 ) ? n1 : n2;
}


double clip(double d){

  return (d> 1.0) ? 1.0 : d;
}



double frac(double d){

  double fraction = d - ((long)d);
  return fraction;
}




//------------------complex numbers -----------------------------------------------------




double c_arg(complex double z)
{
  double arg;
  arg = carg(z);
  if (arg<0.0) arg+= TwoPi ; 
  return arg; 
}

double c_turn(complex double z)
{
  double arg;
  arg = c_arg(z);
  return arg/TwoPi; 
}





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

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


complex double GiveZ( int ix, int iy){
  double Zx = GiveZx(ix);
  double Zy = GiveZy(iy);
	
  return Zx + Zy*I;
	
	


}


// from screen to world coordinate ; linear mapping
// uses global cons
double GiveWx ( int ix)
{
  return (WxMin + ix * wPixelWidth);
}

// uses globaal cons
double GiveWy (int iy) {
  return (WyMax - iy * wPixelHeight);
}				// reverse y axis


complex double GiveW( int ix, int iy){
  double Wx = GiveWx(ix);
  double Wy = GiveWy(iy);
	
  return Wx + Wy*I;
	
	


}




int SetZPlane(complex double center, double radius, double a_ratio){

  ZxMin = creal(center) - radius*a_ratio;	
  ZxMax = creal(center) + radius*a_ratio;	//0.75;
  ZyMin = cimag(center) - radius;	// inv
  ZyMax = cimag(center) + radius;	//0.7;
  return 0;

}







// ****************** DYNAMICS = trap tests ( target sets) ****************************


// compute radius of circle around finite attractor which is independent of the image size ( iWidth/2000.0 )
// input k is a number of pixels ( in case of iWidth = 2000 )
double GiveAR(const double k){

	return k*PixelWidth*iWidth/2000.0 ; 
	
}


/* find such AR for internal LCM/J and LSM that level curves croses critical point and it's preimages
for attracting ( also weakly attracting = parabolic) dynamics

it may fail if one iteration is bigger then smallest distance between periodic point and Julia set
*/
double GiveTunedAR(int i_Max){

  complex double z= zcr0; // criical point
  int i;
  //int i_Max = 1000;
  // critical point escapes very fast here. Higher valus gives infinity
  for (i=0; i< i_Max; ++i ){
  z= z*z*z +c*z; // forward iteration
  }
  double r = cabs(z-za0);
  if ( r > AR_max ) {r = AR_max;}	
	 

	 
 return r;
	
	
}


double GiveMaxFinalRadius(){

  complex double z = ZxMax + ZyMax*I;
  double r = log(cabs(z))/loger - 1.0; // final_z_abs in not in [0,1]

  return r;
}
	
double GiveNormalizedFinalRadius(complex double z){

  double FinalRadius = log(cabs(z))/loger - 1.0; // final_z_abs in not in [0,1]
  return (FinalRadius/ MaxFinalRadius);

}







// bailout test
// z escapes when 
// abs(z)> ER or cabs2(z)> ER2 
// https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Julia_set#Boolean_Escape_time
// this function is not used !!!! dead code 

int Escapes(complex double z){
  // here target set (trap) is the exterior  circle with radsius = ER ( EscapeRadius) 
  // with ceter = origin z= 0
  // on the Riemann sphere it is a circle with point at infinity as a center  
   
  if (cabs(z)>ER) return 1;
  return 0;
}






// =====================
int IsPointInsideTrap0(complex double  z){

	
	if (cabs(z - za0) <AR) {return 1;} // circle around periodic point
	return 0; // outside



}



// =====================
int IsPointInsideTrap1(complex double  z){

	
	if (cabs(z - za1) <AR) {return 1;} // circle around periodic point
	
	return 0; // outside



}





/* -----------  array functions = drawing -------------- */

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


// ***********************************************************************************************
// ********************** edge detection usung Sobel filter ***************************************
// ***************************************************************************************************

// 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)); // 
 
  // 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)
  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  */
 
 
  fprintf(stderr, "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;
}

// =============================  tests ============================================


// Check Orientation of z-plane image : mark first quadrant of complex plane 
// it should be in the upper right position
// uses global var :  ...
int CheckZPlaneOrientation(unsigned char A[] )
{
 
  double Zx, Zy; //  Z= Zx+ZY*i;
  unsigned i; /* index of 1D array */
  unsigned int ix, iy;		// pixel coordinate 
	
  fprintf(stderr, "compute image CheckOrientation\n");
  // for all pixels of image 
#pragma omp parallel for schedule(dynamic) private(ix,iy, i, Zx, Zy) shared(A, ixMax , iyMax) 
  for (iy = iyMin; iy <= iyMax; ++iy){
    fprintf (stderr, " %d from %d \r", iy, iyMax);	//info 
    for (ix = ixMin; ix <= ixMax; ++ix){
      // from screen to world coordinate 
      Zy = GiveZy(iy);
      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;
}




int IsInsideWWindow(complex double w){

  if ( creal(w) < WxMax && creal(w) > WxMin &&
       cimag(w) < WyMax && cimag(w) > WyMin) {return 1;}
	
	
  return 0;
	
		


}


/*

  Array A should have image of z-plane ( not w-plane) 
  compare of image of array A unchanged
  image of w window shows part of z window and outside of z-window
 
  "Note that the flower-shaped hole in the center is originally the edge boundary of the grid."
  http://xahlee.info/SpecialPlaneCurves_dir/Inversion_dir/inversion.html
 
  https://mathworld.wolfram.com/ConformalMapping.html
  http://home.iitk.ac.in/~saiwal/engineering/complex-mappings/
 
 
*/ 
int ShowWWindowOnZWindow(unsigned char A[] )
{
 
  complex double z;
  //double Zx, Zy; //  Z= Zx+ZY*i;
  complex double w;
  unsigned i; /* index of 1D array */
  unsigned int ix, iy;		// pixel coordinate 
	
  fprintf(stderr, "compute image ShowWWindowOnZWindow\n");
  // for all pixels of image 
#pragma omp parallel for schedule(dynamic) private(ix,iy, i, w, z) shared(A, ixMax , iyMax) 
  for (iy = iyMin; iy <= iyMax; ++iy){
    fprintf (stderr, " %d from %d \r", iy, iyMax);	//info 
    for (ix = ixMin; ix <= ixMax; ++ix){
    			
      z = GiveZ(ix,iy); // from screen to world coordinate 
      w = 1/z; // invert complex plane z 
      if (IsInsideWWindow(w)){
	i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */
	A[i]=255-A[i];   // marking w window on z window
      }
    }
  }
   
   
  return 0;
}


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




int IsInsideZWindow(complex double z){

  if ( creal(z) < ZxMax && creal(z) > ZxMin &&
       cimag(z) < ZyMax && cimag(z) > ZyMin) {return 1;}
	
	
  return 0;
	
		


}


/*

  Array A should have image of w-plane ( not z-plane) 
  compare of image of array A unchanged
  image of w window shows part of z window and outside of z-window
 
  "Note that the flower-shaped hole in the center is originally the edge boundary of the grid."
  http://xahlee.info/SpecialPlaneCurves_dir/Inversion_dir/inversion.html
 
  https://mathworld.wolfram.com/ConformalMapping.html
  http://home.iitk.ac.in/~saiwal/engineering/complex-mappings/
 
 
*/ 
int ShowZWindowOnWWindow(unsigned char A[] )
{
 
  complex double z;
  //double Zx, Zy; //  Z= Zx+ZY*i;
  complex double w;
  unsigned i; /* index of 1D array */
  unsigned int ix, iy;		// pixel coordinate 
	
  fprintf(stderr, "compute image ShowZWindowOnWWindow\n");
  // for all pixels of image 
#pragma omp parallel for schedule(dynamic) private(ix,iy, i, w, z) shared(A, ixMax , iyMax) 
  for (iy = iyMin; iy <= iyMax; ++iy){
    fprintf (stderr, " %d from %d \r", iy, iyMax);	//info 
    for (ix = ixMin; ix <= ixMax; ++ix){
    			
      w = GiveW(ix,iy); // from screen to world coordinate 
      z = 1/w; // invert complex plane z 
      if (IsInsideZWindow(z)){
	i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */
	A[i]=255-A[i];   // marking w window on z window
      }
    }
  }
   
   
  return 0;
}



//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// ++++++++++++++++++++++++++++++++++++++++++ color +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

unsigned char ComputeColorOfFatou (complex double z )
{



	int nMax = iterMax;
  	int n;
  	double r;

  	for (n=0; n < nMax; n++){ //forward iteration


		z = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */ 
		r =cabs(z);
		
      		if (r > ER) // esaping = exterior
		{
	  		//uExterior += 1;
	  		return iColorOfExterior;
		}			
	
		if ( IsPointInsideTrap0(z)) {
			//uInterior +=1;
			return iColorOfInterior1;}
	
		if (IsPointInsideTrap1(z)){
			//uInterior +=1;
			return iColorOfInterior2;}




    	}

  	
  	return iColorOfUnknown;


}


// ***************************************************************************************************************************
// ************************** internal LSM/J*****************************************
// ****************************************************************************************************************************

unsigned char ComputeColorOfIntLSM (complex double z)
{



	
	
  double r;
  int nMax = iterMax;
  int n;

  
  for (n=0; n < nMax; n++){ //forward iteration
    


		

     z = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */ 
		r =cabs(z);
		
      if (r > ER) // esaping = exterior
	{
	  //uExterior += 1;
	  return iColorOfExterior;
	}			
	
      // internal level sets 	the same for both sets
      if ( IsPointInsideTrap0(z) || IsPointInsideTrap1(z)) {
	
	if ( n % 2 ) 
		{return (n*10) % 255  ;}
		else {return (n*11) % 255;}
      } 
	
		
	

    }

  //uUnknown += 1;
  return iColorOfUnknown;




}



// ***************************************************************************************************************************
// ************************** external LSM/J*****************************************
// ****************************************************************************************************************************




int GiveEscapeTime(complex double z){


  int nMax = iterMax_LSM;
  double cabsz;
  int n;

  for (n=0; n < nMax; n++){ //forward iteration
    cabsz = cabs(z);
    if (cabsz > ER_LSM) {break;} // escaping = exterior
    z = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */ 
  }
  	
  	 
  return n;

}


unsigned char ComputeColorOfExtLSM(complex double z){

  unsigned char iColor;
  int n; // escape time

	
  n = GiveEscapeTime(z);
	
  // manually udjusted series of ordered colors ( shades of gray )
  if (n==iterMax_LSM) 
  	{ iColor = iColorOfInterior;}
  	else  iColor = 255 - 230.0*((double) n)/18.0; // nMax or lower values in denominator
  
  
  return iColor;


}




// ***************************************************************************************************************************
// ************************** LSM/J = both external and internal *****************************************
// ****************************************************************************************************************************


int GiveEscapeAndAttractionTime(complex double z){


  int nMax = iterMax_LSM;
  double cabsz;
  int n;

  for (n=0; n < nMax; n++){ //forward iteration
    cabsz = cabs(z);
    if (cabsz > ER_LSM) break; // escaping = exterior
    if ( IsPointInsideTrap0(z) || IsPointInsideTrap1(z)) {break;} // attracted to finite attractor = interior
    //if (cabsz< PixelWidth) break; // fails into finite attractor = interior, but not for disconnected Julia sets, then critical point and its preimages  !!!!
    z = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */ 
  }
  	
  	 
  return n;

}


unsigned char ComputeColorOfLSM(complex double z){

  unsigned char iColor;
  int n; // escape time

	
  n = GiveEscapeAndAttractionTime(z);
	
  // manually udjusted series of ordered colors ( shades of gray )
  iColor = 255 - 230.0*((double) n)/18.0; // nMax or lower values in denominator
  
  
  return iColor;


}




// ***************************************************************************************************************************
// ************************** DEM/J*****************************************
// ****************************************************************************************************************************

unsigned char ComputeColorOfDEM(const complex double z0){
  // https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Julia_set#DEM.2FJ


  
  int nMax = iterMax_DEM;
  complex double z = z0;
  complex double dz = 1.0; //  is first derivative with respect to z.
  //double distance;
  
   ER_DEM  = 4.0;
 double cabsz;
	
  int n;

  for (n=0; n < nMax; n++){ //forward iteration
    cabsz = cabs(z);
    if ( IsPointInsideTrap0(z) || IsPointInsideTrap1(z)) {return iColorOfInterior;} // attracted to finite attractor = interior
  			
    if (cabsz > ER_DEM || cabs(dz)> 1e60) {break;} // big values 
    
    //dz = (3.0*z*z + c)*dz; 
    z = z*z*z + c*z; /* forward iteration : complex cubic polynomial */ 
  }
  
  
 // distance = 2.0 * cabsz* log(cabsz)/ cabs(dz);
 // if (distance <distanceMax) return iColorOfBoundary; // distanceMax = BoundaryWidth*PixelWidth;
  // else
  
  return iColorOfExterior;

 
}






// ***************************************************************************************************************************
// ************************** only boundary by  DEM/J*****************************************
// ****************************************************************************************************************************

unsigned char ComputeColorOfDEMJ_boundary(complex double z){
  // https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Julia_set#DEM.2FJ


  
  int nMax = iterMax_DEM;
  complex double dz = 1.0; //  is first derivative with respect to z.
  double distance;
  double cabsz;
	
  int n;

  for (n=0; n < nMax; n++){ //forward iteration
    cabsz = cabs(z);
    if (cabsz > 1e60 || cabs(dz)> 1e60) break; // big values 
    //if (cabsz< PixelWidth) return iColorOfInterior; // falls into finite attractor = interior
  			
    dz = (3.0*z*z + c)*dz; 
    z = z*z*z +c*z ; /* forward iteration : complex  qubic polynomial */ 
  }
  
  
  distance = 2.0 * cabsz* log(cabsz)/ cabs(dz);
  if (distance <distanceMax) return iColorOfBoundary; // distanceMax = BoundaryWidth*PixelWidth;
  // else
  
  return iColorOfExterior;

 
}



// plots raster point (ix,iy) 
int DrawPointOfDEMJ_boundary (unsigned char A[], int PlaneInversion, int ix, int iy, unsigned char iColor0)
{
  int i;			/* index of 1D array */
  unsigned char iColor;
  complex double z;


  i = Give_i (ix, iy);		/* compute index of 1D array from indices of 2D array */
  
  if (PlaneInversion)
    { 
      complex double w;
      w = GiveW(ix,iy);
      z = 1/w;
    }
  else {  z = GiveZ(ix,iy);}
  	
  iColor = ComputeColorOfDEMJ_boundary(z);
  if (iColor == iColorOfBoundary) // check if it is boundary
    { A[i] = iColor0 ;} // draw only boundary without changing other parts using color iColor0		
  
  return 0;
}




// fill array 
// uses global var :  ...
// scanning complex plane 
int DrawImageOfDEMJ_boundary (unsigned char A[], int PlaneInversion, const unsigned char iColor)
{
  unsigned int ix, iy;		// pixel coordinate 

  fprintf(stderr, "compute image DEM boundary\n");
  // for all pixels of image 
#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax)
  for (iy = iyMin; iy <= iyMax; ++iy){
    fprintf (stderr, " %d from %d \r", iy, iyMax);	//info 
    for (ix = ixMin; ix <= ixMax; ++ix)
      DrawPointOfDEMJ_boundary(A, PlaneInversion, ix, iy, iColor);	//  
  }

  return 0;
}










// ***************************************************************************************************************************
// ************************** Unknown: boundary and slow dynamics *****************************************
// ****************************************************************************************************************************

unsigned char ComputeColorOfUnknown(complex double z){



  
  int nMax = 20; // very low value
  
  double cabsz;
	
  int n;

  for (n=0; n < nMax; n++){ //forward iteration
    cabsz = cabs(z);
    if (cabsz > 10000000000*ER )  return iColorOfExterior; // big values
    if (cabsz < (PixelWidth/100)) return iColorOfInterior; // falls into finite attractor = interior
  			
    
    z = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */ 
  }
  
  
  
  
  //printf("found \n");
  return iColorOfUnknown;

 
}






// ***************************************************************************************************************************
// ************************** binary decomposition BD/J*****************************************
// ****************************************************************************************************************************

unsigned char ComputeColorOfBD(complex double z){

  int nMax = iterMax_LSM;
  double cabsz;
  unsigned char iColor;
	
  int n;

  for (n=0; n < nMax; n++){ //forward iteration
    cabsz = cabs(z);
    if (cabsz > ER_LSM) break; // esacping
    //if (cabsz< PixelWidth) break; // fails into finite attractor = interior but not for disconnected Julia sets, then critical point and its preimages  !!!!
    
    // attracted to finite attractor = interior
    if ( IsPointInsideTrap0(z)) { if (cimag(z)>za0im) {return 255;} else return 0;}
    if ( IsPointInsideTrap1(z)) { if (cimag(z)>za1im) {return 255;} else return 0;}			
   
    z = z*z*z +c*z ; /* forward iteration : complex cubic polynomial */ 
  }
  // exterior
  if (cimag(z)>0.0) 
    iColor = 255; 
  else iColor = 0;
  
  
  return iColor;


}





// ***************************************************************************************************************************
// ************************** modified binary decomposition MBD *****************************************
// ****************************************************************************************************************************

unsigned char ComputeColorOfMBD(complex double z){
  // const number of iterations
  int nMax = 7;
  //double cabsz;
  unsigned char iColor;
	
  int n;

  for (n=0; n < nMax; n++){ //forward iteration
    //cabsz = cabs(z);
    //if (cabsz > ER) break; // esacping
    //if (cabsz< PixelWidth) break; // falls into finite attractor = interior
  			
   
    z = z*z*z +c*z  ; /* forward iteration : complex cubic polynomial */ 
  }
  
  //if (cabs(z) > 2.0)
  { // exterior
    if (creal(z)>0.0) 
      iColor = 255; 
    else iColor = 0;
  }
  //	else iColor = iColorOfInterior;
  	
  return iColor;


}



// ***************************************************************************************************************************
// ************************** binary decomposition boundaries with texture mapping  *****************************************
// ****************************************************************************************************************************

// https://fractalforums.org/programming/11/how-many-different-ways-are-there-to-show-such-set/3874

/* 
   
   to add
   https://www.iquilezles.org/www/articles/distfunctions2d/distfunctions2d.htm
   
   
   2D Gray gradient  = 2D gray texture
   input x and y is in [0,1] range
   
*/

double Give2DGrayGradient(double x, double y, const int k ){

	
  double d;  // position of the color in the gradient . It is in [0,1]] range
	
	
  switch(k){
	
  case 0: {d =  max(fabs(x - 0.5) ,fabs(y-0.5)); break;} // 
  	
  case 1: {d =  min(x,y); break;}
  		
  case 2 : {d = fabs(x)+fabs(y) -0.5; break;}
  		
  case 3 : {d = y; break;}
  		
  case 4 : {d = x; break;}
  		
    // gradients 5,6,7 are similar , difference : 1, 1,5, 2.0 
  case 5: {x =x - 0.5; y =y - 0.5; d = cabs(x+y*I); break;} // cabs(z)
  		
  case 6: {x =1.5*(x - 0.5); y =1.5*(y - 0.5); d = cabs(x+y*I); break;} // cabs(z)
  		
  case 7: {x =2.0*(x - 0.5); y =2.0*(y - 0.5); d = cabs(x+y*I); break;} // cabs(z)
  
  default:{ d= 0.0; }
  }
  	
  return d;
}
	





unsigned char ComputeColorOfTexture(complex double z, const int k){

  int nMax = iterMax_LSM;
  double cabsz;
  unsigned char iColor;
	
  int n;

  for (n=0; n < nMax; n++){ //forward iteration
    cabsz = cabs(z);
    if (cabsz > ER_LSM) break; // esacping
    //if (cabsz< PixelWidth) break; // fails into finite attractor = interior but not for disconnected Julia sets, then critical point and its preimages  !!!!
  			
   
    z = z*z*z +c*z  ; /* forward iteration : complex cubic polynomial */ 
  }
  
  // https://gitlab.com/adammajewski/mandelbrot-book_book/-/blob/master/README.md#final-angle
  //if (n < nMax) // exterior 
  // {
  //double et = ((double)n)/nMax; // ok but the same for all points inside level set so segmentation
  double final_angle = c_turn(z); // in [0,1] range
  //double final_radius = GiveNormalizedFinalRadius(z); // = final_absz should be in [0,1]
        
        
  //}

  
  
  double y = frac(n-log(log(cabsz)));
  // inside each cell point has additional coordinate w = (final_angle, final_radius) in [0,1]x[0,1]
  double gray = Give2DGrayGradient(final_angle, y, k);
  iColor = gray*255;
  
  // bd : mark each cell 
  //if (cimag(z)>0.0) iColor =255 -iColor;
  
  return iColor;


}



// plots raster point (ix,iy) 
int DrawPointOfTexture (unsigned char A[], int ix, int iy, const int k)
{
  int i;			/* index of 1D array */
  unsigned char iColor;
  complex double z;


  i = Give_i (ix, iy);		/* compute index of 1D array from indices of 2D array */
  z = GiveZ(ix,iy);
  iColor = ComputeColorOfTexture(z, k);
  A[i] = iColor ;		// interior
  
  return 0;
}




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

  fprintf(stderr, "compute image texture k = %d\n", k);
  // for all pixels of image 
#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax)
  for (iy = iyMin; iy <= iyMax; ++iy){
    fprintf (stderr, " %d from %d \r", iy, iyMax);	//info 
    for (ix = ixMin; ix <= ixMax; ++ix)
      DrawPointOfTexture(A, ix, iy, k);	//  
  }

  return 0;
}








// ***********************************************************************************************
//*************************************** SAC/J **************************************************
// *****************************************************************************************
// https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/stripeAC
// SAC = Stripe Average Coloring

//

// the addend function
// input : complex number z
// output : double number t 
double Give_t(double complex z){

  return 0.5+0.5*sin(s*carg(z));

}

/*
  input :
  - complex number
  - intege
  output = average
 
*/
double Give_Arg(double complex z , int iMax)
{
  int i=0; // iteration 
   
   
  //double complex Z= 0.0; // initial value for iteration Z0
  double A = 0.0; // A(n)
  double prevA = 0.0; // A(n-1)
  double R; // =radius = cabs(Z)
  double d; // smooth iteration count
  double complex dz = 1.0; // first derivative with respect to z
  double de; // Distance Estimation from DEM/J  
   
    
  // iteration = computing the orbit
  for(i=0;i<iMax;i++)
    { 
    
      dz = (3.0*z*z + c)*dz; 
      z = z*z*z +c*z  ; /* forward iteration : complex cubic polynomial */ 
      
      if (i>i_skip) A += Give_t(z); // 
      
      R = cabs(z);
      // if(R > EscapeRadius) break; // exterior of M set
      if (R > 1e60 || cabs(dz)> 1e60) break; // prevent NAN 	 	
      prevA = A; // save value for interpolation
        
    } // for(i=0
   
   
  if (i == iMax) 
    A = -1.0; // interior 
  else { // exterior
    de = 2 * R * log(R) / cabs(dz);
    if (de < distanceMax) A = FP_ZERO; //  boundary
    else {
      // computing interpolated average
      A /= (i - i_skip) ; // A(n)
      prevA /= (i - i_skip - 1) ; // A(n-1) 
      // smooth iteration count
      d = i + 1 + log(lnER/log(R))/M_LN2;
      d = d - (int)d; // only fractional part = interpolation coefficient
      // linear interpolation
      A = d*A + (1.0-d)*prevA;
    }   
  }
    
  return A;  
}
 
 
 
 
 
unsigned char ComputeColorOfSAC(complex double z){

  unsigned char iColor;
  double arg;
  
   
   
  arg = Give_Arg( z, 2500); //   N in wiki 
	
  // color is proportional to arg 
  if (arg < 0.0)
           
    iColor = 0;  // interior                        
    
  else //  
    {if (arg == FP_ZERO) 
	iColor = 255; // boundary     
      else iColor = (unsigned char) (255 - 255*arg );// exterior
    }
      
    
  return iColor;


}








 

// ***************************************************************************************************************************
// ************************** DLD/J*****************************************
// ****************************************************************************************************************************



/* partial pnorm 
   input: z , zn = f(z), p
   output ppn
   
   
*/
double
ppnorm (complex double z, complex double zn, double p)
{

  double s[2][3];		// array for 2 points on the Riemann sphere
  int j;
  double d;			// denominator 
  double x;
  double y;

  double ds;
  double ppn = 0.0;

  // map from complex plane to riemann sphere
  // z
  x = creal (z);
  y = cimag (z);
  d = x * x + y * y + 1.0;

  s[0][0] = (2.0 * x) / d;
  s[0][1] = (2.0 * y) / d;
  s[0][2] = (d - 2.0) / d;	// (x^2 + y^2 - 1)/d

  // zn
  x = creal (zn);
  y = cimag (zn);
  d = x * x + y * y + 1.0;
  s[1][0] = (2.0 * x) / d;
  s[1][1] = (2.0 * y) / d;
  s[1][2] = (d - 2.0) / d;	// (x^2 + y^2 - 1)/d

  // sum 
  for (j = 0; j < 3; ++j)
    {
      ds = fabs (s[1][j] - s[0][j]);
      //  normal:  neither zero, subnormal, infinite, nor NaN
      //if (fpclassify (ds) !=FP_INFINITE)
      //if (isnormal(ds)) 
      // it is solved by if (cabs(z) > 1e60 ) break; procedure in parent function 
      ppn += pow (ds, p);	// |ds|^p
      //      else {ppn = 10000.0; printf("ds = infty\t");} // 

    }


  return ppn;







}

// DLD = Discret Lagrangian Descriptior
double
lagrangian (complex double z0, complex double c, int iMax, double p)
{

  int i;			// number of iteration
  double d = 0.0;		// DLD = sum
  double ppn;			// partial pnorm
  complex double z = z0;
  complex double zn;		// next z

  for (i = 0; i < iMax; ++i)
    {




      zn = z*z*z +c*z  ; /* forward iteration : complex cubic polynomial */ 
      ppn = ppnorm (z, zn, p);
      d += ppn;			// sum
      //
      z = zn;

      //if (! isnormal(d)) { return 0.0; } // not works
      if (cabs (z) > ER_DLD ) //1e6)
	break;			// exterior : big values produces artifacts on the image  



    }





  //if (d<0.0) {// interior
  // d(z1a) - d(z21) = -0.0804163521959989        
  //      d = - d;
  //      d = (db - d) /dd ; // normalize, see test_interior
  //d = d*d;
  //if (d>1.0) {printf("d int > 1.0\n");
  ///     }
  //      else {

  d = d / ((double) i);		// averaging not summation
  //d = d*me;} // exterior

  return d;




}





unsigned char
ComputeColorOfDLD (complex double z)
{


  //double cabsz;
  int iColor;
  double d;
  int N = iterMax_DLD; // N in wiki = fixed number : maximal number of iterations 

  //if (FatouType == 1)
  // {				// interior
  // d = lagrangian (z, c, N, p);
  // modify gradient position

  //{d = d - (int)d;} // only fractional part
  // d = d * d * mi;
  //if ( d< 1.0 ) d = 0.0;

  // }				//  
  //else
  //{
  d = lagrangian (z, c, N, p); //  
  //}

  iColor = (int) (d * 255) % 255;	// nMax or lower walues in denominator



  return (unsigned char) iColor;


}




//=========================================

 
 

// ***************************************************************************************************************************
// ************************** NPM/J = Normal Potential *****************************************
// ****************************************************************************************************************************





/* 
   The dot product of two vectors a = [a1, a2, ..., an] and b = [b1, b2, ..., bn] is defined as:[1]
   d = a1b1 + a2b2
  
*/
double cdot(double complex a, double complex b){
 
 
  return creal(a)*creal(b)+cimag(a)*cimag(b); 


}


// 
// output 
// 
double GiveReflection(double complex z )
{
  int i=0; // iteration 
  int iMax = 2000;
   
  // https://en.wikibooks.org/wiki/Fractals/Mathematics/Derivative
  double complex dz = 1.0; // derivative with respect to z 
  double reflection = 0.0; //  
   
  double h2 = 1.5 ; // height factor of the incoming light
  double angle = 45.0/360.0 ; // incoming direction of light in turns 
  double complex v = cexp(2.0*angle *M_PI* I); // = exp(1j*angle*2*pi/360)  // unit 2D vector in this direction
  // incoming light 3D vector = (v.re,v.im,h2)
  
  // https://en.wikipedia.org/wiki/Lambertian_reflectance

   
  double  complex u;
   
   
  z = z*z*z +c*z  ; /* forward iteration : complex cubic polynomial */ 
  dz = 1.0;
   
  for(i=0;i<iMax;i++)
    {  
      dz = (3.0*z*z + c)*dz; 
      z = z*z*z +c*z  ; /* forward iteration : complex cubic polynomial */ 
      
      
      
      if(cabs(z) > ER_NP) 
	{ // exterior
	  u = z / dz;
	  u = u / cabs(u);
	  reflection =  cdot(u, v) + h2;  /* use the simplest model for the shading: 
					     Lambert, which consists in using the dot product of (x,y,1) with a constant vector indicating the direction of the light. */
	  reflection = reflection/(1.0 + h2);  // rescale so that t does not get bigger than 1
	  if (reflection<0.0) reflection =0.0;
	  break;
      
	}
    }
    
  return reflection;  
}





// Potential to color
unsigned char ComputeColorOfNP(complex double z){
  //https://www.math.univ-toulouse.fr/~cheritat/wiki-draw/index.php/Mandelbrot_set#Normal_map_effect


  
  
  double reflection;
  unsigned char iColor;
   
   
   
  // compute 
  reflection = GiveReflection( z);
  
    
  // 
  //if (reflection <  )
  //{ /*  interior  */
  //  iColor = 0;}
  //else // exterior 
        
  { iColor = 255 * reflection;}
     
  return iColor;   
  
 
}



// https://en.wikipedia.org/wiki/Shading

//  normal = perpendicular 
// shading using Normal map and Potential
// https://en.wikipedia.org/wiki/Lambertian_reflectance
// http://www.math.titech.ac.jp/~kawahira/gallery/movies/movies.html
// see 0_1.avi and image 
//






 
 

// ***************************************************************************************************************************
// ************************** NDM/J = Normal Distance *****************************************
// ****************************************************************************************************************************


//  normal = perpendicular 
// shading using Normal map and Potential
// https://en.wikipedia.org/wiki/Lambertian_reflectance
// http://www.math.titech.ac.jp/~kawahira/gallery/movies/movies.html
// see 0_1.avi and image 
//





// 
// output 
// 
double GiveReflectionD(double complex z )
{
  int i=0; // iteration 
  int iMax = 2000;
   
  // https://en.wikibooks.org/wiki/Fractals/Mathematics/Derivative
  double complex dz = 1.0;   // first derivative with respect to z 
  double complex dz2 = 0.0;  // second derivative with respect to z 
  double reflection = 0.0; //  
  double lo;
   
  double h2 = 1.5 ; // height factor of the incoming light
  double angle = 45.0/360.0 ; // incoming direction of light in turns 
  double complex v = cexp(2.0*angle *M_PI* I); // = exp(1j*angle*2*pi/360)  // unit 2D vector in this direction
  // incoming light 3D vector = (v.re,v.im,h2)
   
  
  // https://en.wikipedia.org/wiki/Lambertian_reflectance

   
  double  complex u;
   
   
  z = z*z*z +c*z  ; /* forward iteration : complex cubic polynomial */ 
  dz = 1.0;
  dz2 = 0.0;
   
  for(i=0;i<iMax;i++)
    {  
      
      dz2 = 2.0* ( dz2*z + dz*dz);//2*(der2*z+der**2)
      dz = (3.0*z*z + c)*dz; 
      z = z*z*z +c*z  ; /* forward iteration : complex cubic polynomial */ 
      
      
      
      if(cabs(z) > ER_NP) 
	{ // exterior
      
	  /*
	    lo = 0.5*log(squared_modulus(z))
	    u = z*der*((1+lo)*conj(der**2)-lo*conj(z*der2))
	    u = u/abs(u)
	  */
      
	  lo = 0.5*log(cabs(z));
	  u = z*dz*((1.0+lo)*conj(dz*dz)-lo*conj(z*dz2));
	  //u = z / dz;
	  u = u / cabs(u);
	  reflection =  cdot(u, v) + h2;  // use the simplest model for the shading: Lambert, which consists in using the dot product of (x,y,1) with a constant vector indicating the direction of the light.
	  reflection = reflection/(1.0 + h2);  // rescale so that t does not get bigger than 1
	  if (reflection<0.0) reflection =0.0;
	  break;
      
	}
    }
    
  return reflection;  
}





// Distance to color
unsigned char ComputeColorOfND(complex double z){
  //https://www.math.univ-toulouse.fr/~cheritat/wiki-draw/index.php/Mandelbrot_set#Variation


  
  
  double reflection;
  unsigned char iColor;
   
   
   
  // compute 
  reflection = GiveReflectionD( z);
  
    
  // 
  //if (reflection <  )
  //{ /*  interior  */
  //  iColor = 0;}
  //else // exterior 
        
  { iColor = 255 * reflection;}
     
  return iColor;   
  
 
}


// -------------------------- potential========


double ComputePotential(const complex double z0){

  double potential = 0.0; // interior
  double s = 0.5;
  complex double z = z0;
  double r;
  int iter;
	
  for (iter = 0; iter < iterMax_pot; ++iter){
		
    z = z*z*z +c*z  ; /* forward iteration : complex cubic polynomial */ 
    s *= 0.5;  // 
    r = cabs(z);
    if (r > ER_pot) {break;}
	
	
  }
	
	
	
	
	
  potential =  s*log2(r); // log(zn)* 2^(-n)
  return potential;
	
}


unsigned char ComputeColorOfPOT(complex double z){


  double potential = ComputePotential(z);
	
  if (PlaneInversion) // usung global var 
    {potential /= 4.0;} // manual normalize
  unsigned char iColor = 255 * sqrt(sqrt(potential));
  return iColor;   
  
 
}

double GiveSmoothEscapeTime(complex double z){


  int nMax = iterMax_LSM;
  double cabsz;
  int n;

  for (n=0; n < nMax; n++){ //forward iteration
    cabsz = cabs(z);
    if (cabsz > ER_LSM) break; // esacping
    //if (cabsz< PixelWidth) break; // fails into finite attractor = interior, but not for disconnected Julia sets, then critical point and its preimages  !!!!
    z = z*z*z +c*z  ; /* forward iteration : complex cubic polynomial */ 
  }
  	
  //float sn = n - log2(log2(dot(z,z))) + 4.0;  // equivalent optimized smooth iteration count
  	
  double sn = ER_LSM/cabsz;
  //n- log2(log2(cdot(z,z))) + 4.0; 
  //sn = sn / nMax; // map to [0,1] range
  	
  return sn;

}



// 
unsigned char ComputeColorOfBlend(complex double z){

	
	
	
  double SET = GiveSmoothEscapeTime(z);
  SET = sqrt(SET);
  SET = 1.0 - SET;
  // 
  double ColorSET = SET*255;
	
  //
  double ColorNP = ComputeColorOfNP(z);
		
  unsigned char iColor = (ColorSET+ ColorNP)/ 2.0; // average blend mode
	
  return iColor;   
  
 
}



/*
  int local_setup(int PlaneInversion){

  if (PlaneInversion)
  { MaxImagePotential =ComputePotential( 1.0/ 0.0);}
  //else {MaxImagePotential}

  return 0;

  }

*/


 
 
/* ==================================================================================================
   ============================= Draw functions ===============================================================
   =====================================================================================================
*/ 
unsigned char ComputeColor(FunctionTypeT FunctionType, complex double z){

  unsigned char iColor;
	
	
	
  switch(FunctionType){
  
  case Fatou :{iColor = ComputeColorOfFatou(z); break;}
  
  case IntLSM :{iColor = ComputeColorOfIntLSM(z); break;}
	
  case ExtLSM :{iColor = ComputeColorOfExtLSM(z); break;}
  
  case LSM :{iColor = ComputeColorOfLSM(z); break;}
		
  case DEM : {iColor = ComputeColorOfDEM(z); break;}
		
  case Unknown : {iColor = ComputeColorOfUnknown(z); break;}
		
  case BD : {iColor = ComputeColorOfBD(z); break;}
		
  case MBD : {iColor = ComputeColorOfMBD(z); break;}
		
  case SAC : {iColor = ComputeColorOfSAC(z); break;}
  
  case DLD : {iColor = ComputeColorOfDLD(z); break;}
		
  case ND : {iColor = ComputeColorOfND(z); break;}
		
  case NP : {iColor = ComputeColorOfNP(z); break;}
		
  case POT : {iColor = ComputeColorOfPOT(z); break;}
		
  case Blend : {iColor = ComputeColorOfBlend(z); break;}
		
	
  default: {}
	
	
  }
	
  return iColor;



}


// plots raster point (ix,iy) 
int DrawPoint (FunctionTypeT FunctionType, int PlaneInversion, unsigned char A[], int ix, int iy)
{
  int i;			/* index of 1D array */
  unsigned char iColor;
  complex double z;


  i = Give_i (ix, iy);		/* compute index of 1D array from indices of 2D array */
  
  if (PlaneInversion)
    { 
      complex double w;
      w = GiveW(ix,iy);
      z = 1/w;
    }
  else {  z = GiveZ(ix,iy);}
  

  iColor = ComputeColor(FunctionType, z);
  A[i] = iColor ;		// 
  
  return 0;
}




// fill array 
// uses global var :  ...
// scanning complex plane 
int DrawImage (FunctionTypeT FunctionType, int PlaneInversion, unsigned char A[])
{
  unsigned int ix, iy;		// pixel coordinate 
  	
  //local_setup(PlaneInversion);

  fprintf(stderr, "compute image FunctionType = %d PlaneInversion = %d\n", FunctionType, PlaneInversion);
  // for all pixels of image 
  #pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax)
  for (iy = iyMin; iy <= iyMax; ++iy){
    fprintf (stderr, " %d from %d \r", iy, iyMax);	//info 
    for (ix = ixMin; ix <= ixMax; ++ix)
      DrawPoint(FunctionType, PlaneInversion, A, ix, iy);	//  
  }

  return 0;
}

 
// ------------------------------------------------------------------------------------------------------------
// ---------------------------Additional drawings : traps, attractors and critical orbits -------------------------
// --------------------------------------------------------------------------------------------------------------
int PointIsInsideZRectangle(const double complex z){

	if (ZxMin < creal(z) && creal(z)<ZxMax && ZyMin < cimag(z) && cimag(z)<ZyMax)
		{return 1;}
	return 0;
	




}

int IsInside (int x, int y, int xcenter, int ycenter, int r){

	
	double dx = x- xcenter;
	double dy = y - ycenter;
	double d = sqrt(dx*dx+dy*dy);
	if (d<r) 
		return 1;
	return 0;
	  

} 



int PlotBigPoint(complex double z, unsigned char A[]){

	
	unsigned int ix_seed = (creal(z)-ZxMin)/PixelWidth;
	unsigned int iy_seed = (ZyMax - cimag(z))/PixelHeight;
	unsigned int i;
	
	
	 /* mark seed point by big pixel */
  	int iSide =3.0*iWidth/2000.0 ; /* half of width or height of big pixel */
  	int iY;
  	int iX;
  	for(iY=iy_seed-iSide;iY<=iy_seed+iSide;++iY){ 
    		for(iX=ix_seed-iSide;iX<=ix_seed+iSide;++iX){ 
    			if (IsInside(iX, iY, ix_seed, iy_seed, iSide)) {
      			
      			i= Give_i(iX,iY); /* index of _data array */
      			
      			//A[i]= 255-A[i]; // inverted color
      			A[i] = 0; // allways black
      			}}}
	
	
	return 0;
	
}


// =====================
int IsPointInsideTraps(unsigned int ix, unsigned int iy){

	
	complex double  z = GiveZ (ix, iy);
	
	if ( IsPointInsideTrap0(z)) {return 1;} // circle with prabolic point on it's boundary
	
	if (IsPointInsideTrap1(z)) {return 1;}
	
	return 0; // outside



}





int MarkTraps(unsigned char A[]){

	unsigned int ix, iy;		// pixel coordinate 
	unsigned int i;


  	fprintf (stderr, "Mark traps \n");
  	// for all pixels of image 
	#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(A, ixMax , iyMax)
  	for (iy = iyMin; iy <= iyMax; ++iy)
    	{
      		fprintf (stderr, " %d from %d \r", iy, iyMax);	//info 
      		for (ix = ixMin; ix <= ixMax; ++ix){
			if (IsPointInsideTraps(ix, iy)) {
      				i= Give_i(ix,iy); /* index of _data array */
      				A[i]= 255-A[i]; // inverse color
      				}}}
  	return 0;
}




int DrawForwardOrbit(const complex double z0,  const unsigned long long int iMax, unsigned char A[] )
{
  
  unsigned long long int i; /* nr of point of critical orbit */
  
  complex double z = z0;
  
 
  PlotBigPoint(z, A);
  
  /* forward orbit of critical point  */
  for (i=1;i<iMax ; ++i)
    {
      z  = z*z*z+ c*z;
      // if (cabs2(z - z2a) > 2.0) {return 1;} // escaping
      if (PointIsInsideZRectangle(z))
      		{PlotBigPoint(z, A);}
      		else fprintf(stderr, "bad point z\n");
    }
  

    
   
  return 0;
 
}


void DrawCriticalOrbits(unsigned char A[]){

	unsigned long long int iMax = iterMax_LSM ;
	
	DrawForwardOrbit(zcr0, iMax, A);
	DrawForwardOrbit(zcr1, iMax, A);

}



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




 
int Test()
{
  unsigned int ix, iy;		// pixel coordinate 
  complex double z;
  double SET;
  int ET;
  	
  //local_setup(PlaneInversion);

  fprintf(stderr, "test\n");
  // for all pixels of image 
  ix = 0;
  for (iy = iyMin; iy <= iyMax; ++iy){
    z = GiveZ(ix,iy);
    ET = GiveEscapeTime(z);
    SET = GiveSmoothEscapeTime(z);
    printf(" %d \t %d \t %f \n", iy, ET, SET);		 
    ix = ix +1;
    //  
  }

  return 0;
}

  
 
 
 
 
 
 
// *******************************************************************************************
// ********************************** save A array to pgm file ****************************
// *********************************************************************************************

int
SaveArray2PGMFile (unsigned char A[], char *shortName , char *comment)
{

  FILE *fp;
  const unsigned int MaxColorComponentValue = 255;	/* color component is coded from 0 to 255 ;  it is 8 bit color file */
  
  
  // https://programmerfish.com/create-output-file-names-using-a-variable-in-c-c/
  char fileName[512];
  const char* fileType = ".pgm";
  sprintf(fileName,"%s%s", shortName, fileType); // 
  
  
  
  char long_comment[200];
  sprintf (long_comment, "f(z) = z*z*z +c*z   where c = %f %+f*i ;  %s", creal(c), cimag(c),comment);





  // save image array 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", long_comment, iWidth, iHeight, MaxColorComponentValue);	// write header to the file
  size_t rSize = fwrite (A, sizeof(A[0]), iSize, fp);	// write whole array with image data bytes to the file in one step 
  fclose (fp);

  // info 
  if ( rSize == iSize) 
    {
      printf ("File %s saved ", fileName);
      if (long_comment == NULL || strlen (long_comment) == 0)
	printf ("\n");
      else { printf (". Comment = %s \n", long_comment); }
    }
  else {printf("wrote %zu elements out of %u requested\n", rSize,  iSize);}
  	
  	
  // 
  NumberOfImages +=1; // count images using global variable

  return 0;
}
















int PrintInfoAboutProgam()
{
  printf("Number of pgm images = %d \n", NumberOfImages);	
  
  // display info messages
  printf ("Numerical approximation of Julia set for fc(z)= z^2 + c \n");
  //printf ("iPeriodParent = %d \n", iPeriodParent);
  //printf ("iPeriodOfChild  = %d \n", iPeriodChild);
  printf ("parameter c = %.16f %+.16f*i  \n", creal(c), cimag(c));
  
  printf ("Image Width = %f in world coordinate\n", ZxMax - ZxMin);
  printf ("PixelWidth = %f \n", PixelWidth);
  
  printf("for DEM/J \n");
  if ( distanceMax<0.0 || distanceMax > ER ) printf("bad distanceMax\n");
  printf("Max distance from exterior to the boundary =  distanceMax = %.16f = %f pixels\n",  distanceMax, BoundaryWidth); 
  printf("\n");
  
  
  // image corners in world coordinate
  // center and radius
  // center and zoom
  // GradientRepetition
  printf ("Maximal number of iterations = iterMax = %ld \n", iterMax);
  
  printf ("For LSM/J \n");
  printf ("Maximal number of iterations = iterMax_LSM = %ld \n", iterMax_LSM);
  printf ("Escape Radius = ER_LSM = %f \n", ER_LSM);
  printf("\n");
  
  
  
  
  printf ("ratio of image  = %f ; it should be 1.000 ...\n", ratio);
  //
  printf("gcc version: %d.%d.%d\n",__GNUC__,__GNUC_MINOR__,__GNUC_PATCHLEVEL__); // https://stackoverflow.com/questions/20389193/how-do-i-check-my-gcc-c-compiler-version-for-my-eclipse
  // OpenMP version is displayed in the console 
  return 0;
}






// *****************************************************************************
//;;;;;;;;;;;;;;;;;;;;;;  setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
// **************************************************************************************

int setup ()
{

  fprintf (stderr, "setup start\n");
  // c= 1.0049542069308062 +0.1008317508132964*i
  c =  1.02*cexp(0.1*I);  ; //   https://web.archive.org/web/20161024194536/http://www.ijon.de/mathe/julia/some_julia_sets_3.html
  
  
  
  
	
  /* 2D array ranges */
  
  iWidth = iHeight* DisplayAspectRatio;
  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].
  //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].
  
  center = za1;
  SetZPlane( center, radius,  DisplayAspectRatio );	

  /* Pixel sizes */
  PixelWidth = (ZxMax - ZxMin) / ixMax;	//  ixMax = (iWidth-1)  step between pixels in world coordinate 
  PixelHeight = (ZyMax - ZyMin) / iyMax;
  ratio = ((ZxMax - ZxMin) / (ZyMax - ZyMin)) / ((double) iWidth / (double) iHeight);	// it should be 1.000 ...
	
  wPixelWidth = (WxMax-WxMin)/ixMax;
  wPixelHeight =(WyMax-WyMin)/iyMax;
	
  
  //ER2 = ER * ER; // for numerical optimisation in iteration
  lnER = log(EscapeRadius); // ln(ER) 
  loger = log(ER_LSM); // for texture
  ER_LSM = 3.0; //GiveER(10); // find such ER for LSM/J that level curves croses critical point and it's preimages
  ER_DLD = 3.0; //GiveER(7);
  
  MaxFinalRadius =  GiveMaxFinalRadius();
  
  AR_max = GiveAR(50.0); //PixelWidth*50.0*iWidth/2000.0 ; // moved to main 
  
  // BD
  za0im = cimag(za0);
  za1im = cimag(za1);
  
   	
  /* create dynamic 1D arrays for colors ( shades of gray ) */
  data = malloc (iSize * sizeof (unsigned char));
  edge = malloc (iSize * sizeof (unsigned char));
  edge2 = malloc (iSize * sizeof (unsigned char));
  //
 
  	
  if (data == NULL || edge == NULL || edge2 == NULL ){
    fprintf (stderr, " Could not allocate memory");
    return 1;
  }

  
 	
  
  BoundaryWidth = 1.0*iWidth/2000.0  ; //  measured in pixels ( when iWidth = 2000) 
  distanceMax = BoundaryWidth*PixelWidth;
  
  
  
  fprintf (stderr," end of setup \n");
	
  return 0;

} // ;;;;;;;;;;;;;;;;;;;;;;;;; end of the setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;




int end(){


  fprintf (stderr," allways free memory (deallocate )  to avoid memory leaks \n"); // https://en.wikipedia.org/wiki/C_dynamic_memory_allocation
  free (data);
  free(edge);
  free(edge2);
 
  PrintInfoAboutProgam();
  return 0;

}

// ********************************************************************************************************************
/* -----------------------------------------  main   -------------------------------------------------------------*/
// ********************************************************************************************************************

int main () {
  
  
  
  setup ();
  
  PlaneInversion = 0; 
 /*
  printf("bigger trap radius for speed \n");
  AR = AR_max; 
  
  printf("AR_max = %.16f  = %f pixels\n",AR_max , AR_max/PixelWidth);
  
    
  DrawImage(Fatou, PlaneInversion, data);
  SaveArray2PGMFile (data, "Fatou", "Fatou set"); 
  
  ComputeBoundaries(data,edge);
  SaveArray2PGMFile (edge, "FatouB", "boundaries of the Fatou set "); 
  
  CopyBoundaries(edge,data);
  SaveArray2PGMFile (data, "FatouAndB", "Fatou set and boundaries");  
  
  
  // mark attractors   
  PlotBigPoint(za0, data);
  PlotBigPoint(za1, data);
  SaveArray2PGMFile (data, "Fatou_p", "Fatou and attracting points"); 
  
  MarkTraps(data);
  SaveArray2PGMFile (data, "Fatou_traps", "Fatou and traps"); 
  
  DrawCriticalOrbits(data);
  SaveArray2PGMFile (data, "Fatou_cr", "Fatou, traps and critical orbits");
  
  AR = GiveAR(10);
  
  DrawImage(BD, PlaneInversion, data);
  SaveArray2PGMFile (data, "BD", "BD"); // name of the file is name.png   
  
  */ 
  printf("smaller trap radius for more detailes \n");
  AR = GiveTunedAR(400); 
  printf("tuned AR = %.16f = %.16f *AR_max = %f pixels\n",AR, AR/AR_max, AR/PixelWidth);
  /* 
  DrawImage(IntLSM, PlaneInversion, data);
  SaveArray2PGMFile (data, "IntLSM", "internal level sets = level sets of attraction time"); // name of the file is name.png 
  
  ComputeBoundaries(data,edge);
  SaveArray2PGMFile (edge, "IntLCM", "internal Level Curves = boundaries of Int Level Sests "); // name of the file is name.png 
  
  CopyBoundaries(edge,data);
  SaveArray2PGMFile (data, "IntLSC", "Int LevelSets and it's boundaries = LevelCurves"); // name of the file is name.png 
  */
  DrawImage(LSM, PlaneInversion, data);
  SaveArray2PGMFile (data, "LSM", "internal and external level sets = level sets of attraction and escaping time"); // name of the file is name.png 
  
  
  ComputeBoundaries(data,edge);
  SaveArray2PGMFile (edge, "LCM_zoom", "internal and external Level Curves = boundaries of both Level Sests "); // name of the file is name.png 
  
  DrawCriticalOrbits(edge);
  SaveArray2PGMFile (edge, "LCM_cr_zoom", "LCM and critical orbits");
  
  /*
  CopyBoundaries(edge,data);
  SaveArray2PGMFile (data, "LSCM", "Int and ext LevelSests and it's boundaries = LevelCurves"); // name of the file is name.png 
  
 
   DrawImage(DEM, PlaneInversion, data);
   SaveArray2PGMFile (data, "DEM", "external Distance Estimation Method = boundary of filled Julia set"); // name of the file is name.png 
   
   */
   
   
  //Test();
  
  //
  end();

  return 0;
}


Makefile

all: 
	chmod +x d.sh
	./d.sh

bash source code

#!/bin/bash 
 
# script file for BASH 
# which bash
# save this file as d.sh
# chmod +x d.sh
# ./d.sh
# checked in https://www.shellcheck.net/


# display OMP info 
export  OMP_DISPLAY_ENV="TRUE"

printf "make pgm and txt files \n"
gcc d.c -lm -Wall -march=native -fopenmp

if [ $? -ne 0 ]
then
    echo ERROR: compilation failed !!!!!!
    exit 1
fi


time ./a.out > a.txt

export  OMP_DISPLAY_ENV="FALSE"

printf "convert all pgm files to png using Image Magic convert \n"
# for all pgm files in this directory
for file in *.pgm ; do
  # b is name of file without extension
  b=$(basename "$file" .pgm)
  # convert  using ImageMagic
  convert "${b}".pgm -resize 2000x2000 "${b}".png
  echo "$file"
done

export  OMP_DISPLAY_ENV="TRUE"
printf "display OMP info \n"

printf "delete all pgm files \n"
rm ./*.pgm

 
echo OK
# end


text output

chmod +x d.sh
./d.sh
make pgm and txt files 

OPENMP DISPLAY ENVIRONMENT BEGIN
  _OPENMP = '201511'
  OMP_DYNAMIC = 'FALSE'
  OMP_NESTED = 'FALSE'
  OMP_NUM_THREADS = '8'
  OMP_SCHEDULE = 'DYNAMIC'
  OMP_PROC_BIND = 'FALSE'
  OMP_PLACES = ''
  OMP_STACKSIZE = '0'
  OMP_WAIT_POLICY = 'PASSIVE'
  OMP_THREAD_LIMIT = '4294967295'
  OMP_MAX_ACTIVE_LEVELS = '2147483647'
  OMP_CANCELLATION = 'FALSE'
  OMP_DEFAULT_DEVICE = '0'
  OMP_MAX_TASK_PRIORITY = '0'
  OMP_DISPLAY_AFFINITY = 'FALSE'
  OMP_AFFINITY_FORMAT = 'level %L thread %i affinity %A'
OPENMP DISPLAY ENVIRONMENT END
setup start
 end of setup 
compute image FunctionType = 3 PlaneInversion = 0
bad point z 
...
bad point z
 allways free memory (deallocate )  to avoid memory leaks 

real	19m13,945s
user	138m33,613s
sys	0m5,662s
convert all pgm files to png using Image Magic convert 
LCM_cr_zoom.pgm
LCM_zoom.pgm
LSM.pgm
display OMP info 
delete all pgm files 
OK







smaller trap radius for more detailes 
tuned AR = 0.0038665412903206 = 0.1546539185302452 *AR_max = 77.326959 pixels
File LSM.pgm saved . Comment = f(z) = z*z*z +c*z   where c = 1.014904 +0.101830*i ;  internal and external level sets = level sets of attraction and escaping time 
File LCM_zoom.pgm saved . Comment = f(z) = z*z*z +c*z   where c = 1.014904 +0.101830*i ;  internal and external Level Curves = boundaries of both Level Sests  
File LCM_cr_zoom.pgm saved . Comment = f(z) = z*z*z +c*z   where c = 1.014904 +0.101830*i ;  LCM and critical orbits 
Number of pgm images = 3 
Numerical approximation of Julia set for fc(z)= z^2 + c 
parameter c = 1.0149042485835864 +0.1018300849797647*i  
Image Width = 1.000000 in world coordinate
PixelWidth = 0.000050 
for DEM/J 
Max distance from exterior to the boundary =  distanceMax = 0.0005000250012501 = 10.000000 pixels

Maximal number of iterations = iterMax = 1000000 
For LSM/J 
Maximal number of iterations = iterMax_LSM = 1000 
Escape Radius = ER_LSM = 3.000000 

ratio of image  = 1.000000 ; it should be 1.000 ...
gcc version: 10.2.0

references

  1. Some Julia sets 3 by Michael Becker

Captions

Julia set p(z)= z^3+(1.0149042485835864102+0.10183008497976470119i)*z. Zoom and critical orbit

Items portrayed in this file

depicts

16 May 2021

File history

Click on a date/time to view the file as it appeared at that time.

Date/TimeThumbnailDimensionsUserComment
current17:35, 17 May 2021Thumbnail for version as of 17:35, 17 May 20212,000 × 2,000 (1.65 MB)Soul windsurferUploaded own work with UploadWizard

Global file usage

The following other wikis use this file:

Metadata