git.mutantstargoat.com Git - freeglut/blob - progs/demos/Lorenz/lorenz.c

   1 /*
   2  * Lorenz Strange Attractor
   3  *
   4  * Written by John F. Fay in honor of the "freeglut" 2.0.0 release in July 2003
   5  *
   6  * What it does:
   7  *  This program starts with two particles right next to each other.  The particles
   8  *  move through a three-dimensional phase space governed by the following equations:
   9  *       dx/dt = sigma * ( y - x )
  10  *       dy/dt = r * x - y + x * z
  11  *       dz/dt = x * y + b * z
  12  *  These are the Lorenz equations and define the "Lorenz Attractor."  Any two particles
  13  *  arbitrarily close together will move apart as time increases, but their tracks are
  14  *  confined within a region of the space.
  15  *
  16  * Commands:
  17  *  Arrow keys:  Rotate the view
  18  *  PgUp, PgDn:  Zoom in and out
  19  *  Mouse click:  Center on the nearest point on a particle trajectory
  20  *
  21  *  'r'/'R':  Reset the simulation
  22  *  'm'/'M':  Modify the Lorenz parameters (in the text window)
  23  *  's'/'S':  Stop (the advancement in time)
  24  *  'g'/'G':  Go
  25  *  <spacebar>:  Single-step
  26  *  <Escape>:  Quit
  27  */
  28
  29 /* Include Files */
  30 #include <stdio.h>
  31 #include <stdlib.h>
  32 #include <string.h>
  33 #include <math.h>
  34 #include <GL/freeglut.h>
  35 #ifdef WIN32
  36 /* DUMP MEMORY LEAKS */
  37 #include <crtdbg.h>
  38 #endif
  39
  40
  41 /************************************** Defined Constants ***************************************/
  42 /* Number of points to draw in the curves */
  43 #define NUM_POINTS    512
  44
  45 /* Angle to rotate when the user presses an arrow key */
  46 #define ROTATION_ANGLE  5.0
  47
  48 /* Amount to scale bu when the user presses PgUp or PgDn */
  49 #define SCALE_FACTOR     0.8
  50
  51
  52 /*************************************** Global Variables ***************************************/
  53 /* Lorenz Attractor variables */
  54 double s0 = 10.0, r0 = 28.0, b0 = 8.0/3.0 ;   /* Default Lorenz attactor parameters */
  55 double time_step = 0.03 ;                     /* Time step in the simulation */
  56 double sigma = 10.0, r = 28.0, b = 8.0/3.0 ;  /* Lorenz attactor parameters */
  57 double red_position[NUM_POINTS][3] ;          /* Path of the red point */
  58 double grn_position[NUM_POINTS][3] ;          /* Path of the green point */
  59 int array_index ;                             /* Position in *_position arrays of most recent point */
  60 double distance = 0.0 ;                       /* Distance between the two points */
  61
  62 /* GLUT variables */
  63 double yaw = 0.0, pit = 0.0 ;                 /* Euler angles of the viewing rotation */
  64 double scale = 1.0 ;                          /* Scale factor */
  65 double xcen = 0.0, ycen = 0.0, zcen = 0.0 ;   /* Coordinates of the point looked at */
  66
  67 int animate = 1 ;                             /* 0 - stop, 1 = go, 2 = single-step */
  68
  69
  70 /******************************************* Functions ******************************************/
  71
  72 /* The Lorenz Attractor */
  73 void calc_deriv ( double position[3], double deriv[3] )
  74 {
  75   /* Calculate the Lorenz attractor derivatives */
  76   deriv[0] = sigma * ( position[1] - position[0] ) ;
  77   deriv[1] = ( r + position[2] ) * position[0] - position[1] ;
  78   deriv[2] = -position[0] * position[1] - b * position[2] ;
  79 }
  80
  81 void advance_in_time ( double time_step, double position[3], double new_position[3] )
  82 {
  83   /* Move a point along the Lorenz attractor */
  84   double deriv0[3], deriv1[3], deriv2[3], deriv3[3] ;
  85   int i ;
  86   memcpy ( new_position, position, 3 * sizeof(double) ) ;  /* Save the present values */
  87
  88   /* First pass in a Fourth-Order Runge-Kutta integration method */
  89   calc_deriv ( position, deriv0 ) ;
  90   for ( i = 0; i < 3; i++ )
  91     new_position[i] = position[i] + 0.5 * time_step * deriv0[i] ;
  92
  93   /* Second pass */
  94   calc_deriv ( new_position, deriv1 ) ;
  95   for ( i = 0; i < 3; i++ )
  96     new_position[i] = position[i] + 0.5 * time_step * deriv1[i] ;
  97
  98   /* Third pass */
  99   calc_deriv ( position, deriv2 ) ;
 100   for ( i = 0; i < 3; i++ )
 101     new_position[i] = position[i] + time_step * deriv2[i] ;
 102
 103   /* Second pass */
 104   calc_deriv ( new_position, deriv3 ) ;
 105   for ( i = 0; i < 3; i++ )
 106     new_position[i] = position[i] + 0.1666666666666666667 * time_step *
 107                       ( deriv0[i] + 2.0 * ( deriv1[i] + deriv2[i] ) + deriv3[i] ) ;
 108 }
 109
 110 static void
 111 checkedFGets ( char *s, int size, FILE *stream )
 112 {
 113   if ( fgets ( s, size, stream ) == NULL ) {
 114     fprintf ( stderr, "fgets failed\n");
 115     exit ( EXIT_FAILURE );
 116   }
 117 }
 118
 119
 120 /* GLUT callbacks */
 121
 122 #define INPUT_LINE_LENGTH 80
 123
 124 void key_cb ( unsigned char key, int x, int y )
 125 {
 126   int i ;
 127   char inputline [ INPUT_LINE_LENGTH ] ;
 128
 129   switch ( key )
 130   {
 131   case 'r' :  case 'R' :  /* Reset the simulation */
 132     /* Reset the Lorenz parameters */
 133     sigma = s0 ;
 134     b = b0 ;
 135     r = r0 ;
 136     /* Set an initial position */
 137     red_position[0][0] = (double)rand() / (double)RAND_MAX ;
 138     red_position[0][1] = (double)rand() / (double)RAND_MAX ;
 139     red_position[0][2] = (double)rand() / (double)RAND_MAX ;
 140     grn_position[0][0] = (double)rand() / (double)RAND_MAX ;
 141     grn_position[0][1] = (double)rand() / (double)RAND_MAX ;
 142     grn_position[0][2] = (double)rand() / (double)RAND_MAX ;
 143     array_index = 0 ;
 144     /* Initialize the arrays */
 145     for ( i = 1; i < NUM_POINTS; i++ )
 146     {
 147       memcpy ( red_position[i], red_position[0], 3 * sizeof(double) ) ;
 148       memcpy ( grn_position[i], grn_position[0], 3 * sizeof(double) ) ;
 149     }
 150
 151     break ;
 152
 153   case 'm' :  case 'M' :  /* Modify the Lorenz parameters */
 154     printf ( "Please enter new value for <sigma> (default %f, currently %f): ", s0, sigma ) ;
 155     checkedFGets ( inputline, sizeof ( inputline ), stdin ) ;
 156     sscanf ( inputline, "%lf", &sigma ) ;
 157
 158     printf ( "Please enter new value for <b> (default %f, currently %f): ", b0, b ) ;
 159     checkedFGets ( inputline, sizeof ( inputline ), stdin ) ;
 160     sscanf ( inputline, "%lf", &b ) ;
 161
 162     printf ( "Please enter new value for <r> (default %f, currently %f): ", r0, r ) ;
 163     checkedFGets ( inputline, sizeof ( inputline ), stdin ) ;
 164     sscanf ( inputline, "%lf", &r ) ;
 165
 166     break ;
 167
 168   case 's' :  case 'S' :  /* Stop the animation */
 169     animate = 0 ;
 170     break ;
 171
 172   case 'g' :  case 'G' :  /* Start the animation */
 173     animate = 1 ;
 174     break ;
 175
 176   case ' ' :  /* Spacebar:  Single step */
 177     animate = 2 ;
 178     break ;
 179
 180   case 27 :  /* Escape key */
 181     glutLeaveMainLoop () ;
 182     break ;
 183   }
 184 }
 185
 186 void special_cb ( int key, int x, int y )
 187 {
 188   switch ( key )
 189   {
 190   case GLUT_KEY_UP :  /* Rotate up a little */
 191     glRotated ( ROTATION_ANGLE, 0.0, 1.0, 0.0 ) ;
 192     break ;
 193
 194   case GLUT_KEY_DOWN :  /* Rotate down a little */
 195     glRotated ( -ROTATION_ANGLE, 0.0, 1.0, 0.0 ) ;
 196     break ;
 197
 198   case GLUT_KEY_LEFT :  /* Rotate left a little */
 199     glRotated ( ROTATION_ANGLE, 0.0, 0.0, 1.0 ) ;
 200     break ;
 201
 202   case GLUT_KEY_RIGHT :  /* Rotate right a little */
 203     glRotated ( -ROTATION_ANGLE, 0.0, 0.0, 1.0 ) ;
 204     break ;
 205
 206   case GLUT_KEY_PAGE_UP :  /* Zoom in a little */
 207     glScaled ( 1.0 / SCALE_FACTOR, 1.0 / SCALE_FACTOR, 1.0 / SCALE_FACTOR ) ;
 208     break ;
 209
 210   case GLUT_KEY_PAGE_DOWN :  /* Zoom out a little */
 211     glScaled ( SCALE_FACTOR, SCALE_FACTOR, SCALE_FACTOR ) ;
 212     break ;
 213   }
 214
 215   glutPostRedisplay () ;
 216 }
 217
 218 void mouse_cb ( int button, int updown, int x, int y )
 219 {
 220   if ( updown == GLUT_DOWN )
 221   {
 222     double dist = 1.0e20 ;  /* A very large number */
 223     dist = 0.0 ;  /* so we don't get "unused variable" compiler warning */
 224     /* The idea here is that we go into "pick" mode and pick the nearest point
 225        to the mouse click position.  Unfortunately I don't have the time to implement
 226        it at the moment. */
 227   }
 228 }
 229
 230 void draw_curve ( int index, double position [ NUM_POINTS ][3] )
 231 {
 232   int i = index ;
 233
 234   glBegin ( GL_LINE_STRIP ) ;
 235   do
 236   {
 237     i = ( i == NUM_POINTS-1 ) ? 0 : i + 1 ;
 238     glVertex3dv ( position[i] ) ;
 239   }
 240   while ( i != index ) ;
 241
 242   glEnd () ;
 243 }
 244
 245 void display_cb ( void )
 246 {
 247   char string [ 80 ] ;
 248
 249   glClear ( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT ) ;
 250
 251   glColor3d ( 1.0, 1.0, 1.0 ) ;  /* White */
 252   /* Draw some axes */
 253   glBegin ( GL_LINES ) ;
 254   glVertex3d ( 0.0, 0.0, 0.0 ) ;
 255   glVertex3d ( 2.0, 0.0, 0.0 ) ;
 256   glVertex3d ( 0.0, 0.0, 0.0 ) ;
 257   glVertex3d ( 0.0, 1.0, 0.0 ) ;
 258   glVertex3d ( 0.0, 0.0, 0.0 ) ;
 259   glVertex3d ( 0.0, 0.0, 1.0 ) ;
 260   glEnd () ;
 261
 262   glColor3d ( 1.0, 0.0, 0.0 ) ;  /* Red */
 263   draw_curve ( array_index, red_position ) ;
 264
 265   glColor3d ( 0.0, 1.0, 0.0 ) ;  /* Green */
 266   draw_curve ( array_index, grn_position ) ;
 267
 268   /* Print the distance between the two points */
 269   glColor3d ( 1.0, 1.0, 1.0 ) ;  /* White */
 270   sprintf ( string, "Distance: %10.6f", distance ) ;
 271   glRasterPos2i ( 10, 10 ) ;
 272   glutBitmapString ( GLUT_BITMAP_HELVETICA_12, (unsigned char*)string ) ;
 273
 274   glutSwapBuffers();
 275 }
 276
 277 void reshape_cb ( int width, int height )
 278 {
 279   float ar;
 280   glViewport ( 0, 0, width, height ) ;
 281   glMatrixMode ( GL_PROJECTION ) ;
 282   glLoadIdentity () ;
 283   ar = (float) width / (float) height ;
 284   glFrustum ( -ar, ar, -1.0, 1.0, 10.0, 100.0 ) ;
 285   glMatrixMode ( GL_MODELVIEW ) ;
 286   glLoadIdentity () ;
 287   xcen = 0.0 ;
 288   ycen = 0.0 ;
 289   zcen = 0.0 ;
 290   glTranslated ( xcen, ycen, zcen - 50.0 ) ;
 291 }
 292
 293
 294 void timer_cb ( int value )
 295 {
 296   /* Function called at regular intervals to update the positions of the points */
 297   double deltax, deltay, deltaz ;
 298   int new_index = array_index + 1 ;
 299
 300   /* Set the next timed callback */
 301   glutTimerFunc ( 30, timer_cb, 0 ) ;
 302
 303   if ( animate > 0 )
 304   {
 305     if ( new_index == NUM_POINTS ) new_index = 0 ;
 306     advance_in_time ( time_step, red_position[array_index], red_position[new_index] ) ;
 307     advance_in_time ( time_step, grn_position[array_index], grn_position[new_index] ) ;
 308     array_index = new_index ;
 309
 310     deltax = red_position[array_index][0] - grn_position[array_index][0] ;
 311     deltay = red_position[array_index][1] - grn_position[array_index][1] ;
 312     deltaz = red_position[array_index][2] - grn_position[array_index][2] ;
 313     distance = sqrt ( deltax * deltax + deltay * deltay + deltaz * deltaz ) ;
 314
 315     if ( animate == 2 ) animate = 0 ;
 316   }
 317
 318   glutPostRedisplay () ;
 319 }
 320
 321
 322
 323 /* The Main Program */
 324
 325 int main ( int argc, char *argv[] )
 326 {
 327   int pargc = argc ;
 328
 329   /* Initialize the random number generator */
 330   srand ( 1023 ) ;
 331
 332   /* Set up the OpenGL parameters */
 333   glEnable ( GL_DEPTH_TEST ) ;
 334   glClearColor ( 0.0, 0.0, 0.0, 0.0 ) ;
 335   glClearDepth ( 1.0 ) ;
 336
 337   /* Initialize GLUT */
 338   glutInitWindowSize ( 600, 600 ) ;
 339   glutInit ( &pargc, argv ) ;
 340   glutInitDisplayMode ( GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH ) ;
 341
 342   /* Create the window */
 343   glutCreateWindow ( "Lorenz Attractor" ) ;
 344   glutKeyboardFunc ( key_cb ) ;
 345   glutMouseFunc ( mouse_cb ) ;
 346   glutSpecialFunc ( special_cb ) ;
 347   glutDisplayFunc ( display_cb ) ;
 348   glutReshapeFunc ( reshape_cb ) ;
 349   glutTimerFunc ( 30, timer_cb, 0 ) ;
 350
 351   /* Initialize the attractor:  The easiest way is to call the keyboard callback with an
 352    * argument of 'r' for Reset.
 353    */
 354   key_cb ( 'r', 0, 0 ) ;
 355
 356   /* Enter the GLUT main loop */
 357   glutMainLoop () ;
 358
 359 #ifdef WIN32
 360   /* DUMP MEMORY LEAK INFORMATION */
 361   _CrtDumpMemoryLeaks () ;
 362 #endif
 363
 364   return 0 ;
 365 }
 366