FREETIME;: Proyecto en OpenGL

Proyecto en OpenGL:

Sistemas Lineales Universidad Nacional Mayor de San Marcos

Facultad de Ciencias Matématicas

Computación Científica

Resumen En el presente trabajo mostramos un programa realizado en OpenGL para vizualizar el Retrato Fase( 2 y 3 dimensiones) de un sistema lineal de EDO's. Mostramos la solución aproximada usando el método de Euler.

1 Conceptos Previos

1.1 OpenGL

OpenGL (ogl en adelante) es la interfaz software de hardware gráfico. Es un motor 3D cuyas rutinas están integradas en tarjetas gráficas 3D. OpenGL posee todas las características necesarias para la representación mediante computadoras de escenas 3D modeladas con polígonos, desde el pintado más básico de triángulos, hasta el mapeado de texturas, iluminación o NURBS.

En OpenGL se proporciona una biblioteca básica de funciones para especificar primitivas gráficas, atribuios, transformaciones geométricas, transformaciones de visualización y muchas otras operaciones. Está diseñada para ser independiente del hardware, por tanto, muchas operaciones, tales como las subrutinas de entrada y salida, no están incluidas en la biblioteca básica. Sin embargo, las subrutinas de entrada y salida y muchas funciones adicionales están disponibles en bibliotecas auxiliares que se han desarrollado para programas OpenGL.

1.1.1 Primitivas Básicas

Una primitiva es simplemente la interpretación de un conjunto de vértices, dibujados de una manera específica en pantalla. Hay diez primitivas distintas en OpenGL, pero en este trabajo expondremos solamente las más importantes: puntos (GL_POINTS), líneas (GL_LINES), . Comentaremos también la primitiva GL_LINES_STRIP utilizadas para definir “tiras” de líneas,

Para crear primitivas en Opengl se utilizan las funciones glBegin y glEnd. La sintaxis de estas funciones sigue el siguiente modelo:

glBegin();
glVertex(...);
glVertex(...);
...
glVertex(...);
glEnd();

Puede observarse que glBegin y glEnd actúan como llaves (“{“ y “}”) de las primitivas, por eso es común añadirle tabulados a las glVertex contenidos entre ellas. El parámetro de glBeing es del tipo glEnum (definido por OpenGL), y será el flag con el nombre de la primitiva .

Dibujo de puntos (GL_POINTS)

glBegin(GL_POINTS);
glVertex3f(0.0f, 0.0f, 0.0f);
glVertex3f(10.0f, 10.0f, 10.0f);
glEnd();

Si en vez de GL_LINES utilizásemos GL_LINE_STRIP, OpenGL ya no trataría los vértices en parejas, si no que el primer vértice y el segundo definirían una línea, y el final de ésta definiría otra línea con el siguiente vértice, y así sucesivamente, definiendo una segmento contínuo. Veamos el siguiente código como ejemplo:

glBegin(GL_LINE_STRIP);
glVertex3f(0.0f, 0.0f, 0.0f);
glVertex3f(2.0f, 1.0f, 0.0f);
glVertex3f(2.0f, 2.0f, 0.0f);
glEnd();

Dibujo de conos(Cone)

glutSolidCone(base, height, slices, stacks)

Rotación

Se define mediante la primitiva glRotatef(alpha, x, y, z). Esta función multiplica la matriz actual por una matriz de rotación de alpha grados respecto al eje (x, y, z).

Traslación

Se define mediante la primitiva glTranslatef(x, y, z). Aplica una translación en (x, y, z) sobre la matriz actual.

Escalado

Se define mediante la primitiva glScalef(sx, sy, sz) .Escalado de cada uno de los ejes.

1.2 EDO's Autónomas Lineales

1.2.1 Sistemas de EDO's de primer orden
$x'=f(x) ,\; donde\; x\in\mathbb{R}^{n} \; y \; f:\mathbb{R}^{n}\; \rightarrow\mathbb{R}^{n}\; ,\; f=(f_{1},\ldots,f_{n})$

