//
// A particle motion on a parabola curve.
//
// To compile
// g++ -Rb -lm one_ball_on_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.10.27 (renamed k/m, from "factor" --> "osq")
//
// 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;
double const GRAVITY_CONST = 9.8; // gravity acceleration
// 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;
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 = 80.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 = -10.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->pos[0] = 5.0; // x
particle->pos[1] = 0.0; // x' = v = (dx/dt)
set_color(particle);
particle->sphere_radius = 0.3;
particle->sphere_stacks = 10;
particle->sphere_slices = 20;
run_mode = run_mode_(Stop);
}
void equation_of_motion(double *pos, double *dpos, double dt)
{
// Since Lagrangian is
// L(x,x') = (m/2)*(x'^2+4*x^2*x'^2) - (k/2)*(s-l0)^2,
// s = sqrt(x^4*3*x^2+1),
// Lagrange's eqs. of motion are:
// x' = v
// v' = -(1/(1+4*x^2)*(4*x*v^2 + (k/m)*(s-l0)/s*(2*x^2-3*x))
//
const double L0 = 2.0;
const double K = 0.2;
const double M = 1.0;
double x = pos[0], xsq = x*x;
double v = pos[1], vsq = v*v;
double s, osq = K/M;
s = sqrt(xsq*xsq+3*xsq+1.0);
dpos[0] = ( v ) * dt;
dpos[1] = ( -(1.0/(1+4*xsq)*(4*x*vsq + osq*(s-L0)/s*(2*xsq*x+3*x))) ) * dt;
}
void runge_kutta4_advance(double *p, double *p1, double *dp, double factor)
{
for (int i=0; i<2; 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[2], work_pos[2];
double dpos1[2], dpos2[2], dpos3[2], dpos4[2];
for (int i=0; i<2; i++) {
pos_before[i] = particle->pos[i];
}
//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
for (int i=0; i<2; i++) {
particle->pos[i] = pos_before[i] + ( ONE_SIXTH*dpos1[i]
+ ONE_THIRD*dpos2[i]
+ ONE_THIRD*dpos3[i]
+ ONE_SIXTH*dpos4[i] );
}
}
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_a_ball(particle_ *p)
{
double x = p->pos[0];
glPushMatrix();
glTranslated(x,x*x+1,0.0);
glMaterialfv(GL_FRONT, GL_DIFFUSE, p->color_diffuse);
glutSolidSphere(p->sphere_radius,p->sphere_slices,p->sphere_stacks);
glPopMatrix();
}
void draw_connector_line(particle_ *p)
{
glLineStipple(1, 0x3F07);
glEnable(GL_LINE_STIPPLE);
glPushMatrix();
glBegin(GL_LINES);
glColor3f(0.6,0.6,0.6);
glVertex3f(0.0, 0.0, 0.0);
double x = p->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);
draw_a_ball(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.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);
glPushMatrix();
glColor3f(0.3,0.3,0.6);
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.5, 0.5, 0.5, 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);
particle = new particle_;
initialize_particle(particle);
step = 0;
t = 0.0;
dt = 0.005; // 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;
}