//
// A particle's motion on a cone
//
// To compile
// g++ -lm particle_on_cone.cpp -framework GLUT -framework OpenGL -framework Foundation
//
// Developed by Akira Kageyama
// on 2010.10.18
// revised on 2011.11.17
// for the lecture "Analytical Mechanics B"
// at Kobe University
//
// 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
const double CONE_SEMI_ANGLE_ALPHA = PI / 12;
// For Calculation-------------------
const int MAX_STEP = 15000;
int step;
double dt, t;
// Particle information--------------
struct particle_ {
double mass; // mass of the particle
double pos[4]; // 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 = 100.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_cone;
//---- 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)
{
// if ( sin(particle->pos[1]) >= 0.0 ) {
for (int i=0; i<4; i++) {
particle->color_diffuse[i] = COLOR_OVER_THE_FLOOR[i];
}
// }
// else {
// for (int i=0; i<4; i++) {
// particle->color_diffuse[i] = COLOR_UNDER_THE_FLOOR[i];
// }
// }
}
void initialize_particle(struct particle_ *particle)
{
particle->pos[0] = 50.0; // r
particle->pos[1] = PI/3; // phi
particle->pos[2] = -10.0; // u (=dr/dt)
particle->pos[3] = 2.0; // v (=dphi/dt)
// particle->pos[2] = 0.0; // u (=dr/dt)
// particle->pos[3] = 0.0; // v (=dphi/dt)
set_color(particle);
particle->sphere_radius = 1.0;
particle->sphere_stacks = 20;
particle->sphere_slices = 40;
run_mode = run_mode_(Stop);
}
void equation_of_motion(double *pos, double *dpos, double dt)
{
double sin_alpha = sin(CONE_SEMI_ANGLE_ALPHA);
double cos_alpha = cos(CONE_SEMI_ANGLE_ALPHA);
double sin_alpha_sq = sin_alpha*sin_alpha;
double r = pos[0];
double phi = pos[1];
double u = pos[2];
double v = pos[3];
dpos[0] = ( u ) * dt;
dpos[1] = ( v ) * dt;
dpos[2] = ( r*sin_alpha_sq*v*v - GRAVITY_CONST*cos_alpha*sin_alpha*(1+cos(phi)) ) * dt;
dpos[3] = ( -2*u*v/r + GRAVITY_CONST*cos_alpha*sin(phi)/ (r*sin_alpha) ) * 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];
for (int i=0; i<4; i++) { // for (x,y,z)-components.
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<4; i++) { // for (r,phi,u,v)
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++;
set_color(particle);
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(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_cone);
glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
glEnable(GL_BLEND);
// ball
glEnable(GL_LIGHTING);
glCallList(gl_display_list_id_axes);
double sin_alpha = sin(CONE_SEMI_ANGLE_ALPHA);
double sin_alpha_sq = sin_alpha*sin_alpha;
double cos_alpha = cos(CONE_SEMI_ANGLE_ALPHA);
double cos_alpha_sq = cos_alpha*cos_alpha;
double r = particle->pos[0];
double phi = particle->pos[1];
double x = r*cos_alpha_sq - r*sin_alpha_sq*cos(phi);
double y = -r*sin_alpha*sin(phi);
double z = r*sin_alpha*(1+cos(phi));
glTranslated(x,y,z);
glMaterialfv(GL_FRONT, GL_DIFFUSE, particle->color_diffuse);
glutSolidSphere(particle->sphere_radius,
particle->sphere_slices,
particle->sphere_stacks);
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_cone(int id)
{
const double R_MAX = far;
int circum_ndiv = 30;
double dphi = (2*PI)/ circum_ndiv;
GLfloat col[] = {0.5, 0.5, 0.5, 1.0};
glNewList(id, GL_COMPILE);
glPushMatrix();
glColor3f(0.5,0.5,0.5);
glBegin(GL_LINES);
for (int i=0; i<circum_ndiv; i++) {
glVertex3f(0.0,0.0,0.0);
float sin_alpha = sin(CONE_SEMI_ANGLE_ALPHA);
float sin_alpha_sq = sin_alpha*sin_alpha;
float cos_alpha = cos(CONE_SEMI_ANGLE_ALPHA);
float cos_alpha_sq = cos_alpha*cos_alpha;
float phi = i*dphi;
float x = R_MAX*cos_alpha_sq - R_MAX*sin_alpha_sq*cos(phi);
float y = -R_MAX*sin_alpha*sin(phi);
float z = R_MAX*sin_alpha*(1+cos(phi));
glVertex3f(x,y,z);
}
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_cone = 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_cone (gl_display_list_id_cone);
}
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();
delete particle;
return 0;
}