/*
 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>
#include <time.h>

int iX, iY;
#define  iXmax  4000
#define iYmax  4000
#define ThreadsNumber 8

int buffor[iXmax][iYmax][3];
int main()
{
    double s1, s2;
    /* screen ( integer) coordinate */

    /* 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 # */
    static unsigned char color[3];
    /* 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*/



// --------------------------------------- SEQUENCE EXECUTION
    clock_t start, stop;
    double time;
    start = clock();
    for (iY = 0; iY < iYmax; iY++)
    {
        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;
            };
            /* compute  pixel color (24 bit = 3 bytes) */
            if (Iteration == IterationMax)
            { /*  interior of Mandelbrot set = black */
                color[0] = 0;
                color[1] = 0;
                color[2] = 0;
            }
            else
            {                   /* exterior of Mandelbrot set = white */
                color[0] = 255; /* Red*/
                color[1] = 255; /* Green */
                color[2] = 255; /* Blue */
            };
            /*write color to the file*/
            buffor[iX][iY][0] = color[0];
            buffor[iX][iY][1] = color[1];
            buffor[iX][iY][2] = color[2];
        }
    }
    stop = clock();
    time = (double)(((stop - start) / (double)CLOCKS_PER_SEC) * 1000);
    printf("timedif_0,  = %.4f ms\n", time);
// --------------------------------------- SEQUENCE EXECUTION

// ---------------------------------------  PARALLEL EXECUTION
int id, sum[ThreadsNumber];
for(int i =0; i < ThreadsNumber; i++)
{
    sum[i] = 0;
}
    s1 = omp_get_wtime();
#pragma omp parallel private(id) num_threads(ThreadsNumber)
    {
        id = omp_get_thread_num();
#pragma omp for private(color, iY, iX, Iteration, Cy, Cx, Zx, Zy, Zx2, Zy2) schedule(static, 10)
        for (iY = 0; iY < iYmax; iY++)
        {
            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;
                };
                sum[id] += Iteration;
                /* compute  pixel color (24 bit = 3 bytes) */
                if (Iteration == IterationMax)
                { /*  interior of Mandelbrot set = black */
                    color[0] = 0;
                    color[1] = 0;
                    color[2] = 0;
                }
                else
                {                                        /* exterior of Mandelbrot set = white */
                    color[0] = (id+3 * 30) % 255; /* Red*/
                    color[1] = (id+3 * 30) % 255; /* Green */
                    color[2] = (id+3 * 30) % 255; /* Blue */
                };
                /*write color to the file*/
                buffor[iX][iY][0] = color[0];
                buffor[iX][iY][1] = color[1];
                buffor[iX][iY][2] = color[2];
            }
            sum[id] += Iteration;
        }

    }
    s2 = omp_get_wtime();


    for (int i = 0; i < ThreadsNumber; i++)
        printf("Thread %d : %d \n", i, sum[i]);

    printf("timedif_1,  = %.4f \n", (s2 - s1) * 1000);


// ---------------------------------------  PARALLEL EXECUTION

    for (int i = 0; i < iYmax; i++)
    {
        for (int j = 0; j < iYmax; j++)
        {
            fwrite(&buffor[j][i][0], 1, 1, fp);
            fwrite(&buffor[j][i][1], 1, 1, fp);
            fwrite(&buffor[j][i][2], 1, 1, fp);
        }
    }
    fclose(fp);
    return 0;
}