1 /* -*- Mode: C; tab-width: 4; c-basic-offset: 4 -*- */
2 /* flow --- flow of strange bees */
3
4 #if 0
5 static const char sccsid[] = "@(#)flow.c 5.00 2000/11/01 xlockmore";
6
7 #endif
8
9 /*-
10 * Copyright (c) 1996 by Tim Auckland <tda10.geo AT yahoo.com>
11 * Incorporating some code from Stephen Davies Copyright (c) 2000
12 *
13 * Search code based on techniques described in "Strange Attractors:
14 * Creating Patterns in Chaos" by Julien C. Sprott
15 *
16 * Permission to use, copy, modify, and distribute this software and its
17 * documentation for any purpose and without fee is hereby granted,
18 * provided that the above copyright notice appear in all copies and that
19 * both that copyright notice and this permission notice appear in
20 * supporting documentation.
21 *
22 * This file is provided AS IS with no warranties of any kind. The author
23 * shall have no liability with respect to the infringement of copyrights,
24 * trade secrets or any patents by this file or any part thereof. In no
25 * event will the author be liable for any lost revenue or profits or
26 * other special, indirect and consequential damages.
27 *
28 * "flow" shows a variety of continuous phase-space flows around strange
29 * attractors. It includes the well-known Lorentz mask (the "Butterfly"
30 * of chaos fame), two forms of Rossler's "Folded Band" and Poincare'
31 * sections of the "Birkhoff Bagel" and Duffing's forced occilator. "flow"
32 * can now discover new attractors.
33 *
34 * Revision History:
35 *
36 * 29-Oct-2004: [TDA] Discover Attractors unknown to science.
37 * Replace 2D rendering of Periodic Attractors with a 3D
38 * 'interrupted' rendering. Replace "-/+allow2d" with "-/+periodic"
39 * Replace all ODE formulae with completely generic forms.
40 * Add '-search' option to perform background high-speed discovery
41 * for completely new attractors without impacting rendering
42 * performance.
43 * Use gaussian distribution for initial point positions and for
44 * parameter search.
45 * Add "+dbuf" option to allow Double-Buffering to be turned off on
46 * slow X servers.
47 * Remove redundant '-zoom' option. Now automatically zooms if both
48 * rotation and riding are permitted.
49 * Replace dynamic bounding box with static one pre-calculated
50 * during discovery phase.
51 * Simplify and fix bounding box clipping code. Should now be safe
52 * to run without double buffer on all XFree86 servers if desired.
53 * 12-Oct-2004: [TDA] Merge Xscreensaver and Xlockmore branches
54 * Added Chalky's orbital camera, but made the zooming work by
55 * flying the camera rather than interpolating the view transforms.
56 * Added Chalky's Bounding Box, but time-averaged the boundaries to
57 * let the lost bees escape.
58 * Added Chalky's 'view-frustrum' clipping, but only applying it to
59 * the Bounding Box. Trails make clipping less useful.
60 * Added Chalky's "-slow" and "-freeze" options for compatibility,
61 * but haven't implemented the features, since the results are ugly
62 * and make no mathematical contribution.
63 * Added Double-Buffering as a work-around for a persistent XFree86
64 * bug that left debris on the screen.
65 * 21-Mar-2003: [TDA] Trails added (XLockmore branch)
66 * 01-Nov-2000: [TDA] Allocation checks (XLockmore branch)
67 * 21-Feb-2000: [Chalky] Major hackage (Stephen Davies, chalky@null.net)
68 * (Xscreensaver branch)
69 * Forced perspective mode, added 3d box around attractor which
70 * involved coding 3d-planar-clipping against the view-frustrum
71 * thingy. Also made view alternate between piggybacking on a 'bee'
72 * to zooming around outside the attractor. Most bees slow down and
73 * stop, to make the structure of the attractor more obvious.
74 * 28-Jan-1999: [TDA] Catch 'lost' bees in flow.c and disable them.
75 * (XLockmore branch)
76 * I chose to disable them rather than reinitialise them because
77 * reinitialising can produce fake attractors.
78 * This has allowed me to relax some of the parameters and initial
79 * conditions slightly to catch some of the more extreme cases. As a
80 * result you may see some bees fly away at the start - these are the ones
81 * that 'missed' the attractor. If the bee with the camera should fly
82 * away the mode will restart :-)
83 * 31-Nov-1998: [TDA] Added Duffing (what a strange day that was :) DAB)
84 * Duffing's forced oscillator has been added to the formula list and
85 * the parameters section has been updated to display it in Poincare'
86 * section.
87 * 30-Nov-1998: [TDA] Added travelling perspective option
88 * A more exciting point-of-view has been added to all autonomous flows.
89 * This views the flow as seen by a particle moving with the flow. In the
90 * metaphor of the original code, I've attached a camera to one of the
91 * trained bees!
92 * 30-Nov-1998: [TDA] Much code cleanup.
93 * 09-Apr-1997: [TDA] Ported to xlockmore-4
94 * 18-Jul-1996: Adapted from swarm.c Copyright (c) 1991 by Patrick J. Naughton.
95 * 31-Aug-1990: Adapted from xswarm by Jeff Butterworth. (butterwo@ncsc.org).
96 */
97
98 #ifdef STANDALONE
99 # define MODE_flow
100 # define DEFAULTS "*delay: 1000 \n" \
101 "*count: 3000 \n" \
102 "*size: -10 \n" \
103 "*cycles: 10000 \n" \
104 "*ncolors: 200 \n" \
105
106 # define free_flow 0
107 # define reshape_flow 0
108 # define flow_handle_event 0
109 # include "xlockmore.h" /* in xscreensaver distribution */
110 # ifndef MI_DEPTH
111 # define MI_DEPTH MI_WIN_DEPTH
112 # endif
113 #else /* STANDALONE */
114 # include "xlock.h" /* in xlockmore distribution */
115 #endif /* STANDALONE */
116
117 #ifdef MODE_flow
118
119 #define DEF_ROTATE "TRUE"
120 #define DEF_RIDE "TRUE"
121 #define DEF_BOX "TRUE"
122 #define DEF_PERIODIC "TRUE"
123 #define DEF_SEARCH "TRUE"
124 #define DEF_DBUF "TRUE"
125
126 static Bool rotatep;
127 static Bool ridep;
128 static Bool boxp;
129 static Bool periodicp;
130 static Bool searchp;
131 static Bool dbufp;
132
133 static XrmOptionDescRec opts[] = {
134 {(char *) "-rotate", (char *) ".flow.rotate", XrmoptionNoArg, (caddr_t) "on"},
135 {(char *) "+rotate", (char *) ".flow.rotate", XrmoptionNoArg, (caddr_t) "off"},
136 {(char *) "-ride", (char *) ".flow.ride", XrmoptionNoArg, (caddr_t) "on"},
137 {(char *) "+ride", (char *) ".flow.ride", XrmoptionNoArg, (caddr_t) "off"},
138 {(char *) "-box", (char *) ".flow.box", XrmoptionNoArg, (caddr_t) "on"},
139 {(char *) "+box", (char *) ".flow.box", XrmoptionNoArg, (caddr_t) "off"},
140 {(char *) "-periodic", (char *) ".flow.periodic", XrmoptionNoArg, (caddr_t) "on"},
141 {(char *) "+periodic", (char *) ".flow.periodic", XrmoptionNoArg, (caddr_t) "off"},
142 {(char *) "-search", (char *) ".flow.search", XrmoptionNoArg, (caddr_t) "on"},
143 {(char *) "+search", (char *) ".flow.search", XrmoptionNoArg, (caddr_t) "off"},
144 {(char *) "-dbuf", (char *) ".flow.dbuf", XrmoptionNoArg, (caddr_t) "on"},
145 {(char *) "+dbuf", (char *) ".flow.dbuf", XrmoptionNoArg, (caddr_t) "off"},
146 };
147
148 static argtype vars[] = {
149 {(void *) &rotatep, (char *) "rotate",
150 (char *) "Rotate", (char *) DEF_ROTATE, t_Bool},
151 {(void *) &ridep, (char *) "ride",
152 (char *) "Ride", (char *) DEF_RIDE, t_Bool},
153 {(void *) &boxp, (char *) "box",
154 (char *) "Box", (char *) DEF_BOX, t_Bool},
155 {(void *) &periodicp, (char *) "periodic",
156 (char *) "Periodic", (char *) DEF_PERIODIC, t_Bool},
157 {(void *) &searchp, (char *) "search",
158 (char *) "Search", (char *) DEF_SEARCH, t_Bool},
159 {(void *) &dbufp, (char *) "dbuf",
160 (char *) "Dbuf", (char *) DEF_DBUF, t_Bool},
161 };
162
163 static OptionStruct desc[] = {
164 {(char *) "-/+rotate", (char *) "turn on/off rotating around attractor."},
165 {(char *) "-/+ride", (char *) "turn on/off ride in the flow."},
166 {(char *) "-/+box", (char *) "turn on/off bounding box."},
167 {(char *) "-/+periodic", (char *) "turn on/off periodic attractors."},
168 {(char *) "-/+search", (char *) "turn on/off search for new attractors."},
169 {(char *) "-/+dbuf", (char *) "turn on/off double buffering."},
170 };
171
172 ENTRYPOINT ModeSpecOpt flow_opts =
173 {sizeof opts / sizeof opts[0], opts, sizeof vars / sizeof vars[0], vars, desc};
174
175 #ifdef USE_MODULES
176 ModStruct flow_description = {
177 "flow", "init_flow", "draw_flow", "release_flow",
178 "refresh_flow", "init_flow", (char *) NULL, &flow_opts,
179 1000, 1024, 10000, -10, 200, 1.0, "",
180 "Shows dynamic strange attractors", 0, NULL
181 };
182
183 #endif
184
185 typedef struct { double x, y, z; } dvector;
186
187 #define N_PARS 20 /* Enough for Full Cubic or Periodic Cubic */
188 typedef dvector Par[N_PARS];
189 enum { /* Name the parameter indices to make it easier to write
190 standard examples */
191 C,
192 X,XX,XXX,XXY,XXZ,XY,XYY,XYZ,XZ,XZZ,
193 Y,YY,YYY,YYZ,YZ,YZZ,
194 Z,ZZ,ZZZ,
195 SINY = XY /* OK to overlap in this case */
196 };
197
198 /* Camera target [TDA] */
199 typedef enum {
200 ORBIT = 0,
201 BEE = 1
202 } Chaseto;
203
204 /* Macros */
205 #define IX(C) ((C) * segindex + sp->cnsegs[(C)])
206 #define B(t,b) (sp->p + (t) + (b) * sp->taillen)
207 #define X(t,b) (B((t),(b))->x)
208 #define Y(t,b) (B((t),(b))->y)
209 #define Z(t,b) (B((t),(b))->z)
210 #define balance_rand(v) ((LRAND()/MAXRAND*(v))-((v)/2)) /* random around 0 */
211 #define LOST_IN_SPACE 2000.0
212 #define INITIALSTEP 0.04
213 #define EYEHEIGHT 0.005
214 #define MINTRAIL 2
215 #define BOX_L 36
216
217 /* Points that make up the box (normalized coordinates) */
218 static const double box[][3] = {
219 {1,1,1}, /* 0 */
220 {1,1,-1}, /* 1 */
221 {1,-1,-1}, /* 2 */
222 {1,-1,1}, /* 3 */
223 {-1,1,1}, /* 4 */
224 {-1,1,-1}, /* 5 */
225 {-1,-1,-1},/* 6 */
226 {-1,-1,1}, /* 7 */
227 {1, .8, .8},
228 {1, .8,-.8},
229 {1,-.8,-.8},
230 {1,-.8, .8},
231 { .8,1, .8},
232 { .8,1,-.8},
233 {-.8,1,-.8},
234 {-.8,1, .8},
235 { .8, .8,1},
236 { .8,-.8,1},
237 {-.8,-.8,1},
238 {-.8, .8,1},
239 {-1, .8, .8},
240 {-1, .8,-.8},
241 {-1,-.8,-.8},
242 {-1,-.8, .8},
243 { .8,-1, .8},
244 { .8,-1,-.8},
245 {-.8,-1,-.8},
246 {-.8,-1, .8},
247 { .8, .8,-1},
248 { .8,-.8,-1},
249 {-.8,-.8,-1},
250 {-.8, .8,-1}
251 };
252
253 /* Lines connecting the box dots */
254 static const double lines[][2] = {
255 {0,1}, {1,2}, {2,3}, {3,0}, /* box */
256 {4,5}, {5,6}, {6,7}, {7,4},
257 {0,4}, {1,5}, {2,6}, {3,7},
258 {4+4,5+4}, {5+4,6+4}, {6+4,7+4}, {7+4,4+4},
259 {4+8,5+8}, {5+8,6+8}, {6+8,7+8}, {7+8,4+8},
260 {4+12,5+12}, {5+12,6+12}, {6+12,7+12}, {7+12,4+12},
261 {4+16,5+16}, {5+16,6+16}, {6+16,7+16}, {7+16,4+16},
262 {4+20,5+20}, {5+20,6+20}, {6+20,7+20}, {7+20,4+20},
263 {4+24,5+24}, {5+24,6+24}, {6+24,7+24}, {7+24,4+24},
264 };
265
266 typedef struct {
267 /* Variables used in rendering */
268 dvector cam[3]; /* camera flight path */
269 int chasetime;
270 Chaseto chaseto;
271 Pixmap buffer; /* Double Buffer */
272 dvector circle[2]; /* POV that circles around the scene */
273 dvector centre; /* centre */
274 int beecount; /* number of bees */
275 XSegment *csegs; /* bee lines */
276 int *cnsegs;
277 XSegment *old_segs; /* old bee lines */
278 int nold_segs;
279 int taillen;
280
281 /* Variables common to iterators */
282 dvector (*ODE) (Par par, double x, double y, double z);
283 dvector range; /* Initial conditions */
284 double yperiod; /* ODE's where Y is periodic. */
285
286 /* Variables used in iterating main flow */
287 Par par;
288 dvector *p; /* bee positions x[time][bee#] */
289 int count;
290 double lyap;
291 double size;
292 dvector mid; /* Effective bounding box */
293 double step;
294
295 /* second set of variables, used for parallel search */
296 Par par2;
297 dvector p2[2];
298 int count2;
299 double lyap2;
300 double size2;
301 dvector mid2;
302 double step2;
303
304 } flowstruct;
305
306 static flowstruct *flows = (flowstruct *) NULL;
307
308 /*
309 * Private functions
310 */
311
312
313 /* ODE functions */
314
315 /* Generic 3D Cubic Polynomial. Includes all the Quadratics (Lorentz,
316 Rossler) and much more! */
317
318 /* I considered offering a separate 'Quadratic' option, since Cubic is
319 clearly overkill for the standard examples, but the performance
320 difference is too small to measure. The compute time is entirely
321 dominated by the XDrawSegments calls anyway. [TDA] */
322 static dvector
Cubic(Par a,double x,double y,double z)323 Cubic(Par a, double x, double y, double z)
324 {
325 dvector d;
326 d.x = a[C].x + a[X].x*x + a[XX].x*x*x + a[XXX].x*x*x*x + a[XXY].x*x*x*y +
327 a[XXZ].x*x*x*z + a[XY].x*x*y + a[XYY].x*x*y*y + a[XYZ].x*x*y*z +
328 a[XZ].x*x*z + a[XZZ].x*x*z*z + a[Y].x*y + a[YY].x*y*y +
329 a[YYY].x*y*y*y + a[YYZ].x*y*y*z + a[YZ].x*y*z + a[YZZ].x*y*z*z +
330 a[Z].x*z + a[ZZ].x*z*z + a[ZZZ].x*z*z*z;
331
332 d.y = a[C].y + a[X].y*x + a[XX].y*x*x + a[XXX].y*x*x*x + a[XXY].y*x*x*y +
333 a[XXZ].y*x*x*z + a[XY].y*x*y + a[XYY].y*x*y*y + a[XYZ].y*x*y*z +
334 a[XZ].y*x*z + a[XZZ].y*x*z*z + a[Y].y*y + a[YY].y*y*y +
335 a[YYY].y*y*y*y + a[YYZ].y*y*y*z + a[YZ].y*y*z + a[YZZ].y*y*z*z +
336 a[Z].y*z + a[ZZ].y*z*z + a[ZZZ].y*z*z*z;
337
338 d.z = a[C].z + a[X].z*x + a[XX].z*x*x + a[XXX].z*x*x*x + a[XXY].z*x*x*y +
339 a[XXZ].z*x*x*z + a[XY].z*x*y + a[XYY].z*x*y*y + a[XYZ].z*x*y*z +
340 a[XZ].z*x*z + a[XZZ].z*x*z*z + a[Y].z*y + a[YY].z*y*y +
341 a[YYY].z*y*y*y + a[YYZ].z*y*y*z + a[YZ].z*y*z + a[YZZ].z*y*z*z +
342 a[Z].z*z + a[ZZ].z*z*z + a[ZZZ].z*z*z*z;
343
344 return d;
345 }
346
347 /* 3D Cubic in (x,z) with periodic sinusoidal forcing term in x. y is
348 the independent periodic (time) axis. This includes Birkhoff's
349 Bagel and Duffing's Attractor */
350 static dvector
Periodic(Par a,double x,double y,double z)351 Periodic(Par a, double x, double y, double z)
352 {
353 dvector d;
354
355 d.x = a[C].x + a[X].x*x + a[XX].x*x*x + a[XXX].x*x*x*x +
356 a[XXZ].x*x*x*z + a[XZ].x*x*z + a[XZZ].x*x*z*z + a[Z].x*z +
357 a[ZZ].x*z*z + a[ZZZ].x*z*z*z + a[SINY].x*sin(y);
358
359 d.y = a[C].y;
360
361 d.z = a[C].z + a[X].z*x + a[XX].z*x*x + a[XXX].z*x*x*x +
362 a[XXZ].z*x*x*z + a[XZ].z*x*z + a[XZZ].z*x*z*z + a[Z].z*z +
363 a[ZZ].z*z*z + a[ZZZ].z*z*z*z;
364
365 return d;
366 }
367
368 /* Numerical integration of the ODE using 2nd order Runge Kutta.
369 Returns length^2 of the update, so that we can detect if the step
370 size needs reducing. */
371 static double
Iterate(dvector * p,dvector (* ODE)(Par par,double x,double y,double z),Par par,double step)372 Iterate(dvector *p, dvector(*ODE)(Par par, double x, double y, double z),
373 Par par, double step)
374 {
375 dvector k1, k2, k3;
376
377 k1 = ODE(par, p->x, p->y, p->z);
378 k1.x *= step;
379 k1.y *= step;
380 k1.z *= step;
381 k2 = ODE(par, p->x + k1.x, p->y + k1.y, p->z + k1.z);
382 k2.x *= step;
383 k2.y *= step;
384 k2.z *= step;
385 k3.x = (k1.x + k2.x) / 2.0;
386 k3.y = (k1.y + k2.y) / 2.0;
387 k3.z = (k1.z + k2.z) / 2.0;
388
389 p->x += k3.x;
390 p->y += k3.y;
391 p->z += k3.z;
392
393 return k3.x*k3.x + k3.y*k3.y + k3.z*k3.z;
394 }
395
396 /* Memory functions */
397
398 #define deallocate(p,t) if (p!=NULL) {free(p); p=(t*)NULL; }
399 #define allocate(p,t,s) if ((p=(t*)malloc(sizeof(t)*s))==NULL)\
400 {free_flow_screen(sp);return;}
401
402 static void
free_flow_screen(flowstruct * sp)403 free_flow_screen(flowstruct *sp)
404 {
405 if (sp == NULL) {
406 return;
407 }
408 deallocate(sp->csegs, XSegment);
409 deallocate(sp->cnsegs, int);
410 deallocate(sp->old_segs, XSegment);
411 deallocate(sp->p, dvector);
412 sp = NULL;
413 }
414
415 /* Generate Gaussian random number: mean 0, "amplitude" A (actually
416 A is 3*standard deviation). */
417
418 /* Note this generates a pair of gaussian variables, so it saves one
419 to give out next time it's called */
420 static double
Gauss_Rand(double A)421 Gauss_Rand(double A)
422 {
423 static double d;
424 static Bool ready = 0;
425 if(ready) {
426 ready = 0;
427 return A/3 * d;
428 } else {
429 double x, y, w;
430 do {
431 x = 2.0 * (double)LRAND() / MAXRAND - 1.0;
432 y = 2.0 * (double)LRAND() / MAXRAND - 1.0;
433 w = x*x + y*y;
434 } while(w >= 1.0);
435
436 w = sqrt((-2 * log(w))/w);
437 ready = 1;
438 d = x * w;
439 return A/3 * y * w;
440 }
441 }
442
443 /* Attempt to discover new atractors by sending a pair of bees on a
444 fast trip through the new flow and computing their Lyapunov
445 exponent. Returns False if the bees fly away.
446
447 If the bees stay bounded, the new bounds and the Lyapunov exponent
448 are stored in sp and the function returns True.
449
450 Repeat invocations continue the flow and improve the accuracy of
451 the bounds and the Lyapunov exponent. Set sp->count2 to zero to
452 start a new test.
453
454 Acts on alternate variable set, so that it can be run in parallel
455 with the main flow */
456
457 static Bool
discover(ModeInfo * mi)458 discover(ModeInfo * mi)
459 {
460 flowstruct *sp;
461 double l = 0;
462 dvector dl;
463 dvector max, min;
464 double dl2, df, rs, lsum = 0, s, maxv2 = 0, v2;
465
466 int N, i, nl = 0;
467
468 if (flows == NULL)
469 return 0;
470 sp = &flows[MI_SCREEN(mi)];
471
472 if(sp->count2 == 0) {
473 /* initial conditions */
474 sp->p2[0].x = Gauss_Rand(sp->range.x);
475 sp->p2[0].y = (sp->yperiod > 0)?
476 balance_rand(sp->range.y) : Gauss_Rand(sp->range.y);
477 sp->p2[0].z = Gauss_Rand(sp->range.z);
478
479 /* 1000 steps to find an attractor */
480 /* Most cases explode out here */
481 for(N=0; N < 1000; N++){
482 (void) Iterate(sp->p2, sp->ODE, sp->par2, sp->step2);
483 if(sp->yperiod > 0 && sp->p2[0].y > sp->yperiod)
484 sp->p2[0].y -= sp->yperiod;
485 if(fabs(sp->p2[0].x) > LOST_IN_SPACE ||
486 fabs(sp->p2[0].y) > LOST_IN_SPACE ||
487 fabs(sp->p2[0].z) > LOST_IN_SPACE) {
488 return 0;
489 }
490 sp->count2++;
491 }
492 /* Small perturbation */
493 sp->p2[1].x = sp->p2[0].x + 0.000001;
494 sp->p2[1].y = sp->p2[0].y;
495 sp->p2[1].z = sp->p2[0].z;
496 }
497
498 /* Reset bounding box */
499 max.x = min.x = sp->p2[0].x;
500 max.y = min.y = sp->p2[0].y;
501 max.z = min.z = sp->p2[0].z;
502
503 /* Compute Lyapunov Exponent */
504
505 /* (Technically, we're only estimating the largest Lyapunov
506 Exponent, but that's all we need to know to determine if we
507 have a strange attractor.) [TDA] */
508
509 /* Fly two bees close together */
510 for(N=0; N < 5000; N++){
511 for(i=0; i< 2; i++) {
512 v2 = Iterate(sp->p2+i, sp->ODE, sp->par2, sp->step2);
513 if(sp->yperiod > 0 && sp->p2[i].y > sp->yperiod)
514 sp->p2[i].y -= sp->yperiod;
515
516 if(fabs(sp->p2[i].x) > LOST_IN_SPACE ||
517 fabs(sp->p2[i].y) > LOST_IN_SPACE ||
518 fabs(sp->p2[i].z) > LOST_IN_SPACE) {
519 return 0;
520 }
521 if(v2 > maxv2) maxv2 = v2; /* Track max v^2 */
522 }
523
524 /* find bounding box */
525 if ( sp->p2[0].x < min.x ) min.x = sp->p2[0].x;
526 else if ( sp->p2[0].x > max.x ) max.x = sp->p2[0].x;
527 if ( sp->p2[0].y < min.y ) min.y = sp->p2[0].y;
528 else if ( sp->p2[0].y > max.y ) max.y = sp->p2[0].y;
529 if ( sp->p2[0].z < min.z ) min.z = sp->p2[0].z;
530 else if ( sp->p2[0].z > max.z ) max.z = sp->p2[0].z;
531
532 /* Measure how much we have to pull the two bees to prevent
533 them diverging. */
534 dl.x = sp->p2[1].x - sp->p2[0].x;
535 dl.y = sp->p2[1].y - sp->p2[0].y;
536 dl.z = sp->p2[1].z - sp->p2[0].z;
537
538 dl2 = dl.x*dl.x + dl.y*dl.y + dl.z*dl.z;
539 if(dl2 > 0) {
540 df = 1e12 * dl2;
541 rs = 1/sqrt(df);
542 sp->p2[1].x = sp->p2[0].x + rs * dl.x;
543 sp->p2[1].y = sp->p2[0].y + rs * dl.y;
544 sp->p2[1].z = sp->p2[0].z + rs * dl.z;
545 lsum = lsum + log(df);
546 nl = nl + 1;
547 l = M_LOG2E / 2 * lsum / nl / sp->step2;
548 }
549 sp->count2++;
550 }
551 /* Anything that didn't explode has a finite attractor */
552 /* If Lyapunov is negative then it probably hit a fixed point or a
553 * limit cycle. Positive Lyapunov indicates a strange attractor. */
554
555 sp->lyap2 = l;
556
557 sp->size2 = max.x - min.x;
558 s = max.y - min.y;
559 if(s > sp->size2) sp->size2 = s;
560 s = max.z - min.z;
561 if(s > sp->size2) sp->size2 = s;
562
563 sp->mid2.x = (max.x + min.x) / 2;
564 sp->mid2.y = (max.y + min.y) / 2;
565 sp->mid2.z = (max.z + min.z) / 2;
566
567 if(sqrt(maxv2) > sp->size2 * 0.2) {
568 /* Flowing too fast, reduce step size. This
569 helps to eliminate high-speed limit cycles,
570 which can show +ve Lyapunov due to integration
571 inaccuracy. */
572 sp->step2 /= 2;
573 }
574 return 1;
575 }
576
577 /* Sets up initial conditions for a flow without all the extra baggage
578 that goes with init_flow */
579 static void
restart_flow(ModeInfo * mi)580 restart_flow(ModeInfo * mi)
581 {
582 flowstruct *sp;
583 int b;
584
585 if (flows == NULL)
586 return;
587 sp = &flows[MI_SCREEN(mi)];
588 sp->count = 0;
589
590 /* Re-Initialize point positions, velocities, etc. */
591 for (b = 0; b < sp->beecount; b++) {
592 X(0, b) = Gauss_Rand(sp->range.x);
593 Y(0, b) = (sp->yperiod > 0)?
594 balance_rand(sp->range.y) : Gauss_Rand(sp->range.y);
595 Z(0, b) = Gauss_Rand(sp->range.z);
596 }
597 }
598
599 /* Returns true if line was behind a clip plane, or it clips the line */
600 /* nx,ny,nz is the normal to the plane. d is the distance from the origin */
601 /* s and e are the end points of the line to be clipped */
602 static int
clip(double nx,double ny,double nz,double d,dvector * s,dvector * e)603 clip(double nx, double ny, double nz, double d, dvector *s, dvector *e)
604 {
605 int front1, front2;
606 dvector w, p;
607 double t;
608
609 front1 = (nx*s->x + ny*s->y + nz*s->z >= -d);
610 front2 = (nx*e->x + ny*e->y + nz*e->z >= -d);
611 if (!front1 && !front2) return 1;
612 if (front1 && front2) return 0;
613 w.x = e->x - s->x;
614 w.y = e->y - s->y;
615 w.z = e->z - s->z;
616
617 /* Find t in line equation */
618 t = ( -d - nx*s->x - ny*s->y - nz*s->z) / ( nx*w.x + ny*w.y + nz*w.z);
619
620 p.x = s->x + w.x * t;
621 p.y = s->y + w.y * t;
622 p.z = s->z + w.z * t;
623
624 /* Move clipped point to the intersection */
625 if (front2) {
626 *s = p;
627 } else {
628 *e = p;
629 }
630 return 0;
631 }
632
633 /*
634 * Public functions
635 */
636
637 ENTRYPOINT void
init_flow(ModeInfo * mi)638 init_flow(ModeInfo * mi)
639 {
640 flowstruct *sp;
641 char *name;
642
643 #ifdef WIN32
644 /* This is needed because we don't have resource management
645 working on Windows yet so all the defaults are being
646 ignored. */
647 rotatep = 1;
648 ridep = 1;
649 boxp = 1;
650 periodicp = 1;
651 searchp = 1;
652 dbufp = 1;
653 #endif
654
655
656 MI_INIT(mi, flows);
657 sp = &flows[MI_SCREEN(mi)];
658
659 sp->count2 = 0;
660
661 sp->taillen = MI_SIZE(mi);
662 if (sp->taillen < -MINTRAIL) {
663 /* Change by sqrt so it seems more variable */
664 sp->taillen = NRAND((int)sqrt((double) (-sp->taillen - MINTRAIL + 1)));
665 sp->taillen = sp->taillen * sp->taillen + MINTRAIL;
666 } else if (sp->taillen < MINTRAIL) {
667 sp->taillen = MINTRAIL;
668 }
669
670 if(!rotatep && !ridep) rotatep = True; /* We need at least one viewpoint */
671
672 /* Start camera at Orbit or Bee */
673 if(rotatep) {
674 sp->chaseto = ORBIT;
675 } else {
676 sp->chaseto = BEE;
677 }
678 sp->chasetime = 1; /* Go directly to target */
679
680 sp->lyap = 0;
681 sp->yperiod = 0;
682 sp->step2 = INITIALSTEP;
683
684 /* Zero parameter set */
685 (void) memset(sp->par2, 0, N_PARS * sizeof(dvector));
686
687 /* Set up standard examples */
688 switch (NRAND((periodicp) ? 5 : 3)) {
689 case 0:
690 /*
691 x' = a(y - x)
692 y' = x(b - z) - y
693 z' = xy - cz
694 */
695 name = (char *) "Lorentz";
696 sp->par2[Y].x = 10 + balance_rand(5*0); /* a */
697 sp->par2[X].x = - sp->par2[Y].x; /* -a */
698 sp->par2[X].y = 28 + balance_rand(5*0); /* b */
699 sp->par2[XZ].y = -1;
700 sp->par2[Y].y = -1;
701 sp->par2[XY].z = 1;
702 sp->par2[Z].z = - 2 + balance_rand(1*0); /* -c */
703 break;
704 case 1:
705 /*
706 x' = -(y + az)
707 y' = x + by
708 z' = c + z(x - 5.7)
709 */
710 name = (char *) "Rossler";
711 sp->par2[Y].x = -1;
712 sp->par2[Z].x = -2 + balance_rand(1); /* a */
713 sp->par2[X].y = 1;
714 sp->par2[Y].y = 0.2 + balance_rand(0.1); /* b */
715 sp->par2[C].z = 0.2 + balance_rand(0.1); /* c */
716 sp->par2[XZ].z = 1;
717 sp->par2[Z].z = -5.7;
718 break;
719 case 2:
720 /*
721 x' = -(y + az)
722 y' = x + by - cz^2
723 z' = 0.2 + z(x - 5.7)
724 */
725 name = (char *) "RosslerCone";
726 sp->par2[Y].x = -1;
727 sp->par2[Z].x = -2; /* a */
728 sp->par2[X].y = 1;
729 sp->par2[Y].y = 0.2; /* b */
730 sp->par2[ZZ].y = -0.331 + balance_rand(0.01); /* c */
731 sp->par2[C].z = 0.2;
732 sp->par2[XZ].z = 1;
733 sp->par2[Z].z = -5.7;
734 break;
735 case 3:
736 /*
737 x' = -z + b sin(y)
738 y' = c
739 z' = 0.7x + az(0.1 - x^2)
740 */
741 name = (char *) "Birkhoff";
742 sp->par2[Z].x = -1;
743 sp->par2[SINY].x = 0.35 + balance_rand(0.25); /* b */
744 sp->par2[C].y = 1.57; /* c */
745 sp->par2[X].z = 0.7;
746 sp->par2[Z].z = 1 + balance_rand(0.5); /* a/10 */
747 sp->par2[XXZ].z = -10 * sp->par2[Z].z; /* -a */
748 sp->yperiod = 2 * M_PI;
749 break;
750 default:
751 /*
752 x' = -ax - z/2 - z^3/8 + b sin(y)
753 y' = c
754 z' = 2x
755 */
756 name = (char *) "Duffing";
757 sp->par2[X].x = -0.2 + balance_rand(0.1); /* a */
758 sp->par2[Z].x = -0.5;
759 sp->par2[ZZZ].x = -0.125;
760 sp->par2[SINY].x = 27.0 + balance_rand(3.0); /* b */
761 sp->par2[C].y = 1.33; /* c */
762 sp->par2[X].z = 2;
763 sp->yperiod = 2 * M_PI;
764 break;
765
766 }
767
768 sp->range.x = 5;
769 sp->range.z = 5;
770
771 if(sp->yperiod > 0) {
772 sp->ODE = Periodic;
773 /* periodic flows show either uniform distribution or a
774 snapshot on the 'time' axis */
775 sp->range.y = NRAND(2)? sp->yperiod : 0;
776 } else {
777 sp->range.y = 5;
778 sp->ODE = Cubic;
779 }
780
781 /* Run discoverer to set up bounding box, etc. Lyapunov will
782 probably be innaccurate, since we're only running it once, but
783 we're using known strange attractors so it should be ok. */
784 (void) discover(mi);
785 if(MI_IS_VERBOSE(mi))
786 (void) fprintf(stdout,
787 "flow: Lyapunov exponent: %g, step: %g, size: %g (%s)\n",
788 sp->lyap2, sp->step2, sp->size2, name);
789 /* Install new params */
790 sp->lyap = sp->lyap2;
791 sp->size = sp->size2;
792 sp->mid = sp->mid2;
793 sp->step = sp->step2;
794 (void) memcpy(sp->par, sp->par2, sizeof(sp->par2));
795
796 sp->count2 = 0; /* Reset search */
797
798 free_flow_screen(sp);
799 sp->beecount = MI_COUNT(mi);
800 if (sp->beecount < 0) { /* random variations */
801 sp->beecount = NRAND(-sp->beecount) + 1; /* Minimum 1 */
802 }
803
804 if(dbufp) { /* Set up double buffer */
805 if (sp->buffer != None)
806 XFreePixmap(MI_DISPLAY(mi), sp->buffer);
807 sp->buffer = XCreatePixmap(MI_DISPLAY(mi), MI_WINDOW(mi),
808 MI_WIDTH(mi), MI_HEIGHT(mi), MI_DEPTH(mi));
809 } else {
810 sp->buffer = (Pixmap) MI_WINDOW(mi);
811 }
812 /* no "NoExpose" events from XCopyArea wanted */
813 XSetGraphicsExposures(MI_DISPLAY(mi), MI_GC(mi), False);
814
815 /* Make sure we're using 'thin' lines */
816 XSetLineAttributes(MI_DISPLAY(mi), MI_GC(mi), 0, LineSolid, CapNotLast,
817 JoinMiter);
818
819 /* Clear the background (may be slow depending on user prefs). */
820 MI_CLEARWINDOW(mi);
821
822 /* Allocate memory. */
823 if (sp->csegs == NULL) {
824 allocate(sp->csegs, XSegment,
825 (sp->beecount + BOX_L) * MI_NPIXELS(mi) * sp->taillen);
826 allocate(sp->cnsegs, int, MI_NPIXELS(mi));
827 allocate(sp->old_segs, XSegment, sp->beecount * sp->taillen);
828 allocate(sp->p, dvector, sp->beecount * sp->taillen);
829 }
830
831 /* Initialize point positions, velocities, etc. */
832 restart_flow(mi);
833
834 /* Set up camera tail */
835 X(1, 0) = sp->cam[1].x = 0;
836 Y(1, 0) = sp->cam[1].y = 0;
837 Z(1, 0) = sp->cam[1].z = 0;
838 }
839
840 ENTRYPOINT void
draw_flow(ModeInfo * mi)841 draw_flow(ModeInfo * mi)
842 {
843 int b, i;
844 int col, begin, end;
845 double M[3][3]; /* transformation matrix */
846 flowstruct *sp = NULL;
847 int swarm = 0;
848 int segindex;
849
850 if (flows == NULL)
851 return;
852 sp = &flows[MI_SCREEN(mi)];
853 if (sp->csegs == NULL)
854 return;
855
856 /* multiplier for indexing segment arrays. Used in IX macro, etc. */
857 segindex = (sp->beecount + BOX_L) * sp->taillen;
858
859 if(searchp){
860 if(sp->count2 == 0) { /* start new search */
861 sp->step2 = INITIALSTEP;
862 /* Pick random parameters. Actual range is irrelevant
863 since parameter scale determines flow speed but not
864 structure. */
865 for(i=0; i< N_PARS; i++) {
866 sp->par2[i].x = Gauss_Rand(1.0);
867 sp->par2[i].y = Gauss_Rand(1.0);
868 sp->par2[i].z = Gauss_Rand(1.0);
869 }
870 }
871 if(!discover(mi)) { /* Flow exploded, reset. */
872 sp->count2 = 0;
873 } else {
874 if(sp->lyap2 < 0) {
875 sp->count2 = 0; /* Attractor found, but it's not strange */
876 }else if(sp->count2 > 1000000) { /* This one will do */
877 sp->count2 = 0; /* Reset search */
878 if(MI_IS_VERBOSE(mi))
879 (void) fprintf(stdout,
880 "flow: Lyapunov exponent: %g, step: %g, size: %g (unnamed)\n",
881 sp->lyap2, sp->step2, sp->size2);
882 /* Install new params */
883 sp->lyap = sp->lyap2;
884 sp->size = sp->size2;
885 sp->mid = sp->mid2;
886 sp->step = sp->step2;
887 (void) memcpy(sp->par, sp->par2, sizeof(sp->par2));
888
889 /* If we're allowed to zoom out, do so now, so that we
890 get a look at the new attractor. */
891 if(sp->chaseto == BEE && rotatep) {
892 sp->chaseto = ORBIT;
893 sp->chasetime = 100;
894 }
895 /* Reset initial conditions, so we don't get
896 misleading artifacts in the particle density. */
897 restart_flow(mi);
898 }
899 }
900 }
901
902 /* Reset segment buffers */
903 for (col = 0; col < MI_NPIXELS(mi); col++)
904 sp->cnsegs[col] = 0;
905
906 MI_IS_DRAWN(mi) = True;
907
908 /* Calculate circling POV [Chalky]*/
909 sp->circle[1] = sp->circle[0];
910 sp->circle[0].x = sp->size * 2 * sin(sp->count / 100.0) *
911 (-0.6 + 0.4 *cos(sp->count / 500.0)) + sp->mid.x;
912 sp->circle[0].y = sp->size * 2 * cos(sp->count / 100.0) *
913 (0.6 + 0.4 *cos(sp->count / 500.0)) + sp->mid.y;
914 sp->circle[0].z = sp->size * 2 * sin(sp->count / 421.0) + sp->mid.z;
915
916 /* Timed chase instead of Chalkie's Bistable oscillator [TDA] */
917 if(rotatep && ridep) {
918 if(sp->chaseto == BEE && NRAND(1000) == 0){
919 sp->chaseto = ORBIT;
920 sp->chasetime = 100;
921 }else if(NRAND(4000) == 0){
922 sp->chaseto = BEE;
923 sp->chasetime = 100;
924 }
925 }
926
927 /* Set up orientation matrix */
928 {
929 double x[3], p[3], x2=0, xp=0;
930 int j;
931
932 /* Chasetime is here to guarantee the camera makes it all the
933 way to the target in a finite number of steps. */
934 if(sp->chasetime > 1)
935 sp->chasetime--;
936
937 if(sp->chaseto == BEE){
938 /* Camera Head targets bee 0 */
939 sp->cam[0].x += (X(0, 0) - sp->cam[0].x)/sp->chasetime;
940 sp->cam[0].y += (Y(0, 0) - sp->cam[0].y)/sp->chasetime;
941 sp->cam[0].z += (Z(0, 0) - sp->cam[0].z)/sp->chasetime;
942
943 /* Camera Tail targets previous position of bee 0 */
944 sp->cam[1].x += (X(1, 0) - sp->cam[1].x)/sp->chasetime;
945 sp->cam[1].y += (Y(1, 0) - sp->cam[1].y)/sp->chasetime;
946 sp->cam[1].z += (Z(1, 0) - sp->cam[1].z)/sp->chasetime;
947
948 /* Camera Wing targets bee 1 */
949 sp->cam[2].x += (X(0, 1) - sp->cam[2].x)/sp->chasetime;
950 sp->cam[2].y += (Y(0, 1) - sp->cam[2].y)/sp->chasetime;
951 sp->cam[2].z += (Z(0, 1) - sp->cam[2].z)/sp->chasetime;
952 } else {
953 /* Camera Head targets Orbiter */
954 sp->cam[0].x += (sp->circle[0].x - sp->cam[0].x)/sp->chasetime;
955 sp->cam[0].y += (sp->circle[0].y - sp->cam[0].y)/sp->chasetime;
956 sp->cam[0].z += (sp->circle[0].z - sp->cam[0].z)/sp->chasetime;
957
958 /* Camera Tail targets diametrically opposite the middle
959 of the bounding box from the Orbiter */
960 sp->cam[1].x +=
961 (2*sp->circle[0].x - sp->mid.x - sp->cam[1].x)/sp->chasetime;
962 sp->cam[1].y +=
963 (2*sp->circle[0].y - sp->mid.y - sp->cam[1].y)/sp->chasetime;
964 sp->cam[1].z +=
965 (2*sp->circle[0].z - sp->mid.z - sp->cam[1].z)/sp->chasetime;
966 /* Camera Wing targets previous position of Orbiter */
967 sp->cam[2].x += (sp->circle[1].x - sp->cam[2].x)/sp->chasetime;
968 sp->cam[2].y += (sp->circle[1].y - sp->cam[2].y)/sp->chasetime;
969 sp->cam[2].z += (sp->circle[1].z - sp->cam[2].z)/sp->chasetime;
970 }
971
972 /* Viewpoint from Tail of camera */
973 sp->centre.x=sp->cam[1].x;
974 sp->centre.y=sp->cam[1].y;
975 sp->centre.z=sp->cam[1].z;
976
977 /* forward vector */
978 x[0] = sp->cam[0].x - sp->cam[1].x;
979 x[1] = sp->cam[0].y - sp->cam[1].y;
980 x[2] = sp->cam[0].z - sp->cam[1].z;
981
982 /* side */
983 p[0] = sp->cam[2].x - sp->cam[1].x;
984 p[1] = sp->cam[2].y - sp->cam[1].y;
985 p[2] = sp->cam[2].z - sp->cam[1].z;
986
987
988 /* So long as X and P don't collide, these can be used to form
989 three mutually othogonal axes: X, (X x P) x X and X x P.
990 After being normalised to unit length, these form the
991 Orientation Matrix. */
992
993 for(i=0; i<3; i++){
994 x2+= x[i]*x[i]; /* X . X */
995 xp+= x[i]*p[i]; /* X . P */
996 M[0][i] = x[i]; /* X */
997 }
998
999 for(i=0; i<3; i++) /* (X x P) x X */
1000 M[1][i] = x2*p[i] - xp*x[i]; /* == (X . X) P - (X . P) X */
1001
1002 M[2][0] = x[1]*p[2] - x[2]*p[1]; /* X x P */
1003 M[2][1] = -x[0]*p[2] + x[2]*p[0];
1004 M[2][2] = x[0]*p[1] - x[1]*p[0];
1005
1006 /* normalise axes */
1007 for(j=0; j<3; j++){
1008 double A=0;
1009 for(i=0; i<3; i++) A+=M[j][i]*M[j][i]; /* sum squares */
1010 A=sqrt(A);
1011 if(A>0)
1012 for(i=0; i<3; i++) M[j][i]/=A;
1013 }
1014
1015 if(sp->chaseto == BEE) {
1016 X(0, 1)=X(0, 0)+M[1][0]*sp->step; /* adjust neighbour */
1017 Y(0, 1)=Y(0, 0)+M[1][1]*sp->step;
1018 Z(0, 1)=Z(0, 0)+M[1][2]*sp->step;
1019 }
1020 }
1021
1022 /* <=- Bounding Box -=> */
1023 if(boxp) {
1024 for (b = 0; b < BOX_L; b++) {
1025
1026 /* Chalky's clipping code, Only used for the box */
1027 /* clipping trails is slow and of little benefit. [TDA] */
1028 int p1 = (int) lines[b][0];
1029 int p2 = (int) lines[b][1];
1030 dvector A1, A2;
1031 double x1=box[p1][0]* sp->size/2 + sp->mid.x - sp->centre.x;
1032 double y1=box[p1][1]* sp->size/2 + sp->mid.y - sp->centre.y;
1033 double z1=box[p1][2]* sp->size/2 + sp->mid.z - sp->centre.z;
1034 double x2=box[p2][0]* sp->size/2 + sp->mid.x - sp->centre.x;
1035 double y2=box[p2][1]* sp->size/2 + sp->mid.y - sp->centre.y;
1036 double z2=box[p2][2]* sp->size/2 + sp->mid.z - sp->centre.z;
1037
1038 A1.x=M[0][0]*x1 + M[0][1]*y1 + M[0][2]*z1;
1039 A1.y=M[1][0]*x1 + M[1][1]*y1 + M[1][2]*z1;
1040 A1.z=M[2][0]*x1 + M[2][1]*y1 + M[2][2]*z1 + EYEHEIGHT * sp->size;
1041 A2.x=M[0][0]*x2 + M[0][1]*y2 + M[0][2]*z2;
1042 A2.y=M[1][0]*x2 + M[1][1]*y2 + M[1][2]*z2;
1043 A2.z=M[2][0]*x2 + M[2][1]*y2 + M[2][2]*z2 + EYEHEIGHT * sp->size;
1044
1045 /* Clip in 3D before projecting down to 2D. A 2D clip
1046 after projection wouldn't be able to handle lines that
1047 cross x=0 */
1048 if (clip(1.0, 0.0, 0.0,-1.0, &A1, &A2) || /* Screen */
1049 clip(1.0, 2.0, 0.0, 0.0, &A1, &A2) || /* Left */
1050 clip(1.0,-2.0, 0.0, 0.0, &A1, &A2) || /* Right */
1051 clip(1.0,0.0, 2.0*MI_WIDTH(mi)/MI_HEIGHT(mi), 0.0, &A1, &A2)||/*UP*/
1052 clip(1.0,0.0,-2.0*MI_WIDTH(mi)/MI_HEIGHT(mi), 0.0, &A1, &A2))/*Down*/
1053 continue;
1054
1055 /* Colour according to bee */
1056 col = b % (MI_NPIXELS(mi) - 1);
1057
1058 sp->csegs[IX(col)].x1 = (short) (MI_WIDTH(mi)/2 + MI_WIDTH(mi) * A1.y/A1.x);
1059 sp->csegs[IX(col)].y1 = (short) (MI_HEIGHT(mi)/2 + MI_WIDTH(mi) * A1.z/A1.x);
1060 sp->csegs[IX(col)].x2 = (short) (MI_WIDTH(mi)/2 + MI_WIDTH(mi) * A2.y/A2.x);
1061 sp->csegs[IX(col)].y2 = (short) (MI_HEIGHT(mi)/2 + MI_WIDTH(mi) * A2.z/A2.x);
1062 sp->cnsegs[col]++;
1063 }
1064 }
1065
1066 /* <=- Bees -=> */
1067 for (b = 0; b < sp->beecount; b++) {
1068 if(fabs(X(0, b)) > LOST_IN_SPACE ||
1069 fabs(Y(0, b)) > LOST_IN_SPACE ||
1070 fabs(Z(0, b)) > LOST_IN_SPACE){
1071 if(sp->chaseto == BEE && b == 0){
1072 /* Lost camera bee. Need to replace it since
1073 rerunning init_flow could lose us a hard-won new
1074 attractor. Try moving it very close to a random
1075 other bee. This way we have a good chance of being
1076 close to the attractor and not forming a false
1077 artifact. If we've lost many bees this may need to
1078 be repeated. */
1079 /* Don't worry about camera wingbee. It stays close
1080 to the main camera bee no matter what happens. */
1081 int newb = 1 + NRAND(sp->beecount - 1);
1082 X(0, 0) = X(0, newb) + 0.001;
1083 Y(0, 0) = Y(0, newb);
1084 Z(0, 0) = Z(0, newb);
1085 if(MI_IS_VERBOSE(mi))
1086 (void) fprintf(stdout,
1087 "flow: resetting lost camera near bee %d\n",
1088 newb);
1089 }
1090 continue;
1091 }
1092
1093 /* Age the tail. It's critical this be fast since
1094 beecount*taillen can be large. */
1095 (void) memmove(B(1, b), B(0, b), (sp->taillen - 1) * sizeof(dvector));
1096
1097 (void) Iterate(B(0,b), sp->ODE, sp->par, sp->step);
1098
1099 /* Don't show wingbee since he's not quite in the flow. */
1100 if(sp->chaseto == BEE && b == 1) continue;
1101
1102 /* Colour according to bee */
1103 col = b % (MI_NPIXELS(mi) - 1);
1104
1105 /* Fill the segment lists. */
1106
1107 begin = 0; /* begin new trail */
1108 end = MIN(sp->taillen, sp->count); /* short trails at first */
1109 for(i=0; i < end; i++){
1110 double x = X(i,b)-sp->centre.x;
1111 double y = Y(i,b)*(sp->yperiod < 0? (sp->size/sp->yperiod) :1)
1112 -sp->centre.y;
1113 double z = Z(i,b)-sp->centre.z;
1114 double XM=M[0][0]*x + M[0][1]*y + M[0][2]*z;
1115 double YM=M[1][0]*x + M[1][1]*y + M[1][2]*z;
1116 double ZM=M[2][0]*x + M[2][1]*y + M[2][2]*z + EYEHEIGHT * sp->size;
1117 short absx, absy;
1118
1119 swarm++; /* count the remaining bees */
1120 if(sp->yperiod > 0 && Y(i,b) > sp->yperiod){
1121 int j;
1122 Y(i,b) -= sp->yperiod;
1123 /* hide tail to prevent streaks in Y. Streaks in X,Z
1124 are ok, they help to outline the Poincare'
1125 slice. */
1126 for(j = i; j < end; j++) Y(j,b) = Y(i,b);
1127 begin = i + 1;
1128 break;
1129 }
1130
1131 if(XM <= 0){ /* off screen - new trail */
1132 begin = i + 1;
1133 continue;
1134 }
1135 absx = (short) (MI_WIDTH(mi)/2 + MI_WIDTH(mi) * YM/XM);
1136 absy = (short) (MI_HEIGHT(mi)/2 + MI_WIDTH(mi) * ZM/XM);
1137 /* Performance bottleneck */
1138 if(absx <= 0 || absx >= MI_WIDTH(mi) ||
1139 absy <= 0 || absy >= MI_HEIGHT(mi)) {
1140 /* off screen - new trail */
1141 begin = i + 1;
1142 continue;
1143 }
1144 if(i > begin) { /* complete previous segment */
1145 sp->csegs[IX(col)].x2 = absx;
1146 sp->csegs[IX(col)].y2 = absy;
1147 sp->cnsegs[col]++;
1148 }
1149
1150 if(i < end -1){ /* start new segment */
1151 sp->csegs[IX(col)].x1 = absx;
1152 sp->csegs[IX(col)].y1 = absy;
1153 }
1154 }
1155 }
1156
1157 /* Erase */
1158 XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_BLACK_PIXEL(mi));
1159 if (dbufp) { /* In Double Buffer case, prepare off-screen copy */
1160 /* For slow systems, this can be the single biggest bottleneck
1161 in the program. These systems may be better of not using
1162 the double buffer. */
1163 XFillRectangle(MI_DISPLAY(mi), (Drawable) sp->buffer,
1164 MI_GC(mi), 0, 0, MI_WIDTH(mi), MI_HEIGHT(mi));
1165 } else { /* Otherwise, erase previous segment list directly */
1166 XDrawSegments(MI_DISPLAY(mi), (Drawable) sp->buffer,
1167 MI_GC(mi), sp->old_segs, sp->nold_segs);
1168 }
1169
1170 /* Render */
1171 if (MI_NPIXELS(mi) > 2){ /* colour */
1172 int mn = 0;
1173 for (col = 0; col < MI_NPIXELS(mi) - 1; col++)
1174 if (sp->cnsegs[col] > 0) {
1175 if(sp->cnsegs[col] > mn) mn = sp->cnsegs[col];
1176 XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_PIXEL(mi, col+1));
1177 /* This is usually the biggest bottleneck on most
1178 systems. The maths load is insignificant compared
1179 to this. */
1180 XDrawSegments(MI_DISPLAY(mi),
1181 (Drawable) sp->buffer, MI_GC(mi),
1182 sp->csegs + col * segindex, sp->cnsegs[col]);
1183 }
1184 } else { /* mono handled separately since xlockmore uses '1' for
1185 mono white! */
1186 XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_WHITE_PIXEL(mi));
1187 XDrawSegments(MI_DISPLAY(mi), (Drawable) sp->buffer, MI_GC(mi),
1188 sp->csegs, sp->cnsegs[0]);
1189 }
1190 if (dbufp) { /* In Double Buffer case, this updates the screen */
1191 XCopyArea(MI_DISPLAY(mi), (Drawable) sp->buffer,
1192 MI_WINDOW(mi), MI_GC(mi), 0, 0,
1193 MI_WIDTH(mi), MI_HEIGHT(mi), 0, 0);
1194 } else { /* Otherwise, screen is already updated. Copy segments
1195 to erase-list to be erased directly next time. */
1196 int c = 0;
1197 for (col = 0; col < MI_NPIXELS(mi) - 1; col++) {
1198 (void) memcpy(sp->old_segs + c, sp->csegs + col * segindex,
1199 sp->cnsegs[col] * sizeof(XSegment));
1200 c += sp->cnsegs[col];
1201 }
1202 sp->nold_segs = c;
1203 }
1204
1205 if(sp->count > 1 && swarm == 0) { /* all gone */
1206 if(MI_IS_VERBOSE(mi))
1207 (void) fprintf(stdout, "flow: all gone at %d\n", sp->count);
1208 init_flow(mi);
1209 }
1210
1211 if(sp->count++ > MI_CYCLES(mi)){ /* Time's up. If we haven't
1212 found anything new by now we
1213 should pick a new standard
1214 flow */
1215 init_flow(mi);
1216 }
1217 }
1218
1219 ENTRYPOINT void
release_flow(ModeInfo * mi)1220 release_flow(ModeInfo * mi)
1221 {
1222 if (flows != NULL) {
1223 int screen;
1224
1225 for (screen = 0; screen < MI_NUM_SCREENS(mi); screen++)
1226 free_flow_screen(&flows[screen]);
1227 free(flows);
1228 flows = (flowstruct *) NULL;
1229 }
1230 }
1231
1232 #ifndef STANDALONE
1233 ENTRYPOINT void
refresh_flow(ModeInfo * mi)1234 refresh_flow(ModeInfo * mi)
1235 {
1236 if(!dbufp) MI_CLEARWINDOW(mi);
1237 }
1238 #endif
1239
1240 XSCREENSAVER_MODULE ("Flow", flow)
1241
1242 #endif /* MODE_flow */
1243