$\; donde \;f_{i}:\mathbb{R}^{n}\rightarrow\mathbb{R} \;\; para\; i=1,\ldots,n.$

1.2.2 Clasificación:

$1. Sist.\: de\; EDO's \; Lineales:\; Si \: f(x)=Ax, donde \; A\in M_{nxn}(\mathbb{R}).$
$2. Sist. \: de\: EDO's \:no\: Lineales:$

$Si\; es\; un\; Sist.\; de\; EDO's\; que \; no \; satisfaga:\; f(x)=Ax.$

1.2.3 Sistema Bidimensional Lineal
$X'=AX$
$donde \; $X(t)=(x(t),y(t))^{t}$ y \; $A\; \in M_{2x2}(\mathbb{R})$.\; Es decir: \begin{eqnarray*}\\ x' & = & f_{1}(x,y)=a_{11}x+a_{12}y\\ y' & = & f_{2}(x,y)=a_{21}x+a_{22}y\end{eqnarray*}$

1.2.3 Sistema Tridimensional Lineal

$X'=AX$
$donde \; $X(t)=(x(t),y(t),z(t))^{t}$ \; y \; $A\in M_{3x3}(\mathbb{R})$.\; Es decir:\\ x' & =f_{1}(x,y,z)= & a_{11}x+a_{12}y+a_{13}z\\ y' & =f_{2}(x,y,z)= & a_{21}x+a_{22}y+a_{23}z\\ z' & =f_{3}(x,y,z)= & a_{31}x+a_{32}y+a_{33}z$

El siguiente paso es graficar las velocidades para vizualizar el comportamiento de las soluciones.

2.1 Contrución del vector velocidad(flecha)

Debemos construir un objeto que represente al vector velocidad. Luego lo
situaremos en su lugar correspondiente en cada punto del espacio por medio de
transformaciones.
void echa(){
GL oat x1=0,y1=0,z1=0,x2=1,y2=0,z2=0;
glLineWidth(1.0f);
glBegin(GL_LINES);
glVertex3f(x1,y1,z1);
glVertex3f(x2,y2,z2);
glEnd();
glTranslatef(0.8,0,0);
glRotatef(90, 0,1, 0);
glutSolidCone(0.2, 0.4, 10, 10);
}
Primero construimos el cuerpo que no es mas que una linea.
Luego usamos glutSolidCone() para la cabeza. Debemos hacer una traslación
y una rotación para situarla en su lugar.
Tambien tenemos que construir los ejes coordenados.
Ahora hay que hacer las traslaciones y rotaciones necesarias para ubicar
nuestro vector en la direción correspondiente.
Como el vector esta inicialmente en el eje x, debemos rotarlo un ángulo
con respecto a un eje que sea perpedicular al plano que forma el vector velocidad
con el eje x. Es decir el eje z.
De esta manera ubicamos nuestro vector en la dirección correspondiente en
cada punto discretizado del espacio.
Tambien mostramos la una curva solución dada una condición inicial. Usamos
el método de Euler para aproximar la curva. Es decir:
$X(t+\Delta t)=X(t)+\Delta tX'(t)$

3 Ideas Básicas y Construcción(3D)

Generalizando a tres dimensiones, ahora tenemos que comenzar discretizando el espacio XYZ. La idea es graficar en cada punto del espacio el vector velocidad.

3.1 Contrución del vector velocidad(flecha)

Debemos construir un objeto que represente al vector velocidad. Luego lo situaremos en su lugar correspondiente en cada punto del espacio por medio de transformaciones.

void echa(){
GL oat x1=0,y1=0,z1=0,x2=1,y2=0,z2=0;
//cuerpo
glBegin(GL_LINES);
glVertex3f(x1,y1,z1);
glVertex3f(x2,y2,z2);
glEnd();
//cabeza
glTranslatef(0.8,0,0);
glRotatef(90, 0,1, 0);
glutSolidCone(0.2, 0.4, 10, 10);
}

