main.cpp

Go to the documentation of this file.

#include <GL/glut.h>
#include <iostream>



#include <MC_particle_engine.hpp>
#include <MC_grid_3d_scalar.hpp>
#include <MC_grid_3d_scalar_marching_cube.hpp>
#include <MC_potential_sphere.hpp>

#include <mesh_converter.hpp>
#include <vector>
#include <algorithm>


using namespace mesh_conv;
using namespace std;

//glut functions
static void display_callback();
static void reshape_callback (int width, int height);
static void keyboard_callback(unsigned char key,int x,int y);
static void idle_callback(int value);
static void motion_callback (int x, int y);
static void mouse_callback (int button, int action, int x, int y);


//main functions
void draw_orientation();
void draw_pointer();
void draw_fps();
void draw_extra();
void draw_sphere(const MC_v3d& center,const double& radius);

int setup_opengl(int argc,char **argv);
int draw_scene(); //draw the model
void init_parameters();
double get_rand();//just a rand
int evolve();



// GLOBAL VARIABLES
// ************************************ //
bool is_shift_on=false;
bool is_draw_grid=false;
bool is_draw_sphere=true;
bool is_run=false;
bool is_draw_force=false;
bool is_draw_isosurface=false;
bool is_plane=false;
bool is_plane_cleared;
bool is_throw_ball=false;
bool is_collision_detection=true;
bool is_potential=false;

unsigned int counter_global=0;

GLuint id_draw_sphere=0;
unsigned int counter=0;
mesh_conv::MC_navigator_tool navigator;

unsigned int thrower_timer=0;

mesh_conv::MC_particle_engine engine;
mesh_conv::MC_v3d_vector COLOR_TABLE;
mesh_conv::MC_grid_3d_scalar potential;
mesh_conv::MC_mesh_index_vector_draw mesh_iso;
double isovalue=2;
unsigned int N_smooth=2;
unsigned int thrower_time=90;
double sphere_radius=0.03;

unsigned int dim_grid_3d=30;


double plane_Lz=0.0;
double plane_L=0.5;
mesh_conv::MC_mesh_index_vector_draw plane_mesh;


void update_potential();
void update_marching_cube();
void draw_plane();
void plane_clear();
void update_plane(const unsigned int& N);
void update_plane_density();
void export_scene(unsigned int number);

void engine_run();
unsigned int engine_timer=0;
void set_light();

