/*
 * Bryan Ward
 * hsv_mandel.c - generate colored mandelbrot set - single processor
 * This currently uses the distance estimator coloring algorithm
 * in conjunction with the escape iteration algorithm 
 */

#include <stdio.h>
#include <complex.h>
#include <math.h>
#include <float.h>
#include <stdlib.h>

#define width 800
#define height 800
#define iter 255
#define xmin -2.0
#define xmax 2.0
#define ymin -2.0
#define ymax 2.0
#define PI 3.141592653589793238462643383279502884197169399375105

typedef struct red_green_blue{
        unsigned char r;
        unsigned char g;
        unsigned char b;
} rgb;

typedef struct hue_saturation_value{
        unsigned char h;
        unsigned char s;
        unsigned char v;
} hsv;

rgb HSVtoRGB(hsv in)
{
        int i;
        float f, p, q, t;
        rgb out;

        if( in.s == 0 ) {
                // achromatic (grey)
                out.r = out.g = out.b = in.v;
                return;
        }

        in.h /= 60;                     // sector 0 to 5
        i = floor( in.h );
        f = in.h - i;                   // factorial part of h
        p = in.v * ( 1 - in.s );
        q = in.v * ( 1 - in.s * f );
        t = in.v * ( 1 - in.s * ( 1 - f ) );

        switch( i ) {
                case 0:
                        out.r = in.v;
                        out.g = t;
                        out.b = p;
                case 1:
                        out.r = q;
                        out.g = in.v;
                        out.b = p;
                case 2:
                        out.r = p;
                        out.g = in.v;
                        out.b = t;
                case 3:
                        out.r = p;
                        out.g = q;
                        out.b = in.v;
                case 4:
                        out.r = t;
                        out.g = p;
                        out.b = in.v;
                default:                // case 5:
                        out.r = in.v;
                        out.g = p;
                        out.b = q;
        }
        return out;

}

unsigned char checkPoint (complex c, double *d)
{
	float fz,fdz;
        double complex z = 0;
	double complex dz = 0;
        unsigned char j = 0;
        while ((creal (z) * creal (z) + cimag (z) * cimag (z)) < 4
               && j < iter)
        {
		dz = 2 * z * dz + 1;
                z = z * z + c;
                j++;
        }
	fz = sqrt(creal(z)*creal(z) + cimag(z)*cimag(z));
	fdz = sqrt(creal(dz)*creal(dz) + cimag(dz)*cimag(dz));
	
	*d = log(fz*fz) * fz / fdz;
//	printf("result: %f\n",log(fz*fz)*fz/fdz);
        return j;
}

void writePixel(rgb in, FILE *fl)
{
	fputc(in.r,fl);
	fputc(in.g,fl);
	fputc(in.b,fl);
}

int main ()
{
        FILE *fl;
        int x, y, i;
        double xpix, ypix;
        // double complex z;
        double complex c;
	unsigned char *escape = malloc(sizeof(unsigned char)*width*height);
	double *distance = malloc(sizeof(double)*width*height);
	hsv tmp;
	if(escape == NULL)
	{
		printf("out of memory\n");
		exit(1);
	}
        if ((fl = fopen ("hsv_mandel.ppm", "w")) == NULL)
                fprintf (stderr, "Cannot open %s\n", "mandelbrot");
        fprintf (fl, "P6\n");
        fprintf (fl, "#Copyright (C) 2007 Bryan Ward\n");
        fprintf (fl, "%d\t%d\n", width, height);
        fprintf (fl, "255\n");

        ypix = (ymax - ymin) / height;
        xpix = (xmax - xmin) / width;
	
	for (i = 0; i < width * height; i++)
	{
		c = xmin + (i%width * xpix) + (ymin + (i/height * ypix))* I;
		escape[i] = checkPoint(c,&distance[i]);
	}
        for (i = 0; i < width * height; i++)
	{
		tmp.h = .125;
		tmp.s = escape[i];
		tmp.v = 0-255*log(distance[i]);//25;
		writePixel(HSVtoRGB(tmp),fl);
		/*fputc(escape[i],fl);
		fputc(distance[i],fl);
		fputc(0,fl);*/
	}
        fclose (fl);
        return 0;
}