/*
	This is a simple application that are able to draw the Mandelbrot Fractal. It is possible to zoom 
	in and out using the mouse. It uses the GNU Multi Precision Library to do the calculations and the
	SDL Library to do the graphics and the mouse-control. Distribute as much as you like.
   
	Using gcc it can be compiled with 'gcc mandel.c -lSDL -lgmp'
	
	Developed by Jonas Lindstrøm Jensen (jonas@imf.au.dk).
	http://home.imf.au.dk/jonas
*/

   
#include <stdio.h>
#include <SDL/SDL.h>
#include <gmp.h>

#define WIDTH 640 // Width of screen
#define HEIGHT 480 // Height of screen
#define BPP 4 // Bytes per pixel
#define DEPTH 32
#define MAX_ITERATIONS 255 // The number of iterations before we do our conclusion
#define ZOOM_RATE 10 // The rate of zoom per click
#define FULLSCREEN 0 // Fullscreen?
#define OUTPUT_DECIMALS 8 // The number of decimals displayed of the position of the mouse cursor

void setpixel(SDL_Surface *screen, int x, int y, Uint8 r, Uint8 g, Uint8 b)
{
	Uint32 *pixmem32;
	Uint32 colour;  
 
	colour = SDL_MapRGB( screen->format, r, g, b );
  
	pixmem32 = (Uint32*) screen->pixels  + y + x;
	*pixmem32 = colour;
}


/* Draws the mandelbrot fractal with upper left corner in (sx,sy) and with each pixel
   having dimensions (dx,dy) */
void DrawScreenGMP(SDL_Surface* screen, mpf_t sx, mpf_t sy, mpf_t dx, mpf_t dy)
{ 
	mpf_t current_x, current_y; // The point currently considered
	mpf_t a, b, u, v, norm; // Computation-variables
	double color_ratio; // = number of colors / MAX_ITERATIONS 
	int x, y, ytimesw, i;
  
	// Lock the surface...
	if(SDL_MUSTLOCK(screen)) {
		if(SDL_LockSurface(screen) < 0) return;
	}

	// There are 255 colours, so what shall we multiply the number of iterations with to get the color
	color_ratio = 255 / MAX_ITERATIONS;
    
	// The point we are currently considering
	mpf_init(current_x); 
	mpf_init(current_y);
  	
	// a+bi is the current value in the iteration
	mpf_init(a); 
	mpf_init(b);
  	
	// Used for temporary storage
	mpf_init(u);
	mpf_init(v);
  	
	// The square norm of a+bi (a^2+b^2)
	mpf_init(norm); 

	for(y = 0; y < screen->h; y++)
	{
		ytimesw = y*screen->pitch/BPP;
		
		// How far are we drawing the image?
		printf("%d%%\r", 100*y / (HEIGHT-1)); 
		fflush(stdout);
		
		for(x = 0; x < screen->w; x++) 
		{
			// current_x = sx + dx*x, current_y = sy + dy*y
			mpf_mul_ui(current_x, dx, x);
			mpf_add(current_x, current_x, sx); 
			mpf_mul_ui(current_y, dy, y);
			mpf_add(current_y, current_y, sy); 
			
			// Set a = b = 0
			mpf_set_ui(a, 0);
			mpf_set_ui(b, 0); 
			
			// norm = a*a + b*b
			mpf_mul(u, a, a);
			mpf_mul(v, b, b);
			mpf_add(norm, u, v); 
			
			// Now iterate using the formula f(a+bi) = (a+bi)^2 + current_x + i current_y
			for(i = 0; (i < MAX_ITERATIONS) && (mpf_cmp_ui(norm, 4) < 0); i++) 
			{
				// norm = a*a + b*b
				mpf_mul(u, a, a);
				mpf_mul(v, b, b);
				mpf_add(norm, u, v); 
				
				// u = a*a - b*b + current_x;
				mpf_sub(u, u, v);
				mpf_add(u, u, current_x); 
				
				// v = 2*a*b + current_y;
				mpf_mul(v, a, b);
				mpf_mul_ui(v, v, 2);
				mpf_add(v, v, current_y); 

				// a = u, b = v
				mpf_set(a,u); 
				mpf_set(b,v);
			}
			setpixel(screen, x, ytimesw, (int) i*color_ratio, (int) i*color_ratio, (int) i*color_ratio);
		}
	}

	// Unlock the surface and flip
	if(SDL_MUSTLOCK(screen)) SDL_UnlockSurface(screen);
	SDL_Flip(screen); 
}


