/* 
 c program:
  1. draws Fatou set for Fc(z)=z*z 
  using  binary decomposition
 -------------------------------         
 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)
 ---------------------------------
 */
 #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.35,Cy=0.0; // Thx to Ingvar Kullberg http://klippan.seths.se/fractals/articles/3.pdf
     /* screen coordinate = coordinate of pixels */      
    int iX, iY, 
        iXmin=0, iXmax=20000,
        iYmin=0, iYmax=20000,
        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=6;     
    /* world ( double) coordinate = parameter plane*/
    const double ZxMin=-2.5;
    const double ZxMax=2.5;
    const double ZyMin=-2.5;
    const double ZyMax=2.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  */
           //NewZx, NewZy,
           DeltaX, DeltaY,
           SqrtDeltaX, SqrtDeltaY,
           AlphaX, AlphaY,
           BetaX,BetaY, /* repelling fixed point Beta */
           AbsLambdaA,AbsLambdaB;
          /*  */
        int Iteration,
            IterationMax=40000000,
            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=AR*AR,
           d,dX,dY; /*  distance from attractor : d=sqrt(dx*dx+dy*dy) */
        /* PPM file */
    FILE * fp;
    char *filename="40c.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);
   /* initial value of orbit Z=Z0 is repelling fixed point */
   Zy=BetaY; /*  */
   Zx=BetaX; 
   /*-------------------------------------------------------------------*/
    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;
             /* */
             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 bytes) */
           if (Iteration==IterationMax)
                    { /*  interior of Filled-in Julia set  =  */
                          /* Z = Z0 */
                          Zx=Z0x;
                          Zy=Z0y;
                          Zx2=Zx*Zx;
                          Zy2=Zy*Zy;
                          dX=Zx-AlphaX;
                            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;
                        };
                         /* */
                           if (Zy>AlphaY) 
                           { 
                             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  */
              /*  */
                       if (Zy>0) 
                       { 
                         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 */    
                       };                      
   /* 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];  
     }
    } 
 
 /* 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 saved");
            //getchar();
            }
    free(array);
    return 0;
    } /* if (array ..  else ... */
 }