Primero construimos el cuerpo que no es mas que una linea.

Luego usamos glutSolidCone() para la cabeza. Debemos De esta manera ubicamos nuestro vector en la dirección correspondiente en
cada punto discretizado del espacio.
Tambien mostramos la una curva solución dada una condición inicial. Usamos
el método de Euler para aproximar la curva. Es decir:hacer una traslación y una rotación para situarla en su lugar.

Tambien tenemos que construir los ejes coordenados.

void ejes(){
// X
glColor3f(0,0,0);
glBegin(GL_LINES);
glVertex3f(0,0,0);
glVertex3f(5,0,0);
glVertex3f(5,0.2,0);
glVertex3f(5.2,0.4,0);
glVertex3f(5,0.4,0);
glVertex3f(5.2,0.2,0);
glEnd();
// Y
glBegin(GL_LINES);
glVertex3f(0,0,0);
glVertex3f(0,5,0);
glVertex3f(-0.1,5,0);
glVertex3f(0.1,5.2,0);
glVertex3f(-0.1,5.2,0);
glVertex3f(0,5.1,0);
glEnd();
// Z
glBegin(GL_LINES);
glVertex3f(0,0,0);
glVertex3f(0,0,5);
glVertex3f(0,-0.1,5);
glVertex3f(0,-0.1,5.2);
glVertex3f(0,0.1,5);
glVertex3f(0,0.1,5.2);
glVertex3f(0,-0.1,5.2);
glVertex3f(0,0.1,5);
glEnd();
}

Ahora hay que hacer las traslaciones y rotaciones necesarias para ubicar
nuestro vector en la direción correspondiente.
Como el vector esta inicialmente en el eje x, debemos rotarlo un ángulo alfa
con respecto a un eje que sea perpedicular al plano que forma el vector velocidad
con el eje x.

