/* 
 c program:
  1. draws Julia set for z^14-z
http://www.ijon.de/mathe/julia/node6.html
  using  escape time 
  aproximates julia set near fixed point z=0 using 13-arm star
 http://www.ijon.de/mathe/julia/node6.html
 -------------------------------         
 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 creating graphic can't be simpler

Maxima CAS session : 
(%i1) z:x+y*%i;
(%o1) %i*y+x
(%i2) z1:z^14-z;
(%o2) (%i*y+x)^14-%i*y-x
(%i3) realpart(z1);
(%o3) -y^14+91*x^2*y^12-1001*x^4*y^10+3003*x^6*y^8-3003*x^8*y^6+1001*x^10*y^4-91*x^12*y^2+x^14-x
(%i4) imagpart(z1);
(%o4) 14*x*y^13-364*x^3*y^11+2002*x^5*y^9-3432*x^7*y^7+2002*x^9*y^5-364*x^11*y^3+14*x^13*y-y



gcc z.c -Wall -lm

author : Adam Majewski
fraktal.republika.pl
 
 */
#include <stdio.h>
#include <stdlib.h> // malloc
#include <math.h>
#include <complex.h>
#include <string.h> // strcat
 
  const int iXmax = 1000; 
 const int iYmax = 1000;
 int iSize ; //= iXmax*iYmax; //
// target set for points falling into alfa fixed point
// is a circle around alfa fixed point
// with radius = AR
double AR ; // radius of target set around alfa fixed point in world coordinate AR = PixelWidth*TargetWidth; 
double AR2; // =AR*AR;
double TargetRadiusInPixels; // radius of target set in pixels ; 

// memmory 1D arrays 
 unsigned char *data;
 unsigned char *edge;
 int i; // index of array
double TwoPi=2.0*M_PI;

/* gives position of 2D point (iX,iY) in 1D array  ; uses also global variable iWidth */
unsigned int Give_i(unsigned int ix, unsigned int iy)
{ return ix + iy*iYmax; }

int NumberOfPetals=26;
 
// save data array to pgm file 
int SaveArray2PGMFile( unsigned char data[], double t)
{
 
  FILE * fp;
  const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
  char name [10]; /* name of file */
  sprintf(name,"%f", t); /*  */
  char *filename =strcat(name,".pgm");
  char *comment="# ";/* comment should start with # */
 
  /* save image to the pgm file  */      
  fp= fopen(filename,"wb"); /*create new file,give it a name and open it in binary mode  */
  fprintf(fp,"P5\n %s\n %u %u\n %u\n",comment,iXmax,iYmax,MaxColorComponentValue);  /*write header to the file*/
  fwrite(data,iSize,1,fp);  /*write image data bytes to the file in one step */
  printf("File %s saved. \n", filename);
  fclose(fp);
 
  return 0;
}

double GiveTurn(double complex z)
{
  double argument;
 
  argument = carg(z); //   argument in radians from -pi to pi
  if (argument<0) argument=argument + TwoPi; //   argument in radians from 0 to 2*pi
  return argument/TwoPi ; // argument in turns from 0.0 to 1.0
}


// all points of interior fall into fixed point z=0
// thru 2*13 petals
unsigned char GiveColorOfInterior(double x, double y)

{
int i;
int iMax=NumberOfPetals; 
double angle;
double  offset = 1.0/(2.0*iMax);


 //  works !!!!!!
  angle=GiveTurn(x+y*I); 
  if (angle< offset  || (1.0-offset)<angle) return 231; // 
for(i=0;i<iMax;++i) if (angle< offset+ (1.0*i)/iMax) return ((unsigned char)(230-i));  

    return (230-iMax - 1); // 
      
 
}

int AddBoundaries(unsigned char data[])
{
 
  unsigned int iX,iY; /* indices of 2D virtual array (image) = integer coordinate */
  unsigned int i; /* index of 1D array  */
  /* sobel filter */
  unsigned char G, Gh, Gv; 
 
 
 
 
  printf(" find boundaries in data array using  Sobel filter\n");   

  for(iY=1;iY<iYmax-1;++iY){ 
    for(iX=1;iX<iXmax-1;++iX){ 
      Gv= data[Give_i(iX-1,iY+1)] + 2*data[Give_i(iX,iY+1)] + data[Give_i(iX-1,iY+1)] - data[Give_i(iX-1,iY-1)] - 2*data[Give_i(iX-1,iY)] - data[Give_i(iX+1,iY-1)];
      Gh= data[Give_i(iX+1,iY+1)] + 2*data[Give_i(iX+1,iY)] + data[Give_i(iX-1,iY-1)] - data[Give_i(iX+1,iY-1)] - 2*data[Give_i(iX-1,iY)] - data[Give_i(iX-1,iY-1)];
      G = sqrt(Gh*Gh + Gv*Gv);
      i= Give_i(iX,iY); /* compute index of 1D array from indices of 2D array */
      if (G==0) {edge[i]=245;} /* background */
      else {edge[i]=0;}  /* boundary */
    }
  }
 
  // copy boundaries from edge array to data array 
  //for(iY=1;iY<iyMax-1;++iY){ 
  //   for(iX=1;iX<ixMax-1;++iX){i= Give_i(iX,iY); if (edge[i]==0) data[i]=0;}}
 
 
 
  return 0;
}
 
