/* 
   c console program
   1. draws  Julia setfor Fc(z)=z*z +c using :
   colors exterior of Julia set using modified decomposition
   dynamic 1D array for 24-bit color values
   -------------------------------         
   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 manual creating graphic can't be simpler
   -------------------------------------------------
   Adam Majewski 
   fraktal.republika.pl 

   ======================
   Linux console : 
   save as n.c
   to compile :
   gcc -lm n.c
   to run :
   ./a.out
   It creates ppm file. 
   To see it use image viewer
   To convert use Image Magic command line tools :
   convert j.ppm -resize 1000x1000 j.png
   Source code was formatted with emacs
*/

#include <stdio.h>
#include <stdlib.h> /* for ISO C Random Number Functions */
#include <math.h>



int main()
{      const double Cx=0.0,Cy=0.0;
     
     
  /* screen coordinate = coordinate of pixels */      
  int iX, iY, 
    iXmin=0, iXmax=10000,
    iYmin=0, iYmax=10000,
    iWidth=iXmax-iXmin+1,
    iHeight=iYmax-iYmin+1,
    /* 3D data : X , Y, color */
    /* number of bytes = number of pixels of image * number of bytes of color */
    iLength=iWidth*iHeight*3,/* 3 bytes of color  */
    index; /* of array */
        
  /*  int iXinc, iYinc,iIncMax=12;     */
   
  /* world ( double) coordinate = parameter plane*/
  const double ZxMin=-5;
  const double ZxMax=5;
  const double ZyMin=-5;
  const double ZyMax=5;
  /* */
  double PixelWidth=(ZxMax-ZxMin)/iWidth;
  double PixelHeight=(ZyMax-ZyMin)/iHeight;
  double Zx,Zy,    /* Z=Zx+Zy*i   */
    Z0x, Z0y,  /* Z0 = Z0x + Z0y*i */
    Zx2,Zy2; 
  /*  */
  int Iteration,
    IterationMax=600 , /*for  loop  */
    iTemp;
  /* bail-out value , radius of circle ;  */
  const int EscapeRadius=100;
    
  /* PPM file */
  FILE * fp;
  char *filename="m.ppm";
  char *comment="# this is julia set for c= 0";/* comment should start with # */
  const int MaxColorComponentValue=255;/* color component ( R or G or B) is coded from 0 to 255 */
  /* dynamic 1D array for 24-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 
    {         
      /* fill the data array with white points */       
      for(index=0;index<iLength-1;++index) array[index]=255;

      /* ---------------------------------------------------------------*/
      for(iY=0;iY<iYmax;++iY)
	{
	  Z0y=ZyMin + 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;
		 
	      /*----------- main loop without checking of abs(zn)  -------------*/
	      Zx=Z0x;
	      Zy=Z0y;
	      Zx2=Zx*Zx;
	      Zy2=Zy*Zy;
	      Iteration=0;
	      while (Iteration<IterationMax && (fabs(Zy)<EscapeRadius) ) 
		{
		  Zy=2*Zx*Zy + Cy;
		  Zx=Zx2-Zy2 +Cx;
		  Zx2=Zx*Zx;
		  Zy2=Zy*Zy;
		  Iteration+=1;
		}
              /* --------------- compute  pixel color (24 bit = 3 bajts) */
          
	      iTemp=((iYmax-iY-1)*iXmax+iX)*3;   
             
	      /* interior of Filled-in Julia set = green */
	      if (Iteration==IterationMax)
		{ 
		  array[iTemp]=0; /* Red*/
		  array[iTemp+1]=255;  /* Green */ 
		  array[iTemp+2]=0;/* Blue */
		} 
	      else
		/* exterior of Filled-in Julia set  */
		if (Zy>0 ) 
		  { 
		    array[iTemp]=255; /* Red*/
		    array[iTemp+1]=255;  /* Green */ 
		    array[iTemp+2]=255;/* Blue */
		  }
		else
		  {
		    array[iTemp]=0; /* Red*/
		    array[iTemp+1]=0;  /* Green */ 
		    array[iTemp+2]=0;/* Blue */    
		  };    
                       
	      /*----------------- check the orientation of Z-plane by marking first quadrant of cartesian plane ----- */
	      /*  if (Z0x>0 && Z0y>0) array[((iYmax-iY-1)*iXmax+iX)*3]=255-array[((iYmax-iY-1)*iXmax+iX)*3];   */
	    }
	} 

       
      /* ---------------------  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,"P6\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 %s saved\n",filename);
	  getchar();
	}
 
      free(array);
      return 0;
    } /* if (array ..  else ... */
}