/*
 c program:
 --------------------------------
  1. draws Mandelbrot set for Fc(z)=z*z +c
  using Mandelbrot algorithm ( boolean escape time )
 -------------------------------
 2. technique 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 <math.h>
#include <omp.h>

const int iXmax = 3000;
const int iYmax = 3000;

unsigned char color[iYmax][iXmax][3];
int stat[8] = { 0 };

 int main()
 {
          /* screen ( integer) coordinate */
        int iX,iY;
        // const int iXmax = 800;
        // const int iYmax = 800;
        /* world ( double) coordinate = parameter plane*/
        double Cx,Cy;
        const double CxMin=-2.5;
        const double CxMax=1.5;
        const double CyMin=-2.0;
        const double CyMax=2.0;
        /* */
        double PixelWidth=(CxMax-CxMin)/iXmax;
        double PixelHeight=(CyMax-CyMin)/iYmax;
        /* color component ( R or G or B) is coded from 0 to 255 */
        /* it is 24 bit color RGB file */
        const int MaxColorComponentValue=255;
        FILE * fp;
        char *filename="new1.ppm";
        char *comment="# ";/* comment should start with # */

        /* Z=Zx+Zy*i  ;   Z0 = 0 */
        double Zx, Zy;
        double Zx2, Zy2; /* Zx2=Zx*Zx;  Zy2=Zy*Zy  */
        /*  */
        int Iteration;
        const int IterationMax=200;
        /* bail-out value , radius of circle ;  */
        const double EscapeRadius=2;
        double ER2=EscapeRadius*EscapeRadius;
        /*create new file,give it a name and open it in binary mode  */
        fp= fopen(filename,"wb"); /* b -  binary mode */
        /*write ASCII header to the file*/
        fprintf(fp,"P6\n %s\n %d\n %d\n %d\n",comment,iXmax,iYmax,MaxColorComponentValue);
        /* compute and write image data bytes to the file*/

    double start_time, end_time;

    start_time = omp_get_wtime();
    int id = 0;

    #pragma omp parallel private (id) num_threads(8)
    {
    #pragma omp for private (iY, iX, Zx, Zy, Zx2, Zy2, Cy, Cx, id) schedule(guided, 50)
        for(iY=0;iY<iYmax;iY++)
        {
         id = omp_get_thread_num();
             Cy=CyMin + iY*PixelHeight;
             if (fabs(Cy)< PixelHeight/2) Cy=0.0; /* Main antenna */

             for(iX=0;iX<iXmax;iX++)
             {
                        Cx=CxMin + iX*PixelWidth;
                        /* initial value of orbit = critical point Z= 0 */
                        Zx=0.0;
                        Zy=0.0;
                        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;
                        };
            stat[id] += Iteration;
                        /* compute  pixel color (24 bit = 3 bytes) */
                        if (Iteration==IterationMax)
                        { /*  interior of Mandelbrot set = black */
                           color[iY][iX][0]=0;
                           color[iY][iX][1]=0;
                           color[iY][iX][2]=0;
                        }
                     else
                        { /* exterior of Mandelbrot set = white */
                             color[iY][iX][0]=(id) * 30; /* Red*/
                             color[iY][iX][1]=(id) * 30; /* Green */
                             color[iY][iX][2]=(id) * 30;/* Blue */
                        };
                        /*write color to the file*/
                }
        }
    }
        fwrite(color,1,3*iXmax*iYmax,fp);
    end_time = omp_get_wtime();
        fclose(fp);
    printf("Work took %f seconds\n", end_time - start_time);
     for (int i=0; i<8; i++)
        printf("stat[%d] = %d\n", i, stat[i]);

        return 0;
 }