Felipe Gonzalez Web Page

#***********************************************************************#
#*                      BOUNCING SPHERES                               *#
#*                                                                     *#
#* This program finds the numerical solution to the Newton equation    *#
#*                 r''(t) = -g                                         *#
#* where g is the the gravity acceleration.                            *#
#***********************************************************************#
# AUTHOR: FELIPE GONZALEZ CATALDO, March 2008.


#ifdef __APPLE__
  #include <OpenGL/gl.h>
  #include <OpenGL/glu.h>
  #include <GLUT/glut.h>
#else
  #include <GL/gl.h>
  #include <GL/glu.h>
 #include <GL/glut.h>
#endif


#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "class-ptcl.h"
#include "class-vect.h"
#include "class-azar.h"

using namespace std;

/***************************************************************/
/************************* WINDOW ******************************/
const double ancho=1200.0;
const double alto=1200.0;
const double profundidad=1200.0;
unsigned char Buttons[3] = {0};
double xmax=ancho/2, ymax=alto/2, xmin=-ancho/2, ymin=-alto/2, zmin=-profundidad/2, zmax=profundidad/2;
double zoom = 2.1*ymax/(2*tan(M_PI/8));
double rotx = -70, roty = 0, rotz = 0, tx = 0, ty = 0;
int lastx=0, lasty=0;
/***************************************************************/
/***************************************************************/
const int num=20;
double dc=1;                    // (valor) [dc] = 1 [segundos]
double mts=100;         // (valor) pixeles representan 1 metro real
double g=9.8*mts/(dc*dc);
double t=0;
double dt=1.5e-2*dc;
Ptcl ball[num];
double x[num][2];
double y[num][2];
double z[num][2];
double fx[num][2];
double fy[num][2];
double fz[num][2];
double k1x[num][2];
double k2x[num][2];
double k3x[num][2];
double k4x[num][2];
double k1y[num][2];
double k2y[num][2];
double k3y[num][2];
double k4y[num][2];
double k1z[num][2];
double k2z[num][2];
double k3z[num][2];
double k4z[num][2];


void esfera(double radio, double x, double y, double z){
        double puntos=20;
        glBegin(GL_POLYGON);
        for(int j=0;j<puntos;j++){
                for (int i=0; i<puntos; i++){
                        glVertex3f(x+radio*sin((M_PI*j)/puntos)*cos((2*M_PI*i)/puntos),y+radio*sin((M_PI*j)/puntos)*sin((2*M_PI*i)/puntos),z+radio*cos((M_PI*j)/puntos));
                }
        }
        glEnd();
}

//------------------------------------------GLUT FUNCTIONS-----------------------------------------------------------//
void Inicializa(void){
        glClearColor(0.0, 0.0, 0.0, 0.0);
        glShadeModel (GL_FLAT);
        glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);       // (3D)
        gluOrtho2D(xmin, xmax, ymin, ymax);
        glEnable(GL_DEPTH_TEST);                                                                        // (3D)
}

void Time(void){ t+= dt;        glutPostRedisplay();}

void Mouse(int b,int s,int x,int y){
        lastx=x;
        lasty=y;
        switch(b){
        case GLUT_LEFT_BUTTON:
                Buttons[0] = ((GLUT_DOWN==s)?1:1);
                break;
        case GLUT_MIDDLE_BUTTON:
                Buttons[1] = ((GLUT_DOWN==s)?1:0);
                break;
        case GLUT_RIGHT_BUTTON:
                Buttons[2] = ((GLUT_DOWN==s)?1:0);
                break;
        default:
                break;
        }
        glutPostRedisplay();
}