$Se\; sabe\: que:\\ $cos(\alpha)=\dfrac{\hat{x}\vec{V}}{\Vert\hat{x}\Vert\Vert\vec{V}\Vert}=\dfrac{(1,0,0).(x',y',z')}{1\Vert\vec{V}\Vert}=\dfrac{x'}{\Vert\vec{V}\Vert}$\\ $\alpha=arccos(\dfrac{x'}{\Vert\vec{V}\Vert})$.$

Un vector perpendicular a x y a V es :

$\left[\begin{array}{ccc} i & j & k\\ 1 & 0 & 0\\ x' & y' & z'\end{array}\right]=(0,-z',y')$

Luego hay que rotar grados respecto a (0;-z'; -y'). Antes debemos trasladar
nuestro vector al punto donde queremos tener el vector velocidad.

v1=A*x+B*y+C*z;
v2=D*x+E*y+F*z;
v3=G*x+H*y+I*z;
mod=sqrt(pow(A*x+B*y+C*z,2)+: : :
pow(D*x+E*y+F*z,2)+pow(G*x+H*y+I*z,2));
alfa=acos(v1/mod)*180/3.14159;
glTranslatef(x,y,z);
glRotatef(alfa,0,-v3,v2);
glScalef(.5, .5, 0);
flecha();

De esta manera ubicamos nuestro vector en la dirección correspondiente en
cada punto discretizado del espacio.
Tambien mostramos la una curva solución dada una condición inicial. Usamos
el método de Euler para aproximar la curva. Es decir:

$X(t+\Delta t)=X(t)+\Delta tX'(t)$

xf=xi+0.05*(A*xi+B*yi+C*zi);
yf=yi+0.05*(D*xi+E*yi+F*zi);
zf=zi+0.05*(G*xi+H*yi+I*zi);

4 Ejecución

El código completo es para el plano fase es:

#include
#include
#include

float A , B, C, D ; //matriz
float alpha , beta ;
float x,y;//condicion inicial
int   xA , yA;

using namespace std;

void init()
{
glClearColor(1,1,1,0);//color de fondo
}

void ejes(){

//       X
glColor3f(0,0,0);

glBegin(GL_LINES);
        glVertex3f(0,0,0);
        glVertex3f(5,0,0);
        glVertex3f(5,0.2,0);
        glVertex3f(5.2,0.4,0);
        glVertex3f(5,0.4,0);
        glVertex3f(5.2,0.2,0);
glEnd();

//      Y
glBegin(GL_LINES);
        glVertex3f(0,0,0);
        glVertex3f(0,5,0);
        glVertex3f(-0.1,5,0);
        glVertex3f(0.1,5.2,0);
        glVertex3f(-0.1,5.2,0);
        glVertex3f(0,5.1,0);
glEnd();

}

void sol(){
int k;
float xf,yf,xi,yi;

glColor3f(1, 0, 0);

glPointSize(8.0);

glBegin(GL_POINTS);
    glVertex3f(x, y, 0);
glEnd();

glColor3f(0.10f, 0.07f, 0.33f);
glLineWidth(2.0f);
xi=x;
yi=y;

glBegin(GL_LINE_STRIP);
    for(k=0;k<200;k++) {

        xf=xi+0.05*(A*xi+B*yi);
        yf=yi+0.05*(C*xi+D*yi);

        if ((yf>4)||(yf<-4)||(xf>4)||(xf<-4)){
            break;
            }
        glVertex3f(xi,yi,0);
        glVertex3f(xf,yf,0);
        xi=xf;
        yi=yf;

    }
glEnd();

xi=x;
yi=y;

glBegin(GL_LINE_STRIP);
    for(k=0;k<200;k++) {
        xf=xi-0.05*(A*xi+B*yi);
        yf=yi-0.05*(C*xi+D*yi);


        if ((yf>4)||(yf<-4)||(xf>4)||(xf<-4)){
            break;
            }
        glVertex3f(xi,yi,0);
        glVertex3f(xf,yf,0);
        xi=xf;
        yi=yf;

    }
glEnd();

}

void flecha(){

    GLfloat x1=0,y1=0,z1=0,x2=1,y2=0,z2=0;
    glLineWidth(1.0f);
    glBegin(GL_LINES);
        glVertex3f(x1,y1,z1);
        glVertex3f(x2,y2,z2);

    glEnd();
    glTranslatef(0.8,0,0);
    glRotatef(90, 0,1, 0);
    glutSolidCone(0.2, 0.4, 10, 10);

    }

void flechas(){

int i, j,k;
float h=1,mod,x,y,z,v1,v2,v3,alfa;
for (i=-3;i<4;i++){
        x=0+i*h;
        glColor3f((x+3)/6,(3-x)/6,1);
        for (j=-3;j<4;j++){
            y=0+j*h;


                    glPushMatrix();


                    v1=A*x+B*y;
                    v2=C*x+D*y;


                    mod=sqrt(pow(v1,2)+pow(v2,2));
                    alfa=acos(v1/mod)*180/3.14159;
                    glTranslatef(x,y,z);
                    glRotatef(alfa,0,0,v2);
                    glScalef(.5, .5, 0);

                    flecha();
                    glPopMatrix();


            }
        }

}


void ingresa_matriz(){
    cout<<"Generador de plano Fase(3D)\n";
    cout<<"Ingrese la matriz 2x2 :\n";
    cout<<"M[1][1] -->";
    cin>>A;
    cout<<"M[1][2] -->";
    cin>>B;
    cout<<"M[2][1] -->";
    cin>>C;
    cout<<"M[2][2] -->";
    cin>>D;


    cout<<"\nIngrese la condicion inicial :\n";
    cout<<"x0 =";
    cin>>x;
    cout<<"y0 =";
    cin>>y;


    }

void display(void) {


init();
glClear(GL_COLOR_BUFFER_BIT);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
glOrtho(-7, 7, -7, 7, -7, 7);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
gluLookAt(0.0f, 0.0f, 5.0f,
0.0f, 0.0f, 0.0f,
0.0f, 1.0f, 0.0f);
glRotatef(alpha, 1.0f, 0.0f, 0.0f);
glRotatef(beta, 0.0f, 1.0f, 0.0f);

ejes();

sol();

flechas();

glutSwapBuffers();
}
void onMouse(int button, int state, int x, int y) {
if ( (button == GLUT_LEFT_BUTTON) & (state == GLUT_DOWN) ) {
xA = x; yA = y;
}
}
void onMotion(int x, int y) {
alpha = (alpha + (y - yA));
beta = (beta + (x - xA));
xA = x; yA = y;
glutPostRedisplay();
}

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

ingresa_matriz();
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
glutInitWindowSize(500, 500);
glutInitWindowPosition(200, 200);
glutCreateWindow("ESPACIO FASE");

glutDisplayFunc(display);
glutMouseFunc(onMouse);
glutMotionFunc(onMotion);
glutMainLoop();
}

Ejecutemos el programa desde el terminal

El código completo es para el espacio fase es:

#include
#include
#include

float A , B, C, D , E, F, G, H , I ; //matriz
float alpha , beta ;
float x,y,z;//condicion inicial
int   xA , yA;

using namespace std;

void init()
{
glClearColor(1,1,1,0);//color de fondo
}

void ejes(){

//       X
glColor3f(0,0,0);

glBegin(GL_LINES);
        glVertex3f(0,0,0);
        glVertex3f(5,0,0);
        glVertex3f(5,0.2,0);
        glVertex3f(5.2,0.4,0);
        glVertex3f(5,0.4,0);
        glVertex3f(5.2,0.2,0);
glEnd();

//      Y
glBegin(GL_LINES);
        glVertex3f(0,0,0);
        glVertex3f(0,5,0);
        glVertex3f(-0.1,5,0);
        glVertex3f(0.1,5.2,0);
        glVertex3f(-0.1,5.2,0);
        glVertex3f(0,5.1,0);
glEnd();

//     Z
glBegin(GL_LINES);
        glVertex3f(0,0,0);
        glVertex3f(0,0,5);
        glVertex3f(0,-0.1,5);
        glVertex3f(0,-0.1,5.2);
        glVertex3f(0,0.1,5);
        glVertex3f(0,0.1,5.2);
        glVertex3f(0,-0.1,5.2);
        glVertex3f(0,0.1,5);
glEnd();

}

void sol(){
int k;
float xf,yf,zf,xi,yi,zi;

glColor3f(1, 0, 0);

glPointSize(8.0);

glBegin(GL_POINTS);
    glVertex3f(x, y, z);
glEnd();

glColor3f(0.10f, 0.07f, 0.33f);
glLineWidth(2.0f);
xi=x;
yi=y;
zi=z;
glBegin(GL_LINE_STRIP);
    for(k=0;k<200;k++) {

        xf=xi+0.05*(A*xi+B*yi+C*zi);
        yf=yi+0.05*(D*xi+E*yi+F*zi);
        zf=zi+0.05*(G*xi+H*yi+I*zi);
        if ((zf>4)||(zf<-4)||(yf>4)||(yf<-4)||(xf>4)||(xf<-4)){
            break;
            }
        glVertex3f(xi,yi,zi);
        glVertex3f(xf,yf,zf);
        xi=xf;
        yi=yf;
        zi=zf;
    }
glEnd();

xi=x;
yi=y;
zi=z;
glBegin(GL_LINE_STRIP);
    for(k=0;k<200;k++) {
        xf=xi-0.05*(A*xi+B*yi+C*zi);
        yf=yi-0.05*(D*xi+E*yi+F*zi);
        zf=zi-0.05*(G*xi+H*yi+I*zi);

        if ((zf>4)||(zf<-4)||(yf>4)||(yf<-4)||(xf>4)||(xf<-4)){
            break;
            }
        glVertex3f(xi,yi,zi);
        glVertex3f(xf,yf,zf);
        xi=xf;
        yi=yf;
        zi=zf;
    }
glEnd();

}

void flecha(){

    GLfloat x1=0,y1=0,z1=0,x2=1,y2=0,z2=0;
    glLineWidth(1.0f);
    glBegin(GL_LINES);
        glVertex3f(x1,y1,z1);
        glVertex3f(x2,y2,z2);

    glEnd();
    glTranslatef(0.8,0,0);
    glRotatef(90, 0,1, 0);
    glutSolidCone(0.2, 0.4, 10, 10);

    }

void flechas(){

int i, j,k;
float h=1,mod,x,y,z,v1,v2,v3,alfa;
for (i=-3;i<4;i++){
        x=0+i*h;
        glColor3f((x+3)/6,(3-x)/6,1);
        for (j=-3;j<4;j++){
            y=0+j*h;
            for (k=-3;k<4;k++){
                z=0+k*h;
                    glPushMatrix();


                    v1=A*x+B*y+C*z;
                    v2=D*x+E*y+F*z;
                    v3=G*x+H*y+I*z;

                    mod=sqrt(pow(v1,2)+pow(v2,2)+pow(v3,2));
                    alfa=acos(v1/mod)*180/3.14159;
                    glTranslatef(x,y,z);
                    glRotatef(alfa,0,-v3,v2);
                    glScalef(.5, .5, 0);

                    flecha();
                    glPopMatrix();

                }
            }
        }

}


void ingresa_matriz(){
    cout<<"Generador de Espacio Fase(3D)\n";
    cout<<"Ingrese la matriz 3x3 :\n";
    cout<<"M[1][1] -->";
    cin>>A;
    cout<<"M[1][2] -->";
    cin>>B;
    cout<<"M[1][3] -->";
    cin>>C;
    cout<<"M[2][1] -->";
    cin>>D;
    cout<<"M[2][2] -->";
    cin>>E;
    cout<<"M[2][3] -->";
    cin>>F;
    cout<<"M[3][1] -->";
    cin>>G;
    cout<<"M[3][2] -->";
    cin>>H;
    cout<<"M[3][3] -->";
    cin>>I;
    cout<<"\nIngrese la condicion inicial :\n";
    cout<<"x0 =";
    cin>>x;
    cout<<"y0 =";
    cin>>y;
    cout<<"z0 =";
    cin>>z;

    }

void display(void) {


init();
glClear(GL_COLOR_BUFFER_BIT);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
glOrtho(-7, 7, -7, 7, -7, 7);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
gluLookAt(0.0f, 0.0f, 5.0f,
0.0f, 0.0f, 0.0f,
0.0f, 1.0f, 0.0f);
glRotatef(alpha, 1.0f, 0.0f, 0.0f);
glRotatef(beta, 0.0f, 1.0f, 0.0f);

ejes();

sol();

flechas();

glutSwapBuffers();
}
void onMouse(int button, int state, int x, int y) {
if ( (button == GLUT_LEFT_BUTTON) & (state == GLUT_DOWN) ) {
xA = x; yA = y;
}
}
void onMotion(int x, int y) {
alpha = (alpha + (y - yA));
beta = (beta + (x - xA));
xA = x; yA = y;
glutPostRedisplay();
}

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

ingresa_matriz();
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
glutInitWindowSize(500, 500);
glutInitWindowPosition(200, 200);
glutCreateWindow("ESPACIO FASE");

glutDisplayFunc(display);
glutMouseFunc(onMouse);
glutMotionFunc(onMotion);
glutMainLoop();
}

Ejecutemos el programa desde el terminal

Referencias
[1] Jaime E. Villate.Introducao aos sistemas dinamicos. Universidad de Porto.
[2] Donald Hearn.Grá cos por Computadora con Opengl. Indiana University.

FREETIME;

domingo, 5 de diciembre de 2010

Proyecto en OpenGL

No hay comentarios:

Publicar un comentario

Visitas