/* 
   c console program
   1. draws  Julia setfor Fc(z)=z*z +c using :
   IIM/J 
   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 e.c -lm -Wall
   to run :
   ./a.out
*/

#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;};
};

/*
  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 
*/

int main()

{      const double Cx=-1.0,Cy=0.0;
  /* screen coordinate = coordinate of pixels */      
  int iX, iY, 
    iXmin=0, iXmax=1000,
    iYmin=0, iYmax=1000,
    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=-1.7;
  const double ZxMax=1.7;
  const double ZyMin=-1.7;
  const double ZyMax=1.7;
  /* */
  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  */
    NewZx, NewZy,
    DeltaX, DeltaY,
    SqrtDeltaX, SqrtDeltaY,
    AlphaX, AlphaY,
    BetaX,BetaY, /* repelling fixed point Beta */
    AbsLambdaA,AbsLambdaB;
  /*  */
  int Iteration,
    IterationMax=6 , /*for modified loop  */
    iTemp;
    
  /* PPM file */
  FILE * fp;
  char *filename="6_basilica_.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 repelling fixed point ---------------------------------*/
  /* 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);

   

   

   

    

    

    

  /*-------------------------------------------------------------------*/

    

  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;

              /* Z = Z0 */

              Zx=Z0x;

              Zy=Z0y;

              

              Zx2=Zx*Zx;

              Zy2=Zy*Zy;

	      /*----------- modified loop without checking of abs(zn)  -------------*/

	      for (Iteration=0;Iteration<IterationMax;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) */

          

             

	      /* exterior of Filled-in Julia set  */

                         

	      /* binary decomposition  */

	      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];  

	    }

	} 

      /*-------------------- draw julia set using IIM/J ------------------------------------------*/ 

      /* initial value of orbit Z=Z0 is repelling fixed point */

      Zy=BetaY; 

      Zx=BetaX; 

              

      for (Iteration=0;Iteration<10000000;Iteration++)

	{

	  /* Zn*Zn=Z(n+1)-c */

	  Zx=Zx-Cx;

	  Zy=Zy-Cy;

	  /* sqrt of complex number algorithm from Peitgen, Jurgens, Saupe: Fractals for the classroom */

	  if (Zx>0)

	    {

	      NewZx=sqrt((Zx+sqrt(Zx*Zx+Zy*Zy))/2);

	      NewZy=Zy/(2*NewZx);        

	    }

	  else /* ZX <= 0 */

	    {

	      if (Zx<0)

		{

		  NewZy=sign(Zy)*sqrt((-Zx+sqrt(Zx*Zx+Zy*Zy))/2);

		  NewZx=Zy/(2*NewZy);        

		}

	      else /* Zx=0 */

		{

		  NewZx=sqrt(fabs(Zy)/2);

		  if (NewZx>0) NewZy=Zy/(2*NewZx);

		  else NewZy=0;    

               

		}

	    };

	  if (rand()<(RAND_MAX/2))

	    {   

	      Zx=NewZx;

	      Zy=NewZy; 

	    }

	  else {Zx=-NewZx;

	    Zy=-NewZy; }

	  /* translate from world to screen coordinate */

	  //  iX=(Zx-ZxMin)/PixelWidth;

	  //  iY=(ZyMax-Zy)/PixelHeight; /* reverse Y  axis */
	  iX=(Zx-ZxMin)/PixelWidth;

	  iY=(Zy-ZyMin)/PixelHeight; /*  */		

	  /* plot  pixel =  boundary of Filled-in Julia set  =  Julia set*/
	  iTemp=((iYmax-iY-1)*iXmax+iX)*3;     

	  array[iTemp]=255; /* Red*/

	  array[iTemp+1]=0;  /* Green */ 

	  array[iTemp+2]=0;/* Blue */

	};            

                       

       

      /* ---------------------  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);

            

	}
      free(array);
      return 0;

    } /* if (array ..  else ... */

  

}