int main(int argc, char *argv[])
{
  cout<<"**hello**\n"<<endl;
  std::cout<<"run "<<argv[0]<<" with "<<argc-1<<" parameters ... \n"<<std::endl;


  //load camera & light
  std::string file_cam("saved_position/saved_camera.dat");
  navigator.read_cam1_file(file_cam);


  //*************//
  //set parameters
  //*************//

  //size of the 3D grid of the potential
  dim_grid_3d=20;
  //frequency of the creation of new sphere
  thrower_time=90;
  //radius of the spheres
  sphere_radius=0.03;
  //isovalue
  isovalue=2;
  //smooth the isosurface N times
  N_smooth=2;
  //does the sphere collide
  is_collision_detection=true;
  //does the potential interaction force computed
  is_potential=false;



  //init two spheres
//  engine.add(MC_sphere_physic(MC_v3d(0.5,0.2,0.5),sphere_radius,MC_v3d(0,0.1,1)));
//  engine.add(MC_sphere_physic(MC_v3d(0.5,0.8,0.5),sphere_radius,MC_v3d(0,-0.1,1)));
  engine.add(MC_sphere_physic(MC_v3d(0.5,0.2,0.5),sphere_radius,MC_v3d(0,1,1)));
  engine.add(MC_sphere_physic(MC_v3d(0.5,0.8,0.5),sphere_radius,MC_v3d(0,-1,1)));


  //*************//
  //start OpenGL
  //*************//
  navigator.x_screen_size_1=800;
  navigator.y_screen_size_1=800;
  setup_opengl(argc,argv);



  //*************//
  //load  models
  //*************//
  init_parameters();

  //*************//
  //run
  //*************//
  std::cout<<"**RUN**\n"<<std::endl;
  std::cout<<"Navigation with keys:\n"<<std::endl;
  std::cout<<"=========================="<<std::endl;
  std::cout<<"Space Displacement"<<std::endl;
  std::cout<<"=========================="<<std::endl;
  std::cout<<"  z    r"<<std::endl;
  std::cout<<"q s d  f\n"<<std::endl;
  std::cout<<"left mouse: rotate "<<std::endl;
  std::cout<<"right mouse: scale "<<std::endl;
  std::cout<<"shift+left mouse: translate in-plane "<<std::endl;
  std::cout<<"shift+right mouse: translate out-of-plane "<<std::endl<<std::endl;
  std::cout<<"=========================="<<std::endl;
  std::cout<<"View:"<<std::endl;
  std::cout<<"=========================="<<std::endl;
  std::cout<<"W: wireframe mode"<<std::endl;
  std::cout<<"F: filled mode"<<std::endl;
  std::cout<<"G: Grid mode"<<std::endl;
  std::cout<<"X: Print cursor coordinates"<<std::endl;
  std::cout<<"S: Draw spheres"<<std::endl;
  std::cout<<"V: Draw forces"<<std::endl;
  std::cout<<"Q: exit"<<std::endl<<std::endl;
  std::cout<<"=========================="<<std::endl;
  std::cout<<"Run:"<<std::endl;
  std::cout<<"=========================="<<std::endl;
  std::cout<<"R: Start/Stop simulation"<<std::endl;
  std::cout<<"T: Start/Stop add spheres"<<std::endl;
  std::cout<<"I: Compute isosurface"<<std::endl;
  std::cout<<"U: Create plane (+shift mouse left/shift mouse right)"<<std::endl;
  std::cout<<"==========================\n"<<std::endl;
  glutMainLoop ();


  return 0;
}



//update the 3D grid
void update_potential()
{
    //initial value of the grid initialized at -1
    potential.fill(-1.0);

    //get every sphere
    std::map<unsigned int,MC_sphere_physic> particle=engine.particle_system();
    for(std::map<unsigned int,MC_sphere_physic>::const_iterator it=particle.begin(),it_end=particle.end();it!=it_end;++it)
    {
        const MC_sphere_physic& s=it->second;
        implicit::MC_potential_sphere pot(s);

        //add the potential on the sphere
        potential.add_potential(pot);
    }
}

//update the marching cube
void update_marching_cube()
{
    //run the marching cube algorithm
    mesh_iso=MC_grid_3d_scalar_marching_cube::marching_cube(potential,isovalue);

    //smooth the mesh
    //mesh_iso.laplacian_smoothing(0.3,N_smooth,false);
    mesh_iso.normal_vertex_update();
}




//evole through time
void engine_run()
{

    //run particle system
    engine.evolve(is_collision_detection,is_potential);

    //update isosurface
    if(is_draw_isosurface)
    {
        update_potential();
        update_marching_cube();
    }

    //export_scene(counter_global++);

}


//particles evolve through time and add new spheres
int evolve()
{
    //set-up timer
    int ut=glutGet(GLUT_ELAPSED_TIME);


    if(ut-engine_timer>25)
    {
        if(is_run)
        {
            //add new balls
            if(is_throw_ball)
            {
                if(ut-thrower_timer>thrower_time)
                {
                    engine.add(MC_sphere_physic(MC_v3d(0.5,0.5,0.5),sphere_radius,MC_v3d(1,1,2)*1+MC_v3d(1+0.1*cos(navigator.current_time*34232),1+0.2*sin(5245.34*navigator.current_time*34232),1)));
                    thrower_timer=ut;
                }
            }

            //evolve particles
            engine_run();
            engine_timer=ut;
        }
    }

    return 0;
}


//init some parameters
void init_parameters()
{
    COLOR_TABLE=MC_v3d_vector(MC_v3d(1,0,0))<<MC_v3d(0,1,0)<<MC_v3d(0,0,1)<<MC_v3d(0,1,1)<<MC_v3d(1,0,1)<<MC_v3d(1,1,0);
    id_draw_sphere=MC_opengl_drawer::send_gpu(MC_mesh_fast_draw(MC_mesh_index_vector::build_sphere(2)));
    potential=MC_grid_3d_scalar(MC_int_vector(dim_grid_3d,dim_grid_3d,dim_grid_3d),MC_v3d(1,1,1),MC_v3d(0,0,0));

}








