/*

  c console program, for CPU, one thread. numbers type : long double 
  It can be compiled and run under Linux, windows, Mac 
  It needs gcc

  draw :
  * check 2 algorithms :
  ** binary escape time
  ** add boundary computed by DEM/M
  * save it to the pgm file 

 
 
  

  -----------------------------------------
  1.pgm file code is  based on the code of Claudio Rocchini
  http://en.wikipedia.org/wiki/Image:Color_complex_plot.jpg
  create 8 bit color graphic file ,  portable gray map file = pgm 
  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
  ---------------------------
  2. first it creates data array which is used to store color values of pixels,
  fills tha array with data and after that writes the data from array to pgm file.
  It alows free ( non sequential) acces to "pixels"
    
  -------------------------------------------
  Adam Majewski   fraktal.republika.pl 
 
  

  to compile : 
  gcc m.c  -lm -Wall -march=native -fopenmp -O3
  to run ( Linux console) :
  time ./a.out

  save  data array to the pgm file 
  File big2er2000.pgm saved. 
  save  graphic data to the text file 
  File big2er2000.txt saved. 
  allways free memory 
 
  real	12m22.161s
  user	94m14.431s
  sys	0m3.460s

   convert a.pgm -resize 2000x2000 a.png
 
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <omp.h>

long double m; // DistanceMax = PixelWidth*m

/* iXmax/iYmax = 1 */
unsigned int iSide = 10000; /* side of rectangle in pixels */
unsigned int iXmax; // ((int)(m*iSide)) /* height of image in pixels */
unsigned int iYmax ; //= iSide;
unsigned int iLength; // = (iXmax*iYmax) /* number of pixels */
/*  world ( long double) coordinate */
long double dSide =  1.5;
long double CxMin ; // = 0.0;
long double CxMax ; // =(m*dSide);
long double CyMin ; //= dSide;
long double CyMax ; // = dSide; 
long double dWidth, dHeight;
long double CenterX, CenterY;

/* (CxMax-CxMin)/(CyMax-CyMin)==iXmax/iYmax = 1 */

unsigned long long int IterationMax; // = (iXmax*100) /* proportional to resolution of picture */

long double PixelWidth; //=  ((CxMax-CxMin)/iXmax)
long double PixelHeight ;//= ((CyMax-CyMin)/iYmax)

long double DistanceMax ; //= PixelWidth; /* proportional to pixel size */

/* fc(z) = z*z + c */

long double EscapeRadius = 2000.0;  /* radius of circle around origin; its complement is a target set for escaping points */
long double ER2 ; //= (EscapeRadius*EscapeRadius)

/* colors = shades of gray from 0=black  to 255=white */
unsigned char iExterior = 255; /* exterior of Julia set  */
unsigned char iBoundary = 0; /* border , boundary*/
unsigned char iInterior = 100;

static inline unsigned int f(unsigned int _iX, unsigned int _iY)
/* 
   gives position of point (iX,iY) in 1D array  ; uses also global variables 
   it does not check if index is good  so memory error is possible 
*/
{return (_iX + (iYmax-_iY-1)*iXmax );}

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

  Hyunsuk Kim  : 
  For Julia sets, z is the variable and c is a constant. Therefore df[n+1](z)/dz = 2*f[n]*f'[n] -- you don't add 1.

  For the Mandelbrot set on the parameter plane, you start at z=0 and c becomes the variable. df[n+1](c)/dc = 2*f[n]*f'[n] + 1. 
  http://iquilezles.org/www/articles/distancefractals/distancefractals.htm

*/

