/* Antoniu Melica Computer Graphics 1 California State University, Stanislaus 6 November, 2000 Project 3 - The Interactive Pendulum This program simulates the natural motion of a pendulum. It allows the user to interact with the pendulum by using the mouse to place the bob at arbitrary initial angles and string lengths. It also allows the user to vary the gravitational field in real-time for values between 0.0 and 20.0 m/s^2. The maximum angle, current angle, angular velocity, angular acceleration, string length, acceleration due to gravity, period, and animation frame rate are displayed on the screen. Derivation of the pendulum motion equation ------------------------------------------------------------------ The period of a pendulum is dependent only on its initial dropoff angle, the string length, and the acceleration due to gravity; mass plays no role. A heavier pendulum does not swing faster than a lighter one for the same reason that a heavier stone does not fall faster than a lighter one: both its weight and its inertia are proportional to its mass. This was demonstrated by Galileo four hundred years ago. During the derivation of the motion equation, the mass drops out. We start with the restoring force on the pendulum (the component of the force of gravity acting on the bob along the tangent of the pendulum arc) F = -mg sin(theta) (1) where m is the mass of the bob, g is the acceleration due to gravity, and theta is the angle the string makes with the down vertical. According to Newton's Second Law, F = ma = m dd(L*theta)/ddt (2) where a is the tangential acceleration, and L is the string length. Equating equations (1) and (2), we see that the mass cancels out: m dd(L*theta)/ddt = -mg sin(theta) dd(L*theta)/ddt = -g sin(theta) Furthermore, L is a constant, so we can take it outside of the derivative and divide both sides by it, giving us the final formula for the pendulum motion (angular acceleration in this case): ..................................... . . . dd(theta)/ddt = -g/L sin(theta) . (3) . . ..................................... There is no analytical solution for equation (3), but good approximations can be calculated using numerical methods. This program uses the 4th order Runge-Kutta numerical integration method to calculate the position of the pendulum with respect to time. The period equation ------------------------------------------------------------------ The formula for the period of the pendulum involves elliptic integrals: ................................................. . . . T = 4 * sqrt(L/g) * K(sin(1/2 * theta_max)) . (4) . . ................................................. where K is the complete elliptic integral of the first kind: pi/2 K(m) = Int [(1/sqrt(1-m^2 * sin^2(theta)))] dtheta 0 This program uses Simpson's Rule to numerically integrate the period equation (4). It is accurate to four significant digits up to 178 degrees. Note that if the pendulum crosses the separatrix (and goes into rotation), the period value displayed is meaningless; it refers to the period of the last theta_max angle recorded. Small angle approximations not used ------------------------------------------------------------------ For small angles of theta (where sin(theta) ~ theta), equation (3) can be approximated by dd(theta)/ddt = -g/L theta (5) and the period T becomes (approximately, for small angles only) T = 2pi sqrt(L/g) (6) This program does not use approximations (5) and (6) since the pendulum swing is not limited to small angles. Conclusion ------------------------------------------------------------------ Therefore, the period of the pendulum is proportional to L and the initial dropoff angle, and inversely proportional to g. Only when the initial dropoff angle is sufficiently small, does its significance diminish in the exact period formula, and the period becomes largely dependent on L and g only (as in equation (6)). For more information, see the following links: http://www.britannica.com/bcom/eb/article/2/0,5716,119062+6+110307,00.html http://www.treasure-troves.com/physics/Pendulum.html ------------------------------------------------------------------ Platform developed on: Win2000 / VC++ 6.0 / Intel PIII 733 / 256 MB RAM / No hardware OpenGL acceleration. OpenGL libraries used: opengl32.lib, glu32.lib, glut32.lib */ #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(); }