/* 
 /* 
   c program:
   1. draws Filled-in Julia set for Fc(z)=z*z +c
   using  level sets of attraction 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)
   ---------------------------------
   3. interior of set is coloured using Level Set Method 
   but iteration is a measure of attraction of Zn to fixed 
   attractin point Z=Alpha
   ( second fixed point Z=Beta is used in IIM , not here )
   it works only if abs(2*Alpha)<1 
   it means when Alpha is attracting 
   ----------------------------------------------------
   I think that creating graphic can't be simpler
*/
#include <stdio.h>
#include <stdlib.h> /* for ISO C Random Number Functions */
#include <math.h>
/*  gives sign of number */
double sign(double d)
{
  if (d<0)
    {return -1.0;}
  else {return 1.0;};
};
/* ----------------------*/
int main()
{   
  const double Cx=0.0,Cy=0.5;
  /* 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=24;     
  /* world ( double) coordinate = parameter plane*/
  const double ZxMin=-1.5;
  const double ZxMax=1.5;
  const double ZyMin=-1.5;
  const double ZyMax=1.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, /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
           
    DeltaX, DeltaY,
    SqrtDeltaX, SqrtDeltaY,
    AlphaX, AlphaY,
    BetaX,BetaY, /* repelling fixed point Beta */
    AbsLambdaA,AbsLambdaB,
    Z_cr_x=0.0, Z_cr_y=0.0; /* critical point */
  /*  */
  int Iteration,
    IterationMax=5000,
    iTemp;
  /* bail-out value , radius of circle ;  */
  const int EscapeRadius=40;
  int ER2=EscapeRadius*EscapeRadius;
  double  /* AR=PixelWidth, minimal distance from attractor = Attractor Radius */
    AR2=1.0e-19, /* > AR*AR */
    d,dX,dY; /*  distance from attractor : d=sqrt(dx*dx+dy*dy) */
  /* PPM file */
  FILE * fp;
  char *filename="c.ppm";
  char *comment="# this is julia set for c= ";/* 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;
  /*  ---------  find fixed points ---------------------------------*/
  /* Delta=1-4*c */
  DeltaX=1-4*Cx;
  DeltaY=-4*Cy;
  /* SqrtDelta = sqrt(Delta) */
  /* sqrt of complex number algorithm from Peitgen, Jurgens, Saupe: Fractals for the classroom */
  if (DeltaX>0)
    {
      SqrtDeltaX=sqrt((DeltaX+sqrt(DeltaX*DeltaX+DeltaY*DeltaY))/2);
      SqrtDeltaY=DeltaY/(2*SqrtDeltaX);        
    }
  else /* DeltaX <= 0 */
    {
      if (DeltaX<0)
	{
          SqrtDeltaY=sign(DeltaY)*sqrt((-DeltaX+sqrt(DeltaX*DeltaX+DeltaY*DeltaY))/2);
          SqrtDeltaX=DeltaY/(2*SqrtDeltaY);        
	}
      else /* DeltaX=0 */
	{
	  SqrtDeltaX=sqrt(fabs(DeltaY)/2);
	  if (SqrtDeltaX>0) SqrtDeltaY=DeltaY/(2*SqrtDeltaX);
	  else SqrtDeltaY=0;    
	}
    };
  /* Beta=(1-sqrt(delta))/2 */
  BetaX=0.5+SqrtDeltaX/2;
  BetaY=SqrtDeltaY/2;
  /* Alpha=(1+sqrt(delta))/2 */
  AlphaX=0.5-SqrtDeltaX/2;
  AlphaY=-SqrtDeltaY/2;
  AbsLambdaA=2*sqrt(AlphaX*AlphaX+AlphaY*AlphaY);
  AbsLambdaB=2*sqrt(BetaX*BetaX+BetaY*BetaY);
  printf(" Cx= %f\n",Cx);
  printf(" Cy= %f\n",Cy); 
  printf(" Beta= %f , %f\n",BetaX,BetaY);
  //printf(" BetaY= %f\n",BetaY);
  printf(" Alpha= %f, %f\n",AlphaX,AlphaY);
  //printf(" AlphaY= %f\n",AlphaY);
  printf(" abs(Lambda (Alpha))= %f\n",AbsLambdaA);
  printf(" abs(lambda(Beta))= %f\n",AbsLambdaB);
  /*------ compute radius for finite attractor ---------------------------*/
  /* based on the code by Evgeny Demidov http://www.ibiblio.org/e-notes/MSet/Contents.htm */
  /* http://www.ibiblio.org/e-notes/MSet/Inside.htm */
  /* Z0 = Z_critical */
  Zx=Z_cr_x;
  Zy=Z_cr_y;
  Zx2=Zx*Zx;
  Zy2=Zy*Zy;
  dX=Zx-AlphaX;
  dY=Zy-AlphaY;
  d=dX*dX+dY*dY;
  for (Iteration=0;Iteration<IterationMax && (d>1e-6);Iteration++)
    {
      Zy=2*Zx*Zy + Cy;
      Zx=Zx2-Zy2 +Cx;
      Zx2=Zx*Zx;
      Zy2=Zy*Zy;
      dX=Zx-AlphaX;
      dY=Zy-AlphaY;
      d=dX*dX+dY*dY;
    };
  AR2=d; /* last distance from attractor of iteration of critical point*/
  printf(" AR2= %f\n",AR2);
  /*--------------------------------------------------------*/
  array = malloc( iLength * sizeof(unsigned char) );
  if (array == NULL)
    {
      fprintf(stderr,"Could not allocate memory");
      getchar();
      return 1;
    }
  else 
    {         
      printf(" I'm working. 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)
	{
          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;
              /* 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;
		};
	      iTemp=((iYmax-iY-1)*iXmax+iX)*3;        
	      /* compute  pixel color (24 bit = 3 bajts) */
	      if (Iteration==IterationMax)
		{ /*  interior of Filled-in Julia set  =  */
		  /* Z = Z0 */
		  Zx=Z0x;
		  Zy=Z0y;
		  Zx2=Zx*Zx;
		  Zy2=Zy*Zy;
		  dX=Zx-AlphaX; /* it uses fixed point as an attractor so it works only when period =1 */
		  dY=Zy-AlphaY;
		  d=dX*dX+dY*dY;
		  for (Iteration=0;Iteration<IterationMax && (d>AR2);Iteration++)
		    {
		      Zy=2*Zx*Zy + Cy;
		      Zx=Zx2-Zy2 +Cx;
		      Zx2=Zx*Zx;
		      Zy2=Zy*Zy;
		      dX=Zx-AlphaX;
		      dY=Zy-AlphaY;
		      d=dX*dX+dY*dY;
		    };
		  /* LSM */
		  if ((Iteration%2)==0) /* even or odd number */
		    { 
		      array[iTemp]=0; /* Red*/
		      array[iTemp+1]=0;  /* Green */ 
		      array[iTemp+2]=0;/* Blue */
		    }
		  else 
		    {
		      array[iTemp]=255; /* Red*/
		      array[iTemp+1]=255;  /* Green */ 
		      array[iTemp+2]=255;/* Blue */    
		    };                      
		}
              else 
		/* exterior of Filled-in Julia set  */
		/*  black solid color */
		{ 
		  array[iTemp]=0; /* Red*/
		  array[iTemp+1]=0;  /* Green */ 
		  array[iTemp+2]=0;/* Blue */
		}
              /* check the orientation of Z-plane */
	      /* mark 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];  
	    }
	} 
      /* draw fixed points ----------------------------------------------------*/             
      /* translate from world to screen coordinate */
      iX=(AlphaX-ZxMin)/PixelWidth;
      iY=(AlphaY-ZxMin)/PixelHeight; /*  */
      /* plot  big green pixel = 6 pixel wide */
      for(iYinc=-iIncMax;iYinc<iIncMax;++iYinc){
	for(iXinc=-iIncMax;iXinc<iIncMax;++iXinc)
	  { 
	    iTemp=((iYmax-iY-1+iYinc)*iXmax+iX+iXinc)*3;                                          
	    array[iTemp]=0;
	    array[iTemp+1]=255;
	    array[iTemp+2]=0;  
	  }
      }
      /* translate from world to screen coordinate */
      iX=(BetaX-ZxMin)/PixelWidth;
      iY=(BetaY-ZyMin)/PixelHeight; /*  */
      /* plot  big red pixel = 6 pixel wide */
      for(iYinc=-iIncMax;iYinc<iIncMax;++iYinc){
	for(iXinc=-iIncMax;iXinc<iIncMax;++iXinc)
	  {  
	    iTemp=((iYmax-iY-1+iYinc)*iXmax+iX+iXinc)*3;
	    array[iTemp]=255;
	    array[iTemp+1]=0;
	    array[iTemp+2]=0;  
	  }
      }          
      /* draw critical point */
      /* translate from world to screen coordinate */
      iX=(Z_cr_x-ZxMin)/PixelWidth;
      iY=(Z_cr_y-ZyMin)/PixelHeight; /*  */
      /* plot  big blue pixel = 6 pixel wide */
      for(iYinc=-iIncMax;iYinc<iIncMax;++iYinc){
	for(iXinc=-iIncMax;iXinc<iIncMax;++iXinc)
	  {  
	    iTemp=((iYmax-iY-1+iYinc)*iXmax+iX+iXinc)*3;
	    array[iTemp]=0;
	    array[iTemp+1]=0;
	    array[iTemp+2]=255;  
	  }
      }   
      /* 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 ... */
}