#include <stdio.h>
#include <math.h>
#include <GLUT/glut.h>
#include <OpenGL/gl.h>
#include <OpenGL/glu.h>
#include <OpenGL/OpenGL.h>

#define X 500
#define Y 500


float image[X][Y];		    //store image color data here
typedef GLfloat point2d[2];
point2d point = {0.0, 0.0};

void display(void)
{
    //clear background to glClearColor
    glClear ( GL_COLOR_BUFFER_BIT );
    
    //point color
    float red=0.0, green=0.0, blue=0.0;
    glColor3f ( red, green, blue );
    
    int x, y;
    
    for ( x=0; x<X; x++){
	for ( y=0; y<Y; y++){
	    
	    point[0] = x;
	    point[1] = y;
	    
	    //use image to set color
	    if ( image[x][y] == 1 ){
		red = 0.0;
		blue = 0.0;
		green =0.0;
	    }else{
		red = image[x][y]; 
		blue = image[x][y]; 
		green = image[x][y];
	    }
	    
	    //set color
	    glColor3f( red, green, blue );
	    
	    //draw points
	    glBegin (GL_POINTS);
		glColor3f ( red, green, blue );
		glVertex2fv(point);
	    glEnd();
	}
    }
   
    
    
    
    //start processing
    glFlush();

}

void mandelbrot()
{
    int i, j, k;
    float min  = -1.0;
    float max = 1.0;
    int max_iterations = 200;
    float complex0, real0, complex1, real1, complex2, real2, magnitude;
    
    //scale min/max onto image size (X, Y)
    float distance = max - min;
    float length = X;
    
    float scale = distance/length;
    
  
    for ( i=0; i<X; i++){
	
	for ( j=0; j<Y; j++){
	    
	    k=0;
	    real0 = 0.0; 
	    complex0 = 0.0; 

	    magnitude = 0.0;
	   
	    //scale our position onto the image
	    real2 = min + scale * i;
	    complex2 = min + scale * j;
	   
	    
	    while ( k<max_iterations){
		
		//loop count
		k++;
		
		//multiply( real0 + complex0 )  * ( real0 + complex0 ) and return ( real1 + complex1 )
		//( a + bi ) * ( a + bi ) = a^2 + 2*a*b - b^2
		real1 = real0 * real0 - complex0 * complex0;
		complex1 = 2 * real0 * complex0;
		
		//add our result to our position
		real0 = real1 + real2;
		complex0 = complex1 + complex2;
		
		//check the magnitude of this number 
		//magnitude = square root of ( real^2 + complex^2 );
		magnitude = sqrt ( real1 * real1 + complex1 * complex1 );
		
		//see if part of mandelbrot set
		if ( magnitude > 2.0 ) break;
		
	    }
	    
	    
	    //set up image data
	    
	    
	    
	    if ( magnitude > 1.0 ) magnitude = 1.0;
	    image[i][j] = magnitude;
	    
	    
	}// end i loop
    }// end j loop
    
    
    
}


void init (void) 
{
    
    //select clearing color
    glClearColor (0.0, 0.0, 0.0, 0.0);
   
    //remove hidden stuff for faster rendering
    glEnable(GL_DEPTH_TEST);

    //initialize viewing values  
    glMatrixMode(GL_PROJECTION);
    gluOrtho2D( 0.0, X, 0.0, Y);
    glMatrixMode( GL_MODELVIEW );
    
    

}




int main( int argc, char** argv )
{
   
    mandelbrot(); 
    
   glutInit(&argc, argv);
   glutInitDisplayMode (GLUT_SINGLE | GLUT_RGB);
   glutInitWindowSize ( X, Y); 
   glutInitWindowPosition (100, 100);
   glutCreateWindow ("test");
   init ();
   glutDisplayFunc(display); 
   glutMainLoop();
   
   return 0;  
}