void draw_extra()
{
    glEnable(GL_LIGHTING);
    glColor3d(0.8,0.8,0.8);

    //draw what you want

}











//****************************************************************************//
//****************************************************************************//
//****************************************************************************//
//****************************************************************************//









//*******************//
//GLUT STUFF
//*******************//


//glut functions
int setup_opengl(int argc,char **argv)
{

    //init window
    glutInit (&argc, argv);
    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
    glutInitWindowPosition(0,0);
    glutInitWindowSize(navigator.x_screen_size_1,navigator.y_screen_size_1);
    glutCreateWindow("win_1");


    //OpenGL set
    glClearColor(1.0f,1.0f,1.0f,1.0f);
    glPointSize(5);
    glLineWidth(2);

    glShadeModel(GL_SMOOTH);
    glEnable(GL_COLOR_MATERIAL);
    //glColorMaterial(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE );
    glEnable(GL_DEPTH_TEST);
    glDepthMask(GL_TRUE);


    glFrontFace(GL_CCW);//clockwise CW or CCW
    glDisable(GL_CULL_FACE);
    //glEnable(GL_CULL_FACE);

    glDisable(GL_BLEND);
    glDisable(GL_TEXTURE_2D);
    glDisable(GL_FOG);

    glEnable(GL_LIGHTING);
    glEnable(GL_LIGHT0);
    glFrontFace( GL_CW ) ;
    glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, GL_TRUE);

    set_light();




    glEnable(GL_NORMALIZE);





    glEnable(GL_LINE_SMOOTH);
    glHint(GL_LINE_SMOOTH_HINT,GL_NICEST);

    //callback
    glutDisplayFunc (display_callback);
    glutReshapeFunc (reshape_callback);
    glutKeyboardFunc(keyboard_callback);
    glutTimerFunc(33,idle_callback,1);
    glutMotionFunc (motion_callback);
    glutMouseFunc (mouse_callback);

    return 0;
}



