/* 
  c console program
  Adam Majewski 
   fraktal.republika.pl 
   1. draws  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)
   ---------------------------------
   http://fraktal.republika.pl/cpp_argphi.html

*/
 
#include <stdio.h>

// this function is based on cpp function   
// in the file mandelXXsrc.zip .
// from program mandel by Wolf Jung
// http://www.mndynamics.com/indexp.html

int GiveLastIteration(double Zx, double Zy,double Cx ,double Cy , int iter_max,  double bailout2)
{
  // z= x+ y*i   
  long double x ,y ;
  //
  long double  u,  v;
  int iter; // iteration
 
  iter=0;//there was no iteration
 
  y =Zy; 
  x =Zx; 
  u = x*x ; 
  v = y*y;
  if ( u + v > bailout2 ) return  iter;  // point is in target set =no need to iterate 
  //
  do
    {
      // Fc(z)= z*z +c
      y = 2 * x * y + Cy;
      x = u - v + Cx; 
      u = x*x; 
      v = y*y;
      iter+=1;     
    } 
  while (( u + v <= bailout2 ) && iter<iter_max) ;
  return iter;
 
}





/* this function is based on function by W Jung */

double jlogphi(double zx0, double zy0, double cx, double cy)
{ 
  int j; 
  double 
    zx=zx0,
    zy=zy0,
    s = 0.5, 
    zx2=zx*zx,
    zy2=zy*zy,
    t;
  //
  for (j = 1; j < 400; j++)
    { s *= 0.5; 
      zy = 2 * zx * zy + cy;
      zx = zx2 - zy2 + cx; 
      zx2 = zx*zx; 
      zy2 = zy*zy;
      t = fabs(zx2 + zy2); // abs(z)
      if ( t > 1e24) break; 
    } 
  return s*log2(t);  // log(zn)* 2^(-n)
}//jlogphi



int main()
{
  const double Cx=-0.12,Cy=0.665;
        
       
  /* screen ( integer) coordinate */
  int iX,iY;
  const int iXmax = 1000, iXmin=0; 
  const int iYmax = 1000, 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 = parameter plane*/
  const double length=2.0;
  const double ZxMin=-length;
  const double ZxMax=length;
  const double ZyMin=-length;
  const double ZyMax=length;
  /* */
  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="p__10.pgm";
  char *comment="# Cx=-0.12, Cy=0.665; EscapeRadius=3 IterationMax=4000;";/* comment should start with # */
       
  /* Z=Zx+Zy*i  ;   Z0 = 0 */
  double Zx, Zy;
  double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
  /*  */
  double potential,temp;
  int Level, PreviousLevel;
  int LastIteration;
  const int IterationMax=4000;
  /* bail-out value , radius of circle ;  */
  const double EscapeRadius=3;
  double ER2=EscapeRadius*EscapeRadius;
  /* dynamic 1D array for 8-bit color values */    
  unsigned char *array;
  unsigned char color;
  /*-------------------------------------------------------------------*/
  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;
    }
  /* ---------------------------------------------------------------*/
  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;
	  //
	  LastIteration=GiveLastIteration(Zx, Zy, Cx , Cy , IterationMax, ER2);
	  if (LastIteration!=IterationMax)
	    /*   */
	    {    potential=jlogphi(Zx, Zy,Cx,Cy);// GivePotential(Zx, Zy,Cx,Cy,ER2,IterationMax );
	      // printf(" potential= %f\n",potential);
	      temp=255*potential;
	      if (temp>255) 
		color=255;
	      else color= temp; 
	      array[((iYmax-iY-1)*iXmax+iX)]= color ;
	    }
	}
    }
        
  /* 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,"OK : file  "); fprintf(stderr,filename); fprintf(stderr,"  saved.\nDone !!!");
      getchar();
    }
  free(array);
       
  getchar();
  return 0;
}