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