#include <GL/gl.h>            /* gl include */
#include <stdlib.h>           
#include <GL/glut.h>              /* glut.h */

typedef GLfloat point2[2];
point2 pp={-1.0,0.0};
point2 qq={1.0,0.0};
double sqrt();

int times;
float sd;

float gauss(float s,float a)
{
float x=0;
int i;
for(i=0;i<12;i++)
   {
   x+=drand48();
   }
return(s*(x-6.0)+a);
}


void split(point2 p, point2 q, float d, int n)
{
int i;
point2 mid, nor, new;
float s;
if(n==0) 
{
glBegin(GL_LINES);
   glVertex2fv(p);
   glVertex2fv(q);
glEnd();
glFlush();
}
else
{
for(i=0;i<2;i++) mid[i]=(p[i]+q[i])/2.0;
nor[0]=p[1]-q[1];
nor[1]=q[0]-p[0];
s=sqrt(nor[0]*nor[0]+nor[1]*nor[1]);
for(i=0;i<2;i++) nor[i]/=s;
s=gauss(d,0.0);
for(i=0;i<2;i++)
    new[i]= mid[i]+ s*nor[i];
s=sqrt((new[0]-p[0])*(new[0]-p[0])+(new[1]-p[1])*(new[1]-p[1]));
split(p,new,sd*s,n-1);
s=sqrt((new[0]-q[0])*(new[0]-q[0])+(new[1]-q[1])*(new[1]-q[1]));
split(new,q,sd*s,n-1);
}
}

void display()
{
float s;
	glClear(GL_COLOR_BUFFER_BIT);
	glColor3f(0.0, 0.0, 0.0);

glFlush();
s=sqrt((pp[0]-qq[0])*(pp[0]-qq[0])+(pp[1]-qq[1])*(pp[1]-qq[1]));
split(pp ,qq,sd*s,times);
}

void myinit(void)

{

	glClearColor (1.0, 1.0, 1.0, 1.0);


	glMatrixMode (GL_PROJECTION);
	glLoadIdentity ();
	gluOrtho2D(-1.0, 1.0, -1.0, 1.0); 
	glMatrixMode (GL_MODELVIEW);
	glLoadIdentity ();


}

int main(int argc, char** argv)
{

        times = atoi(argv[1]);
        sd= atof(argv[2]);

	glutInit(&argc,argv);
	glutInitDisplayMode (GLUT_SINGLE | GLUT_RGB); 
	glutInitWindowSize(500,500);
	glutInitWindowPosition(0,0);
	glutCreateWindow("Fractal Brownian Motion");
        myinit();
	glutDisplayFunc(display);
	glutMainLoop();

}