/* c console program which creates ppm file */
#include <stdio.h>
#include <stdlib.h> /* for ISO C Random Number Functions */
#include <math.h>

#define PI 3.14159265


 
/*  gives sign of number */
double sign(double d)
{
       if (d<0)
       {return -1.0;}
       else {return 1.0;};
};

 



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

 int main()
{   /* 
    epsilon:0.001;
    radius:1+epsilon;
    angle:(3-sqrt(5))/2; 
    l:radius * exp(%i*angle*2*%pi);
    c:l/2 - (l*l)/4;
    c =0.58762398718188*%i-0.39095328938716
    -------------------------
    epsilon:-0.001;
    0.58595232559771*%i-0.39012849476251
    --------------------------------
    epsilon:0.0;
    0.58678790734697*%i-0.3905408702184
    -------------------------
    trifurcation angle:1/3;
    
    
    -------------------
    epsilon:0.0;
    radius:1+epsilon;
    angle:1/4;  tetrafurcation
    -0.4<y <0.4
    --------------------
    1/5 -0.2 <Zx  <0.4;
    0.0 < ZyM< 0.8
    -----------------------
    

*/
    double angle=1.0/6.0, 
           epsilon=-0.0005,
           radius=1.0+epsilon,
           /* L=Lx+Ly*i=lambda=radius * exp(i*angle*2*PI); */
           Lx=cos(2*PI*angle)*radius, 
           Ly=sin(2*PI*angle)*radius,
           /* C= Cx+Cy*i = L/2 - L*L/4  */
           Cx=-(cos(4*PI*angle)*radius*radius)/4+(cos(2*PI*angle)*radius)/2,
           Cy=((sin(2*PI*angle)*(epsilon+1))/2-(sin(4*PI*angle)*radius*radius)/4);
           
    
    /* screen coordinate = coordinate of pixels */      
    int iX, iY, 
        iXmin=0, iXmax=500,
        iYmin=0, iYmax=500,
        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  */
        period,
        index; /* of array */
        
   int iXinc, iYinc,iIncMax=1, /* width of big pixel;  for iXmax=2000 chose 6, for 500 choose 2 */    
       iIncMax2;
   /* world ( double) coordinate = parameter plane*/
    const double ZxMin=-0.2;
    const double ZxMax=0.6;
    const double ZyMin=0.0;
    const double ZyMax=0.8;
    /* */
    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,
           Z_cr_x=0.0, Z_cr_y=0.0; /* critical point */
     
        
        /*  */
        int Iteration,
            IterationMax=100000,
            iTemp;
     /* bail-out value , radius of circle ;  */
    const int EscapeRadius=40;
    int ER2=EscapeRadius*EscapeRadius;
    double  /*  AR= minimal distance from attractor = Attractor Radius */
           AR2=1.0e-6, /* AR2=AR*AR */
           d,dX,dY; /*  distance from attractor : d=sqrt(dx*dx+dy*dy) */

    
    /* PPM file */
    FILE * fp;
    char *filename="co_a_1_6_m00005_5.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;
    
    
    
    
    

   printf(" angle= %f\n",angle);
   printf(" epsilon= %f\n",epsilon);   
   printf(" radius= %f\n",radius);
     
   printf(" Lx= %f\n",Lx);
   printf(" Ly= %f\n",Ly);
   printf(" Cx= %f\n",Cx);  
   printf(" Cy= %f\n",Cy);  
   printf(" ZxMax-ZxMin= %f\n",ZxMax-ZxMin);
   printf(" ZyMax-ZyMin= %f\n",ZyMax-ZyMin);

   /*-----------------*/
    
    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 black points */       
      for(index=0;index<iLength-1;++index) array[index]=0;
      /* ---------------------------------------------------------------*/
      
  /*------ draw orbit of critical point  ---------------------------*/
  /* Z0 = Z_critical */
   Zx=Z_cr_x;
   Zy=Z_cr_y;
  
   Zx2=Zx*Zx;
   Zy2=Zy*Zy;
  
   //dX=Zx-AlphaX;
//    dY=Zy-AlphaY;
    d=Zx2+Zy2;
  //iIncMax2=0;
   for (Iteration=0;Iteration<IterationMax ;Iteration++) /* && (d>AR2) */
   {
    Zy=2*Zx*Zy + Cy;
    Zx=Zx2-Zy2 +Cx;
    
    /* translate from world to screen coordinate */
    iX=(Zx-ZxMin)/PixelWidth;
    iY=(Zy-ZyMin)/PixelHeight; /*  */
    //for(iYinc=-iIncMax;iYinc<iIncMax2;++iYinc)
//     for(iXinc=-iIncMax;iXinc<iIncMax2;++iXinc)
     { 
      iTemp=((iYmax-iY-1)*iXmax+iX)*3; //((iYmax-iY-1+iYinc)*iXmax+iX+iXinc)*3;
      if(iTemp<=iLength && iTemp>-1) // safety                                         
      {array[iTemp]=255; /* white point */
      array[iTemp+1]=255;
      array[iTemp+2]=255; 
      } 
      else    printf(" out\n");
     }
    /* ---------*/
    Zx2=Zx*Zx;
    Zy2=Zy*Zy;
    //dX=Zx-AlphaX;
//    dY=Zy-AlphaY;
//    d=dX*dX+dY*dY;
      d=Zx2+Zy2;
   };  
    
    
   
 
    
    

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