/* 
c console program
 
 comments : Adam Majewski 
 fraktal.republika.pl 
*/
/* 
c console program:
 1. draws Level curves of escape time  for Fc(z)=z*z +c
 
 
 -------------------------------         
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 8 bit color ( gray scale ) graphic file ,  portable pixmap file = PGM  (P5)
see http://en.wikipedia.org/wiki/Portable_pixmap
to see the file use external application ( graphic viewer)
---------------------------------
3. 
Algorithm of drawing level curves is based on paper :
Drawing the Mandelbrot set by the method of escape lines. M. Romera et al.
http://www.iec.csic.es/~miguel/Preprint3.ps
this is translations of  BASIC program of M. Romera)
 */
#include <stdio.h>
#include <math.h>
int main()
{
        const double Cx=-0.12256,Cy=0.74486;
        
       
         /* screen ( integer) coordinate */
       int iX,iY;
       const int iXmax = 10000, iXmin=0; 
       const int iYmax = 10000, iYmin=0;
       int iWidth=iXmax-iXmin+1,
       iHeight=iYmax-iYmin+1,
       /* number of bytes = number of pixels of image * number of bytes of color */
       iLength=iWidth*iHeight*1,/* 1 bytes of color  */
       index; /* of array */
       
       /* world ( double) coordinate = dynamic plane ( z-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 8 bit color RGB file */
       const int MaxColorComponentValue=255; 
       FILE * fp;
       char *filename="LCMJ_rabbit.pgm";
       char *comment="# Cx=-0.12256, Cy=0.74486; EscapeRadius=1000 IterationMax=200;";/* comment should start with # */
       
       /* Z=Zx+Zy*i  ;   Z0 = 0 */
       double Zx, Zy;
       double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
       /*  */
       int Iteration, PreviousIter;
       const int IterationMax=200;
       /* bail-out value , radius of circle ;  */
       const double EscapeRadius=1000;
       double ER2=EscapeRadius*EscapeRadius;
       /* dynamic 1D array for 8-bit color values */    
       unsigned char *array;
       /*-------------------------------------------------------------------*/
       array = malloc( iLength * sizeof(unsigned char) );
       if (array == NULL)
       {
       fprintf(stderr,"Could not allocate memory");
       getchar();
       return 1;
       }
       else 
            {   fprintf(stderr,"I'm working. Please wait (:-))\n ");
                /* fill the data array with white points */       
            for(index=0;index<iLength-1;++index) array[index]=255;
            }
       /* ---------------------------------------------------------------*/
       
       
       
       
       /* first coat of paint */
       for(iY=0;iY<iYmax;iY++)
       {
            
            
            for(iX=0;iX<iXmax;iX++)
            {         /* compute Zx and Zy for each point */
                       Zy=ZyMin + iY*PixelHeight;
                       if (fabs(Zy)< PixelHeight/2) Zy=0.0; /*  */
                       Zx=ZxMin + iX*PixelWidth;
                       /* initial value of orbit  */
                       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;
                       };
                       /* plot point of Level Curve */
                       if (iX!=0 && Iteration!=PreviousIter)
                       { 
                          array[((iYmax-iY-1)*iXmax+iX)]=0;
                          PreviousIter=Iteration;                       
                       }
                 
                        
               }
       }
       
       /* second coat of paint */
       for(iX=0;iX<iXmax;iX++)
            {         
            
            for(iY=0;iY<iYmax;iY++)
              {/* compute Zx and Zy for each point */
               Zx=ZxMin + iX*PixelWidth;
               Zy=ZyMin + iY*PixelHeight;
               if (fabs(Zy)< PixelHeight/2) Zy=0.0; /*  */
            
            
                       /* initial value of orbit = Z */
                       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;
                       };
                       /* plot point of Level Curve */
                       if (iX!=0 && Iteration!=PreviousIter)
                       { 
                          array[((iYmax-iY-1)*iXmax+iX)]=0;
                          PreviousIter=Iteration;                       
                       }
                 
                        
               }
       }
       
       
       /* write the whole data array to ppm file in one step */      
     /*create new file,give it a name and open it in binary mode  */
     fp= fopen(filename,"wb"); /* b -  binary mode */
     if (fp == NULL){ fprintf(stderr,"file error"); }
           else
           {
           /*write ASCII header to the file*/
           fprintf(fp,"P5\n %s\n %d\n %d\n %d\n",comment,iXmax,iYmax,MaxColorComponentValue);
           /*write image data bytes to the file*/
           fwrite(array,iLength ,1,fp);
           fclose(fp);
           fprintf(stderr,"file saved\n");
           getchar();
           }
     free(array);
       
       getchar();
       return 0;
}