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