/* 
   c program:
   m_int
   = Integer Escape Time - Palette gradient  Mandelbrot set 
   = LSM/M = Level Set Method for Mandelbrot set
   
   detection of interior :
   "if c is a hyperbolic parameter then its orbit tends to an attracting cycle 
   and therefore the derivative ( with respect to z)  will tend to 0. 
   The converse does not exactly hold, but the counterexample values of c are sparse."
   
   
   
   
   https://en.wikibooks.org/wiki/Color_Theory/Color_gradient#Palette
   
   --------------------------------
   1. draws Mandelbrot set for Fc(z)=z*z +c
   using Mandelbrot algorithm ( boolean escape time )
   -------------------------------         
   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)
-----
 it is example  for : 
 https://www.math.univ-toulouse.fr/~cheritat/wiki-draw/index.php/Mandelbrot_set
 
 -------------
 compile : 

 
 
   gcc m_int.c -lm -Wall
 
 
   ./a.out
   
   
   -------- git --------
   
   
   cd existing_folder
git init
git remote add origin git@gitlab.com:adammajewski/mandelbrot_wiki_ACh.git
git add m_int.c
git commit -m "detection of interior"
git push -u origin master



 ==============
time ./a.out
with detection versus without detection
 real	0m0.383s versus 0m8.692s

so it is 22,6945169713 faster !!!!
 
*/
#include <stdio.h>
#include <math.h>
#include <complex.h> // https://stackoverflow.com/questions/6418807/how-to-work-with-complex-numbers-in-c
 
 

 
 
 
 
/* screen ( integer) coordinate */
int iX,iY;
const int iXmax = 8000; 
const int iYmax = 6000;
/* world ( double) coordinate = parameter plane*/
double Cx,Cy;
const double CxMin=-2.5;
const double CxMax=1.5;
const double CyMin=-1.5;
const double CyMax=1.5;
/* */
double PixelWidth; //=(CxMax-CxMin)/iXmax;
double PixelHeight; // =(CyMax-CyMin)/iYmax;
/* color component ( R or G or B) is coded from 0 to 255 */
/* it is 24 bit color RGB file */
const int MaxColorComponentValue=255; 
FILE * fp;
char *filename="M_int.ppm"; 
char *comment="# ";/* comment should start with # */
        
static unsigned char color[3]; // 24-bit rgb color

// 1D array = palette 
// default.map from fractint = VGA
unsigned char palette[3*8]= 
{0,0,255, // first color p=0
111,31,202, // second color p=1 is 
222,62,149, // green
0,188,188,
77,93,96,
188,124,43,
43,155,246,
154,186,193
};

/* Z=Zx+Zy*i  ;   Z0 = 0 */
//double Zx, Zy;
//double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
/*  */
int Iteration;
const int IterationMax=10000;
/* bail-out value , radius of circle ;  */
const double EscapeRadius=10.0;
const double eps=0.001;
        
        
        
        
 double complex give_c(int iX, int iY){
  double Cx,Cy;
  Cy=CyMin + iY*PixelHeight;
  if (fabs(Cy)< PixelHeight/2) Cy=0.0; /* Main antenna */
  Cx=CxMin + iX*PixelWidth;
   
  return Cx+Cy*I;
 
 
}
 

// print 24-bit color from 1D array = palette 
void test_palette(unsigned char p){
printf("p = %d  color :  R= %d G= %d B = %d  \n", p, palette[p*3], palette[p*3+1], palette[p*3+2]);


}


// gives last iterate = escape time
// output 0< i < iMax
 int iterate(double complex C , int iMax)
  {
   int i=0;
   double complex Z= C; // initial value for iteration Z0
   complex double D = 1.0; // derivative with respect to z 
   
   for(i=0;i<iMax;i++)
    { if(cabs(Z)>EscapeRadius) break; // exterior
      if(cabs(D)<eps) break; // interior
      D = 2.0*D*Z;
      Z=Z*Z+C; // complex quadratic polynomial
      
    }
   return i; 
 }
 
 
 
 
 
int compute_color(complex double c, unsigned char color[3]){
 
   int i; // last iteration
   int p=6; // index of palette   
   
   
   // compute escape time = last iteration:  1<i <= IterationMax
   i = iterate( c, IterationMax);
  
    
      
     i= i-1;
     if (i<3) p = 0;
       else {if (i <10) p = 1;
               else {if (i < 31) p = 2; 
                      else if (i<100) p=3; 
                             else { if (i<300) p=4; 
                                    else { if (i<900) p = 5; }
                             } 
                      }
                      
             }
     
     
     
     p = 3*p; // map to [0 : length(paleette) = 3*255 ] range, 
        
     
     color[0]=palette[p];  /* Red*/
     color[1]=palette[p+1];  /* Green */ 
     color[2]=palette[p+2];  /* Blue */
                            
     
 
   
  return 0;
}
 
 
 
 void setup(){
 
  //
  PixelWidth=(CxMax-CxMin)/iXmax;
  PixelHeight=(CyMax-CyMin)/iYmax;
        
         
  /*create new file,give it a name and open it in binary mode  */
  fp= fopen(filename,"wb"); /* b -  binary mode */
  /*write ASCII header to the file*/
  fprintf(fp,"P6\n %s\n %d\n %d\n %d\n",comment,iXmax,iYmax,MaxColorComponentValue);
 
 }
 





void info(){




 double distortion;
 // widt/height
 double PixelsAspectRatio = (double)iXmax/iYmax;  // https://en.wikipedia.org/wiki/Aspect_ratio_(image) 
 double WorldAspectRatio = (CxMax-CxMin)/(CyMax-CyMin);

 distortion = PixelsAspectRatio - WorldAspectRatio;
 printf("distortion = %.16f ( it should be zero !)\n", distortion );

 // file  
 printf("file %s saved. Compare M_int.svg \n", filename);

 
}
 













 
 void close(){
 
 
 
 fclose(fp);
 info(); 
 
 
 
 
 }
 
 
 
 
 
// ************************************* main ************************* 
int main()
{
        
  complex double c;
        
        
 
  setup();      
        
        
  printf(" render = compute and write image data bytes to the file \n");
 
  for(iY=0;iY<iYmax;iY++)
    for(iX=0;iX<iXmax;iX++)
      { // compute pixel coordinate        
	c = give_c(iX, iY);  
	/* compute  pixel color (24 bit = 3 bytes) */
	compute_color(c,color);         
	/*write color to the file*/
	fwrite(color,1,3,fp);
      }
        
  
  
  close();
  
         
  return 0;
}