// boolean Escape time and DEM/M in one loop  
static inline unsigned char GiveColor( long double Cx, long double Cy, int iMax, long double DistanceMax)
{ 

  // C = Cx + Cy* I = point of parameter c-plane

  long long int i; /* iteration number */

  long double Zx, Zy; /* Z = Zx + Zy*I  point of dynamical plane */
  long double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
  long double temp;
  long double absZ2;

  // http://en.wikipedia.org/wiki/Complex_quadratic_polynomial
  // first derivative  of fc(zcr) with respect to c =  dZ = dfc(zcr)/dc = 2*Z*dZ  = dZx + dZy*I 
  long double dZx = 0.0;
  long double dZy = 0.0; 
  long double absdZ2; // = abs(dZ)* abs(dZ) = dZx*dZx + dZy*dZy

  long double distance;
  unsigned char color = iInterior;

  

  /* initial value of orbit  = critical point zcr = 0 */
  Zx=0.0; 
  Zy=0.0;
  //
  Zx2 = Zx*Zx;
  Zy2 = Zy*Zy;
  absZ2= Zx*Zx + Zy*Zy;
  absdZ2= dZx*dZx + dZy*dZy;

  // iteration of critical point z= 0 on the dynamical z-plane       
  for (i=0; i<iMax; i++)
    { // check if not escaping : abs(z)>ER
      if (absZ2  > ER2 ) { /*color=  iExterior;*/ break ; } // exterior when escapes
      //if (absdZ2 > 1e60) { return iInterior; } //  interior when derivative explodes 
     

      // in the loop, the derivative should be calculated before the new z
      /* first derivative   zp = 2*z*zp  = xp + yp*i; */
      temp = 2*(Zx*dZx - Zy*dZy) + 1.0 ; /*  */ 
      dZy = 2*(Zx*dZy + Zy*dZx); 
      dZx = temp;

      // z = fc(z) = z*z + c 
      Zy=2*Zx*Zy + Cy;
      Zx=Zx2-Zy2 +Cx;
      
     

      // abs 
      Zx2 = Zx*Zx;      
      Zy2 = Zy*Zy;
      absZ2= Zx2 + Zy2; // nz = x*x + y*y;
      absdZ2= dZx*dZx + dZy*dZy; // 
    };

  if (i<iMax) // exterior 
    { 
     

      // long double distance = sqrt(z2) * log(z2) / (cabs(dc) * pixel_spacing);
      distance = sqrt(absZ2/absdZ2)*log(absZ2); //
      // printf("color =  %d \n", (int)(255.0*distance/DistanceMax));

      // distance = pow( 4.0*distance, 0.25 );
      //if (distance > 0.0 && distance<DistanceMax*) 
      //   { //color = ((int)(255.0*distance)) % 255;

      // d = sqrt( dot(z,z)/dot(dz,dz) )*log(dot(z,z));
      // do some soft coloring based on distance
      //d = pow( 4.0*d/zoom, 0.25 );

      //if (distance < (DistanceMax)) color = iBoundary ; // boundary = black
      // else color=iExterior;

      if (distance< DistanceMax) 
	//color = (int)(255.0*sin(M_PI_2*distance)); // near boundary = shades of gray
	color = (int)(255.0*distance/DistanceMax) ; 
          
      //printf("color =  %d \n", color);
      //   }
      else 
	color=iExterior; // rest of the exterior 
      //}    
    }
  // if (nz < bailout) return 1; //still not escaping after iteration , rare
  // if (absdZ2 < absZ2) color= iExterior;  //includes escaping through 0 // som eiterior points are coded as exterior = error
  else color = iInterior;          

 
  return color; 
}  

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

