/* 
Campo Direccional de una ecuacion Diferencial
dy/dx = x^2 - y^2
*/
#include <GL/glut.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

double f(double x, double y)
{
 return (x*x - y*y);
}

void Vector(double xc, double yc)
{
 double m, b, r;
 double a0, a1, a2, tmp;
 double x1, x2;
 double y1, y2;

 r=0.1;
 //y=m*x+b
 m=f(xc,yc);
 b=yc-m*xc;
 //(x-xc)^2 + (y-yc)^2 = r^2
 a0=(b-yc)*(b-yc)+xc*xc-r*r;
 a1=2.0*m*(b-yc)-2.0*xc;
 a2=1.0+m*m;
 tmp=sqrt(a1*a1-4.0*a2*a0);
 x1=(-a1+tmp)/(2.0*a2);
 x2=(-a1-tmp)/(2.0*a2);
 y1=m*x1+b;
 y2=m*x2+b;
 glBegin(GL_LINES);
 glColor3f (1.0, 1.0, 0.5);  
 glVertex2f(x1, y1);
 glColor3f (1.0, 0.0, 0.0);  
 glVertex2f(x2, y2);
 glEnd();
}

void Campo_Direcciones(double xinicial, double xfinal,
					   double yinicial, double yfinal)
{
 double x, y;

 y=yinicial;
 while (y<=yfinal)
	{
	x=xinicial;
	while (x<=xfinal)
		{	
		Vector(x, y);
		x=x+0.2;
		}
	y=y+0.2;
	}
}

void display(void)
{
 glClear(GL_COLOR_BUFFER_BIT);
 glMatrixMode( GL_MODELVIEW_MATRIX );
 glLoadIdentity();

 Campo_Direcciones(-2.0, 2.0, -2.0, 2.0);
 
 glFlush ();
}

void reshape (int w, int h)
{
 if (!h)
	return;
 glViewport(0, 0, w, h);
 glMatrixMode(GL_PROJECTION);
 glLoadIdentity();
 gluOrtho2D(-2.2, 2.2, -2.2, 2.2);
 glMatrixMode(GL_MODELVIEW);
 glLoadIdentity();
}

void keyboard(unsigned char key, int x, int y)
{
 switch (key) 
   {
   case 27: exit(0);
             break;
   }
}    

int main(int argc, char** argv)
{
 glutInit(&argc, argv);
 glutInitDisplayMode (GLUT_SINGLE | GLUT_RGB);
 glutInitWindowSize (400, 400); 
 glutInitWindowPosition (0, 0);
 glutCreateWindow ("Campo de Direcciones");
 glutDisplayFunc(display); 
 glutReshapeFunc(reshape);
 glutKeyboardFunc(keyboard);
 glutMainLoop();
 return 0;
}