void Motion(int x,int y){
        int diffx=x-lastx;
        int diffy=y-lasty;
        lastx=x;
        lasty=y;

        if( Buttons[0] && Buttons[1] ){ zoom -= (double) 10.1f * diffx;}
        else{
                if( Buttons[0] ){
                        rotx += (double) 0.5f * diffy;
                        roty += (double) 0.5f * diffx;
                }
                else{
                        if( Buttons[1] ){
                                tx += (double) 0.8f * diffx;
                                ty -= (double) 0.8f * diffy;
                        }
                }
        }
        glutPostRedisplay();
}

void reshape3D(int w, int h){
        // prevent divide by 0 error when minimised
        if(w==0)
                h = 1;

        glViewport(0,0,w,h);
        glMatrixMode(GL_PROJECTION);
        glLoadIdentity();
        gluPerspective(45,(double)w/h,0.1,10000);
        glMatrixMode(GL_MODELVIEW);
        glLoadIdentity();
        gluLookAt (0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);
}

void Up_Down(int key, int x, int y){
        switch(key) {
                case GLUT_KEY_UP:{}
                break;
                case GLUT_KEY_DOWN:{}
                break;
                case GLUT_KEY_RIGHT: dt+=1e-3/5*dc;
                break;
                case GLUT_KEY_LEFT: dt-=1e-3/5*dc;
                break;
                default:break;
        }
        glutPostRedisplay();
}


void Keyboard(unsigned char key, int x, int y){
        switch (key){
                case 27: exit(0);
                break;
                case 'q': exit(0);
                break;
                case 'p':{dt=0;}
                break;
                case 'm':{dt=1.5e-2*dc;}
                break;
                case 'r':{dt=-dt;}
                break;
                case 'a':{zoom=2.1*ymax/(2*tan(M_PI/8)); rotx=0; roty=0;}
                break;
                case 'z':{zoom+=5;}
                break;
                case 'Z':{zoom-=5;}
        }
}

void drawString(char *s){
     unsigned int i;
     for (i = 0; i < strlen(s); i++){   glutBitmapCharacter (GLUT_BITMAP_HELVETICA_12, s[i]);}
}


void drawStringBig(char *s){
     unsigned int i;
     for (i = 0; i < strlen(s); i++){   glutBitmapCharacter (GLUT_BITMAP_HELVETICA_18, s[i]);}
}
//-------------------------------------------------------------------------------------------------------------------//

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////// BOUNDARY CONDITIONS ///////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void    Boundary_Conditions(void){
        randomize();
        // PLANET'S MASSES
        for (int i=0; i<num; i++){ball[i].set_mass(1);}

        // PLANET'S RATIO
        for (int i=0; i<num; i++){      ball[i].set_size(0.2*mts);}

        // INITIAL POSITIONS
        for (int i=0; i<num; i++){      ball[i].set_r(dazar(xmin/2,xmax/2)*e1+dazar(ymin/2,ymax/2)*e2+dazar(0,zmax)*e3);}
        ball[0].set_r(2.4*mts*e3);

        // INICIAL VELOCITIES
        for (int i=0; i<num; i++){      ball[i].set_v(dazar(-5*mts/dc,5*mts/dc)*e1+dazar(-5*mts/dc,5*mts/dc)*e2+dazar(-5*mts/dc,5*mts/dc)*e3);}
        ball[0].set_v(0*e1);

        // COLORS
        for (int i=0; i<num; i++){      ball[i].set_red(dazar(0,1)); ball[i].set_green(dazar(0,1)); ball[i].set_blue(dazar(0,1));}

        // NAMES
        for (int i=0; i<num; i++){      ball[i].set_name("pelota");}

        /**************************Initialize x[][], y[][] & z[][]**************************/
        for (int i=0; i<num; i++){
                x[i][0]=ball[i].r().GetX();             y[i][0]=ball[i].r().GetY();             z[i][0]=ball[i].r().GetZ();
                x[i][1]=ball[i].v().GetX();             y[i][1]=ball[i].v().GetY();             z[i][1]=ball[i].v().GetZ();
        }
        /***********************************************************************************/
}