int CopyBoundaries()
{
 
  unsigned int iX,iY; /* indices of 2D virtual array (image) = integer coordinate */
  unsigned int i; /* index of 1D array  */
      

 printf("copy boundaries from edge array to data array \n");
 for(iY=1;iY<iYmax-1;++iY)
    for(iX=1;iX<iXmax-1;++iX)
       {i= Give_i(iX,iY); if (edge[i]==0) data[i]=0;}
 
 
 
  return 0;
}

 int main()
 {
 
  /* screen ( integer) coordinate */
  int iX,iY;
 
  
  iSize = iXmax*iYmax; //
 /* world ( double) coordinate = parameter plane*/
 const double xMin=-1.2;
 const double xMax=1.2;
 const double yMin=-1.2;
 const double yMax=1.2;
 /* */
 double PixelWidth=(xMax-xMin)/iXmax;
 double PixelHeight=(yMax-yMin)/iYmax;
 //
 /* color component ( R or G or B) is coded from 0 to 255 */
 /* it is 24 bit color RGB file */
 
 unsigned char color;

//

 double x, y, tempx,    /* Z=x+y*i   */
 x0, y0,  /* Z0 = x0 + y0*i */
 x2, y2; /* x2=x*x;  y2=y*y  */
 double x2y2; //, y4, x4, x2my2;
 double x10,y10;
 /*  */
 int Iteration;
 const int IterationMax=2000;

  
 /* bail-out value , radius of circle ;  */
 const int EscapeRadius=2;
 int ER2=EscapeRadius*EscapeRadius;
 TargetRadiusInPixels = 280.0; // = 15.0; // radius of target set in pixels ; Maybe increase to 20 = 1000/50
 AR = PixelWidth*TargetRadiusInPixels; // !!!! important value  ( and precision and time of creation of the pgm image ) 
 AR2= AR*AR;

/* create dynamic 1D arrays for colors ( shades of gray ) */
  data = malloc( iSize * sizeof(unsigned char) );
  edge = malloc( iSize * sizeof(unsigned char) );
  if (edge == NULL || edge == NULL )
    {
      fprintf(stderr," Could not allocate memory\n");
      return 1;
    }
  else fprintf(stderr," memory is OK \n");


 

 /* compute and write image data bytes to the file*/
 for(iY=0;iY<iYmax;++iY)
    {
     y0=yMax - iY*PixelHeight; /* reverse Y  axis */
     if (fabs(y0)<PixelHeight/2) y0=0.0; /*  */  
     printf(" iy = %d from %d\r", iY, iYmax); //info   
     for(iX=0;iX<iXmax;++iX)
       {    /* initial value of orbit Z0 */
        x0=xMin + iX*PixelWidth;
        /* Z = Z0 */
       x=x0;
       y=y0;
            
       x2=x*x;
       y2=y*y;
       x2y2= x2+y2;
      // color=100;
       /* */
       for (Iteration=0;Iteration<IterationMax ;Iteration++)
          { 
             if (x2y2>ER2) { color=255; break;}/* exterior of Filled-in Julia set   */
             // check test every 26 iteration 
             if (!(Iteration % NumberOfPetals) && (x2y2<AR2)) { color=GiveColorOfInterior(x,y); break;}/* interior of Filled-in Julia set   */
             // z^14- z
             y10= y2*y2*y2*y2*y2;
             x10= x2*x2*x2*x2*x2;


//realpart(z1);
// -y^14           +91*x^2*y^12     -1001*x^4*y^10    +3003*x^6*y^8                -3003*x^8*y^6                +1001*x^10*y^4    -91*x^12*y^2    +x^14      -x
tempx = -y10*y2*y2 + 91.0*x2*y10*y2 -1001.0*x2*x2*y10 +3003.0*x2*x2*x2*y2*y2*y2*y2 -3003.0*x2*x2*x2*x2*y2*y2*y2 +1001.0*x10*y2*y2 -91.0*x10*x2*y2 +x10*x2*x2 -x ; 
//  14*x*y^13       -364*x^3*y^11     +2002*x^5*y^9                 -3432*x^7*y^7                 +2002*x^9*y^5                  -364*x^11*y^3     +14*x^13*y       -y 
y = 14.0*x*y10*y2*y -364.0*x2*x*y10*y +2002.0*x*x2*x2*y2*y2*y2*y2*y -3432.0*x2*x2*x2*x*y2*y2*y2*y +2002.0*x2*x2*x2*x2*x*y2*y2*y  -364.0*x10*x*y*y2 +14.0*x10*x2*x*y -y;

	     x=tempx;
             x2=x*x;
             y2=y*y; 
              x2y2= x2+y2;     
             
          };

         i=Give_i(iX,iY);
     data[i]=color;
                         
    
  
  }
 }
SaveArray2PGMFile(data,TargetRadiusInPixels+0.0100); 
AddBoundaries(data); 
SaveArray2PGMFile(edge,TargetRadiusInPixels+0.0200); 
CopyBoundaries();
SaveArray2PGMFile(data,TargetRadiusInPixels+0.0300);  

//
  free(data);
  free(edge);

 return 0;
 }