int main(){

  unsigned int iX,iY, /* indices of 2D virtual array (image) = integer coordinate */
    i; /* index of 1D array  */
  long double Cx, Cy; 
  
  //int  ExtLastIteration; //
  unsigned char color;

  printf(" Setup \n");
  iXmax= iSide; /* height of image in pixels */
  iYmax  = iSide;
  iLength = iXmax*iYmax; /* number of pixels */

  // http://www.pbase.com/duncanc/image/120385219
  CxMin = -1.7492046336;    
  CxMax = -1.7492046333;
  CyMin = -0.00028684675;
  CyMax = -0.00028684645; 
  dWidth =  CxMax-CxMin;
  dHeight = CyMax-CyMin;
  CenterX = CxMax -dWidth/2.0;
  CenterY = CyMax -dHeight/2.0;
  IterationMax = 1000000000;
  PixelWidth=  dWidth/iXmax;
  PixelHeight = dHeight/iYmax;
  ER2 = EscapeRadius*EscapeRadius;
  m= 1.5; //1.0/PixelWidth; // for computing DistanceMax after zoom ; without zoom it should be 1.???
  DistanceMax = PixelWidth/m;

  /* dynamic 1D array for colors ( shades of gray ) */
  unsigned char *data;
  data = malloc( iLength * sizeof(unsigned char) );
  if (data == NULL )
    {
      fprintf(stderr," Could not allocate memory");
      return 1;
    }
    
  

  printf(" compute color \n");
  /* 
     You should make color and i thread-private instead of thread-shared to 
     avoid possible image corruption due to a race condition between your 
     openmp threads (one thread could set color in between another thread 
     setting color and using it).

     Generally it's better (more performance, less overhead) to parallelize 
     the outer loop instead of the inner loop.  Harder to give progress 
     reports then though it's still possible (need atomic ops or critical 
     sections).
     Claude 
     http://mathr.co.uk/blog/

  */
#pragma omp  parallel for schedule(dynamic) shared(data)  private(iY,iX,Cy, Cx, i, color)
  for(iY=0;iY<iYmax;++iY){ 
    Cy=CyMin + iY*PixelHeight; /*  */
    printf("row %u from %u \n",iY, iYmax); 

   
    for(iX=0;iX<iXmax;++iX){ 
      Cx=CxMin + iX*PixelWidth;
      color=GiveColor(Cx, Cy, IterationMax,  DistanceMax);
      
      i= f(iX,iY); /* compute index of 1D array from indices of 2D array */
      data[i]=color;  /* change the pixel color */
	     
    }
  }
      

   

  printf(" save  data array to the pgm file \n");
  unsigned int length = 30;
  char filename[length] ;
  snprintf(filename, length,  "big%.0Lfer%.0Lf%s", m ,EscapeRadius, ".pgm");
  char *comment="#  ";/* comment should start with # */
  const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
  FILE * fp;     
  fp = fopen(filename,"wb"); /*create new file,give it a name and open it in binary mode  */
  fprintf(fp,"P5\n %s\n %u\n %u\n %u\n",comment,iXmax,iYmax,MaxColorComponentValue);  /*write header to the file*/
  fwrite(data,iLength,1,fp);  /*write image data bytes to the file in one step */
  printf("File %s saved. \n", filename);
  fclose(fp);

  printf(" save  graphic data to the text file \n");
  char tfilename[length] ;
  snprintf(tfilename, length, "big%.0Lfer%.0Lf%s", m ,EscapeRadius, ".txt");
  fp = fopen(tfilename,"wb"); /*create new file,give it a name and open it in binary mode  */
  fprintf(fp,"IterationMax = %llu \n", IterationMax);
  fprintf(fp,"EscapeRadius = %Lf ER2 = %Lf \n", EscapeRadius, ER2);
  fprintf(fp,"\n" );  /* */
  fprintf(fp,"C plane : \n" );  /* */
  fprintf(fp,"PixelWidth = %.18Lf ; PixelHeight = %.18Lf \n", PixelHeight, PixelWidth);
  fprintf(fp," dWidth  = %.18Lf \n dHeight = %.18Lf \n", CxMax-CxMin, CyMax-CyMin);  /* */
  fprintf(fp," CxMin = %.18Lf \n", CxMin);  /* */
  fprintf(fp," CxMax = %.18Lf \n", CxMax);  /* */
  fprintf(fp," CyMin = %.18Lf \n", CyMin);  /* */
  fprintf(fp," CyMax = %.18Lf \n", CyMax);  /* */
  fprintf(fp," center C = %.18Lf ; %.18Lf \n",CenterX,CenterY );  /* */
  
  fprintf(fp,"\n" );
  fprintf(fp,"DistanceMax = PixelWidth/m = %.18Lf = %.18Lf * %.18Lf \n", DistanceMax, PixelWidth, m);
  fprintf(fp,"\n" );

  printf("File %s saved. \n", tfilename);
  fclose(fp);

  printf(" allways free memory \n ");
  free(data);
  
  

  return 0;
}