//display
static void display_callback()
{


    //vieport (all screen)
    glViewport (0, 0, navigator.x_screen_size_1,navigator.y_screen_size_1);
    glClear (GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    //Projection
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity ();
    mesh_conv::MC_double_vector proj_param=navigator.projection_param();
    glFrustum(proj_param[0], -proj_param[0], proj_param[1], -proj_param[1], proj_param[2], proj_param[3]);


    //Modelview
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();

    glTranslated(0,0,static_cast<GLdouble>(navigator.dist_cam1)); //zoom
    MC_matrix temp_rotation = navigator.trackball_cam1.quaternion().conjugated().matrix().to_matrix4();
    glMultMatrixd(temp_rotation.pointer()); //rotation
    glTranslated(static_cast<GLdouble>(navigator.tr_cam1[0]),static_cast<GLdouble>(navigator.tr_cam1[1]),static_cast<GLdouble>(navigator.tr_cam1[2])); //translation


    //draw
    draw_scene();

    //draw_pointer
    draw_pointer();

    //draw orientation (clear DEPTH-BUFFER)
    draw_orientation();

    //update fps
    navigator.increase_frame_number();
    navigator.current_time=glutGet(GLUT_ELAPSED_TIME);
    if(navigator.delta_time()>1000)
    {navigator.update_fps();}

    //draw fps
    draw_fps();

    //error OpenGL_drawer_2
    GLenum errCode;
    const GLubyte *errString;
    if((errCode=glGetError()) != GL_NO_ERROR)
    {
        errString = gluErrorString(errCode);
        std::cout<<"OPENGL ERROR : "<<errString<<std::endl;
        exit(-2);
    }

    counter++;


  //redraw
  glFlush ();
  glutSwapBuffers ();

}


static void reshape_callback (int width, int height)
{
    navigator.x_screen_size_1=width;navigator.y_screen_size_1=height;
    glViewport(0,0, width, height);
}


//keyboard
static void keyboard_callback(unsigned char key,int x,int y)
{

    x=0;y=0;
    double dL=0.001*(1+100*abs(navigator.dist_cam1));

#define GLUT_SPACE 32

    int mod = glutGetModifiers();
    if((mod&GLUT_ACTIVE_CTRL)!=0)
    {
        std::cout<<"ctrl on"<<std::endl;

    }


    switch(key)
    {
    case 'Q'://quit
        {
            std::cout<<"\n[EXIT OK]\n\n"<<std::endl;

            std::string file_cam("saved_position/saved_camera.dat");
            navigator.write_cam1_file(file_cam);

            exit(0);
        }
        break;


    case 'z':
        navigator.go_forward_trackball_cam1(dL);
        break;
    case 's':
        navigator.go_forward_trackball_cam1(-dL);
        break;
    case 'q':
        navigator.go_right_trackball_cam1(-dL);
        break;
    case 'd':
        navigator.go_right_trackball_cam1(dL);
        break;
    case 'r':
        navigator.go_up_trackball_cam1(dL);
        break;
    case 'f':
        navigator.go_up_trackball_cam1(-dL);
        break;

    case 'G':
        is_draw_grid==1?is_draw_grid=0:is_draw_grid=1;
        break;

    case 'W'://wireframe
        glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
        break;
    case 'F'://Filled
        glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
        break;

    case 'X'://cursor position
        std::cout<<"cursor position: "<<-navigator.tr_cam1<<std::endl;
        break;

    case 'R'://run
        is_run=!is_run;
        std::cout<<"run: "<<is_run<<std::endl;
        break;

    case 'V': //draw forces
        is_draw_force=!is_draw_force;
        break;

    case 'S': //draw sphere
        is_draw_sphere=!is_draw_sphere;
        break;

    case 'I': //isosurface
        is_draw_isosurface=!is_draw_isosurface;
        if(is_draw_isosurface==true)
        {
            update_potential();
            update_marching_cube();
        }
        break;

    case 'U':
        {
            is_plane==true?is_plane=false:is_plane=true;
            is_plane_cleared=false;
            update_plane(2);
            if(is_plane==false)
            {
                update_plane(180);
                update_plane_density();
            }
        }
        break;

            case 'T': //throw ball
                {
                    is_throw_ball=!is_throw_ball;

                }
                break;

                case 'E': //save
                export_scene(counter_global++);
            }



    glutPostRedisplay ();
}


static void idle_callback(int value)
{
  value=1;//useless

  glutPostRedisplay ();
  glutTimerFunc( 10, idle_callback, 1);


  evolve();
}



static void motion_callback (int x, int y)
{
    double x_old=navigator.x_previous,y_old=navigator.y_previous;
    double w=navigator.x_screen_size_1,h=navigator.y_screen_size_1;



    if ( is_shift_on==false && (navigator.left_button) )
    {
        //trackball
        double epsilon=1e-6;
        double x0=(2.0*x_old-w)/w,y0=(2.0*y_old-h)/h,x1=(2*x-w)/w,y1=(2*y-h)/h;
        if(std::sqrt((x0-x1)*(x0-x1)+(y0-y1)*(y0-y1))>epsilon)
        {
            navigator.trackball_cam1.set_2d_coords((2.0*x_old-w)/w,(2.0*y_old-h)/h,(2*x-w)/w,(2*y-h)/h);
            navigator.trackball_cam1.apply_rotation();
        }
        navigator.x_previous=int(x);navigator.y_previous=int(y);

    }
    else if ( is_shift_on==false && (navigator.right_button))
    {
        double u = y-y_old;
        double fu = u/500.0; //factor of magnification

        navigator.dist_cam1 += (fabs(navigator.dist_cam1)+1)*fu; navigator.dist_cam1 = min(navigator.dist_cam1,0.0);
        navigator.y_previous=int(y);

    }
    if(is_plane==false && is_shift_on==true && (navigator.left_button) )
    {
        //PAN
        double dL=0.0001*(1+100*abs(navigator.dist_cam1));//move camera

        MC_v3d tr=MC_v3d(x-x_old,y-y_old,0.0);
        navigator.tr_cam1 += dL*navigator.trackball_cam1.quaternion().conjugated()*tr;
        navigator.x_previous=int(x);navigator.y_previous=int(y);
    }
    else if( is_plane==false && is_shift_on && (navigator.right_button) )
    {

        //PAN
        double dL=0.0002*(1+100*abs(navigator.dist_cam1));//move camera

        MC_v3d tr=MC_v3d(0,0,y-y_old);
        navigator.tr_cam1 += dL*navigator.trackball_cam1.quaternion().conjugated()*tr;
        navigator.x_previous=int(x);navigator.y_previous=int(y);
    }


    if(is_plane==true && is_shift_on && navigator.left_button )
    {
        plane_Lz += 0.01*static_cast<double>(y-y_old);
        navigator.x_previous=int(x);navigator.y_previous=int(y);
        update_plane(2);
    }
    else if(is_plane==true && is_shift_on && navigator.right_button )
    {
        double d1=x-x_old;
        double d2=y-y_old;
        double d;

        plane_L += 0.01*(d1+d2);
        navigator.x_previous=int(x);navigator.y_previous=int(y);
        update_plane(2);
    }

    glutPostRedisplay ();
}
static void mouse_callback (int button, int action, int x, int y)
{
    navigator.x_previous=x;
    navigator.y_previous=y;

    int mod = glutGetModifiers();
    if((mod&GLUT_ACTIVE_SHIFT)!=0)
    {
        is_shift_on=true;
    }
    else
        is_shift_on=false;

  if (button == GLUT_LEFT_BUTTON && action == GLUT_DOWN) {
    navigator.left_button=true;
  }
  if (button == GLUT_LEFT_BUTTON && action == GLUT_UP)
    navigator.left_button=false;

  if (button == GLUT_RIGHT_BUTTON && action == GLUT_DOWN) {
    navigator.right_button=true;
  }
  if (button == GLUT_RIGHT_BUTTON && action == GLUT_UP)
    navigator.right_button=false;

  glutPostRedisplay ();
}

//just a rand number \in [0,1]
double get_rand()
{
  int precision=50000;
  return double(rand()%precision)/double(precision);
}


void draw_orientation()
{
    glDisable(GL_LIGHTING);
    //rotation
    glClear(GL_DEPTH_BUFFER_BIT);
    glMatrixMode(GL_PROJECTION);
    glPushMatrix();
    glLoadIdentity ();
    glOrtho(-1,1,1,-1,0.1,10.0);
    glMatrixMode(GL_MODELVIEW);
    glPushMatrix();
    glLoadIdentity();
    glTranslated(0.78,0.78,-1);
    // get pointer
    MC_matrix M=navigator.trackball_cam1.quaternion().conjugated().matrix().to_matrix4();
    // rotation
    glMultMatrixd(M.pointer());

    double L=0.2;
    glBegin(GL_LINES);
    glColor3d(1,0,0);
    glVertex3d(0,0,0);
    glVertex3d(L,0,0);
    glColor3d(0,1,0);
    glVertex3d(0,0,0);
    glVertex3d(0,L,0);
    glColor3d(0,0,1);
    glVertex3d(0,0,0);
    glVertex3d(0,0,L);
    glEnd();
    //pop matrices
    glPopMatrix();
    glMatrixMode(GL_PROJECTION);
    glPopMatrix();
    glMatrixMode(GL_MODELVIEW);

}


void draw_pointer()
{

    glPushMatrix();
    glLoadIdentity();

    glTranslatef(0,0,GLfloat(navigator.dist_cam1)); //zoom
    MC_matrix temp_m=navigator.trackball_cam1.quaternion().conjugated().matrix().to_matrix4();
    glMultMatrixd(temp_m.pointer()); //rotation

    //local basis
    double L=0.01;

    glBegin(GL_LINES);
    glColor3d(1,0,0);
    glVertex3d(-L,0,0);
    glVertex3d(L,0,0);
    glColor3d(0,1,0);
    glVertex3d(0,-L,0);
    glVertex3d(0,L,0);
    glColor3d(0,0,1);
    glVertex3d(0,0,-L);
    glVertex3d(0,0,L);
    glEnd();

    glPopMatrix();

}


void draw_fps()
{
    glDisable(GL_DEPTH_TEST);

    int int_fps=int(navigator.current_fps);
    int deci_fps=int(10*(navigator.current_fps-int_fps));
    double truncated_fps=int_fps+double(deci_fps)/10.0;

    std::string fps=(mesh_conv::MC_string_converter::to_string(truncated_fps)+" fps");

    glPushMatrix();
    glLoadIdentity();
    glRasterPos3f(-0.5f,-0.5f,-1.4f);
    for (unsigned int k=0,N=fps.length();k<N;++k)
        glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_24,fps[k]);
    glPopMatrix();

    std::string particles=(mesh_conv::MC_string_converter::to_string(engine.particle_system().size())+" particles");
    glPushMatrix();
    glLoadIdentity();
    glRasterPos3f(-0.5f,+0.5f,-1.4f);
    for (unsigned int k=0,N=particles.length();k<N;++k)
        glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_24,particles[k]);
    glPopMatrix();

    glEnable(GL_DEPTH_TEST);
}



