/* 
 c console program:

  1. draws Filled-in Julia setfor Fc(z)=z*z +c
  using basic algorithm ( boolean escape time )
 -------------------------------         
 2. technic 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)
 ---------------------------------
 I think that creating graphic can't be simpler

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


 
 */
 #include <stdio.h>
#include <math.h>


 int main()
 {
  const double Cx=-1.0;
  const double Cy=0.1;
  /* screen ( integer) coordinate */
  int iX,iY;
  const int iXmax = 2000; 
 const int iYmax = 2000;
 /* world ( double) coordinate = parameter 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;
 /* 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="f1.ppm";
 char *comment="# ";/* comment should start with # */
 static unsigned char color[3];
 double Zx, Zy,    /* Z=Zx+Zy*i   */
 Z0x, Z0y,  /* Z0 = Z0x + Z0y*i */
 Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
 /*  */
 int Iteration;
 const int IterationMax=2000;
 /* bail-out value , radius of circle ;  */
 const int EscapeRadius=400;
 int ER2=EscapeRadius*EscapeRadius;
 /*---------------------------------------------------------------*/
 /*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);
 /* compute and write image data bytes to the file*/
 for(iY=0;iY<iYmax;++iY)
    {
     Z0y=ZyMax - iY*PixelHeight; /* reverse Y  axis */
     if (fabs(Z0y)<PixelHeight/2) Z0y=0.0; /*  */    
     for(iX=0;iX<iXmax;++iX)
       {    /* initial value of orbit Z0 */
        Z0x=ZxMin + iX*PixelWidth;
        /* Z = Z0 */
       Zx=Z0x;
       Zy=Z0y;
       /* */         
       Zx2=Zx*Zx;
       Zy2=Zy*Zy;
       /* */
       for (Iteration=0;Iteration<IterationMax && ((Zx2+Zy2)<ER2);Iteration++)
          {
             Zy=2*Zx*Zy + Cy;
             Zx=Zx2-Zy2 +Cx;
             Zx2=Zx*Zx;
             Zy2=Zy*Zy;
          };
     /* compute  pixel color (24 bit = 3 bajts) */
    if (Iteration==IterationMax)
    { /*  interior of Filled-in Julia set  = black */
      color[0]=0;
      color[1]=0;
      color[2]=0;                           
    }
    else 
    { /* exterior of Filled-in Julia set  = white */
      color[0]=255; /* Red*/
      color[1]=255;  /* Green */ 
      color[2]=255;/* Blue */
    };
  /*write color to the file*/
  fwrite(color,1,3,fp);
  }
 }
 fclose(fp);
 return 0;
 }