//
// Two particles on two parabola curves.
//
// To compile
// g++ -Rb -lm two_particles_parabola.cpp -framework GLUT -framework OpenGL -framework Foundation
//
// Developed by Akira Kageyama
// on 2010.10.25
// for the lecture "Analytical Mechanics B"
// at Kobe University
// Revised on 2011.11.03
// for SGPESS outreach talk at Kobe Univ.
//
// Methods summary
// Numerical integration: 4th-order Runge-Kutta method.
// Visualization: OpenGL + GLUT
//
// Usage
// Right mouse button: menu
// Drag: rotation of the view
// Ctrl+drag: zoom
// Shift+drag: translation
// Escape: exit (or choose it from the menu)
//
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
// #ifdef MACOSX
#include <GLUT/glut.h>
// #else
// #include <GL/glut.h>
// #endif
// Math and Physics------------------
const double PI = 3.14159265358979323846264338;
// For Calculation-------------------
const int MAX_STEP = 150000;
int step;
double dt, t;
// Particle information--------------
struct particle_ {
double mass; // mass of the particle
double pos[6]; // r, phi, u(=\dot{r}), v(=\dot{phi}).
GLfloat color_diffuse[4];
GLdouble sphere_radius;
GLint sphere_slices, sphere_stacks;
};
struct particle_ * particle[2];
const GLfloat COLOR_OVER_THE_FLOOR[] = {1.0, 1.0, 0.0, 1.0};
const GLfloat COLOR_UNDER_THE_FLOOR[] = {0.5, 0.5, 0.0, 1.0};
// For OpenGL------------------------
float dist = 30.0; // distance
float fovy = 30.0; // field of view
float near = 0.5;
float far = 1000.0;
int height, width;
bool rotation = false;
bool zoom = false;
bool translate = false;
float rotx_angle = 0.0;
float roty_angle = 0.0;
double translate_x = 0.0;
double translate_y = 0.0;
int mouse_first = 0;
GLuint gl_display_list_id_floor;
GLuint gl_display_list_id_axes;
GLuint gl_display_list_id_parabola;
//---- for Manupulation
enum run_mode_ {Stop, Running};
run_mode_ run_mode;
enum mouse_drag_modifier_ {ROTATION, ZOOM, TRANSLATE};
mouse_drag_modifier_ mouse_drag_modifier;
void set_color(particle_ *particle)
{
for (int i=0; i<4; i++) {
particle->color_diffuse[i] = COLOR_OVER_THE_FLOOR[i];
}
}
void initialize_particle(struct particle_ **particle)
{
particle[0]->pos[0] = 2.0; // x
particle[0]->pos[1] = 0.0; // x' = v = (dx/dt)
particle[1]->pos[0] = -3.0; // chaos
particle[1]->pos[1] = 0.0; // x' = v = (dx/dt)
for (int p=0; p<2; p++) {
set_color(particle[p]);
particle[p]->sphere_radius = 0.5;
particle[p]->sphere_stacks = 10;
particle[p]->sphere_slices = 20;
}
run_mode = run_mode_(Stop);
}
void equation_of_motion(double *pos, double *dpos, double dt)
{
// Lagrangian
// L(x0,x1,x0',x1') = (m/2)*(x0'^2+4*x0^2*x0'^2)
// + (m/2)*(x1'^2+4*x1^2*x1'^2)
// - (k/2)*(s-l0)^2
// with
// s = sqrt(dx^2+dy^2), dx=x0-x1, dy=x0^2+x1^2+2
//
const double L0 = 2.0;
const double K = 5.0;
const double M = 0.2;
double x0 = pos[0];
double v0 = pos[1];
double x1 = pos[2];
double v1 = pos[3];
double dx = x0 - x1;
double x0sq = x0*x0, v0sq = v0*v0;
double x1sq = x1*x1, v1sq = v1*v1;
double dy = x0sq + x1sq + 2;
double s = sqrt(dx*dx + dy*dy);
double f0 = (K/M)*(s-L0)/s*( dx+2*x0*dy);
double f1 = (K/M)*(s-L0)/s*(-dx+2*x1*dy);
dpos[0] = ( v0 ) * dt;
dpos[1] = ( -1.0/(1+4*x0sq)*(4*x0*v0sq + f0) ) * dt;
dpos[2] = ( v1 ) * dt;
dpos[3] = ( -1.0/(1+4*x1sq)*(4*x1*v1sq + f1) ) * dt;
}
void runge_kutta4_advance(double *p, double *p1, double *dp, double factor)
{
for (int i=0; i<4; i++) {
p[i] = p1[i] + factor*dp[i];
}
}
void runge_kutta4(struct particle_ **particle, double dt)
{
const double ONE_SIXTH = 1.0/6.0;
const double ONE_THIRD = 1.0/3.0;
double pos_before[4], work_pos[4];
double dpos1[4], dpos2[4], dpos3[4], dpos4[4];
pos_before[0] = particle[0]->pos[0];
pos_before[1] = particle[0]->pos[1];
pos_before[2] = particle[1]->pos[0];
pos_before[3] = particle[1]->pos[1];
//step 1
equation_of_motion(pos_before, dpos1, dt);
runge_kutta4_advance(work_pos, pos_before, dpos1, 0.5);
//step 2
equation_of_motion(work_pos, dpos2, dt);
runge_kutta4_advance(work_pos, pos_before, dpos2, 0.5);
//step 3
equation_of_motion(work_pos, dpos3, dt);
runge_kutta4_advance(work_pos, pos_before, dpos3, 1.0);
//step 4
equation_of_motion(work_pos, dpos4, dt);
//the result
particle[0]->pos[0] = pos_before[0] + ( ONE_SIXTH*dpos1[0]
+ ONE_THIRD*dpos2[0]
+ ONE_THIRD*dpos3[0]
+ ONE_SIXTH*dpos4[0] );
particle[0]->pos[1] = pos_before[1] + ( ONE_SIXTH*dpos1[1]
+ ONE_THIRD*dpos2[1]
+ ONE_THIRD*dpos3[1]
+ ONE_SIXTH*dpos4[1] );
particle[1]->pos[0] = pos_before[2] + ( ONE_SIXTH*dpos1[2]
+ ONE_THIRD*dpos2[2]
+ ONE_THIRD*dpos3[2]
+ ONE_SIXTH*dpos4[2] );
particle[1]->pos[1] = pos_before[3] + ( ONE_SIXTH*dpos1[3]
+ ONE_THIRD*dpos2[3]
+ ONE_THIRD*dpos3[3]
+ ONE_SIXTH*dpos4[3] );
}
void solver()
{
if(step > MAX_STEP || run_mode == Stop) return;
runge_kutta4(particle, dt);
t += dt;
step++;
if(step > MAX_STEP) run_mode = Stop;
glutPostRedisplay();
}
void keyboard(unsigned char key, int x, int y)
{
switch (key) {
case 27:
exit(0);
break;
case ' ':
glutPostRedisplay();
break;
}
}
void resize(int w,int h)
{
glViewport(0, 0, w, h);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
gluPerspective(fovy, (double)w/(double)h, near, far);
// Modelview matrix
glMatrixMode(GL_MODELVIEW);
width=w;
height=h;
}
void draw_balls(particle_ **p)
{
double x;
glPushMatrix();
x = p[0]->pos[0];
glTranslated(x,x*x+1,0.0);
glMaterialfv(GL_FRONT, GL_DIFFUSE, p[0]->color_diffuse);
glutSolidSphere(p[0]->sphere_radius,p[0]->sphere_slices,p[0]->sphere_stacks);
glPopMatrix();
glPushMatrix();
x = p[1]->pos[0];
glTranslated(x,-x*x-1,0.0);
glMaterialfv(GL_FRONT, GL_DIFFUSE, p[1]->color_diffuse);
glutSolidSphere(p[1]->sphere_radius,p[1]->sphere_slices,p[1]->sphere_stacks);
glPopMatrix();
}
void draw_connector_line(particle_ **p)
{
glLineWidth(2.0);
glEnable(GL_LINE_STIPPLE);
glLineStipple (1, 0x0F0F); /* dotted */
glPushMatrix();
glBegin(GL_LINES);
glColor3f(1.0,1.0,1.0);
double x = p[0]->pos[0];
glVertex3f(x,x*x+1,0.0);
x = p[1]->pos[0];
glVertex3f(x,-x*x-1,0.0);
glEnd();
glPopMatrix();
glDisable(GL_LINE_STIPPLE);
}
void draw(void)
{
glClear (GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
gluLookAt(0.0, 0.0, dist*2,
0.0, 0.0, 0.0,
0.0, 1.0, 0.0);
char *s;
glColor3f(1.0,1.0,1.0);
glRasterPos2i(-28, 24);
char title[128];
if (run_mode == Running) sprintf(title,"running");
else if (run_mode == Stop) sprintf(title,"stop");
s = title;
while (*s) {
glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_24, *s);
s++;
}
glLoadIdentity();
gluLookAt(0.0, 0.0, dist,
0.0, 0.0, 0.0,
0.0, 1.0, 0.0);
glTranslated(translate_x, 0.0, 0.0);
glTranslated(0.0, translate_y, 0.0);
glRotatef(roty_angle, 1,0,0);
glRotatef(rotx_angle, 0,0,1);
glTranslated(0.0, 0.0, -10.0);
glPushMatrix();
glCallList(gl_display_list_id_floor);
glCallList(gl_display_list_id_parabola);
glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
glEnable(GL_BLEND);
glEnable(GL_LINE_SMOOTH);
glHint(GL_LINE_SMOOTH_HINT,GL_DONT_CARE);
draw_connector_line(particle);
glEnable(GL_LIGHTING);
// glCallList(gl_display_list_id_axes); do not show axes 111107
draw_balls(particle);
glDisable(GL_BLEND);
glDisable(GL_LIGHTING);
glPopMatrix();
glutSwapBuffers();
}
void initialize_gl_settings_make_floor(int id)
{
const int NX = 101;
const int NY = 101;
const double X_MIN = -far;
const double X_MAX = +far;
const double Y_MIN = -far;
const double Y_MAX = +far;
double dx = (X_MAX-X_MIN)/(NX-1);
double dy = (Y_MAX-Y_MIN)/(NY-1);
double grid_x[NX], grid_y[NY];
for (int i=0; i<NX; i++) {
grid_x[i] = X_MIN + dx*i;
}
for (int j=0; j<NY; j++){
grid_y[j] = Y_MIN + dy*j;
}
glNewList(id, GL_COMPILE);
glPushMatrix();
// Mesh on the floor (z=0 plane)
glColor3f(0.5,0.5,0.5);
glBegin(GL_LINES);
for (int i=0; i<NX; i++) {
glVertex3f(grid_x[i], grid_y[0], 0.0);
glVertex3f(grid_x[i], grid_y[NY-1], 0.0);
}
for (int j=0; j<NY; j++) {
glVertex3f(grid_x[0], grid_y[j], 0.0);
glVertex3f(grid_x[NX-1], grid_y[j], 0.0);
}
glEnd();
glPopMatrix();
glEndList();
}
void initialize_gl_settings_make_axes(int id)
{
const double X_MAX = far;
const double Y_MAX = X_MAX;
int sides = 6;
double dphi = (2*PI)/ sides;
double radius = 0.1;
GLfloat red[] = {1.0, 0.0, 0.0, 1.0};
GLfloat green[] = {0.0, 1.0, 0.0, 1.0};
glNewList(id, GL_COMPILE);
// x-axis
glPushMatrix();
glMaterialfv(GL_FRONT, GL_DIFFUSE, red);
glBegin(GL_QUAD_STRIP);
for (int i = 0; i <= sides; i++) {
double phi = dphi * i;
double y = cos(phi);
double z = sin(phi);
glNormal3d(0, y, z);
glVertex3f(0.0, radius*y, radius*z);
glVertex3f(X_MAX, radius*y, radius*z);
}
glEnd();
glPopMatrix();
// y-axis
glPushMatrix();
glMaterialfv(GL_FRONT, GL_DIFFUSE, green);
glBegin(GL_QUAD_STRIP);
for (int i = 0; i <= sides; i++) {
double phi = dphi * i;
double x = cos(phi);
double z = sin(phi);
glNormal3d(x, 0, z);
glVertex3f(radius*x, 0.0, radius*z);
glVertex3f(radius*x, Y_MAX, radius*z);
}
glEnd();
glPopMatrix();
glEndList();
}
void initialize_gl_settings_make_parabola(int id)
{
const double X_MAX = far;
const double Y_MAX = X_MAX;
int sides = 6;
double dphi = (2*PI)/ sides;
double radius = 0.5;
GLfloat red[] = {1.0, 0.0, 0.0, 1.0};
GLfloat green[] = {0.0, 1.0, 0.0, 1.0};
glNewList(id, GL_COMPILE);
glLineWidth(2.5);
glColor3f(1.0,1.0,1.0);
glPushMatrix();
glBegin(GL_LINE_STRIP);
for (double x = -10.0; x<=10.0; x+=0.1) {
glVertex3f(x,x*x+1,0.0);
}
glEnd();
glPopMatrix();
glPushMatrix();
glBegin(GL_LINE_STRIP);
for (double x = -10.0; x<=10.0; x+=0.1) {
glVertex3f(x,-x*x-1,0.0);
}
glEnd();
glPopMatrix();
glEndList();
}
void initialize_gl_settings(void)
{
GLfloat light_position0[] = { 0.0, 0.0, 1.0, 0.0 };
GLfloat light_ambient[] = { 0.7, 0.7, 0.7, 0.0 };
GLfloat light_diffuse[] = { 1.0, 1.0, 1.0, 0.0 };
GLfloat light_specular[] = { 1.0, 1.0, 1.0, 0.0 };
glClearColor(0.0, 0.0, 0.0, 0.0);
glLightfv(GL_LIGHT0, GL_POSITION, light_position0);
glLightfv(GL_LIGHT0, GL_AMBIENT, light_ambient);
glLightfv(GL_LIGHT0, GL_DIFFUSE, light_diffuse);
glLightfv(GL_LIGHT0, GL_SPECULAR, light_specular);
glEnable(GL_LIGHT0);
GLfloat ambient[4]={0.15,0.15,0.25,0.9};
GLfloat diffuse[4]={0.45,0.45,0.85,0.9};
GLfloat specular[4]={1.0,1.0,1.0,0.9};
glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT, ambient);
glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, diffuse);
glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, specular);
glMaterialf (GL_FRONT, GL_SHININESS, 70.0);
glEnable(GL_DEPTH_TEST);
gl_display_list_id_floor = glGenLists(1);
gl_display_list_id_axes = glGenLists(1);
gl_display_list_id_parabola = glGenLists(1);
initialize_gl_settings_make_floor (gl_display_list_id_floor);
initialize_gl_settings_make_axes (gl_display_list_id_axes);
initialize_gl_settings_make_parabola(gl_display_list_id_parabola);
}
void mouse(int button, int state, int x, int y)
{
int m = glutGetModifiers();
if (m==GLUT_ACTIVE_CTRL) {
mouse_drag_modifier = mouse_drag_modifier_(ZOOM);
} else {
mouse_drag_modifier = mouse_drag_modifier_(ROTATION);
}
if (m==GLUT_ACTIVE_SHIFT) {
mouse_drag_modifier = mouse_drag_modifier_(TRANSLATE);
}
switch (button) {
case GLUT_LEFT_BUTTON:
if (state == GLUT_UP) {
rotation = false;
zoom = false;
translate = false;
} else if(state == GLUT_DOWN){
mouse_first = 1;
switch (mouse_drag_modifier) {
case ROTATION:
rotation = true;
zoom = false;
translate = false;
break;
case ZOOM:
zoom = true;
rotation = false;
translate = false;
break;
case TRANSLATE:
zoom = false;
rotation = false;
translate = true;
break;
}
}
break;
}
glutPostRedisplay();
}
void motion (int x, int y)
{
static int prev_x, prev_y;
double translate_factor = 0.1;
if ( !(rotation | zoom | translate ) ) return;
if (mouse_first == 1) {
prev_x = x; prev_y = y;
mouse_first = 0;
return;
}
if (rotation) {
rotx_angle += (float)(x - prev_x)/width*180.0;
roty_angle += (float)(y - prev_y)/height*180.0;
} else if (zoom) {
dist -= (float)(y - prev_y)/10.0;
if(dist <= 0.0) dist = 0.0;
else if(dist >= far*0.8) dist = far*0.8;
} else if (translate) {
translate_x += (x - prev_x)*translate_factor;
translate_y -= (y - prev_y)*translate_factor;
}
prev_x = x; prev_y = y;
glutPostRedisplay();
}
void reset(struct particle_ **particle)
{
initialize_particle(particle);
step = 0;
t = 0.0;
}
void menu(int value){
switch(value) {
case 1: // initial state
run_mode = run_mode_(Stop);
step = 0;
t = 0.0;
break;
case 2:
run_mode = run_mode_(Running);
break;
case 3:
reset(particle);
break;
case 4:
exit(0);
break;
}
glutPostRedisplay();
}
int main (int argc, char **argv)
{
width = 800;
height = 600;
run_mode = run_mode_(Stop);
mouse_drag_modifier = mouse_drag_modifier_(ROTATION);
for (int p=0; p<2; p++) {
particle[p] = new particle_;
}
initialize_particle(particle);
step = 0;
t = 0.0;
dt = 0.0005; // Find an appropriate value by trial and error.
glutInit (&argc, argv);
glutInitDisplayMode (GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
glutInitWindowSize (width, height);
glutInitWindowPosition (0, 0);
glutCreateWindow ("Particle Motion");
glutKeyboardFunc(keyboard);
glutReshapeFunc(resize);
glutMouseFunc(mouse);
glutMotionFunc(motion);
int Mmenu;
Mmenu = glutCreateMenu(menu);
glutAddMenuEntry("Stop", 1);
glutAddMenuEntry("Start", 2);
glutAddMenuEntry("Reset", 3);
glutAddMenuEntry("Exit", 4);
glutAttachMenu(GLUT_RIGHT_BUTTON);
initialize_gl_settings();
glutDisplayFunc(draw);
glutIdleFunc(solver);
glutMainLoop();
return 0;
}