/* Fractal - Mandelbrot set */

#include "frak.h"

pixmap * pixMain=NULL; /* big pixmap for Mandelbrot or Julia set  */
pixmap * pixPreview=NULL; /* small pixmap for Julia preview mode
	(maybe also for orbits mode in future...) */
	
unsigned char needsRecalcMain = 1; /* pixMain needs to be recalculated */
unsigned char needsRecalcPreview = 0; /* pixPreview needs to be recalculated */
unsigned char needsCompleteRecalcMain = 1; /* force recalc, even when pixmap can be moved */
unsigned char enablePixmove = 1; /* enable pixmove optimization */
unsigned char juliaPreviewMode = 0; /* small preview of Julia set */
unsigned char juliaMode = 0; /* Julia set instead of Mandelbrot set */
unsigned char calculatingMain = 0; /* pixMain is being calculated - keys must wait */

unsigned int maxiter = DEFAULT_MAXITER;

/* center of current scene - complex coordinates */
double centerx = DEFAULT_CENTERX, centery  = DEFAULT_CENTERY;
double scale = DEFAULT_SCALE_MANDELBROT;

/* previous scene parameters - used by the pixmove optimization */
double prevcenterx = DEFAULT_CENTERX, prevcentery  = DEFAULT_CENTERY;
double prevscale = 0;

double jx = 0, jy = 0.9; /* Julia parameters */

unsigned char palette[PALETTE_SIZE][3]; /* pallete - dynamicly generate */

/* ------------------------- FRACTAL FUNCTIONS ------------------------- */ 

void generatePalette(void)
{
	int i;
	for (i=0;i<PALETTE_SIZE;i++)
	{
		palette[i][0] = (int) ((sin(i/16.0-7)+1)*128);
		palette[i][1] = (int) ((sin(i/4.0)+1)*128);
		palette[i][2] = (int) ((sin(i/8.0+4)+1)*128);

	}
}

void colorize(const int iter, const int maxiter,
		const double zx, const double zy,
		unsigned char *r, unsigned char *g, unsigned char *b)
{
	int i = iter;

	if (i == maxiter)
	{	/* in coloring */
		*r=*g=*b=0;
	} else
	{	/* out coloring */
		
		/* r = (int) (sin(((double)atan2((double)zx,(double)zy)*10)+1)*127);
		 * g = (int) ((atan2(zx2,zy2))*255/M_PI);
		 */

		/* a "dynamic" palette - works ok */
		
		i = iter % PALETTE_SIZE;
		*r = palette[i][0];
		*g = palette[i][1];
		*b = palette[i][2];

		/*
		*r=*g=*b=(iter % 32)*8;  / * grayscale * /
		*r = (int) ((sin(zx*zy)+1)*127);
		*g = (int) ((cos(zx*zx)/zx+1)*127); / * really nice :) * /
		*/
	}
}

void calcMandelbrotSet(pixmap *p,
			double centerx, double centery,
			double scale,
			int maxiter,
			int fromx, int fromy, /* corners in pixmap */
			int tox, int toy)
{
	/* z_(n+1) = z_n^2 + c
	 * z0 = 0, c = x + yi
	 *
	 * if |z| >= 2, then z diverges
	 * eg. convergence test: sqrt(zx^2 + zy^2) < 2  --> zx^2 + zy^2 < 4
	 */
	int iter, x, y;
	double zx, zy, zx2, zy2;
	double cx, cy;
	double minx, maxx, miny, maxy;
	unsigned char r=0, g=0, b=0;

	if (!p || !p->pixels) return; /* no pixmap, or void pixmap */

	if (fromx<0 || fromy<0 || tox<0 || toy<0 ||
		fromx>p->width ||fromy>p->height ||
		tox>p->width || toy>p->height)
		return;

	/* corners - in complex plane */
	minx = centerx - 1/scale;
	maxx = centerx + 1/scale;
	miny = centery - p->height/(scale*p->width);
	maxy = centery + p->height/(scale*p->width);

	
	cy = miny + fromy*2/(scale*p->width);
	for (y=fromy; y<toy; y++) { /* iterate over pixmap columns... */
		cx = minx + fromx*2/(scale*p->width);
		for (x=fromx; x<tox; x++) { /* ... and rows */
			zx = zy = 0.0;

			/* it seems, that iteration loop using */
			/* while cycle doesn't work ok */
			
			/* main loop */
			iter = 0;
			do {
				/* z^2 */
				/* compute these squares only once */
				zx2 = zx*zx;
				zy2 = zy*zy;
				zy = 2*zx*zy + cy;
				zx = zx2 - zy2 + cx; /* z^2 = (zx + zy)^2 */

				/* z^3
				zx2 = zx*zx;
				zy2 = zy*zy;
				zy = zy*(3*zx2 - zy2) + cy;
				zx = zx*(zx2 - 3*zy2) + cx;
				*/
				
				iter++;
			} while (((zx2 + zy2) < 4) && (iter<maxiter));
			
			colorize(iter, maxiter, zx, zy, &r, &g, &b);	

			putpixel(p, x, y, r, g, b);
			
			/* next row */
			cx += 2/(scale*p->width);/* cx += (maxx-minx)/p->width; */
		}
		/* next column */
		cy += 2/(scale*p->width); /* cy += (maxy-miny)/p->height; */
	}
}