int main(int argc, char* argv[])
{
	// The keyboard map and a variable indicating if we should exit or not
	Uint8 *keyDown;
	int exit = 0;

	// Pointer to the screen and event handler
	SDL_Surface *screen;
	SDL_Event event;

	mpf_t u, v; // Temporary vars
	mpf_init(u);
	mpf_init(v);
   	
	mpf_t dx, dy, sx, sy, w, h;
	mpf_init_set_d(w, 4); // (w,h) = (4,3) - width and height of the displayed area
	mpf_init_set_d(h, 3);
	mpf_init_set_d(sx, -2); // (sx,sy) = (-2.0,1.5) - upper right corner
	mpf_init_set_d(sy, -1.5); 
	mpf_init(dx); // (dx,dy) - width and height of a pixel
	mpf_init(dy);
	mpf_div_ui(dx, w, WIDTH); // dx = w/WIDTH
	mpf_div_ui(dy, h, HEIGHT); // dy = h/HEIGHT

	if(SDL_Init(SDL_INIT_VIDEO) < 0 ) return 1;

	if(!(screen = SDL_SetVideoMode(WIDTH, HEIGHT, DEPTH, SDL_HWSURFACE))) {
		SDL_Quit();
		return 1;
	}
    
	printf("The Mandelbrot Fractal. By Jonas Lindstrøm Jensen (jonas@imf.au.dk).\n");
    
	/* Draw the fractal with upper left corner in (sx,sy) and where each pixel has dimensions (dx,dy) */
	DrawScreenGMP(screen, sx, sy, dx, dy);

	while(!exit)
	{
		while(SDL_PollEvent(&event)) 
		{
			switch(event.type) 
			{
				case SDL_MOUSEMOTION:
				mpf_mul_ui(u, dx, event.motion.x);
				mpf_mul_ui(v, dy, event.motion.y);
				mpf_add(u, u, sx);
				mpf_add(v, v, sy);
				gmp_printf("(%.*Ff,%.*Ff)  \r",OUTPUT_DECIMALS,u,OUTPUT_DECIMALS,v);
				fflush(stdout);
				break;

				case SDL_MOUSEBUTTONDOWN:
				if(event.button.button == SDL_BUTTON_LEFT) { // Left mouse button zooms in
					mpf_div_ui(w, w, ZOOM_RATE); // w = w/ZOOM_RATE
					mpf_div_ui(h, h, ZOOM_RATE); // h = h/ZOOM_RATE
				}

				if(event.button.button == SDL_BUTTON_RIGHT) { // Right mouse button zooms out
					mpf_mul_ui(w, w, ZOOM_RATE); // w = w*ZOOM_RATE
					mpf_mul_ui(h, h, ZOOM_RATE); // h = h*ZOOM_RATE
				}
				
				// sx = event.button.x*dx + sx - w/2
				mpf_mul_ui(u, dx, event.button.x);
				mpf_add(sx, u, sx);
				mpf_div_ui(u, w, 2);
				mpf_sub(sx, sx, u);

				// sy = event.button.y*dy + sy - h/2
				mpf_mul_ui(u, dy, event.button.y);
				mpf_add(sy, u, sy);
				mpf_div_ui(u, h, 2);
				mpf_sub(sy, sy, u);

				// dx = w/WIDTH, dy = h/HEIGHT
				mpf_div_ui(dx, w, WIDTH); 
				mpf_div_ui(dy, h, HEIGHT);

				// Print the coordinate of the upper left corner and the dimensions of the picture
				gmp_printf("(%.Ff,%.Ff) + %.Ff x %.Ff\n", sx,sy,w,h);

				// Draw the stuff!      	
				DrawScreenGMP(screen, sx, sy, dx, dy);
				break;

				// If a escape or the exit button is pressed exit
				case SDL_QUIT:
				exit = 1;
				break;

				case SDL_KEYDOWN:
				keyDown = SDL_GetKeyState(NULL);
				if(keyDown[SDLK_ESCAPE] == 1) exit = 1;
				break;
			}
		}
	}
    
	SDL_Quit();
	return 0;
}