File:Douady rabbit - final phase after fixed number of iterations without bailout test.gif

/*

   https://commons.wikimedia.org/wiki/File:Julia_IIM_6_basilica.png
  
   c console program:
   
   
   --------------------------------
   1. draws Julia set for Fc(z)=z*z +c
   using  IIM
   -------------------------------         
   2. technique of creating ppm file is  based on the code of Claudio Rocchini
   http://en.wikipedia.org/wiki/Image:Color_complex_plot.jpg
   create 24 bit color graphic file ,  portable pixmap file = PPM 
   see http://en.wikipedia.org/wiki/Portable_pixmap
   to see the file use external application ( graphic viewer)
   
   
   complex point c -> virtual 2D array -> memory 1D array -> ppm file on the disc -> png file 
   
   Z -> pixel (iX,iY)  -> index k  -> 24bit color 
   
   -----
   https://stackoverflow.com/questions/6418807/how-to-work-with-complex-numbers-in-c
   complex numbers are built in type 
 
   --------------
   formated with emacs
   -------------
   to compile :

 
 
   gcc phase_f.c -lm -Wall 
 
 
   ./a.out
   
   
   to convert to png using ImageMagic

   convert phase.ppm phase.png  
   
   use f.sh to create animated gif

   ----------------------
   cd existing_folder
   git init
   git remote add origin git@gitlab.com:adammajewski/dynamic_external_angle.git
   git add phase_f.c
   git commit -m "Initial commit"
   git push -u origin master

  -------------------
 formated with emacs
 
 
*/
#include <stdio.h>
#include <stdlib.h> // malloc and  ISO C Random Number Functions */
#include <math.h>
#include <complex.h> // https://stackoverflow.com/questions/6418807/how-to-work-with-complex-numbers-in-c
#include <string.h> // strncat

complex double C = -0.122561166876654  +0.744861766619744*I; //   period = 3 center
double complex beta; // repelling fixed point
 
/* screen ( integer) coordinate */

const int iWidth  = 1000; 
const int iHeight = 1000;

/* world ( double) coordinate = parameter plane*/
// double complex Z =  Zx + Zy*I ;
double ZxMin;
double ZxMax;
double ZyMin;
double ZyMax;

/* */
double PixelWidth; //=(ZxMax-ZxMin)/iWidth;
double PixelHeight; // =(ZyMax-ZyMin)/iHeight;

/* color component ( R or G or B) is coded from 0 to 255 */
/* it is 24 bit color RGB file */
int ColorBytes = 3; // 3*8 = 24 bit color 
#define iRed 1 
#define iWhite  2

/* iterations  */
int  IterationMax = 10 ;

/* bail-out value , radius of circle ;  */
const double EscapeRadius=2.0;
double m = 2.0; //  multiplier

double TwoPi=2.0*M_PI;

double eps = 0.0000000001; // maximal error in numerical  computing 
       
// memmory virtual 1D array 
unsigned char *data;       
size_t MemmorySize;

/* 
   gives position ( index) in 1D virtual array  of 2D point (iX,iY) ;
   uses also global variable iWidth 
   without bounds check !!
*/
int f(int ix, int iy) // Give_k
{ return ColorBytes*(ix + iy*iWidth); }

/* 
   gives position ( index) in 1D virtual array  of 2D point Z 
   without bounds check !!
*/
int fd(complex double Z){ // double version of Give_k
  /* translate from world to screen coordinate */

  //  iY=(ZyMax-Zy)/PixelHeight; /*  */
  int ix=(creal(Z)-ZxMin)/PixelWidth;
  int iy=(ZyMax - cimag(Z))/PixelHeight; /* reverse Y  axis */		

	
  return f(ix,iy);

}

       
void GiveGrayColor(double position , int k, unsigned char A[]) 
{
  unsigned char b = 255*position;
  A[k]   = b;
  A[k+1] = b;
  A[k+2] = b;

}

void ColorPixel(int iColor, int k, unsigned char A[])
{
  switch (iColor)
    {
    case iRed:    A[k]   = 255; A[k+1] = 0;   A[k+2] = 0;   break;
    case iWhite : A[k]   = 255; A[k+1] = 255; A[k+2] = 255; break;
    }
}

void ColorPixel_d(complex double Z, int iColor,  unsigned char A[]){
  int k = fd(Z); // compute index of 1D array
  ColorPixel(iColor, k, A);
	

}

        
        
        
double complex give_z(int iX, int iY){
  double Zx, Zy;
  Zy = ZyMax - iY*PixelHeight; // inverse y axis
  
  Zx = ZxMin + iX*PixelWidth;
   
  return Zx + Zy*I;
 
 
}

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

