/* 
 * c program:
 *   1. draws Filled-in Julia setfor Fc(z)=z+z^5
 *      using :
 *      - boolean escape time  ( exterior / interior)
 *      - attraction time and checking angle ( components of interior , see GiveColorOfInterior)
 * -------------------------------         
 *   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
 * 
 * gcc z1.c -Wall -lm
 * 
 * time ./a.out
 * 
 * ============= const int iXmax = 1000; const int iYmax = 1000;
 * TargetRadiusInPixels = iXmax/15.0;
 * 
 * real      0m9.259s user      0m7.680s sys       0m0.060s
 * 
 */
#include <stdio.h>
#include <stdlib.h>             // malloc
#include <math.h>
#include <string.h>             // strcat

const int       iXmax = 2000;
const int       iYmax = 2000;
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 1 D arrays
unsigned char  *data;
unsigned char  *edge;
int             i;
//index of array

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

//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;
}

//all points of interior fall into fixed point z = 0
// check at what angle = which quadrant
unsigned char
GiveColorOfInterior(double x, double y)
{
        if (x > 0.0 && y > 0.0)
                return 235;
        if (x > 0.0 && y < 0.0)
                return 225;
        if (x < 0.0 && y > 0.0)
                return 215;
        //if (x < 0 && y < 0)
                return 185;
}

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] = 255;
                        }
                         /* background */ 
                        else {
                                edge[i] = 0;
                        }       /* boundary */
                }
        }

        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 */
                //const int     MaxColorComponentValue = 255;
        //FILE * fp;
        //char         *filename = "f2.pgm";
        //char         *comment = "# "; /* comment should start with # */
        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;
        int             Iteration;
        const int       IterationMax = 2000;

        /* bail-out value , radius of circle ;  */
        double          EscapeRadius = 2.0;
        double          ER2 = EscapeRadius * EscapeRadius;
        TargetRadiusInPixels = iXmax / 15.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;
                                //x2my2 = x2 * y2;

                                for (Iteration = 0; Iteration < IterationMax - 1; Iteration++) {
                                        if ((x2y2) > ER2) {
                                                color = 245;
                                                break;
                                        }       /* exterior of Filled-in Julia set   */
                                        if ((x2y2) < AR2) {
                                                color = GiveColorOfInterior(x, y);
                                                break;
                                        }       /* interior of Filled-in Julia set   */
                                        //z + z ^ 5
                                        tempx = 5.0 * x * y * y * y * y - 10.0 * x * x * x * y * y + x * x * x * x * x + x;
                                        //5 * x * y ^ 4 - 10 * x ^ 3 * y ^ 2 + x ^ 5 + x
                                        y = y * y * y * y * y - 10.0 * x * x * y * y * y + 5.0 * x * x * x * x * y + y;
                                        //y ^ 5 - 10 * x ^ 2 * y ^ 3 + 5 * x ^ 4 * y + y
                                        x = tempx;
                                        //
                                        x2 = x * x;
                                        y2 = y * y;
                                        //y4 = y2 * y2;
                                        //x4 = x2 * x2;
                                        x2y2 = x2 + y2;
                                        //x2my2 = x2 * y2;
                                };

                                //if (Iteration == IterationMax) {
                                printf(" Iteration==IterationMax\r");
                                color = GiveColorOfInterior(x, y);
                        }
                i = Give_i(iX, iY);
                data[i] = color;
                }
        }
        SaveArray2PGMFile(data, 0.22);
        AddBoundaries(data);
        SaveArray2PGMFile(edge, 0.32);
        SaveArray2PGMFile(data, 0.42);

        free(data);
        free(edge);

        return 0;
}