//Draw the scene
int draw_scene()
{

    set_light();


    //draw spheres
    glColor3d(1,0,0);
    if(is_draw_sphere)
    {
        const std::map<unsigned int,MC_sphere_physic>& particles=engine.particle_system();
        for(std::map<unsigned int,MC_sphere_physic>::const_iterator it=particles.begin(),it_end=particles.end();it!=it_end;++it)
        {
            const unsigned int& id=it->first;
            glColor3dv(COLOR_TABLE[id%COLOR_TABLE.size()].pointer());


            const MC_sphere_physic& current_particle=it->second;

            const MC_v3d& c=current_particle.center();
            const double& r=current_particle.radius();
            glPushMatrix();
            glTranslated(c[0],c[1],c[2]);
            glScaled(r,r,r);
            glCallList(id_draw_sphere);
            glPopMatrix();

            //draw forces
            glColor3d(0,0,1);
            if(is_draw_force)
            {
                glDisable(GL_LIGHTING);
                MC_segment f(c,c+0.01*current_particle.force());
                MC_opengl_drawer::draw(f);
            }

        }
    }

    //draw bounding cube
    glColor3d(0,0,0);
    glDisable(GL_LIGHTING);
    MC_opengl_drawer::draw_grid(MC_mesh_index_vector::build_cube());



    //draw isosurface
    glEnable(GL_LIGHTING);
    if(is_draw_isosurface)
    {
        glEnable(GL_POLYGON_OFFSET_FILL);
        glPolygonOffset(2,2);
        glColor3d(0.5,0.5,0.8);
        glEnable(GL_LIGHTING);
        MC_opengl_drawer::draw(mesh_iso);

        glDisable(GL_POLYGON_OFFSET_FILL);

        if(is_draw_grid)
        {
            glColor3d(0.7,0.1,0.1);
            glDisable(GL_LIGHTING);
            MC_opengl_drawer::draw_grid(mesh_iso);
        }


    }


    //draw plane
    {

        if(is_plane==true || is_plane_cleared==false)
        {
            glDisable(GL_LIGHTING);
            glColor3d(0.8,0.8,0.8);
            glEnable(GL_POLYGON_OFFSET_FILL);
            glPolygonOffset(2,2);
            MC_opengl_drawer::draw(plane_mesh);
            glDisable(GL_POLYGON_OFFSET_FILL);

            glDisable(GL_LIGHTING);
            glColor3d(0.2,0.2,1);
            if(is_draw_grid)
                MC_opengl_drawer::draw_grid(plane_mesh);
        }
    }



    draw_extra();

  return 0;
}