void calcJuliaSet(pixmap *p,
			double centerx, double centery,
			double scale,
			int maxiter, double jx, double jy)
{
	/* z_(n+1) = z_n^2 + d
	 * z0 = x + yi
	 */
	int iter, x, y;
	double zx, zy, zx2, zy2;
	double cx, cy;
	double minx, maxx, miny, maxy;
	unsigned char r=0, g=0, b=0;

	if (!p || !p->pixels) return; /* no pixmap, or void pixmap */

	/* corners */
	minx = centerx - 1/scale;
	maxx = centerx + 1/scale;
	miny = centery - p->height/(scale*p->width);
	maxy = centery + p->height/(scale*p->width);

	cy = miny;
	for (y=0; y<p->height; y++) { /* iterate over pixmap columns... */
		cx = minx;
		for (x=0; x<p->width; x++) { /* ... and rows */
			zx = cx;
			zy = cy;

			iter = 0;
			do {
				/* compute these squares only once */
				zx2 = zx*zx;
				zy2 = zy*zy;
				zy = 2*zx*zy + jy;
				zx = zx2 - zy2 + jx; /* z^2 = (zx + zy)^2 */
				
				iter++;
			} while (((zx2 + zy2) < 4) && (iter<maxiter));
			
			colorize(iter, maxiter, zx, zy, &r, &g, &b);	

			putpixel(p, x, y, r, g, b);

			/* next row */
			cx += 2/(scale*p->width);/* cx += (maxx-minx)/p->width; */
		}
		/* next column */
		cy += 2/(scale*p->width); /* cy += (maxy-miny)/p->height; */
	}
}

void calc(void)
{
	int xoff, yoff;
	if (needsRecalcPreview) {
		calcJuliaSet(pixPreview, 0, 0, DEFAULT_SCALE,
			maxiter, jx, jy);
		needsRecalcPreview = 0;
	}

	if (needsRecalcMain) {
		if (juliaMode) {
			calcJuliaSet(pixMain, centerx, centery, scale,
				maxiter, jx, jy);
		} else if (enablePixmove && !needsCompleteRecalcMain && scale == prevscale) {
			/* -- optimization --
			 * When moving around, copy the common part of previous
			 * and this frame and calculate only newly viewed pixels.
			 * It is not used when scale differs (eg. there is
			 * a change in zoom).
			 */

			/* (centerx-prevcenterx) < 2/scale*/

			if (centery == prevcentery) {
				xoff = (int) round((double)(centerx-prevcenterx)*scale*pixMain->width/2);

				/* horizontal movement */
				pixmove(pixMain, xoff, 0);
				if (centerx < prevcenterx) /* left */
					calcMandelbrotSet(pixMain, centerx, centery, scale,
						maxiter, 0, 0, -xoff, pixMain->height);
				else /* right */
					calcMandelbrotSet(pixMain, centerx, centery, scale,
						maxiter, pixMain->width-xoff, 0,
						pixMain->width, pixMain->height);
			} else if (centerx == prevcenterx) {
				yoff = (int) round((double)(centery-prevcentery)*scale*pixMain->width/2);

				/* vertical movement */
				pixmove(pixMain, 0, yoff);
				if (centery < prevcentery) /* down */
					calcMandelbrotSet(pixMain, centerx, centery, scale,
						maxiter, 0, 0, pixMain->width, -yoff);
				else /* up */
					calcMandelbrotSet(pixMain, centerx, centery, scale,
						maxiter, 0, pixMain->height-yoff,
						pixMain->width, pixMain->height);
			}
		} else { /* no optimization */
			calcMandelbrotSet(pixMain, centerx, centery, scale,
				maxiter, 0, 0, pixMain->width,
				pixMain->height);
		}
		needsRecalcMain = 0;
		calculatingMain = 0;
	}
}