void Solver(){
/**************************Initialize x[][], y[][] & z[][]**************************/
        for (int i=0; i<num; i++){
                x[i][0]=ball[i].r().GetX();             y[i][0]=ball[i].r().GetY();             z[i][0]=ball[i].r().GetZ();
                x[i][1]=ball[i].v().GetX();             y[i][1]=ball[i].v().GetY();             z[i][1]=ball[i].v().GetZ();
        }
/***********************************************************************************/
//-----------------------------------------------Fixing f------------------------------------------------------------//
        for (int i=0; i<num; i++){

                fx[i][0]=x[i][1];
                fy[i][0]=y[i][1];
                fz[i][0]=z[i][1];

                fx[i][1]=0;     fy[i][1]=0;     fz[i][1]=0;
                fx[i][1]=0;
                fz[i][1]=-g;
        }
//-------------------------------------------------------------------------------------------------------------------//

//-----------------------------------------------Fixing k1-----------------------------------------------------------//
        for (int i=0; i<num; i++){
                k1x[i][0]=dt*fx[i][0];
                k1y[i][0]=dt*fy[i][0];
                k1z[i][0]=dt*fz[i][0];

                k1x[i][1]=dt*fx[i][1];
                k1y[i][1]=dt*fy[i][1];
                k1z[i][1]=dt*fz[i][1];
        }
//-------------------------------------------------------------------------------------------------------------------//

//-----------------------------------------------Fixing k2-----------------------------------------------------------//
        //**************Next is used to evaluate fxij in xij+k1xij/2************//
        for (int i=0; i<num; i++){
                x[i][0]+=k1x[i][0]/2;   y[i][0]+=k1y[i][0]/2;   z[i][0]+=k1z[i][0]/2;
                x[i][1]+=k1x[i][1]/2;   y[i][1]+=k1y[i][1]/2;   z[i][1]+=k1z[i][1]/2;
                ball[i].set_r(x[i][0]*e1+y[i][0]*e2+z[i][0]*e3);
                ball[i].set_v(x[i][1]*e1+y[i][1]*e2+z[i][1]*e3);
        }
        //**********************************************************************//
        //*******************************Fixing*********************************//
        for (int i=0; i<num; i++){
                k2x[i][0]=dt*x[i][1];
                k2y[i][0]=dt*y[i][1];
                k2z[i][0]=dt*z[i][1];

                k2x[i][1]=0; k2y[i][1]=0; k2z[i][1]=0;
                k2x[i][1]+=dt*0;
                k2y[i][1]+=dt*0;
                k2z[i][1]+=dt*(-g);
        }
        //***********************************************************************//
        //***********Next is used to undo the evaluation of fxij above***********//
        for (int i=0; i<num; i++){
                x[i][0]-=k1x[i][0]/2;   y[i][0]-=k1y[i][0]/2;   z[i][0]-=k1z[i][0]/2;
                x[i][1]-=k1x[i][1]/2;   y[i][1]-=k1y[i][1]/2;   z[i][1]-=k1z[i][1]/2;
                ball[i].set_r(x[i][0]*e1+y[i][0]*e2+z[i][0]*e3);
                ball[i].set_v(x[i][1]*e1+y[i][1]*e2+z[i][1]*e3);
        }
        //**********************************************************************//
//-------------------------------------------------------------------------------------------------------------------//

//-----------------------------------------------Fixing k3-----------------------------------------------------------//
        //**************Next is used to evaluate fxij in xij+k2xij/2************//
        for (int i=0; i<num; i++){
                x[i][0]+=k2x[i][0]/2;   y[i][0]+=k2y[i][0]/2;   z[i][0]+=k2z[i][0]/2;
                x[i][1]+=k2x[i][1]/2;   y[i][1]+=k2y[i][1]/2;   z[i][1]+=k2z[i][1]/2;
                ball[i].set_r(x[i][0]*e1+y[i][0]*e2+z[i][0]*e3);
                ball[i].set_v(x[i][1]*e1+y[i][1]*e2+z[i][1]*e3);
        }
        //**********************************************************************//
        //*******************************Fixing*********************************//
        for (int i=0; i<num; i++){
                k3x[i][0]=dt*x[i][1];
                k3y[i][0]=dt*y[i][1];
                k3z[i][0]=dt*z[i][1];

                k3x[i][1]=0; k3y[i][1]=0; k3z[i][1]=0;
                k3x[i][1]+=dt*0;
                k3y[i][1]+=dt*0;
                k3z[i][1]+=dt*(-g);
        }
        //**********************************************************************//
        //**********Next is used to undo the evaluation of fxij above***********//
        for (int i=0; i<num; i++){
                x[i][0]-=k2x[i][0]/2;   y[i][0]-=k2y[i][0]/2;   z[i][0]-=k2z[i][0]/2;
                x[i][1]-=k2x[i][1]/2;   y[i][1]-=k2y[i][1]/2;   z[i][1]-=k2z[i][1]/2;
                ball[i].set_r(x[i][0]*e1+y[i][0]*e2+z[i][0]*e3);
                ball[i].set_v(x[i][1]*e1+y[i][1]*e2+z[i][1]*e3);
        }
        //**********************************************************************//
//-------------------------------------------------------------------------------------------------------------------//

//-----------------------------------------------Fixing k4-----------------------------------------------------------//
        //**************Next is used to evaluate fxij in xij+k3xij**************//
        for (int i=0; i<num; i++){
                x[i][0]+=k3x[i][0];     y[i][0]+=k3y[i][0];     z[i][0]+=k3z[i][0];
                x[i][1]+=k3x[i][1];     y[i][1]+=k3y[i][1];     z[i][1]+=k3z[i][1];
                ball[i].set_r(x[i][0]*e1+y[i][0]*e2+z[i][0]*e3);
                ball[i].set_v(x[i][1]*e1+y[i][1]*e2+z[i][1]*e3);
        }
        //**********************************************************************//
        //*******************************Fixing*********************************//
        for (int i=0; i<num; i++){
                k4x[i][0]=dt*x[i][1];
                k4y[i][0]=dt*y[i][1];
                k4z[i][0]=dt*z[i][1];

                k4x[i][1]=0; k4y[i][1]=0; k4z[i][1]=0;
                k4x[i][1]+=dt*0;
                k4y[i][1]+=dt*0;
                k4z[i][1]+=dt*(-g);
        }
        //***********************************************************************//
        //***********Next is used to undo the evaluation of fxij above***********//
        for (int i=0; i<num; i++){
                x[i][0]-=k3x[i][0];     y[i][0]-=k3y[i][0];     z[i][0]-=k3z[i][0];
                x[i][1]-=k3x[i][1];     y[i][1]-=k3y[i][1];     z[i][1]-=k3z[i][1];
                ball[i].set_r(x[i][0]*e1+y[i][0]*e2+z[i][0]*e3);
                ball[i].set_v(x[i][1]*e1+y[i][1]*e2+z[i][1]*e3);
        }
        //**********************************************************************//
//-------------------------------------------------------------------------------------------------------------------//

//------------------------------------------Computing the Solutions--------------------------------------------------//
        for (int i=0; i<num; i++){
                x[i][0]+=(k1x[i][0]+2.0*k2x[i][0]+2.0*k3x[i][0]+k4x[i][0])/6.0;
                x[i][1]+=(k1x[i][1]+2.0*k2x[i][1]+2.0*k3x[i][1]+k4x[i][1])/6.0;
                y[i][0]+=(k1y[i][0]+2.0*k2y[i][0]+2.0*k3y[i][0]+k4y[i][0])/6.0;
                y[i][1]+=(k1y[i][1]+2.0*k2y[i][1]+2.0*k3y[i][1]+k4y[i][1])/6.0;
                z[i][0]+=(k1z[i][0]+2.0*k2z[i][0]+2.0*k3z[i][0]+k4z[i][0])/6.0;
                z[i][1]+=(k1z[i][1]+2.0*k2z[i][1]+2.0*k3z[i][1]+k4z[i][1])/6.0;
        }
//-------------------------------------------------------------------------------------------------------------------//

//----------------------------------------New Positions and Velocities-----------------------------------------------//
        for (int i=0; i<num; i++){
                ball[i].set_r(x[i][0]*e1+y[i][0]*e2+z[i][0]*e3);
                ball[i].set_v(x[i][1]*e1+y[i][1]*e2+z[i][1]*e3);
        }

//-------------------------------------------------------------------------------------------------------------------//

//----------------------------------------------- Rebotes -----------------------------------------------------------//
        // La bola se mueve en el plano X-Y, sintiendo gravedad en Z:
        for (int i=0; i<num; i++){
                // Choques con el piso:
                if (ball[i].r().GetZ()<ball[i].size()){
                        ball[i].set_v(ball[i].v().GetX()*e1+ball[i].v().GetY()*e2-ball[i].v().GetZ()*e3);
                }
                // Choques con la pared derecha (plano Y-Z)
                if (0.5*xmax-ball[i].r().GetX()<ball[i].size()){
                        ball[i].set_v(-ball[i].v().GetX()*e1+ball[i].v().GetY()*e2+ball[i].v().GetZ()*e3);
                }
                // Choques con la pared izquierda (plano Y-Z)
                if (fabs(0.5*xmin-ball[i].r().GetX())<ball[i].size()){
                        ball[i].set_v(-ball[i].v().GetX()*e1+ball[i].v().GetY()*e2+ball[i].v().GetZ()*e3);
                }
                // Choques con la pared trasera (plano X-Z)
                if (0.5*ymax-ball[i].r().GetY()<ball[i].size()){
                        ball[i].set_v(ball[i].v().GetX()*e1-ball[i].v().GetY()*e2+ball[i].v().GetZ()*e3);
                }
                // Choques con la pared delantera (plano X-Z)
                if (fabs(0.5*ymin-ball[i].r().GetY())<ball[i].size()){
                        ball[i].set_v(ball[i].v().GetX()*e1-ball[i].v().GetY()*e2+ball[i].v().GetZ()*e3);
                }
        }
//-------------------------------------------------------------------------------------------------------------------//

//----------------------------------------------- Rebotes -----------------------------------------------------------//
//      for (int i=0; i<num; i++){for (int j=0; j<num; j++){
//              if (i!=j && crushing(p[i],p[j])){

//              }
//      }}
//-------------------------------------------------------------------------------------------------------------------//
}

