1 write a chtlopenGL program simulating the double pendulum

1. write a chtlopenGL program simulating the double pendulum of strings of unit length and unit masses.

Solution

#define OSX #ifdef MAC #include \"glut.h\" #endif #ifdef OSX #include #endif #ifdef UNIX #include #endif #ifdef PC #include #endif #include #include #include // Global variables. int main_window; // The mainframe window. int pendulum_window; // The pendulum window. double pi; // pi double g; // The acceleration due to gravity. double length; // The length of the pendulum string. double theta_max[2]; // The maximum angle of swing (rad). double theta; // The current angle (rad). double omega; // The current angular speed (rad/sec). double alpha; // The current angular acceleration (rad/sec^2). double x; // The x position of the bob. double y; // The y position of the bob. clock_t t; // The time of last display. double period; // The period of the pendulum. int frame; // Frame counter. int frame_rate; // Frame rate (frames per second). clock_t frame_t; // Time used to calculate display frame rate. bool pendulum_adjust; // User is currently adjusting bob position. bool gravity_adjust; // User is currently adjusting gravity. // Function prototypes. int min(int a, int b); void printString_8x13(double x, double y, double z, char *s); void printString_9x15(double x, double y, double z, char *s); void printString_TimesRoman24(double x, double y, double z, char *s); void init(void); void main_display(void); void main_reshape(int width, int height); void main_mouseButton(int button, int state, int x, int y); void main_mouseMotion(int x, int y); void pendulum_display(void); void pendulum_reshape(int width, int height); void pendulum_mouseButton(int button, int state, int x, int y); void pendulum_mouseMotion(int x, int y); void idle(void); void integrate_motion(double & theta, double & omega, double & alpha, double length, double g, double delta_t); void integrate_period(double & period, double theta_max, double length, double g); int min(int a, int b) { if (b < a) return b; else return a; } void printString_8x13(double x, double y, double z, char *s) { glRasterPos3d(x, y, z); while (*s != NULL) glutBitmapCharacter(GLUT_BITMAP_8_BY_13, *s++); } void printString_9x15(double x, double y, double z, char *s) { glRasterPos3d(x, y, z); while (*s != NULL) glutBitmapCharacter(GLUT_BITMAP_9_BY_15, *s++); } void printString_TimesRoman24(double x, double y, double z, char *s) { glRasterPos3d(x, y, z); while (*s != NULL) glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_24, *s++); } void init(void) { // Initialize lighting and material properties. float light_position[] = {150.0, 150.0, 250.0, 1.0}; float light_ambient[] = {0.05, 0.05, 0.05, 1.00}; float light_diffuse[] = {0.65, 0.65, 0.65, 1.00}; float light_specular[] = {1.00, 1.00, 1.00, 1.00}; float material_specular[] = {0.6, 0.6, 0.6, 1.0}; float material_specular_exponent[] = {80.0}; glEnable(GL_LIGHT0); glEnable(GL_COLOR_MATERIAL); glEnable(GL_CULL_FACE); glFrontFace(GL_CCW); 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); glColorMaterial(GL_FRONT, GL_AMBIENT_AND_DIFFUSE); glMaterialfv(GL_FRONT, GL_SPECULAR, material_specular); glMaterialfv(GL_FRONT, GL_SHININESS, material_specular_exponent); glClearColor(0.0, 0.0, 0.0, 1.0); // Initialize variables. pi = 4.0 * atan(1.0); g = 9.8; // 9.8 m/s^2 length = 3.0; // 3 m theta_max[0] = theta_max[1] = 0.0; theta = 0.0; omega = 0.0; period = 0.0; frame = 0; frame_rate = 0; pendulum_adjust = true; gravity_adjust = false; // Calculate initial x y bob coordinates. x = length * sin(theta); y = -length * cos(theta); t = frame_t = clock(); } void main_display(void) { glClearColor(0.831, 0.816, 0.784, 1.0); glClear(GL_COLOR_BUFFER_BIT); glColor3ub(0, 0, 0); // Print title. printString_TimesRoman24(554, 50, 0, \"Pendulum Simulator\"); printString_8x13(554, 80, 0, \"Motion equation numerically\"); printString_8x13(554, 93, 0, \"integrated using 4th order\"); printString_8x13(554, 106, 0, \"Runge-Kutta method.\"); printString_8x13(554, 132, 0, \"Period equation numerically\"); printString_8x13(554, 145, 0, \"integrated using Simpson\'s\"); printString_8x13(554, 158, 0, \"Rule (accurate to four\"); printString_8x13(554, 171, 0, \"significant digits up to\"); printString_8x13(554, 184, 0, \"178 degrees).\"); // Print slider description. printString_8x13(554, 250, 0, \"GRAVITY ADJUSTMENT\"); // Print credits. printString_8x13(554, 500, 0, \"Antoniu Melica, 2000\"); printString_8x13(554, 513, 0, \"amelica@csustan.edu\"); // The slider track. glBegin(GL_LINES); glColor3ub(255, 255, 255); glVertex2f(554, 272); glVertex2f(766, 272); glVertex2f(766, 272); glVertex2f(766, 270); glColor3ub(0, 0, 0); glVertex2f(766, 270); glVertex2f(554, 270); glVertex2f(554, 270); glVertex2f(554, 272); glEnd(); // The slider. glBegin(GL_QUADS); glColor3f(0.831, 0.816, 0.784); glVertex2f(554+g*10, 282); glVertex2f(554+g*10+12, 282); glVertex2f(554+g*10+12, 260); glVertex2f(554+g*10, 260); glEnd(); glBegin(GL_LINES); glColor3ub(0, 0, 0); glVertex2f(554+g*10, 282); glVertex2f(554+g*10+12, 282); glVertex2f(554+g*10+12, 282); glVertex2f(554+g*10+12, 260); glColor3ub(255, 255, 255); glVertex2f(554+g*10+12, 260); glVertex2f(554+g*10, 260); glVertex2f(554+g*10, 260); glVertex2f(554+g*10, 282); glEnd(); glutSwapBuffers(); } void main_reshape(int width, int height) { glViewport(0, 0, width, height); glMatrixMode(GL_PROJECTION); glLoadIdentity(); gluOrtho2D(0, width, height, 0); glMatrixMode(GL_MODELVIEW); glLoadIdentity(); glutPostRedisplay(); } void main_mouseButton(int button, int state, int x, int y) { if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && x >= 554 && x <= 766 && y >= 260 && y <= 282) { gravity_adjust = true; main_mouseMotion(x, y); } else if (button == GLUT_LEFT_BUTTON && state == GLUT_UP) { gravity_adjust = false; } } void main_mouseMotion(int x, int y) { if (gravity_adjust == true) { if (x <= 560) g = 0.0; else if (x >= 760) g = 20.0; else g = (double)(x - 560)/10; glutPostRedisplay(); glutSetWindow(pendulum_window); glutPostRedisplay(); } } void pendulum_display(void) { char status[64]; if (pendulum_adjust == false) ++frame; // Calculate frame rate. if (((clock() - frame_t) / (double)CLOCKS_PER_SEC) >= 1.0) { frame_t = clock(); frame_rate = frame; frame = 0; } glClear(GL_COLOR_BUFFER_BIT); glEnable(GL_LIGHTING); // The pendulum support block. glColor3f(0.8, 0.8, 0.8); glBegin(GL_QUADS); glNormal3f( 0.0, 0.0, 1.0); glVertex3d(-0.3, 0.0, 0.0); glVertex3d( 0.3, 0.0, 0.0); glVertex3d( 0.3, 0.3, 0.0); glVertex3d(-0.3, 0.3, 0.0); glEnd(); // The pendulum line. glColor3f(0.8, 0.8, 0.8); glBegin(GL_LINES); glVertex3d(0.0, 0.0, 0.0); glVertex3d(x, y, 0.0); glEnd(); // The pendulum bob. glColor3f(0.2, 0.65, 0.9); glPushMatrix(); glTranslated(x, y, 0.0); glutSolidSphere(0.32, 30, 30); glPopMatrix(); glDisable(GL_LIGHTING); // Status messages: glColor3f(0.9, 0.9, 0.9); // Print theta max, theta, omega, and alpha. sprintf(status, \"THETA MAX %7.2f degrees\", theta_max[0]*180.0/pi); printString_8x13(-3.8, 3.70, 0, status); sprintf(status, \"THETA %7.2f degrees\", theta*180.0/pi); printString_8x13(-3.8, 3.45, 0, status); sprintf(status, \"OMEGA %10.2f degrees/s\", omega*180.0/pi); printString_8x13(-3.8, 3.20, 0, status); sprintf(status, \"ALPHA %10.2f degrees/s^2\", alpha*180.0/pi); printString_8x13(-3.8, 2.95, 0, status); // Print length, acceleration due to gravity, period, and frame rate. sprintf(status, \"LENGTH %7.2f meters\", length); printString_8x13(0.6, 3.70, 0, status); sprintf(status, \"GRAVITY %7.2f m/s^2\", g); printString_8x13(0.6, 3.45, 0, status); sprintf(status, \"PERIOD %7.2f seconds\", period); printString_8x13(0.6, 3.20, 0, status); sprintf(status, \"FRAME RATE %4d fps\", frame_rate); printString_8x13(0.6, 2.95, 0, status); // Print instructions. sprintf(status, \"Use the mouse to position the pendulum bob.\"); printString_9x15(-3.8, -3.85, 0, status); glutSwapBuffers(); } void pendulum_reshape(int width, int height) { glViewport(0, 0, width, height); glMatrixMode(GL_PROJECTION); glLoadIdentity(); // Scale projection view when window is resized. if (width <= height) { glOrtho(-4.0, 4.0, -4.0 * (GLfloat)height / (GLfloat)width, 4.0 * (GLfloat)height / (GLfloat)width, -10.0, 10.0); } else { glOrtho(-4.0 * (GLfloat)width / (GLfloat)height, 4.0 * (GLfloat)width / (GLfloat)height, -4.0, 4.0, -10.0, 10.0); } glMatrixMode(GL_MODELVIEW); glLoadIdentity(); glutPostRedisplay(); } void pendulum_mouseButton(int button, int state, int x, int y) { if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN) { pendulum_adjust = true; omega = 0.0; alpha = 0.0; period = 0.0; frame = 0; frame_rate = 0; pendulum_mouseMotion(x, y); } else if (button == GLUT_LEFT_BUTTON && state == GLUT_UP) { pendulum_adjust = false; t = frame_t = clock(); } } void pendulum_mouseMotion(int win_x, int win_y) { double adj_x; double adj_y; GLint viewport[4]; if (pendulum_adjust == true) { // Get the current viewport settings. glGetIntegerv(GL_VIEWPORT, viewport); // Convert from top-left origin to center origin. adj_x = (double)(win_x - viewport[2]/2); adj_y = (double)(viewport[3]/2 - win_y); // Calculate theta. if (adj_x == fabs(adj_x)) theta_max[0] = theta_max[1] = theta = atan(adj_y/adj_x) + pi/2.0; else theta_max[0] = theta_max[1] = theta = atan(adj_y/adj_x) - pi/2.0; // Calculate (model view) pendulum length. length = sqrt(adj_x * adj_x + adj_y * adj_y) * 8.0 / min(viewport[2], viewport[3]); // Calculate period. integrate_period(period, theta_max[0], length, g); // Calculate (model view) x y bob coordinates. x = length * sin(theta); y = -length * cos(theta); // Request scene redraw. glutPostRedisplay(); } } void idle(void) { double delta_t; // The time span since last display. // If user is in the process of positioning the bob, there is no // need to calculate the pendulum motion. if (pendulum_adjust == true) return; // Calculate time span since last display and reset time variable. delta_t = (clock() - t) / (double)CLOCKS_PER_SEC; t = clock(); // Calculate theta, omega, and alpha of pendulum. integrate_motion(theta, omega, alpha, length, g, delta_t); // Normalize theta. if (theta > pi) theta = theta - 2*pi; else if (theta < -pi) theta = theta + 2*pi; // Measure theta_max for each swing direction and calculate period. if (theta == fabs(theta) && omega == fabs(omega) && theta > theta_max[1]) { theta_max[1] = theta; } else if (theta == -fabs(theta) && omega == -fabs(omega) && theta < theta_max[1]) { theta_max[1] = theta; } else if (theta_max[0] != theta_max[1]) { theta_max[0] = theta_max[1]; integrate_period(period, theta_max[0], length, g); } // Calculate x and y coordinates of the bob based on theta. x = length * sin(theta); y = -length * cos(theta); // Request scene to be redrawn. glutSetWindow(pendulum_window); glutPostRedisplay(); } void integrate_motion(double & theta, double & omega, double & alpha, double length, double g, double delta_t) { double h = 0.00001; // Step size (1/100,000 sec). double k1, k2, k3, k4; // Intermediate terms for Runge-Kutta method. // Calculate theta and omega for the current time using fourth-order // Runge-Kutta method (RK4) with \'h\' as the stepsize and with theta // and omega as the last calculated (or given) values. // // The equation that is being solved is the second-order, non-linear // differential equation: // // ddtheta/ddt = -(g/l) * sin(theta) // // where theta is the angle the pendulum string makes with the \"down\" // vertical, t is the time, g is the acceleration due to gravity, and // l is the length of the string. // // To solve the equation, it is split into two first-order differential // equations: // // dtheta/dt = omega // domega/dt = -(g/l) sin(theta) // // where omega is the angular speed. // // These equations are solved numerically using RK4 as shown below, // for the time interval delta_t with h as the step size. // while (delta_t >= 0.0) { k1 = h * omega; k2 = h * (omega + h/2.0); k3 = h * (omega + h/2.0); k4 = h * (omega + h); theta += (k1 + 2.0*k2 + 2.0*k3 + k4) / 6.0; k1 = h * (-g/length * sin(theta)); k2 = h * (-g/length * sin(theta + h/2.0)); k3 = h * (-g/length * sin(theta + h/2.0)); k4 = h * (-g/length * sin(theta + h)); omega += (k1 + 2.0*k2 + 2.0*k3 + k4) / 6.0; delta_t -= h; } // Calculate angular acceleration. alpha = -g/length * sin(theta); } void integrate_period(double & period, double theta_max, double length, double g) { double pi; double m; // Period equation parameter for elliptic integral. double M, L, R; // Midpoint, left, and right evaluation terms. double theta; double delta_theta; period = 0.0; pi = 4.0 * atan(1.0); m = sin(0.5*theta_max); // Evaluating at 100 points between 0 and pi/2. delta_theta = (pi/2.0)/100.0; for (theta = delta_theta; theta <= pi/2.0; theta += delta_theta) { // Midpoint rule. M = 1.0 / sqrt(1.0 - pow(m * sin(theta - 0.5*delta_theta), 2.0)); // Left rectangle rule. L = 1.0 / sqrt(1.0 - pow(m * sin(theta - delta_theta), 2.0)); // Right rectangle rule. R = 1.0 / sqrt(1.0 - pow(m * sin(theta), 2.0)); // Simpson\'s Rule. period += (4.0*M + L + R); } // The rest of Simpson\'s Rule. period *= delta_theta / 6.0; // Period equation proportionality constant. period *= 4.0 * sqrt(length/g); } int main(int argc, char *argv[]) { // Initialize GLUT. glutInit(&argc, argv); glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA); glutInitWindowPosition(50, 50); glutInitWindowSize(800, 540); // Set up the main window. main_window = glutCreateWindow(\"Pendulum Simulator\"); glutDisplayFunc(main_display); glutReshapeFunc(main_reshape); glutMouseFunc(main_mouseButton); glutMotionFunc(main_mouseMotion); // Set up the pendulum window. pendulum_window = glutCreateSubWindow(main_window, 20, 20, 500, 500); glutDisplayFunc(pendulum_display); glutReshapeFunc(pendulum_reshape); glutMouseFunc(pendulum_mouseButton); glutMotionFunc(pendulum_mouseMotion); // Register idle callback. glutIdleFunc(idle); // Miscellaneous initializations. init(); // Call the main event loop. glutMainLoop(); }