void set_light()
{
    MC_v3d_vector axis_cam=navigator.axis_cam1();

    MC_v3d Light_pos=MC_v3d(navigator.current_cam1().to_matrix3().first*(axis_cam[0]))+axis_cam[1]*10;//+eye_n*(+15.0)+eye_up*10.0;
    //MC_v3d Light_pos=axis_cam[0]+axis_cam[1]*0.5+axis_cam[2]*0.5;

    //light
    glEnable(GL_LIGHTING);
    glEnable(GL_LIGHT0);

    float light_ambient[4]={0.4f,0.4f,0.4f,1.0f};
    float light_diffuse[4]={0.4f,0.4f,0.4f,1.0f};
    float light_specular[4]={1.6f,1.6f,1.6f,1.0f};
    float light_position[4]={Light_pos[0],Light_pos[1],Light_pos[2],1.0};

    float material_ambient [4] = {1.0f,1.0f,1.0f,1.0f};
    float material_diffuse [4] = {1.0f,1.0f,1.0f,1.0f};
    float material_specular [4] = {1.0f,1.0f,1.0f,1.0f};
    float material_exponent [1] = {100.0f};



    glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, GL_TRUE);
    glLightfv(GL_LIGHT0, GL_POSITION, light_position);
    glLightfv(GL_LIGHT0, GL_AMBIENT, light_ambient);
    glLightfv(GL_LIGHT0, GL_DIFFUSE, light_diffuse);
    glLightfv(GL_LIGHT0, GL_SPECULAR, light_specular);

    glEnable(GL_COLOR_MATERIAL);
    glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT, material_ambient);
    glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, material_diffuse);
    glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, material_specular);
    glMaterialfv(GL_FRONT_AND_BACK, GL_SHININESS, material_exponent);
}






