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

/* Inverted (1/z) or standard Mandelbrot set */
#define INVERTED 0

/* Color map name (see cmap_list below) */
#define CMAP "cyan"

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

/* Circular cyan ramp */
color_t cyan_ramp[] = {
    { 0,            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 },
    { 1,            1, 0, 1 }
};

/* Custom palette */
color_t isomage1[] = {
    { 0,            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 },
    { 1,            1.000, 1.000, 1.000 }
};

typedef struct {
    char    *name;
    color_t *step;
    int      size;
}
cmap_t;

cmap_t cmap_list[] = {
    { "cyan",     cyan_ramp, sizeof(cyan_ramp) / sizeof(cyan_ramp[0]) },
    { "isomage1", isomage1,  sizeof(isomage1)  / sizeof(isomage1[0])  }
};

cmap_t *get_cmap(char *name) {
    int cmaps, i;

    cmaps = sizeof(cmap_list) / sizeof(cmap_list[0]);

    for (i = 0; i < cmaps; ++i)
        if (strcasecmp(name, cmap_list[i].name) == 0)
            return &(cmap_list[i]);

    return NULL;
}

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

void get_color(double x, cmap_t *cmap, double *r, double *g, double *b) {
    int i;

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

    for (i = 0; i < cmap->size; ++i) {
        int i1 = (i + 1) % cmap->size;
        if (x >= cmap->step[i].t && x < cmap->step[i1].t) {
            *r = interpolate(x, cmap->step[i].t, cmap->step[i1].t, cmap->step[i].r, cmap->step[i1].r);
            *g = interpolate(x, cmap->step[i].t, cmap->step[i1].t, cmap->step[i].g, cmap->step[i1].g);
            *b = interpolate(x, cmap->step[i].t, cmap->step[i1].t, cmap->step[i].b, cmap->step[i1].b);
            break;
        }
    }
}

void write_image(int width, int height, double **data, cmap_t *cmap, 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;
            get_color(data[i][j], cmap, &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, x, y, zx, zy;
    cmap_t *cmap;

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

    if ((cmap = get_cmap(CMAP)) == NULL) {
        fprintf(stderr, "%s: unknown color map\n", CMAP);
        exit(1);
    }

    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) {
        x = xmin + i * dx;
            
        for (j = 0; j < height; ++j) {
            y = ymin + j * dy;

            if (INVERTED) {
                double r2 = x*x + y*y;

                if (r2 == 0.0) {
                    data[i][j] = 1.0;
                    continue;
                }

                cx =  x / r2;
                cy = -y / r2;
            }
            else {
                cx = x;
                cy = y;
            }

            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, cmap, "ms.png");

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