1 /* -*- Mode: C; tab-width: 4; c-basic-offset: 4 -*- */
2 /* flow --- flow of strange bees */
5 static const char sccsid[] = "@(#)flow.c 5.00 2000/11/01 xlockmore";
9 * Copyright (c) 1996 by Tim Auckland <tda10.geo@yahoo.com>
10 * Incorporating some code from Stephen Davies Copyright (c) 2000
12 * Search code based on techniques described in "Strange Attractors:
13 * Creating Patterns in Chaos" by Julien C. Sprott
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.
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.
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.
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
42 * Use gaussian distribution for initial point positions and for
44 * Add "+dbuf" option to allow Double-Buffering to be turned off on
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.
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'
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
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).
99 # define DEFAULTS "*delay: 10000 \n" \
102 "*cycles: 10000 \n" \
105 # define release_flow 0
106 # include "xlockmore.h" /* in xscreensaver distribution */
107 #else /* STANDALONE */
108 # include "xlock.h" /* in xlockmore distribution */
109 #endif /* STANDALONE */
113 #define DEF_ROTATE "TRUE"
114 #define DEF_RIDE "TRUE"
115 #define DEF_BOX "TRUE"
116 #define DEF_PERIODIC "TRUE"
117 #define DEF_SEARCH "TRUE"
118 #define DEF_DBUF "TRUE"
123 static Bool periodicp;
127 static XrmOptionDescRec opts[] = {
128 {"-rotate", ".flow.rotate", XrmoptionNoArg, "on"},
129 {"+rotate", ".flow.rotate", XrmoptionNoArg, "off"},
130 {"-ride", ".flow.ride", XrmoptionNoArg, "on"},
131 {"+ride", ".flow.ride", XrmoptionNoArg, "off"},
132 {"-box", ".flow.box", XrmoptionNoArg, "on"},
133 {"+box", ".flow.box", XrmoptionNoArg, "off"},
134 {"-periodic", ".flow.periodic", XrmoptionNoArg, "on"},
135 {"+periodic", ".flow.periodic", XrmoptionNoArg, "off"},
136 {"-search", ".flow.search", XrmoptionNoArg, "on"},
137 {"+search", ".flow.search", XrmoptionNoArg, "off"},
138 {"-dbuf", ".flow.dbuf", XrmoptionNoArg, "on"},
139 {"+dbuf", ".flow.dbuf", XrmoptionNoArg, "off"},
142 static argtype vars[] = {
143 {&rotatep, "rotate", "Rotate", DEF_ROTATE, t_Bool},
144 {&ridep, "ride", "Ride", DEF_RIDE, t_Bool},
145 {&boxp, "box", "Box", DEF_BOX, t_Bool},
146 {&periodicp, "periodic", "Periodic", DEF_PERIODIC, t_Bool},
147 {&searchp, "search", "Search", DEF_SEARCH, t_Bool},
148 {&dbufp, "dbuf", "Dbuf", DEF_DBUF, t_Bool},
151 static OptionStruct desc[] = {
152 {"-/+rotate", "turn on/off rotating around attractor."},
153 {"-/+ride", "turn on/off ride in the flow."},
154 {"-/+box", "turn on/off bounding box."},
155 {"-/+periodic", "turn on/off periodic attractors."},
156 {"-/+search", "turn on/off search for new attractors."},
157 {"-/+dbuf", "turn on/off double buffering."},
160 ENTRYPOINT ModeSpecOpt flow_opts =
161 {sizeof opts / sizeof opts[0], opts,
162 sizeof vars / sizeof vars[0], vars, desc};
165 ModStruct flow_description = {
166 "flow", "init_flow", "draw_flow", "release_flow",
167 "refresh_flow", "init_flow", NULL, &flow_opts,
168 1000, 1024, 10000, -10, 200, 1.0, "",
169 "Shows dynamic strange attractors", 0, NULL
174 typedef struct { double x, y, z; } dvector;
176 #define N_PARS 20 /* Enough for Full Cubic or Periodic Cubic */
177 typedef dvector Par[N_PARS];
178 enum { /* Name the parameter indices to make it easier to write
181 X,XX,XXX,XXY,XXZ,XY,XYY,XYZ,XZ,XZZ,
184 SINY = XY /* OK to overlap in this case */
187 /* Camera target [TDA] */
194 #define IX(C) ((C) * segindex + sp->cnsegs[(C)])
195 #define B(t,b) (sp->p + (t) + (b) * sp->taillen)
196 #define X(t,b) (B((t),(b))->x)
197 #define Y(t,b) (B((t),(b))->y)
198 #define Z(t,b) (B((t),(b))->z)
199 #define balance_rand(v) ((LRAND()/MAXRAND*(v))-((v)/2)) /* random around 0 */
200 #define LOST_IN_SPACE 2000.0
201 #define INITIALSTEP 0.04
202 #define EYEHEIGHT 0.005
206 /* Points that make up the box (normalized coordinates) */
207 static const double box[][3] = {
242 /* Lines connecting the box dots */
243 static const double lines[][2] = {
244 {0,1}, {1,2}, {2,3}, {3,0}, /* box */
245 {4,5}, {5,6}, {6,7}, {7,4},
246 {0,4}, {1,5}, {2,6}, {3,7},
247 {4+4,5+4}, {5+4,6+4}, {6+4,7+4}, {7+4,4+4},
248 {4+8,5+8}, {5+8,6+8}, {6+8,7+8}, {7+8,4+8},
249 {4+12,5+12}, {5+12,6+12}, {6+12,7+12}, {7+12,4+12},
250 {4+16,5+16}, {5+16,6+16}, {6+16,7+16}, {7+16,4+16},
251 {4+20,5+20}, {5+20,6+20}, {6+20,7+20}, {7+20,4+20},
252 {4+24,5+24}, {5+24,6+24}, {6+24,7+24}, {7+24,4+24},
256 /* Variables used in rendering */
257 dvector cam[3]; /* camera flight path */
260 Pixmap buffer; /* Double Buffer */
261 dvector circle[2]; /* POV that circles around the scene */
262 dvector centre; /* centre */
263 int beecount; /* number of bees */
264 XSegment *csegs; /* bee lines */
266 XSegment *old_segs; /* old bee lines */
270 /* Variables common to iterators */
271 dvector (*ODE) (Par par, double x, double y, double z);
272 dvector range; /* Initial conditions */
273 double yperiod; /* ODE's where Y is periodic. */
275 /* Variables used in iterating main flow */
277 dvector *p; /* bee positions x[time][bee#] */
281 dvector mid; /* Effective bounding box */
284 /* second set of variables, used for parallel search */
295 static flowstruct *flows = (flowstruct *) NULL;
304 /* Generic 3D Cubic Polynomial. Includes all the Quadratics (Lorentz,
305 Rossler) and much more! */
307 /* I considered offering a seperate 'Quadratic' option, since Cubic is
308 clearly overkill for the standard examples, but the performance
309 difference is too small to measure. The compute time is entirely
310 dominated by the XDrawSegments calls anyway. [TDA] */
312 Cubic(Par a, double x, double y, double z)
315 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 +
316 a[XXZ].x*x*x*z + a[XY].x*x*y + a[XYY].x*x*y*y + a[XYZ].x*x*y*z +
317 a[XZ].x*x*z + a[XZZ].x*x*z*z + a[Y].x*y + a[YY].x*y*y +
318 a[YYY].x*y*y*y + a[YYZ].x*y*y*z + a[YZ].x*y*z + a[YZZ].x*y*z*z +
319 a[Z].x*z + a[ZZ].x*z*z + a[ZZZ].x*z*z*z;
321 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 +
322 a[XXZ].y*x*x*z + a[XY].y*x*y + a[XYY].y*x*y*y + a[XYZ].y*x*y*z +
323 a[XZ].y*x*z + a[XZZ].y*x*z*z + a[Y].y*y + a[YY].y*y*y +
324 a[YYY].y*y*y*y + a[YYZ].y*y*y*z + a[YZ].y*y*z + a[YZZ].y*y*z*z +
325 a[Z].y*z + a[ZZ].y*z*z + a[ZZZ].y*z*z*z;
327 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 +
328 a[XXZ].z*x*x*z + a[XY].z*x*y + a[XYY].z*x*y*y + a[XYZ].z*x*y*z +
329 a[XZ].z*x*z + a[XZZ].z*x*z*z + a[Y].z*y + a[YY].z*y*y +
330 a[YYY].z*y*y*y + a[YYZ].z*y*y*z + a[YZ].z*y*z + a[YZZ].z*y*z*z +
331 a[Z].z*z + a[ZZ].z*z*z + a[ZZZ].z*z*z*z;
336 /* 3D Cubic in (x,z) with periodic sinusoidal forcing term in x. y is
337 the independent periodic (time) axis. This includes Birkhoff's
338 Bagel and Duffing's Attractor */
340 Periodic(Par a, double x, double y, double z)
344 d.x = a[C].x + a[X].x*x + a[XX].x*x*x + a[XXX].x*x*x*x +
345 a[XXZ].x*x*x*z + a[XZ].x*x*z + a[XZZ].x*x*z*z + a[Z].x*z +
346 a[ZZ].x*z*z + a[ZZZ].x*z*z*z + a[SINY].x*sin(y);
350 d.z = a[C].z + a[X].z*x + a[XX].z*x*x + a[XXX].z*x*x*x +
351 a[XXZ].z*x*x*z + a[XZ].z*x*z + a[XZZ].z*x*z*z + a[Z].z*z +
352 a[ZZ].z*z*z + a[ZZZ].z*z*z*z;
357 /* Numerical integration of the ODE using 2nd order Runge Kutta.
358 Returns length^2 of the update, so that we can detect if the step
359 size needs reducing. */
361 Iterate(dvector *p, dvector(*ODE)(Par par, double x, double y, double z),
362 Par par, double step)
366 k1 = ODE(par, p->x, p->y, p->z);
370 k2 = ODE(par, p->x + k1.x, p->y + k1.y, p->z + k1.z);
374 k3.x = (k1.x + k2.x) / 2.0;
375 k3.y = (k1.y + k2.y) / 2.0;
376 k3.z = (k1.z + k2.z) / 2.0;
382 return k3.x*k3.x + k3.y*k3.y + k3.z*k3.z;
385 /* Memory functions */
387 #define deallocate(p,t) if (p!=NULL) {free(p); p=(t*)NULL; }
388 #define allocate(p,t,s) if ((p=(t*)malloc(sizeof(t)*s))==NULL)\
389 {free_flow(mi);return;}
392 free_flow(ModeInfo * mi)
394 flowstruct *sp = &flows[MI_SCREEN(mi)];
395 deallocate(sp->csegs, XSegment);
396 deallocate(sp->cnsegs, int);
397 deallocate(sp->old_segs, XSegment);
398 deallocate(sp->p, dvector);
401 /* Generate Gaussian random number: mean 0, "amplitude" A (actually
402 A is 3*standard deviation). */
404 /* Note this generates a pair of gaussian variables, so it saves one
405 to give out next time it's called */
410 static Bool ready = 0;
417 x = 2.0 * (double)LRAND() / MAXRAND - 1.0;
418 y = 2.0 * (double)LRAND() / MAXRAND - 1.0;
422 w = sqrt((-2 * log(w))/w);
429 /* Attempt to discover new atractors by sending a pair of bees on a
430 fast trip through the new flow and computing their Lyapunov
431 exponent. Returns False if the bees fly away.
433 If the bees stay bounded, the new bounds and the Lyapunov exponent
434 are stored in sp and the function returns True.
436 Repeat invocations continue the flow and improve the accuracy of
437 the bounds and the Lyapunov exponent. Set sp->count2 to zero to
440 Acts on alternate variable set, so that it can be run in parallel
441 with the main flow */
444 discover(ModeInfo * mi)
450 double dl2, df, rs, lsum = 0, s, maxv2 = 0, v2;
456 sp = &flows[MI_SCREEN(mi)];
458 if(sp->count2 == 0) {
459 /* initial conditions */
460 sp->p2[0].x = Gauss_Rand(sp->range.x);
461 sp->p2[0].y = (sp->yperiod > 0)?
462 balance_rand(sp->range.y) : Gauss_Rand(sp->range.y);
463 sp->p2[0].z = Gauss_Rand(sp->range.z);
465 /* 1000 steps to find an attractor */
466 /* Most cases explode out here */
467 for(N=0; N < 1000; N++){
468 Iterate(sp->p2, sp->ODE, sp->par2, sp->step2);
469 if(sp->yperiod > 0 && sp->p2[0].y > sp->yperiod)
470 sp->p2[0].y -= sp->yperiod;
471 if(fabs(sp->p2[0].x) > LOST_IN_SPACE ||
472 fabs(sp->p2[0].y) > LOST_IN_SPACE ||
473 fabs(sp->p2[0].z) > LOST_IN_SPACE) {
478 /* Small perturbation */
479 sp->p2[1].x = sp->p2[0].x + 0.000001;
480 sp->p2[1].y = sp->p2[0].y;
481 sp->p2[1].z = sp->p2[0].z;
484 /* Reset bounding box */
485 max.x = min.x = sp->p2[0].x;
486 max.y = min.y = sp->p2[0].y;
487 max.z = min.z = sp->p2[0].z;
489 /* Compute Lyapunov Exponent */
491 /* (Technically, we're only estimating the largest Lyapunov
492 Exponent, but that's all we need to know to determine if we
493 have a strange attractor.) [TDA] */
495 /* Fly two bees close together */
496 for(N=0; N < 5000; N++){
497 for(i=0; i< 2; i++) {
498 v2 = Iterate(sp->p2+i, sp->ODE, sp->par2, sp->step2);
499 if(sp->yperiod > 0 && sp->p2[i].y > sp->yperiod)
500 sp->p2[i].y -= sp->yperiod;
502 if(fabs(sp->p2[i].x) > LOST_IN_SPACE ||
503 fabs(sp->p2[i].y) > LOST_IN_SPACE ||
504 fabs(sp->p2[i].z) > LOST_IN_SPACE) {
507 if(v2 > maxv2) maxv2 = v2; /* Track max v^2 */
510 /* find bounding box */
511 if ( sp->p2[0].x < min.x ) min.x = sp->p2[0].x;
512 else if ( sp->p2[0].x > max.x ) max.x = sp->p2[0].x;
513 if ( sp->p2[0].y < min.y ) min.y = sp->p2[0].y;
514 else if ( sp->p2[0].y > max.y ) max.y = sp->p2[0].y;
515 if ( sp->p2[0].z < min.z ) min.z = sp->p2[0].z;
516 else if ( sp->p2[0].z > max.z ) max.z = sp->p2[0].z;
518 /* Measure how much we have to pull the two bees to prevent
520 dl.x = sp->p2[1].x - sp->p2[0].x;
521 dl.y = sp->p2[1].y - sp->p2[0].y;
522 dl.z = sp->p2[1].z - sp->p2[0].z;
524 dl2 = dl.x*dl.x + dl.y*dl.y + dl.z*dl.z;
528 sp->p2[1].x = sp->p2[0].x + rs * dl.x;
529 sp->p2[1].y = sp->p2[0].y + rs * dl.y;
530 sp->p2[1].z = sp->p2[0].z + rs * dl.z;
531 lsum = lsum + log(df);
533 l = M_LOG2E / 2 * lsum / nl / sp->step2;
537 /* Anything that didn't explode has a finite attractor */
538 /* If Lyapunov is negative then it probably hit a fixed point or a
539 * limit cycle. Positive Lyapunov indicates a strange attractor. */
543 sp->size2 = max.x - min.x;
545 if(s > sp->size2) sp->size2 = s;
547 if(s > sp->size2) sp->size2 = s;
549 sp->mid2.x = (max.x + min.x) / 2;
550 sp->mid2.y = (max.y + min.y) / 2;
551 sp->mid2.z = (max.z + min.z) / 2;
553 if(sqrt(maxv2) > sp->size2 * 0.2) {
554 /* Flowing too fast, reduce step size. This
555 helps to eliminate high-speed limit cycles,
556 which can show +ve Lyapunov due to integration
563 /* Sets up initial conditions for a flow without all the extra baggage
564 that goes with init_flow */
566 restart_flow(ModeInfo * mi)
573 sp = &flows[MI_SCREEN(mi)];
576 /* Re-Initialize point positions, velocities, etc. */
577 for (b = 0; b < sp->beecount; b++) {
578 X(0, b) = Gauss_Rand(sp->range.x);
579 Y(0, b) = (sp->yperiod > 0)?
580 balance_rand(sp->range.y) : Gauss_Rand(sp->range.y);
581 Z(0, b) = Gauss_Rand(sp->range.z);
585 /* Returns true if line was behind a clip plane, or it clips the line */
586 /* nx,ny,nz is the normal to the plane. d is the distance from the origin */
587 /* s and e are the end points of the line to be clipped */
589 clip(double nx, double ny, double nz, double d, dvector *s, dvector *e)
595 front1 = (nx*s->x + ny*s->y + nz*s->z >= -d);
596 front2 = (nx*e->x + ny*e->y + nz*e->z >= -d);
597 if (!front1 && !front2) return 1;
598 if (front1 && front2) return 0;
603 /* Find t in line equation */
604 t = ( -d - nx*s->x - ny*s->y - nz*s->z) / ( nx*w.x + ny*w.y + nz*w.z);
606 p.x = s->x + w.x * t;
607 p.y = s->y + w.y * t;
608 p.z = s->z + w.z * t;
610 /* Move clipped point to the intersection */
624 init_flow (ModeInfo * mi)
629 MI_INIT (mi, flows, free_flow);
630 sp = &flows[MI_SCREEN(mi)];
634 sp->taillen = MI_SIZE(mi);
635 if (sp->taillen < -MINTRAIL) {
636 /* Change by sqrt so it seems more variable */
637 sp->taillen = NRAND((int)sqrt((double) (-sp->taillen - MINTRAIL + 1)));
638 sp->taillen = sp->taillen * sp->taillen + MINTRAIL;
639 } else if (sp->taillen < MINTRAIL) {
640 sp->taillen = MINTRAIL;
643 if(!rotatep && !ridep) rotatep = True; /* We need at least one viewpoint */
645 /* Start camera at Orbit or Bee */
651 sp->chasetime = 1; /* Go directly to target */
655 sp->step2 = INITIALSTEP;
657 /* Zero parameter set */
658 memset(sp->par2, 0, N_PARS * sizeof(dvector));
660 /* Set up standard examples */
661 switch (NRAND((periodicp) ? 5 : 3)) {
669 sp->par2[Y].x = 10 + balance_rand(5*0); /* a */
670 sp->par2[X].x = - sp->par2[Y].x; /* -a */
671 sp->par2[X].y = 28 + balance_rand(5*0); /* b */
675 sp->par2[Z].z = - 2 + balance_rand(1*0); /* -c */
685 sp->par2[Z].x = -2 + balance_rand(1); /* a */
687 sp->par2[Y].y = 0.2 + balance_rand(0.1); /* b */
688 sp->par2[C].z = 0.2 + balance_rand(0.1); /* c */
690 sp->par2[Z].z = -5.7;
696 z' = 0.2 + z(x - 5.7)
698 name = "RosslerCone";
700 sp->par2[Z].x = -2; /* a */
702 sp->par2[Y].y = 0.2; /* b */
703 sp->par2[ZZ].y = -0.331 + balance_rand(0.01); /* c */
706 sp->par2[Z].z = -5.7;
712 z' = 0.7x + az(0.1 - x^2)
716 sp->par2[SINY].x = 0.35 + balance_rand(0.25); /* b */
717 sp->par2[C].y = 1.57; /* c */
719 sp->par2[Z].z = 1 + balance_rand(0.5); /* a/10 */
720 sp->par2[XXZ].z = -10 * sp->par2[Z].z; /* -a */
721 sp->yperiod = 2 * M_PI;
725 x' = -ax - z/2 - z^3/8 + b sin(y)
730 sp->par2[X].x = -0.2 + balance_rand(0.1); /* a */
731 sp->par2[Z].x = -0.5;
732 sp->par2[ZZZ].x = -0.125;
733 sp->par2[SINY].x = 27.0 + balance_rand(3.0); /* b */
734 sp->par2[C].y = 1.33; /* c */
736 sp->yperiod = 2 * M_PI;
744 if(sp->yperiod > 0) {
746 /* periodic flows show either uniform distribution or a
747 snapshot on the 'time' axis */
748 sp->range.y = NRAND(2)? sp->yperiod : 0;
754 /* Run discoverer to set up bounding box, etc. Lyapunov will
755 probably be innaccurate, since we're only running it once, but
756 we're using known strange attractors so it should be ok. */
758 if(MI_IS_VERBOSE(mi))
760 "flow: Lyapunov exponent: %g, step: %g, size: %g (%s)\n",
761 sp->lyap2, sp->step2, sp->size2, name);
762 /* Install new params */
763 sp->lyap = sp->lyap2;
764 sp->size = sp->size2;
766 sp->step = sp->step2;
767 memcpy(sp->par, sp->par2, sizeof(sp->par2));
769 sp->count2 = 0; /* Reset search */
771 sp->beecount = MI_COUNT(mi);
773 sp->beecount = 1; /* The camera requires 1 or more */
774 } else if (sp->beecount < 0) { /* random variations */
775 sp->beecount = NRAND(-sp->beecount) + 1; /* Minimum 1 */
778 # ifdef HAVE_JWXYZ /* Don't second-guess Quartz's double-buffering */
782 if(dbufp) { /* Set up double buffer */
783 if (sp->buffer != None)
784 XFreePixmap(MI_DISPLAY(mi), sp->buffer);
785 sp->buffer = XCreatePixmap(MI_DISPLAY(mi), MI_WINDOW(mi),
786 MI_WIDTH(mi), MI_HEIGHT(mi), MI_DEPTH(mi));
788 sp->buffer = MI_WINDOW(mi);
790 /* no "NoExpose" events from XCopyArea wanted */
791 XSetGraphicsExposures(MI_DISPLAY(mi), MI_GC(mi), False);
793 /* Make sure we're using 'thin' lines */
794 XSetLineAttributes(MI_DISPLAY(mi), MI_GC(mi), 0, LineSolid, CapNotLast,
797 /* Clear the background (may be slow depending on user prefs). */
800 /* Allocate memory. */
801 if (sp->csegs == NULL) {
802 allocate(sp->csegs, XSegment,
803 (sp->beecount + BOX_L) * MI_NPIXELS(mi) * sp->taillen);
804 allocate(sp->cnsegs, int, MI_NPIXELS(mi));
805 allocate(sp->old_segs, XSegment, (sp->beecount + BOX_L) * sp->taillen);
806 allocate(sp->p, dvector, sp->beecount * sp->taillen);
809 /* Initialize point positions, velocities, etc. */
812 /* Set up camera tail */
813 X(1, 0) = sp->cam[1].x = 0;
814 Y(1, 0) = sp->cam[1].y = 0;
815 Z(1, 0) = sp->cam[1].z = 0;
819 draw_flow (ModeInfo * mi)
823 double M[3][3]; /* transformation matrix */
824 flowstruct *sp = NULL;
830 sp = &flows[MI_SCREEN(mi)];
831 if (sp->csegs == NULL)
834 #ifdef HAVE_JWXYZ /* Don't second-guess Quartz's double-buffering */
835 XClearWindow (MI_DISPLAY(mi), MI_WINDOW(mi));
838 /* multiplier for indexing segment arrays. Used in IX macro, etc. */
839 segindex = (sp->beecount + BOX_L) * sp->taillen;
842 if(sp->count2 == 0) { /* start new search */
843 sp->step2 = INITIALSTEP;
844 /* Pick random parameters. Actual range is irrelevant
845 since parameter scale determines flow speed but not
847 for(i=0; i< N_PARS; i++) {
848 sp->par2[i].x = Gauss_Rand(1.0);
849 sp->par2[i].y = Gauss_Rand(1.0);
850 sp->par2[i].z = Gauss_Rand(1.0);
853 if(!discover(mi)) { /* Flow exploded, reset. */
857 sp->count2 = 0; /* Attractor found, but it's not strange */
858 }else if(sp->count2 > 1000000) { /* This one will do */
859 sp->count2 = 0; /* Reset search */
860 if(MI_IS_VERBOSE(mi))
862 "flow: Lyapunov exponent: %g, step: %g, size: %g (unnamed)\n",
863 sp->lyap2, sp->step2, sp->size2);
864 /* Install new params */
865 sp->lyap = sp->lyap2;
866 sp->size = sp->size2;
868 sp->step = sp->step2;
869 memcpy(sp->par, sp->par2, sizeof(sp->par2));
871 /* If we're allowed to zoom out, do so now, so that we
872 get a look at the new attractor. */
873 if(sp->chaseto == BEE && rotatep) {
877 /* Reset initial conditions, so we don't get
878 misleading artifacts in the particle density. */
884 /* Reset segment buffers */
885 for (col = 0; col < MI_NPIXELS(mi); col++)
888 MI_IS_DRAWN(mi) = True;
890 /* Calculate circling POV [Chalky]*/
891 sp->circle[1] = sp->circle[0];
892 sp->circle[0].x = sp->size * 2 * sin(sp->count / 100.0) *
893 (-0.6 + 0.4 *cos(sp->count / 500.0)) + sp->mid.x;
894 sp->circle[0].y = sp->size * 2 * cos(sp->count / 100.0) *
895 (0.6 + 0.4 *cos(sp->count / 500.0)) + sp->mid.y;
896 sp->circle[0].z = sp->size * 2 * sin(sp->count / 421.0) + sp->mid.z;
898 /* Timed chase instead of Chalkie's Bistable oscillator [TDA] */
899 if(rotatep && ridep) {
900 if(sp->chaseto == BEE && NRAND(1000) == 0){
903 }else if(NRAND(4000) == 0){
909 /* Set up orientation matrix */
911 double x[3], p[3], x2=0, xp=0;
914 /* Chasetime is here to guarantee the camera makes it all the
915 way to the target in a finite number of steps. */
916 if(sp->chasetime > 1)
919 if(sp->chaseto == BEE){
920 /* Camera Head targets bee 0 */
921 sp->cam[0].x += (X(0, 0) - sp->cam[0].x)/sp->chasetime;
922 sp->cam[0].y += (Y(0, 0) - sp->cam[0].y)/sp->chasetime;
923 sp->cam[0].z += (Z(0, 0) - sp->cam[0].z)/sp->chasetime;
925 /* Camera Tail targets previous position of bee 0 */
926 sp->cam[1].x += (X(1, 0) - sp->cam[1].x)/sp->chasetime;
927 sp->cam[1].y += (Y(1, 0) - sp->cam[1].y)/sp->chasetime;
928 sp->cam[1].z += (Z(1, 0) - sp->cam[1].z)/sp->chasetime;
930 /* Camera Wing targets bee 1 */
931 sp->cam[2].x += (X(0, 1) - sp->cam[2].x)/sp->chasetime;
932 sp->cam[2].y += (Y(0, 1) - sp->cam[2].y)/sp->chasetime;
933 sp->cam[2].z += (Z(0, 1) - sp->cam[2].z)/sp->chasetime;
935 /* Camera Head targets Orbiter */
936 sp->cam[0].x += (sp->circle[0].x - sp->cam[0].x)/sp->chasetime;
937 sp->cam[0].y += (sp->circle[0].y - sp->cam[0].y)/sp->chasetime;
938 sp->cam[0].z += (sp->circle[0].z - sp->cam[0].z)/sp->chasetime;
940 /* Camera Tail targets diametrically opposite the middle
941 of the bounding box from the Orbiter */
943 (2*sp->circle[0].x - sp->mid.x - sp->cam[1].x)/sp->chasetime;
945 (2*sp->circle[0].y - sp->mid.y - sp->cam[1].y)/sp->chasetime;
947 (2*sp->circle[0].z - sp->mid.z - sp->cam[1].z)/sp->chasetime;
948 /* Camera Wing targets previous position of Orbiter */
949 sp->cam[2].x += (sp->circle[1].x - sp->cam[2].x)/sp->chasetime;
950 sp->cam[2].y += (sp->circle[1].y - sp->cam[2].y)/sp->chasetime;
951 sp->cam[2].z += (sp->circle[1].z - sp->cam[2].z)/sp->chasetime;
954 /* Viewpoint from Tail of camera */
955 sp->centre.x=sp->cam[1].x;
956 sp->centre.y=sp->cam[1].y;
957 sp->centre.z=sp->cam[1].z;
960 x[0] = sp->cam[0].x - sp->cam[1].x;
961 x[1] = sp->cam[0].y - sp->cam[1].y;
962 x[2] = sp->cam[0].z - sp->cam[1].z;
965 p[0] = sp->cam[2].x - sp->cam[1].x;
966 p[1] = sp->cam[2].y - sp->cam[1].y;
967 p[2] = sp->cam[2].z - sp->cam[1].z;
970 /* So long as X and P don't collide, these can be used to form
971 three mutually othogonal axes: X, (X x P) x X and X x P.
972 After being normalised to unit length, these form the
973 Orientation Matrix. */
976 x2+= x[i]*x[i]; /* X . X */
977 xp+= x[i]*p[i]; /* X . P */
978 M[0][i] = x[i]; /* X */
981 for(i=0; i<3; i++) /* (X x P) x X */
982 M[1][i] = x2*p[i] - xp*x[i]; /* == (X . X) P - (X . P) X */
984 M[2][0] = x[1]*p[2] - x[2]*p[1]; /* X x P */
985 M[2][1] = -x[0]*p[2] + x[2]*p[0];
986 M[2][2] = x[0]*p[1] - x[1]*p[0];
991 for(i=0; i<3; i++) A+=M[j][i]*M[j][i]; /* sum squares */
994 for(i=0; i<3; i++) M[j][i]/=A;
997 if(sp->chaseto == BEE) {
998 X(0, 1)=X(0, 0)+M[1][0]*sp->step; /* adjust neighbour */
999 Y(0, 1)=Y(0, 0)+M[1][1]*sp->step;
1000 Z(0, 1)=Z(0, 0)+M[1][2]*sp->step;
1004 /* <=- Bounding Box -=> */
1006 for (b = 0; b < BOX_L; b++) {
1008 /* Chalky's clipping code, Only used for the box */
1009 /* clipping trails is slow and of little benefit. [TDA] */
1010 int p1 = lines[b][0];
1011 int p2 = lines[b][1];
1013 double x1=box[p1][0]* sp->size/2 + sp->mid.x - sp->centre.x;
1014 double y1=box[p1][1]* sp->size/2 + sp->mid.y - sp->centre.y;
1015 double z1=box[p1][2]* sp->size/2 + sp->mid.z - sp->centre.z;
1016 double x2=box[p2][0]* sp->size/2 + sp->mid.x - sp->centre.x;
1017 double y2=box[p2][1]* sp->size/2 + sp->mid.y - sp->centre.y;
1018 double z2=box[p2][2]* sp->size/2 + sp->mid.z - sp->centre.z;
1020 A1.x=M[0][0]*x1 + M[0][1]*y1 + M[0][2]*z1;
1021 A1.y=M[1][0]*x1 + M[1][1]*y1 + M[1][2]*z1;
1022 A1.z=M[2][0]*x1 + M[2][1]*y1 + M[2][2]*z1 + EYEHEIGHT * sp->size;
1023 A2.x=M[0][0]*x2 + M[0][1]*y2 + M[0][2]*z2;
1024 A2.y=M[1][0]*x2 + M[1][1]*y2 + M[1][2]*z2;
1025 A2.z=M[2][0]*x2 + M[2][1]*y2 + M[2][2]*z2 + EYEHEIGHT * sp->size;
1027 /* Clip in 3D before projecting down to 2D. A 2D clip
1028 after projection wouldn't be able to handle lines that
1030 if (clip(1, 0, 0,-1, &A1, &A2) || /* Screen */
1031 clip(1, 2, 0, 0, &A1, &A2) || /* Left */
1032 clip(1,-2, 0, 0, &A1, &A2) || /* Right */
1033 clip(1,0, 2.0*MI_WIDTH(mi)/MI_HEIGHT(mi), 0, &A1, &A2)||/*UP*/
1034 clip(1,0,-2.0*MI_WIDTH(mi)/MI_HEIGHT(mi), 0, &A1, &A2))/*Down*/
1037 /* Colour according to bee */
1038 col = b % (MI_NPIXELS(mi) - 1);
1040 sp->csegs[IX(col)].x1 = MI_WIDTH(mi)/2 + MI_WIDTH(mi) * A1.y/A1.x;
1041 sp->csegs[IX(col)].y1 = MI_HEIGHT(mi)/2 + MI_WIDTH(mi) * A1.z/A1.x;
1042 sp->csegs[IX(col)].x2 = MI_WIDTH(mi)/2 + MI_WIDTH(mi) * A2.y/A2.x;
1043 sp->csegs[IX(col)].y2 = MI_HEIGHT(mi)/2 + MI_WIDTH(mi) * A2.z/A2.x;
1049 for (b = 0; b < sp->beecount; b++) {
1050 if(fabs(X(0, b)) > LOST_IN_SPACE ||
1051 fabs(Y(0, b)) > LOST_IN_SPACE ||
1052 fabs(Z(0, b)) > LOST_IN_SPACE){
1053 if(sp->chaseto == BEE && b == 0){
1054 /* Lost camera bee. Need to replace it since
1055 rerunning init_flow could lose us a hard-won new
1056 attractor. Try moving it very close to a random
1057 other bee. This way we have a good chance of being
1058 close to the attractor and not forming a false
1059 artifact. If we've lost many bees this may need to
1061 /* Don't worry about camera wingbee. It stays close
1062 to the main camera bee no matter what happens. */
1063 int newb = 1 + NRAND(sp->beecount - 1);
1064 X(0, 0) = X(0, newb) + 0.001;
1065 Y(0, 0) = Y(0, newb);
1066 Z(0, 0) = Z(0, newb);
1067 if(MI_IS_VERBOSE(mi))
1069 "flow: resetting lost camera near bee %d\n",
1075 /* Age the tail. It's critical this be fast since
1076 beecount*taillen can be large. */
1077 memmove(B(1, b), B(0, b), (sp->taillen - 1) * sizeof(dvector));
1079 Iterate(B(0,b), sp->ODE, sp->par, sp->step);
1081 /* Don't show wingbee since he's not quite in the flow. */
1082 if(sp->chaseto == BEE && b == 1) continue;
1084 /* Colour according to bee */
1085 col = b % (MI_NPIXELS(mi) - 1);
1087 /* Fill the segment lists. */
1089 begin = 0; /* begin new trail */
1090 end = MIN(sp->taillen, sp->count); /* short trails at first */
1091 for(i=0; i < end; i++){
1092 double x = X(i,b)-sp->centre.x;
1093 double y = Y(i,b)*(sp->yperiod < 0? (sp->size/sp->yperiod) :1)
1095 double z = Z(i,b)-sp->centre.z;
1096 double XM=M[0][0]*x + M[0][1]*y + M[0][2]*z;
1097 double YM=M[1][0]*x + M[1][1]*y + M[1][2]*z;
1098 double ZM=M[2][0]*x + M[2][1]*y + M[2][2]*z + EYEHEIGHT * sp->size;
1101 swarm++; /* count the remaining bees */
1102 if(sp->yperiod > 0 && Y(i,b) > sp->yperiod){
1104 Y(i,b) -= sp->yperiod;
1105 /* hide tail to prevent streaks in Y. Streaks in X,Z
1106 are ok, they help to outline the Poincare'
1108 for(j = i; j < end; j++) Y(j,b) = Y(i,b);
1113 if(XM <= 0){ /* off screen - new trail */
1117 absx = MI_WIDTH(mi)/2 + MI_WIDTH(mi) * YM/XM;
1118 absy = MI_HEIGHT(mi)/2 + MI_WIDTH(mi) * ZM/XM;
1119 /* Performance bottleneck */
1120 if(absx <= 0 || absx >= MI_WIDTH(mi) ||
1121 absy <= 0 || absy >= MI_HEIGHT(mi)) {
1122 /* off screen - new trail */
1126 if(i > begin) { /* complete previous segment */
1127 sp->csegs[IX(col)].x2 = absx;
1128 sp->csegs[IX(col)].y2 = absy;
1132 if(i < end -1){ /* start new segment */
1133 sp->csegs[IX(col)].x1 = absx;
1134 sp->csegs[IX(col)].y1 = absy;
1140 XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_BLACK_PIXEL(mi));
1141 if (dbufp) { /* In Double Buffer case, prepare off-screen copy */
1142 /* For slow systems, this can be the single biggest bottleneck
1143 in the program. These systems may be better of not using
1144 the double buffer. */
1145 XFillRectangle(MI_DISPLAY(mi), sp->buffer, MI_GC(mi), 0, 0,
1146 MI_WIDTH(mi), MI_HEIGHT(mi));
1147 } else { /* Otherwise, erase previous segment list directly */
1148 XDrawSegments(MI_DISPLAY(mi), sp->buffer, MI_GC(mi),
1149 sp->old_segs, sp->nold_segs);
1153 if (MI_NPIXELS(mi) > 2){ /* colour */
1155 for (col = 0; col < MI_NPIXELS(mi) - 1; col++)
1156 if (sp->cnsegs[col] > 0) {
1157 if(sp->cnsegs[col] > mn) mn = sp->cnsegs[col];
1158 XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_PIXEL(mi, col+1));
1159 /* This is usually the biggest bottleneck on most
1160 systems. The maths load is insignificant compared
1162 XDrawSegments(MI_DISPLAY(mi), sp->buffer, MI_GC(mi),
1163 sp->csegs + col * segindex, sp->cnsegs[col]);
1165 } else { /* mono handled seperately since xlockmore uses '1' for
1167 XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_WHITE_PIXEL(mi));
1168 XDrawSegments(MI_DISPLAY(mi), sp->buffer, MI_GC(mi),
1169 sp->csegs, sp->cnsegs[0]);
1171 if (dbufp) { /* In Double Buffer case, this updates the screen */
1172 XCopyArea(MI_DISPLAY(mi), sp->buffer, MI_WINDOW(mi), MI_GC(mi), 0, 0,
1173 MI_WIDTH(mi), MI_HEIGHT(mi), 0, 0);
1174 } else { /* Otherwise, screen is already updated. Copy segments
1175 to erase-list to be erased directly next time. */
1177 for (col = 0; col < MI_NPIXELS(mi) - 1; col++) {
1178 memcpy(sp->old_segs + c, sp->csegs + col * segindex,
1179 sp->cnsegs[col] * sizeof(XSegment));
1180 c += sp->cnsegs[col];
1185 if(sp->count > 1 && swarm == 0) { /* all gone */
1186 if(MI_IS_VERBOSE(mi))
1187 fprintf(stdout, "flow: all gone at %d\n", sp->count);
1191 if(sp->count++ > MI_CYCLES(mi)){ /* Time's up. If we haven't
1192 found anything new by now we
1193 should pick a new standard
1200 reshape_flow(ModeInfo * mi, int width, int height)
1207 refresh_flow (ModeInfo * mi)
1209 if(!dbufp) MI_CLEARWINDOW(mi);
1213 flow_handle_event (ModeInfo *mi, XEvent *event)
1215 if (screenhack_event_helper (MI_DISPLAY(mi), MI_WINDOW(mi), event))
1224 XSCREENSAVER_MODULE ("Flow", flow)
1226 #endif /* MODE_flow */