void Dibuja(void){
        glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);
        glLoadIdentity();

        glPushMatrix();
        glTranslatef(0,-0.5*ymax,-zoom);
        glTranslatef(tx,ty,0);
        glRotatef(rotx,1,0,0);
        glRotatef(roty,0,1,0);
        glRotatef(rotz,0,0,1);

        Solver();

        //***********************************BALLS**************************************//
        for (int i=0; i<num; i++){ //if (ball[i].is_alive()){
                glColor3d(ball[i].R(),ball[i].G(),ball[i].B());
                glPushMatrix();
                glTranslatef(ball[i].r().GetX(), ball[i].r().GetY(),ball[i].r().GetZ());
                glRotatef(360*ball[i].R(), 0.0, 1.0, 0.0);
                glRotatef(100*t, 0.0, 0.0, 1.0);
                glColor3d(ball[i].R(),ball[i].G(),ball[i].B()); glutWireSphere(ball[i].size(), 20, 16);
                glPopMatrix();
        }//}
        //******************************************************************************//
//              glPushMatrix();
//              glRotatef(0.0, 0.0,0.0, 0.0);
//              glTranslatef(100.0, 0.0,0.0);
//              glColor3d(0,1,1); glutSolidTorus(10, 120, 20, 16);
//              glColor3d(1,0,0); glutSolidSphere(100.0, 20, 16);
//              glPopMatrix();

        glColor3d(0,1,0);
        glBegin(GL_LINES);
        for(double i=xmin/2;i<=xmax/2;i+=50){   glVertex3f(i,ymin/2,0); glVertex3f(i,ymax/2,0);}
        for(double i=ymin/2;i<=ymax/2;i+=50){   glVertex3f(xmin/2,i,0); glVertex3f(xmax/2,i,0);}

        for(double i=xmin/2;i<=xmax/2;i+=50){   glVertex3f(i,ymax/2,0); glVertex3f(i,ymax/2,zmax/2);}
        for(double i=0;i<=zmax/2;i+=50){        glVertex3f(xmin/2,ymax/2,i);    glVertex3f(xmax/2,ymax/2,i);}

        for(double i=ymin/2;i<=ymax/2;i+=50){   glVertex3f(xmin/2,i,0); glVertex3f(xmin/2,i,zmax/2);}
        for(double i=0;i<=zmax/2;i+=50){        glVertex3f(xmin/2,ymin/2,i);    glVertex3f(xmin/2,ymax/2,i);}

        for(double i=ymin/2;i<=ymax/2;i+=50){   glVertex3f(xmax/2,i,0); glVertex3f(xmax/2,i,zmax/2);}
        for(double i=0;i<=zmax/2;i+=50){        glVertex3f(xmax/2,ymin/2,i);    glVertex3f(xmax/2,ymax/2,i);}

        glColor3d(1,0,1);
        for(double i=0;i<=zmax/2;i+=50){        glVertex3f(xmin/2,ymin/2,i);    glVertex3f(xmax/2,ymin/2,i);}
        glEnd();
        glPopMatrix();

        glPushMatrix();
        /**************** Writing in the window ****************/
        glColor3f (0.0, 1.0, 0.0);
        glTranslatef(0,0,-20);
        char *label1="ball.r()";
        char *label2="ball.v()";
        char *label3="Time";
        char cad[30];
        snprintf(cad,15,"%s %4.2f %s %4.2f %s","= (",ball[0].r().GetX(),",",ball[0].r().GetZ(),")");
        glRasterPos3f (-625.0,0, 400.0);        drawString(label1);
        glRasterPos3f (-500.0,0, 400.0);        drawString(cad);

        snprintf(cad,15,"%s %4.2f %s %4.2f %s","= (",ball[0].v().GetX(),",",ball[0].v().GetZ(),")");
        glRasterPos3f (-625.0,0, 360.0);        drawString(label2);
        glRasterPos3f (-500.0,0, 360.0);        drawString(cad);

        snprintf(cad,15,"%s %4.2f %s",": ",t," [dc]");
        glRasterPos3f (-625.0,0, 320.0);        drawString(label3);
        glRasterPos3f (-500.0,0, 320.0);        drawString(cad);
        /********************************************************/
        glPopMatrix();


        glutSwapBuffers();
}
//------------------------------------------------------------------------------------------------//

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

        // Principal Window
        glutInit(&argc, argv);
        glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
        glutInitWindowSize(int(ancho),int(alto));
        glutInitWindowPosition(0,0);
        glutCreateWindow("Free Fall");
        // Program's Functions:
        Inicializa();
        Boundary_Conditions();

        glutIdleFunc(Time);
        glutDisplayFunc(Dibuja);
        glutKeyboardFunc(Keyboard);
        glutSpecialFunc(Up_Down);
        glutMouseFunc(Mouse);
        glutReshapeFunc(reshape3D);
        glutMotionFunc(Motion);

        glutMainLoop();
        return 0;
}
//********************************************************************************//