/* 
   c console  program:
   1. draws Julia setfor Fc(z)=z*z +c
   using DEM/J algorithm ( Distance Esthimation Method for Julia set )
   -------------------------------         
   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
   comments : Adam Majewski 

   gcc d.c -lm

   it creates a.out file. Then run it :

   ./a.out




*/



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



int GiveEscapeTime(double C_x, double C_y, int iMax, double _ER2)
{ 
  int i;
  double Zx, Zy;
  double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
 
  Zx=0.0; /* initial value of orbit = critical point Z= 0 */
  Zy=0.0;
  Zx2=Zx*Zx;
  Zy2=Zy*Zy;
 
  for (i=0;i<iMax && ((Zx2+Zy2)<_ER2);i++)
    {
      Zy=2*Zx*Zy + C_y;
      Zx=Zx2-Zy2 +C_x;
      Zx2=Zx*Zx;
      Zy2=Zy*Zy;
    };
  return i;
}

 


/*
  estimates distance from point c to nearest point in Julia  set 
  for Fc(z)= z*z + c
  z(n+1) = Fc(zn)  
  this function is based on function  mndlbrot::dist  from  mndlbrot.cpp
  from program mandel by Wolf Jung (GNU GPL )
  http://www.mndynamics.com/indexp.html 

  Hyunsuk Kim  : 
  For Julia sets, z is the variable and c is a constant. Therefore df[n+1](z)/dz = 2*f[n]*f'[n] -- you don't add 1.

  For the Mandelbrot set on the parameter plane, you start at z=0 and c becomes the variable. df[n+1](c)/dc = 2*f[n]*f'[n] + 1. 


*/
double jdist(double Zx, double Zy, double Cx, double Cy ,  int iter_max)
{ 
  int i;
  double x = Zx, /* Z = x+y*i */
    y = Zy, 
    /* Zp = xp+yp*1 = 1  */
    xp = 1, 
    yp = 0, 
    /* temporary */
    nz,  
    nzp,
    /* a = abs(z) */
    a; 
  for (i = 1; i <= iter_max; i++)
    { /* first derivative   zp = 2*z*zp  = xp + yp*i; */
      nz = 2*(x*xp - y*yp) ; 
      yp = 2*(x*yp + y*xp); 
      xp = nz;
      /* z = z*z + c = x+y*i */
      nz = x*x - y*y + Cx; 
      y = 2*x*y + Cy; 
      x = nz; 
      /* */
      nz = x*x + y*y; 
      nzp = xp*xp + yp*yp;
      if (nzp > 1e60 || nz > 1e60) break;
    }
  a=sqrt(nz);
  /* distance = 2 * |Zn| * log|Zn| / |dZn| */
  return 2* a*log(a)/sqrt(nzp); 
}
/* ------------------------------------------------------*/
int main(void)
{
  const double Cx=-0.75;
  const double Cy=0.14;
  /* screen ( integer) coordinate */
  int iX,iY;
  const int iXmax = 10000; 
  const int iYmax = 10000;
  /* world ( double) coordinate = parameter plane*/
  const double ZxMin=-2.;
  const double ZxMax=2.0;
  const double ZyMin=-2.0;
  const double ZyMax=2.0;
  /* */
  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="demj-.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;
  double distanceMax=PixelWidth/15; /*jdist( 0,0,Cx,Cy, IterationMax);*/
  /*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 bytes) */
	  if (Iteration==IterationMax)
	    { /*  interior of Julia set  = black */
	      color[0]=0;
	      color[1]=0;
	      color[2]=0;                           
	    }
	  else /* exterior of Filled-in Julia set  =  */
	    {  double distance=jdist(Z0x,Z0y,Cx,Cy,IterationMax);
	      if (distance<distanceMax)
		{ /*  Julia set  = white */
		  color[0]=255; /* Red*/
		  color[1]=255;  /* Green */ 
		  color[2]=255;/* Blue */
		}
	      else 
		{ /*  exterior of Julia set  = black */
		  color[0]=0;
		  color[1]=0;
		  color[2]=0;                           
		};
	    }
	  /* check the orientation of Z-plane */
	  /* mark first quadrant of cartesian plane*/     
	  /*     if (Z0x>0 && Z0y>0) color[0]=255-color[0];  */
	  /*write color to the file*/
	  fwrite(color,1,3,fp);
	}
    }
  fclose(fp);
  printf("file saved ");
  getchar();
  return 0;
}