379 lines
11 KiB
C
379 lines
11 KiB
C
/*
|
|
* Lorenz Strange Attractor
|
|
*
|
|
* Written by John F. Fay in honor of the "freeglut" 2.0.0 release in July 2003
|
|
*
|
|
* What it does:
|
|
* This program starts with two particles right next to each other. The particles
|
|
* move through a three-dimensional phase space governed by the following equations:
|
|
* dx/dt = sigma * ( y - x )
|
|
* dy/dt = r * x - y + x * z
|
|
* dz/dt = x * y + b * z
|
|
* These are the Lorenz equations and define the "Lorenz Attractor." Any two particles
|
|
* arbitrarily close together will move apart as time increases, but their tracks are
|
|
* confined within a region of the space.
|
|
*
|
|
* Commands:
|
|
* Arrow keys: Rotate the view
|
|
* PgUp, PgDn: Zoom in and out
|
|
* Mouse click: Center on the nearest point on a particle trajectory
|
|
*
|
|
* 'r'/'R': Reset the simulation
|
|
* 'm'/'M': Modify the Lorenz parameters (in the text window)
|
|
* 's'/'S': Stop (the advancement in time)
|
|
* 'g'/'G': Go
|
|
* <spacebar>: Single-step
|
|
* <Escape>: Quit
|
|
*/
|
|
|
|
/* Include Files */
|
|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <stdarg.h>
|
|
#include <string.h>
|
|
#include <math.h>
|
|
#include <GL/freeglut.h>
|
|
#ifdef _MSC_VER
|
|
/* DUMP MEMORY LEAKS */
|
|
#include <crtdbg.h>
|
|
#endif
|
|
|
|
|
|
/************************************** Defined Constants ***************************************/
|
|
/* Number of points to draw in the curves */
|
|
#define NUM_POINTS 512
|
|
|
|
/* Angle to rotate when the user presses an arrow key */
|
|
#define ROTATION_ANGLE 5.0
|
|
|
|
/* Amount to scale bu when the user presses PgUp or PgDn */
|
|
#define SCALE_FACTOR 0.8
|
|
|
|
|
|
/*************************************** Global Variables ***************************************/
|
|
/* Lorenz Attractor variables */
|
|
double s0 = 10.0, r0 = 28.0, b0 = 8.0/3.0 ; /* Default Lorenz attactor parameters */
|
|
double time_step = 0.03 ; /* Time step in the simulation */
|
|
double sigma = 10.0, r = 28.0, b = 8.0/3.0 ; /* Lorenz attactor parameters */
|
|
double red_position[NUM_POINTS][3] ; /* Path of the red point */
|
|
double grn_position[NUM_POINTS][3] ; /* Path of the green point */
|
|
int array_index ; /* Position in *_position arrays of most recent point */
|
|
double distance = 0.0 ; /* Distance between the two points */
|
|
|
|
/* GLUT variables */
|
|
double yaw = 0.0, pit = 0.0 ; /* Euler angles of the viewing rotation */
|
|
double scale = 1.0 ; /* Scale factor */
|
|
double xcen = 0.0, ycen = 0.0, zcen = 0.0 ; /* Coordinates of the point looked at */
|
|
|
|
int animate = 1 ; /* 0 - stop, 1 = go, 2 = single-step */
|
|
|
|
|
|
/******************************************* Functions ******************************************/
|
|
|
|
/* The Lorenz Attractor */
|
|
void calc_deriv ( double position[3], double deriv[3] )
|
|
{
|
|
/* Calculate the Lorenz attractor derivatives */
|
|
deriv[0] = sigma * ( position[1] - position[0] ) ;
|
|
deriv[1] = ( r + position[2] ) * position[0] - position[1] ;
|
|
deriv[2] = -position[0] * position[1] - b * position[2] ;
|
|
}
|
|
|
|
void advance_in_time ( double time_step, double position[3], double new_position[3] )
|
|
{
|
|
/* Move a point along the Lorenz attractor */
|
|
double deriv0[3], deriv1[3], deriv2[3], deriv3[3] ;
|
|
int i ;
|
|
memcpy ( new_position, position, 3 * sizeof(double) ) ; /* Save the present values */
|
|
|
|
/* First pass in a Fourth-Order Runge-Kutta integration method */
|
|
calc_deriv ( position, deriv0 ) ;
|
|
for ( i = 0; i < 3; i++ )
|
|
new_position[i] = position[i] + 0.5 * time_step * deriv0[i] ;
|
|
|
|
/* Second pass */
|
|
calc_deriv ( new_position, deriv1 ) ;
|
|
for ( i = 0; i < 3; i++ )
|
|
new_position[i] = position[i] + 0.5 * time_step * deriv1[i] ;
|
|
|
|
/* Third pass */
|
|
calc_deriv ( position, deriv2 ) ;
|
|
for ( i = 0; i < 3; i++ )
|
|
new_position[i] = position[i] + time_step * deriv2[i] ;
|
|
|
|
/* Second pass */
|
|
calc_deriv ( new_position, deriv3 ) ;
|
|
for ( i = 0; i < 3; i++ )
|
|
new_position[i] = position[i] + 0.1666666666666666667 * time_step *
|
|
( deriv0[i] + 2.0 * ( deriv1[i] + deriv2[i] ) + deriv3[i] ) ;
|
|
}
|
|
|
|
static void
|
|
checkedFGets ( char *s, int size, FILE *stream )
|
|
{
|
|
if ( fgets ( s, size, stream ) == NULL ) {
|
|
fprintf ( stderr, "fgets failed\n");
|
|
exit ( EXIT_FAILURE );
|
|
}
|
|
}
|
|
|
|
|
|
/* GLUT callbacks */
|
|
|
|
#define INPUT_LINE_LENGTH 80
|
|
|
|
void key_cb ( unsigned char key, int x, int y )
|
|
{
|
|
int i ;
|
|
char inputline [ INPUT_LINE_LENGTH ] ;
|
|
|
|
switch ( key )
|
|
{
|
|
case 'r' : case 'R' : /* Reset the simulation */
|
|
/* Reset the Lorenz parameters */
|
|
sigma = s0 ;
|
|
b = b0 ;
|
|
r = r0 ;
|
|
/* Set an initial position */
|
|
red_position[0][0] = (double)rand() / (double)RAND_MAX ;
|
|
red_position[0][1] = (double)rand() / (double)RAND_MAX ;
|
|
red_position[0][2] = (double)rand() / (double)RAND_MAX ;
|
|
grn_position[0][0] = (double)rand() / (double)RAND_MAX ;
|
|
grn_position[0][1] = (double)rand() / (double)RAND_MAX ;
|
|
grn_position[0][2] = (double)rand() / (double)RAND_MAX ;
|
|
array_index = 0 ;
|
|
/* Initialize the arrays */
|
|
for ( i = 1; i < NUM_POINTS; i++ )
|
|
{
|
|
memcpy ( red_position[i], red_position[0], 3 * sizeof(double) ) ;
|
|
memcpy ( grn_position[i], grn_position[0], 3 * sizeof(double) ) ;
|
|
}
|
|
|
|
break ;
|
|
|
|
case 'm' : case 'M' : /* Modify the Lorenz parameters */
|
|
printf ( "Please enter new value for <sigma> (default %f, currently %f): ", s0, sigma ) ;
|
|
checkedFGets ( inputline, sizeof ( inputline ), stdin ) ;
|
|
sscanf ( inputline, "%lf", &sigma ) ;
|
|
|
|
printf ( "Please enter new value for <b> (default %f, currently %f): ", b0, b ) ;
|
|
checkedFGets ( inputline, sizeof ( inputline ), stdin ) ;
|
|
sscanf ( inputline, "%lf", &b ) ;
|
|
|
|
printf ( "Please enter new value for <r> (default %f, currently %f): ", r0, r ) ;
|
|
checkedFGets ( inputline, sizeof ( inputline ), stdin ) ;
|
|
sscanf ( inputline, "%lf", &r ) ;
|
|
|
|
break ;
|
|
|
|
case 's' : case 'S' : /* Stop the animation */
|
|
animate = 0 ;
|
|
break ;
|
|
|
|
case 'g' : case 'G' : /* Start the animation */
|
|
animate = 1 ;
|
|
break ;
|
|
|
|
case ' ' : /* Spacebar: Single step */
|
|
animate = 2 ;
|
|
break ;
|
|
|
|
case 27 : /* Escape key */
|
|
glutLeaveMainLoop () ;
|
|
break ;
|
|
}
|
|
}
|
|
|
|
void special_cb ( int key, int x, int y )
|
|
{
|
|
switch ( key )
|
|
{
|
|
case GLUT_KEY_UP : /* Rotate up a little */
|
|
glRotated ( ROTATION_ANGLE, 0.0, 1.0, 0.0 ) ;
|
|
break ;
|
|
|
|
case GLUT_KEY_DOWN : /* Rotate down a little */
|
|
glRotated ( -ROTATION_ANGLE, 0.0, 1.0, 0.0 ) ;
|
|
break ;
|
|
|
|
case GLUT_KEY_LEFT : /* Rotate left a little */
|
|
glRotated ( ROTATION_ANGLE, 0.0, 0.0, 1.0 ) ;
|
|
break ;
|
|
|
|
case GLUT_KEY_RIGHT : /* Rotate right a little */
|
|
glRotated ( -ROTATION_ANGLE, 0.0, 0.0, 1.0 ) ;
|
|
break ;
|
|
|
|
case GLUT_KEY_PAGE_UP : /* Zoom in a little */
|
|
glScaled ( 1.0 / SCALE_FACTOR, 1.0 / SCALE_FACTOR, 1.0 / SCALE_FACTOR ) ;
|
|
break ;
|
|
|
|
case GLUT_KEY_PAGE_DOWN : /* Zoom out a little */
|
|
glScaled ( SCALE_FACTOR, SCALE_FACTOR, SCALE_FACTOR ) ;
|
|
break ;
|
|
}
|
|
|
|
glutPostRedisplay () ;
|
|
}
|
|
|
|
void mouse_cb ( int button, int updown, int x, int y )
|
|
{
|
|
if ( updown == GLUT_DOWN )
|
|
{
|
|
/*double dist = 1.0e20 ; A very large number */
|
|
/* The idea here is that we go into "pick" mode and pick the nearest point
|
|
to the mouse click position. Unfortunately I don't have the time to implement
|
|
it at the moment. */
|
|
}
|
|
}
|
|
|
|
void draw_curve ( int index, double position [ NUM_POINTS ][3] )
|
|
{
|
|
int i = index ;
|
|
|
|
glBegin ( GL_LINE_STRIP ) ;
|
|
do
|
|
{
|
|
i = ( i == NUM_POINTS-1 ) ? 0 : i + 1 ;
|
|
glVertex3dv ( position[i] ) ;
|
|
}
|
|
while ( i != index ) ;
|
|
|
|
glEnd () ;
|
|
}
|
|
|
|
void bitmapPrintf (const char *fmt, ...)
|
|
{
|
|
static char buf[256];
|
|
va_list args;
|
|
|
|
va_start(args, fmt);
|
|
#if defined(WIN32) && !defined(__CYGWIN__)
|
|
(void) _vsnprintf (buf, sizeof(buf), fmt, args);
|
|
#else
|
|
(void) vsnprintf (buf, sizeof(buf), fmt, args);
|
|
#endif
|
|
va_end(args);
|
|
glutBitmapString ( GLUT_BITMAP_HELVETICA_12, (unsigned char*)buf ) ;
|
|
}
|
|
|
|
void display_cb ( void )
|
|
{
|
|
glClear ( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT ) ;
|
|
|
|
glColor3d ( 1.0, 1.0, 1.0 ) ; /* White */
|
|
/* Draw some axes */
|
|
glBegin ( GL_LINES ) ;
|
|
glVertex3d ( 0.0, 0.0, 0.0 ) ;
|
|
glVertex3d ( 2.0, 0.0, 0.0 ) ;
|
|
glVertex3d ( 0.0, 0.0, 0.0 ) ;
|
|
glVertex3d ( 0.0, 1.0, 0.0 ) ;
|
|
glVertex3d ( 0.0, 0.0, 0.0 ) ;
|
|
glVertex3d ( 0.0, 0.0, 1.0 ) ;
|
|
glEnd () ;
|
|
|
|
glColor3d ( 1.0, 0.0, 0.0 ) ; /* Red */
|
|
draw_curve ( array_index, red_position ) ;
|
|
|
|
glColor3d ( 0.0, 1.0, 0.0 ) ; /* Green */
|
|
draw_curve ( array_index, grn_position ) ;
|
|
|
|
/* Print the distance between the two points */
|
|
glColor3d ( 1.0, 1.0, 1.0 ) ; /* White */
|
|
glRasterPos2i ( 1, 1 ) ;
|
|
bitmapPrintf ( "Distance: %10.6f", distance ) ;
|
|
|
|
glutSwapBuffers();
|
|
}
|
|
|
|
void reshape_cb ( int width, int height )
|
|
{
|
|
float ar;
|
|
glViewport ( 0, 0, width, height ) ;
|
|
glMatrixMode ( GL_PROJECTION ) ;
|
|
glLoadIdentity () ;
|
|
ar = (float) width / (float) height ;
|
|
glFrustum ( -ar, ar, -1.0, 1.0, 10.0, 100.0 ) ;
|
|
glMatrixMode ( GL_MODELVIEW ) ;
|
|
glLoadIdentity () ;
|
|
xcen = 0.0 ;
|
|
ycen = 0.0 ;
|
|
zcen = 0.0 ;
|
|
glTranslated ( xcen, ycen, zcen - 50.0 ) ;
|
|
}
|
|
|
|
|
|
void timer_cb ( int value )
|
|
{
|
|
/* Function called at regular intervals to update the positions of the points */
|
|
double deltax, deltay, deltaz ;
|
|
int new_index = array_index + 1 ;
|
|
|
|
/* Set the next timed callback */
|
|
glutTimerFunc ( 30, timer_cb, 0 ) ;
|
|
|
|
if ( animate > 0 )
|
|
{
|
|
if ( new_index == NUM_POINTS ) new_index = 0 ;
|
|
advance_in_time ( time_step, red_position[array_index], red_position[new_index] ) ;
|
|
advance_in_time ( time_step, grn_position[array_index], grn_position[new_index] ) ;
|
|
array_index = new_index ;
|
|
|
|
deltax = red_position[array_index][0] - grn_position[array_index][0] ;
|
|
deltay = red_position[array_index][1] - grn_position[array_index][1] ;
|
|
deltaz = red_position[array_index][2] - grn_position[array_index][2] ;
|
|
distance = sqrt ( deltax * deltax + deltay * deltay + deltaz * deltaz ) ;
|
|
|
|
if ( animate == 2 ) animate = 0 ;
|
|
}
|
|
|
|
glutPostRedisplay () ;
|
|
}
|
|
|
|
|
|
|
|
/* The Main Program */
|
|
|
|
int main ( int argc, char *argv[] )
|
|
{
|
|
int pargc = argc ;
|
|
|
|
/* Initialize the random number generator */
|
|
srand ( 1023 ) ;
|
|
|
|
/* Set up the OpenGL parameters */
|
|
glEnable ( GL_DEPTH_TEST ) ;
|
|
glClearColor ( 0.0, 0.0, 0.0, 0.0 ) ;
|
|
glClearDepth ( 1.0 ) ;
|
|
|
|
/* Initialize GLUT */
|
|
glutInitWindowSize ( 600, 600 ) ;
|
|
glutInit ( &pargc, argv ) ;
|
|
glutInitDisplayMode ( GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH ) ;
|
|
|
|
/* Create the window */
|
|
glutCreateWindow ( "Lorenz Attractor" ) ;
|
|
glutKeyboardFunc ( key_cb ) ;
|
|
glutMouseFunc ( mouse_cb ) ;
|
|
glutSpecialFunc ( special_cb ) ;
|
|
glutDisplayFunc ( display_cb ) ;
|
|
glutReshapeFunc ( reshape_cb ) ;
|
|
glutTimerFunc ( 30, timer_cb, 0 ) ;
|
|
|
|
/* Initialize the attractor: The easiest way is to call the keyboard callback with an
|
|
* argument of 'r' for Reset.
|
|
*/
|
|
key_cb ( 'r', 0, 0 ) ;
|
|
|
|
/* Enter the GLUT main loop */
|
|
glutMainLoop () ;
|
|
|
|
#ifdef _MSC_VER
|
|
/* DUMP MEMORY LEAK INFORMATION */
|
|
_CrtDumpMemoryLeaks () ;
|
|
#endif
|
|
|
|
return 0 ;
|
|
}
|
|
|