/* 
 c console program:
 1. draws components of Filled-in Julia set for Fc(z)=z*z +c
 using   escape time
 it works except GivePeriod
 -------------------------------         
 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;};
 }
 void GiveColorFromList6(int i ,unsigned char c[])
 { /*  */
  switch( i){
   case 0: {
      c[0] = 210;
      c[1] = 0;
      c[2] = 0;
       break;
    };
   case 1: {
      c[0] = 0 ;
      c[1] = 0;
      c[2] = 255;
       break;
    };
   case 2: {
      c[0] = 0;
      c[1] = 255;
      c[2] = 0;
       break;
    };
   case 3: {
      c[0] = 0;
      c[1] = 255 ;
      c[2] = 255;
       break;
    };
   case 4: {
      c[0] = 255;
      c[1] = 0;
      c[2] = 255;
       break;
    };
   case 5: {
      c[0] = 255;
      c[1] = 0;
      c[2] = 255 ;
       break;
    };
   case 6: {
      c[0] = 164;
      c[1] = 136;
      c[2] = 88 ;
       break;
    }; 
   default: {
      c[0] = 100;
      c[1] = 100;
      c[2] = 0;
       break;
     }
   }
 } 
 /* ----------------------*/
 int main()
 {   // 4: -0.15652016683376e+000	-1.03224710892283;  0.28227139076691e+000	-0.53006061757853
 // 4: c=0.281+0.533, AR2=0.001, ER2=8*8, 
 //3: -0.12256116687665e+000	-0.74486176661974e
 // 2: -1
 // 1: 0.4*i
 const double Cx=0.281,
              Cy=0.533;
 /* screen coordinate = coordinate of pixels */      
 int period=4,
     p,
     iX, iY, 
     iXmin=0, iXmax=2000,
     iYmin=0, iYmax=2000,
     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   */
          tempZx,
          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,
          ZAx,ZAy, /* attractor  ZA = ZAx+ZAy*i */
          Z_cr_x=0.0, Z_cr_y=0.0; /* critical point */
  /*  */
  int Iteration,iIteration,
  IterationMax=80,
  IterationMaxBig=1000,
  iTemp;
  /* bail-out value , radius of circle ;  */
  const int EscapeRadius=8;
 int ER2=EscapeRadius*EscapeRadius;
 double //AR=PixelWidth, /* minimal distance from attractor = Attractor Radius */
 AR2=0.001, /* proportional to period  */
 d,dX,dY; /*  distance from attractor : d=sqrt(dx*dx+dy*dy) */
 /* PPM file */
 FILE * fp;
 char *filename="fatou_4b.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;
 static unsigned char color[3];
 /*  ---------  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.0;//0.5-SqrtDeltaX/2;
  AlphaY=0.0;//-SqrtDeltaY/2;
  AbsLambdaA=2*sqrt(AlphaX*AlphaX+AlphaY*AlphaY);
  AbsLambdaB=2*sqrt(BetaX*BetaX+BetaY*BetaY);
  /* -- find attractor  ZA = ZAx+ZAy*i ---------------*/
  /* Z = 0 = critical point */
  Zx=0.0;
  Zy=0.0;
  Zx2=Zx*Zx;
  Zy2=Zy*Zy;
  /* */
  for (Iteration=0;Iteration<IterationMaxBig && ((Zx2+Zy2)<ER2);Iteration++)
                   {
                       Zy=2*Zx*Zy + Cy;
                       Zx=Zx2-Zy2 +Cx;
                       Zx2=Zx*Zx;
                       Zy2=Zy*Zy;
                   };
  
  ZAx=Zx;
  ZAy=Zy;
  /*----*/
  //period= GivePeriod( Cx,Cy,ZAx, ZAy,ER2,AR2 ); 
  /* --------------- info -------------------------*/
  printf(" C= (%f , %f) \n",Cx,Cy);
  printf(" Fixed points : \n");
  printf(" Beta= %f , %f\n",BetaX,BetaY);
  printf(" Alpha= %f, %f\n",AlphaX,AlphaY);
  printf(" abs(Lambda (Alpha))= %f\n",AbsLambdaA);
  printf(" abs(lambda(Beta))= %f\n",AbsLambdaB);
  printf(" Limit distances  : \n");
  printf(" AR2= %f\n",AR2);
  printf(" ER2= %d\n",ER2);
  printf(" Attractor : \n");
  printf(" ZA= %f , %f\n",ZAx,ZAy);
  printf(" period= %d\n",period);
  /*--------------------------------------------------------*/
  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;
                           dY=Zy-AlphaY;
                           d=dX*dX+dY*dY;
                           for (iIteration=0;iIteration<IterationMax && (d>AR2);iIteration++)
                       {
                           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;
                       };
                       /* iLSM */
                     if (iIteration==IterationMax) 
                          { // points of interior which did not reached attractor = solid red
                            array[iTemp]=255; /* Red*/
                            array[iTemp+1]=0;  /* Green */ 
                            array[iTemp+2]=0;/* Blue */
                          }
                          else
                          { /* components of filled Julia set */
                           iIteration%=period;// modulo period
                           GiveColorFromList6(iIteration ,color);
                           array[iTemp]=color[0]; /* Red*/
                           array[iTemp+1]=color[1];  /* Green */ 
                           array[iTemp+2]=color[2];/* Blue */
                          }                      
                       }
             else /* exterior of Filled-in Julia set  */
             { /*  LSM */
               //GiveColorFromList6(Iteration ,color);
                           int i=Iteration*5;
                          {
                           array[iTemp]=(255-i)%255; /* Red*/
                           array[iTemp+1]=(215-i)%255;  /* Green */ 
                           array[iTemp+2]=(200-i)%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];  
    }
   } 
 /* -------------------  draw special points ----------------------------*/             
 /* ---  fixed points : Alpha and Beta ---------*/
 iX=(AlphaX-ZxMin)/PixelWidth; /*translate from world to screen coordinate */
 iY=(AlphaY-ZxMin)/PixelHeight; /*  */
 /* plot  big green pixel = 2*iIncMax  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 = 2*iIncMax 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;  
   }
   }          
 /* -------------- critical point ------------*/
 /* translate from world to screen coordinate */
 iX=(Z_cr_x-ZxMin)/PixelWidth;
 iY=(Z_cr_y-ZyMin)/PixelHeight; /*  */
 /* plot  big blue pixel = 2*iIncMax  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;  
   }
   } 
 /*----------- attracting cycle ---------------------*/
 Zx=ZAx;
 Zy=ZAy; 
 for(p=0;p<period;p++)
 {
  /* translate from world to screen coordinate */
 iX=(Zx-ZxMin)/PixelWidth;
 iY=(Zy-ZyMin)/PixelHeight; /*  */
 /* plot  big white pixel = 2*iIncMax  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]=0;  
   }
   }  
   /* compute new point */
   tempZx= Zx*Zx - Zy*Zy + Cx;
   Zy = 2*Zx*Zy + Cy; 
   Zx = tempZx;
 } 
 /* ----- 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 ... */
 }