/*
Sistema de Primer Orden

Resolver la siguiente ecuacion diferencial de segundo orden:
  y''=-2y'-4y
con las condiciones iniciales y(0)=2, y'(0)=0 en el intervalo [0,5].

nuevas variables
u1=y  u2=y'
convertimos la ecuacion original en un sistema de primer orden
u1'=u2=f(t,u1,u2)
u2'=-2*u2-4*u1=g(t,u1,u2)
*/
#include <GL/glut.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define	MAXPUNTOS	100

double t[MAXPUNTOS+1], u1[MAXPUNTOS+1], u2[MAXPUNTOS+1];

double f(double t, double u1, double u2)
{
 return (u2);
}

double g(double t, double u1, double u2)
{
 return (-2.0*u2 - 4.0*u1);
}

void RK4(double tinicial, double tfinal, 
		 double u1_inicial, double u2_inicial,
		 int n)
{
 double k1, k2, k3, k4;
 double l1, l2, l3, l4;
 double h;
 int i;

 h=(tfinal-tinicial)/n;
 t[0]=tinicial;
 u1[0]=u1_inicial;
 u2[0]=u2_inicial;
 for (i=0; i<n; i++)
	{
	k1=h*f(t[i], u1[i], u2[i]);
	l1=h*g(t[i], u1[i], u2[i]);
	k2=h*f(t[i]+h/2.0, u1[i]+k1/2.0, u2[i]+l1/2.0);
	l2=h*g(t[i]+h/2.0, u1[i]+k1/2.0, u2[i]+l1/2.0);
	k3=h*f(t[i]+h/2.0, u1[i]+k2/2.0, u2[i]+l2/2.0);
	l3=h*g(t[i]+h/2.0, u1[i]+k2/2.0, u2[i]+l2/2.0);
	k4=h*f(t[i]+h, u1[i]+k3, u2[i]+l3);
	l4=h*g(t[i]+h, u1[i]+k3, u2[i]+l3);
	u1[i+1]=u1[i]+(1.0/6.0)*(k1 + 2.0*k2 + 2.0*k3+ k4);
	u2[i+1]=u2[i]+(1.0/6.0)*(l1 + 2.0*l2 + 2.0*l3+ l4);
	t[i+1]=t[i] + h;
	}
}

void init(void)
{
 RK4(0.0, 5.0, 2.0, 0.0, MAXPUNTOS);
}

void display(void)
{
 int i;

 glClear(GL_COLOR_BUFFER_BIT);
 glMatrixMode( GL_MODELVIEW_MATRIX );
 glLoadIdentity();

 //ejes
 glColor3f(0.0, 0.5, 0.8);  
 glBegin(GL_LINES);
  glVertex2f(0.0, 0.0);
  glVertex2f(5.0, 0.0);
  glVertex2f(0.0, -2.5);
  glVertex2f(0.0, 2.5);
 glEnd();

 //grafica de u1 en color Verde ( y(t) )
 glColor3f(0.0, 1.0, 0.0);  
 glBegin(GL_LINE_STRIP);
 for (i=0; i<MAXPUNTOS; i++)
 	 glVertex2f(t[i], u1[i]);
 glEnd();

 //grafica de u2 en color Rojo ( y'(t) )
 glColor3f(1.0, 0.0, 0.0);
 glBegin(GL_LINE_STRIP);
 for (i=0; i<MAXPUNTOS; i++)
 	 glVertex2f(t[i], u2[i]);
 glEnd();

 
 glFlush ();
}

void reshape (int w, int h)
{
 if (!h)
	return;
 glViewport(0, 0, w, h);
 glMatrixMode(GL_PROJECTION);
 glLoadIdentity();
 gluOrtho2D(-1.0, 6.0, -3.0, 3.0);
 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 (350, 350); 
 glutInitWindowPosition (0, 0);
 glutCreateWindow ("Ecuación de Segundo Orden");
 init ();
 glutDisplayFunc(display); 
 glutReshapeFunc(reshape);
 glutKeyboardFunc(keyboard);
 glutMainLoop();
 return 0;
}