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

Conversion to png is made with convert from ImageMagic
convert j.ppm -resize 2000x2000 j.png


*/

#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=0.0,Cy=0.0;
   /* 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=12;     */
/* 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  */
       NewZx, NewZy,
       DeltaX, DeltaY,
       SqrtDeltaX, SqrtDeltaY,
       AlphaX, AlphaY,
       BetaX,BetaY, /* repelling fixed point Beta */
       AbsLambdaA,AbsLambdaB;
   /*  */
     int Iteration,
     IterationMax=6 , /*for modified loop  */
     iTemp;
    /* bail-out value , radius of circle ;  */
  // const int EscapeRadius=100;
   // 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="j.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 */
                       }
                       if (Zy<0 )
                       {
                         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);
            //getchar();
            }
      free(array);
      return 0;
    } /* if (array ..  else ... */
}