void update_plane(const unsigned int& N)
{
    //int N=20;
    plane_mesh=MC_mesh_index_vector::build_square(N,N);
    plane_mesh+=MC_v3d(-0.5,-0.5,0);

    MC_v3d_vector axes=navigator.axis_cam1();
    MC_v3d normal=axes[1];
    MC_v3d right=axes[2];
    MC_v3d up=axes[3];
    MC_v3d center=-navigator.tr_cam1;

    MC_matrix R1=MC_matrix::rotation_axis_to_axis(MC_v3d(0,0,1),normal);
    MC_matrix R2=MC_matrix::rotation_axis_to_axis(R1*MC_v3d(0,1,0),up);

    plane_mesh=R2*(R1*plane_L)*plane_mesh+center+plane_Lz*normal;
}
void plane_clear()
{
    plane_L=0.5;
    plane_Lz=0.0;
    is_plane_cleared=true;
}

void update_plane_density()
{

    plane_mesh.color().resize(plane_mesh.point_set().size());
    for(int k=0,N=plane_mesh.vertex_number();k<N;++k)
    {

        MC_v3d x=plane_mesh.point_set()(k)-MC_v3d(0.0,0.0,0);

        double pot=potential(x);

        double b=0.0;
        double r=0.0;
        if(std::fabs(pot-isovalue)/isovalue<0.2)
        {
            r=0.0;
            b=1.0;
        }
        else
        {
            double u=std::fabs(pot)/(isovalue);
            if(u>1)
                u=-log(u)+1;
            r=pow(sin(M_PI/2.0*u),2.0);
        }

        plane_mesh.color()(k)=MC_v3d(r,0,b);
    }
}

void export_scene(unsigned int number)
{
    std::string dir="exported_data/";
    std::string number_str=MC_string_converter::zero_padding(MC_string_converter::to_string(number),5);

    //export mesh
    MC_io_obj::write_obj_file(dir+"mesh_"+number_str+".obj",mesh_iso);

    //export sphere
    std::string sphere_file=dir+"sphere_"+number_str+".sph";
    std::ofstream file(sphere_file.c_str());
    const std::map<unsigned int,MC_sphere_physic>& particles=engine.particle_system();
    for(std::map<unsigned int,MC_sphere_physic>::const_iterator it=particles.begin(),it_end=particles.end();it!=it_end;++it)
        file<<it->first<<" "<<it->second.center()[0]<<" "<<it->second.center()[1]<<" "<<it->second.center()[2]<<" "<<it->second.radius()<<std::endl;
    file.close();
}