int ComputeAndSavePixelColor(int iX, int iY){
 
  
  
  int i=0; // iteration
  double complex Z; // initial value for iteration Z0
  int k; // index of the 1D array
  
  // argument of complex numer in turns
  double tf; //final turn = z after IterationMax forward iterations   
  
  
  // index of 1D memory array
  k = f(iX, iY);   
  
  Z = give_z(iX, iY);
  
      
        
  // forward iteration with fixed number of iterations
  // without bailout test
  for  (i=0;  i<IterationMax; i++)
    Z=Z*Z+C; // forward iteration 
  
  tf =  GiveTurn(Z); // compute angle in turns of final z = final turn = tf
  GiveGrayColor(tf,k,data);  // color pixel proportionaly to tf
  	
  return 0;
}
 
 
 
void DrawFatouSet(){

  int iX,iY; // screen = integer coordinate in pixels

  // fill the array = render image = scanline 2D  of virtual 2D array 
  for(iY=0;iY<iHeight;iY++)
    for(iX=0;iX<iWidth;iX++)
      ComputeAndSavePixelColor(iX, iY); 
      	

}

double complex GiveBeta(double complex c){

  double complex delta=1-4*c;
  return (1.0-csqrt(delta))/2.0;

}

 
// IIM/J 
// https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/q-iterations#Julia_set_by_IIM.2FJ
void DrawJuliaSet(){

  complex double Z;
  int i; 
  beta = GiveBeta(C);
 
  Z= beta;
 
  for (i=0; i<100000000; i++){
	        
    ColorPixel_d(Z, iRed, data);
    // inverse iteration 
    Z = csqrt(Z-C);
    // random choose of preimage 
    if (rand()<(RAND_MAX/2)) 
      Z= -Z; 
 
 
 
  }

}
 
 
 
int setup(){

  ZxMin= -1.5*EscapeRadius;
  ZxMax=  1.5*EscapeRadius;
  ZyMin= -1.5*EscapeRadius;
  ZyMax=  1.5*EscapeRadius;
 
  //
  PixelWidth=(ZxMax-ZxMin)/iWidth;
  PixelHeight=(ZyMax-ZyMin)/iHeight;
  //
    
  //
  MemmorySize = iWidth * iHeight * ColorBytes * sizeof (unsigned char);	// https://stackoverflow.com/questions/492384/how-to-find-the-sizeof-a-pointer-pointing-to-an-array
        
  /* create dynamic 1D arrays for colors   */
  data = malloc (MemmorySize);
  if (data == NULL )
    { fprintf (stderr, " Error: Could not allocate memory !!! \n"); return 1; }

  printf (" No errors. End of setup \n");
  return 0;

}
 
 
 
 
 
// save dynamic "A" array to pgm file 
int SaveArray_2_PPM_file (unsigned char A[], int iMax)
{

  FILE *fp;
  const unsigned int MaxColorComponentValue = 255;	/* color component is coded from 0 to 255 ;  it is 8 bit color file */
  
  char name [100]; /* name of file */
  snprintf(name, sizeof name, "%d", iMax); /*  */
  char *filename =strncat(name,".ppm", 4);
  
  
  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, "P6\n %s\n %u %u\n %u\n", comment, iWidth, iHeight, MaxColorComponentValue);	/*write header to the file */
  fwrite (A, MemmorySize, 1, fp);	/*write A array to the file in one step */
  printf ("File %s saved. \n", filename);
  fclose (fp);

  return 0;
}

 
 
void CreateImage(int iMax){

  IterationMax = iMax;
  DrawFatouSet();
  DrawJuliaSet();
  
  
  
      	
  SaveArray_2_PPM_file (data, iMax);     	  
} 
 
 
 
void info(){

  printf(" Dynamic plane ( c plane) with Julia and Fatou sets for complex quadratic polynomial fc(z) = z^2 + c\n ");
  printf(" Rectangle part of 2D dynamic plane: corners: \n ZxMin = %f;   ZxMax = %f;  ZyMin = %f; ZyMax = %f \n ", ZxMin, ZxMax, ZyMin, ZyMax);
  printf(" center and radius: \n CenterX = %f;   CenterY = %f;  radius = %f\n ", (ZxMax+ZxMin)/2.0, (ZyMax+ZyMin)/2.0, fabs(ZyMax-ZyMin)/2.0);
  printf(" Mag = zoom = %f\n ",  2.0/fabs(ZyMax-ZyMin));
  printf(" PixelWidth = %f and PixelHeight =%f\n", PixelWidth, PixelHeight);
  printf(" parameter c = %f ; %f \n ", creal(C), cimag(C));
  printf(" beta : repelling fixed point z = %f ; %f \n ", creal(beta), cimag(beta));
  printf(" Escape Radius = %f\n ", EscapeRadius);
  printf(" Iteration Max = %d\n ", IterationMax);

} 
 
 
 
void close(){
 
  info(); 
  free(data); 
}
 
 
 
int main()
{

  int i;
 
  setup();   
     
  for (i=0; i< 11; i++)
    CreateImage(i); 
  	
  	
  	
  close();
  
        
  return 0;
}