#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <math.h>
#include <gd.h>

double interpolate(double x, double a, double b, double a1, double b1) {
    return a1 + (b1 - a1) * (x - a) / (b - a);
}

typedef struct {
    double t, r, g, b;
}
color_t;

/*
// Circular cyan ramp
color_t color_map[] = {
    { (double) 0/6, 1, 0, 1 },
    { (double) 1/6, 0, 0, 1 },
    { (double) 2/6, 0, 1, 1 },
    { (double) 3/6, 0, 1, 0 },
    { (double) 4/6, 1, 1, 0 },
    { (double) 5/6, 1, 0, 0 },
    { (double) 6/6, 1, 0, 1 }
};
*/

color_t color_map[] = {
    { (double) 0/4, 1.000, 1.000, 1.000 },
    { (double) 1/4, 0.011, 0.019, 0.320 },
    { (double) 2/4, 0.390, 0.074, 0.109 },
    { (double) 3/4, 0.996, 0.792, 0.125 },
    { (double) 4/4, 1.000, 1.000, 1.000 }
};

const int colors = sizeof(color_map) / sizeof(color_map[0]);

void getcolor(double x, double *r, double *g, double *b) {
    int c;

    if (x == 1) {
        *r = *g = *b = 0;
        return;
    }

    for (c = 0; c < colors; ++c) {
        int c1 = (c + 1) % colors;
        if (x >= color_map[c].t && x < color_map[c1].t) {
            *r = interpolate(x, color_map[c].t, color_map[c1].t, color_map[c].r, color_map[c1].r);
            *g = interpolate(x, color_map[c].t, color_map[c1].t, color_map[c].g, color_map[c1].g);
            *b = interpolate(x, color_map[c].t, color_map[c1].t, color_map[c].b, color_map[c1].b);
            break;
        }
    }
}

void write_image(int width, int height, double **data, char *filename) {
    gdImagePtr image;
    FILE *file;
    int i, j;

    image = gdImageCreateTrueColor(width, height);

    for (i = 0; i < width; ++i)
        for (j = 0; j < height; ++j) {
            double r, g, b;
            getcolor(data[i][j], &r, &g, &b);
            gdImageSetPixel(image, i, (height - 1) - j, gdTrueColor((int) (255 * r), (int) (255 * g), (int) (255 * b)));
        }

    if ((file = fopen(filename, "w")) == NULL) {
        perror(filename);
        exit(errno);
    }

    gdImagePng(image, file);

    fclose(file);

    gdImageDestroy(image);
}

int main(int argc, char **argv) {
    int width, height, iterations, i, j, n;
    double **data, cx, cy, xwindow, ywindow, xmin, ymin, dx, dy, zx, zy;

    if (argc != 7) {
        fprintf(stderr, "usage: %s <width> <height> <iterations> <center x> <center y> <window size>\n", argv[0]);
        exit(0);
    }

    width      = atoi(argv[1]);
    height     = atoi(argv[2]);
    iterations = atoi(argv[3]);
    cx         = atof(argv[4]);
    cy         = atof(argv[5]);
    xwindow    = atof(argv[6]);

    ywindow = xwindow * (double) height / width;

    xmin = cx - xwindow/2;
    ymin = cy - ywindow/2;

    dx = xwindow / width;
    dy = ywindow / height;

    data = (double **) malloc(width * sizeof(double *));
    for (i = 0; i < width; ++i)
        data[i] = (double *) malloc(height * sizeof(double));

    for (i = 0; i < width; ++i) {
        cx = xmin + i * dx;
            
        for (j = 0; j < height; ++j) {
            cy = ymin + j * dy;

            for (zx = zy = 0.0, n = 0; n < iterations && zx*zx + zy*zy <= 4.0; ++n) {
                double temp = cx + zx*zx - zy*zy;
                zy = cy + 2 * zx * zy;
                zx = temp;
            }

            data[i][j] = (n == iterations) ? 1.0 : (double) (n % 100) / 100;
        }
    }

    write_image(width, height, data, "ms.png");

    for (i = 0; i < width; ++i)
        free(data[i]);
    free(data);
}