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 # include "xlockmore.h" /* in xscreensaver distribution */
106 #else /* STANDALONE */
107 # include "xlock.h" /* in xlockmore distribution */
108 #endif /* STANDALONE */
112 #define DEF_ROTATE "TRUE"
113 #define DEF_RIDE "TRUE"
114 #define DEF_BOX "TRUE"
115 #define DEF_PERIODIC "TRUE"
116 #define DEF_SEARCH "TRUE"
117 #define DEF_DBUF "TRUE"
122 static Bool periodicp;
126 static XrmOptionDescRec opts[] = {
127 {"-rotate", ".flow.rotate", XrmoptionNoArg, "on"},
128 {"+rotate", ".flow.rotate", XrmoptionNoArg, "off"},
129 {"-ride", ".flow.ride", XrmoptionNoArg, "on"},
130 {"+ride", ".flow.ride", XrmoptionNoArg, "off"},
131 {"-box", ".flow.box", XrmoptionNoArg, "on"},
132 {"+box", ".flow.box", XrmoptionNoArg, "off"},
133 {"-periodic", ".flow.periodic", XrmoptionNoArg, "on"},
134 {"+periodic", ".flow.periodic", XrmoptionNoArg, "off"},
135 {"-search", ".flow.search", XrmoptionNoArg, "on"},
136 {"+search", ".flow.search", XrmoptionNoArg, "off"},
137 {"-dbuf", ".flow.dbuf", XrmoptionNoArg, "on"},
138 {"+dbuf", ".flow.dbuf", XrmoptionNoArg, "off"},
141 static argtype vars[] = {
142 {&rotatep, "rotate", "Rotate", DEF_ROTATE, t_Bool},
143 {&ridep, "ride", "Ride", DEF_RIDE, t_Bool},
144 {&boxp, "box", "Box", DEF_BOX, t_Bool},
145 {&periodicp, "periodic", "Periodic", DEF_PERIODIC, t_Bool},
146 {&searchp, "search", "Search", DEF_SEARCH, t_Bool},
147 {&dbufp, "dbuf", "Dbuf", DEF_DBUF, t_Bool},
150 static OptionStruct desc[] = {
151 {"-/+rotate", "turn on/off rotating around attractor."},
152 {"-/+ride", "turn on/off ride in the flow."},
153 {"-/+box", "turn on/off bounding box."},
154 {"-/+periodic", "turn on/off periodic attractors."},
155 {"-/+search", "turn on/off search for new attractors."},
156 {"-/+dbuf", "turn on/off double buffering."},
159 ENTRYPOINT ModeSpecOpt flow_opts =
160 {sizeof opts / sizeof opts[0], opts,
161 sizeof vars / sizeof vars[0], vars, desc};
164 ModStruct flow_description = {
165 "flow", "init_flow", "draw_flow", "release_flow",
166 "refresh_flow", "init_flow", NULL, &flow_opts,
167 1000, 1024, 10000, -10, 200, 1.0, "",
168 "Shows dynamic strange attractors", 0, NULL
173 typedef struct { double x, y, z; } dvector;
175 #define N_PARS 20 /* Enough for Full Cubic or Periodic Cubic */
176 typedef dvector Par[N_PARS];
177 enum { /* Name the parameter indices to make it easier to write
180 X,XX,XXX,XXY,XXZ,XY,XYY,XYZ,XZ,XZZ,
183 SINY = XY /* OK to overlap in this case */
186 /* Camera target [TDA] */
193 #define IX(C) ((C) * segindex + sp->cnsegs[(C)])
194 #define B(t,b) (sp->p + (t) + (b) * sp->taillen)
195 #define X(t,b) (B((t),(b))->x)
196 #define Y(t,b) (B((t),(b))->y)
197 #define Z(t,b) (B((t),(b))->z)
198 #define balance_rand(v) ((LRAND()/MAXRAND*(v))-((v)/2)) /* random around 0 */
199 #define LOST_IN_SPACE 2000.0
200 #define INITIALSTEP 0.04
201 #define EYEHEIGHT 0.005
205 /* Points that make up the box (normalized coordinates) */
206 static const double box[][3] = {
241 /* Lines connecting the box dots */
242 static const double lines[][2] = {
243 {0,1}, {1,2}, {2,3}, {3,0}, /* box */
244 {4,5}, {5,6}, {6,7}, {7,4},
245 {0,4}, {1,5}, {2,6}, {3,7},
246 {4+4,5+4}, {5+4,6+4}, {6+4,7+4}, {7+4,4+4},
247 {4+8,5+8}, {5+8,6+8}, {6+8,7+8}, {7+8,4+8},
248 {4+12,5+12}, {5+12,6+12}, {6+12,7+12}, {7+12,4+12},
249 {4+16,5+16}, {5+16,6+16}, {6+16,7+16}, {7+16,4+16},
250 {4+20,5+20}, {5+20,6+20}, {6+20,7+20}, {7+20,4+20},
251 {4+24,5+24}, {5+24,6+24}, {6+24,7+24}, {7+24,4+24},
255 /* Variables used in rendering */
256 dvector cam[3]; /* camera flight path */
259 Pixmap buffer; /* Double Buffer */
260 dvector circle[2]; /* POV that circles around the scene */
261 dvector centre; /* centre */
262 int beecount; /* number of bees */
263 XSegment *csegs; /* bee lines */
265 XSegment *old_segs; /* old bee lines */
269 /* Variables common to iterators */
270 dvector (*ODE) (Par par, double x, double y, double z);
271 dvector range; /* Initial conditions */
272 double yperiod; /* ODE's where Y is periodic. */
274 /* Variables used in iterating main flow */
276 dvector *p; /* bee positions x[time][bee#] */
280 dvector mid; /* Effective bounding box */
283 /* second set of variables, used for parallel search */
294 static flowstruct *flows = (flowstruct *) NULL;
303 /* Generic 3D Cubic Polynomial. Includes all the Quadratics (Lorentz,
304 Rossler) and much more! */
306 /* I considered offering a seperate 'Quadratic' option, since Cubic is
307 clearly overkill for the standard examples, but the performance
308 difference is too small to measure. The compute time is entirely
309 dominated by the XDrawSegments calls anyway. [TDA] */
311 Cubic(Par a, double x, double y, double z)
314 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 +
315 a[XXZ].x*x*x*z + a[XY].x*x*y + a[XYY].x*x*y*y + a[XYZ].x*x*y*z +
316 a[XZ].x*x*z + a[XZZ].x*x*z*z + a[Y].x*y + a[YY].x*y*y +
317 a[YYY].x*y*y*y + a[YYZ].x*y*y*z + a[YZ].x*y*z + a[YZZ].x*y*z*z +
318 a[Z].x*z + a[ZZ].x*z*z + a[ZZZ].x*z*z*z;
320 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 +
321 a[XXZ].y*x*x*z + a[XY].y*x*y + a[XYY].y*x*y*y + a[XYZ].y*x*y*z +
322 a[XZ].y*x*z + a[XZZ].y*x*z*z + a[Y].y*y + a[YY].y*y*y +
323 a[YYY].y*y*y*y + a[YYZ].y*y*y*z + a[YZ].y*y*z + a[YZZ].y*y*z*z +
324 a[Z].y*z + a[ZZ].y*z*z + a[ZZZ].y*z*z*z;
326 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 +
327 a[XXZ].z*x*x*z + a[XY].z*x*y + a[XYY].z*x*y*y + a[XYZ].z*x*y*z +
328 a[XZ].z*x*z + a[XZZ].z*x*z*z + a[Y].z*y + a[YY].z*y*y +
329 a[YYY].z*y*y*y + a[YYZ].z*y*y*z + a[YZ].z*y*z + a[YZZ].z*y*z*z +
330 a[Z].z*z + a[ZZ].z*z*z + a[ZZZ].z*z*z*z;
335 /* 3D Cubic in (x,z) with periodic sinusoidal forcing term in x. y is
336 the independent periodic (time) axis. This includes Birkhoff's
337 Bagel and Duffing's Attractor */
339 Periodic(Par a, double x, double y, double z)
343 d.x = a[C].x + a[X].x*x + a[XX].x*x*x + a[XXX].x*x*x*x +
344 a[XXZ].x*x*x*z + a[XZ].x*x*z + a[XZZ].x*x*z*z + a[Z].x*z +
345 a[ZZ].x*z*z + a[ZZZ].x*z*z*z + a[SINY].x*sin(y);
349 d.z = a[C].z + a[X].z*x + a[XX].z*x*x + a[XXX].z*x*x*x +
350 a[XXZ].z*x*x*z + a[XZ].z*x*z + a[XZZ].z*x*z*z + a[Z].z*z +
351 a[ZZ].z*z*z + a[ZZZ].z*z*z*z;
356 /* Numerical integration of the ODE using 2nd order Runge Kutta.
357 Returns length^2 of the update, so that we can detect if the step
358 size needs reducing. */
360 Iterate(dvector *p, dvector(*ODE)(Par par, double x, double y, double z),
361 Par par, double step)
365 k1 = ODE(par, p->x, p->y, p->z);
369 k2 = ODE(par, p->x + k1.x, p->y + k1.y, p->z + k1.z);
373 k3.x = (k1.x + k2.x) / 2.0;
374 k3.y = (k1.y + k2.y) / 2.0;
375 k3.z = (k1.z + k2.z) / 2.0;
381 return k3.x*k3.x + k3.y*k3.y + k3.z*k3.z;
384 /* Memory functions */
386 #define deallocate(p,t) if (p!=NULL) {free(p); p=(t*)NULL; }
387 #define allocate(p,t,s) if ((p=(t*)malloc(sizeof(t)*s))==NULL)\
388 {free_flow(sp);return;}
391 free_flow(flowstruct *sp)
393 deallocate(sp->csegs, XSegment);
394 deallocate(sp->cnsegs, int);
395 deallocate(sp->old_segs, XSegment);
396 deallocate(sp->p, dvector);
399 /* Generate Gaussian random number: mean 0, "amplitude" A (actually
400 A is 3*standard deviation). */
402 /* Note this generates a pair of gaussian variables, so it saves one
403 to give out next time it's called */
408 static Bool ready = 0;
415 x = 2.0 * (double)LRAND() / MAXRAND - 1.0;
416 y = 2.0 * (double)LRAND() / MAXRAND - 1.0;
420 w = sqrt((-2 * log(w))/w);
427 /* Attempt to discover new atractors by sending a pair of bees on a
428 fast trip through the new flow and computing their Lyapunov
429 exponent. Returns False if the bees fly away.
431 If the bees stay bounded, the new bounds and the Lyapunov exponent
432 are stored in sp and the function returns True.
434 Repeat invocations continue the flow and improve the accuracy of
435 the bounds and the Lyapunov exponent. Set sp->count2 to zero to
438 Acts on alternate variable set, so that it can be run in parallel
439 with the main flow */
442 discover(ModeInfo * mi)
448 double dl2, df, rs, lsum = 0, s, maxv2 = 0, v2;
454 sp = &flows[MI_SCREEN(mi)];
456 if(sp->count2 == 0) {
457 /* initial conditions */
458 sp->p2[0].x = Gauss_Rand(sp->range.x);
459 sp->p2[0].y = (sp->yperiod > 0)?
460 balance_rand(sp->range.y) : Gauss_Rand(sp->range.y);
461 sp->p2[0].z = Gauss_Rand(sp->range.z);
463 /* 1000 steps to find an attractor */
464 /* Most cases explode out here */
465 for(N=0; N < 1000; N++){
466 Iterate(sp->p2, sp->ODE, sp->par2, sp->step2);
467 if(sp->yperiod > 0 && sp->p2[0].y > sp->yperiod)
468 sp->p2[0].y -= sp->yperiod;
469 if(fabs(sp->p2[0].x) > LOST_IN_SPACE ||
470 fabs(sp->p2[0].y) > LOST_IN_SPACE ||
471 fabs(sp->p2[0].z) > LOST_IN_SPACE) {
476 /* Small perturbation */
477 sp->p2[1].x = sp->p2[0].x + 0.000001;
478 sp->p2[1].y = sp->p2[0].y;
479 sp->p2[1].z = sp->p2[0].z;
482 /* Reset bounding box */
483 max.x = min.x = sp->p2[0].x;
484 max.y = min.y = sp->p2[0].y;
485 max.z = min.z = sp->p2[0].z;
487 /* Compute Lyapunov Exponent */
489 /* (Technically, we're only estimating the largest Lyapunov
490 Exponent, but that's all we need to know to determine if we
491 have a strange attractor.) [TDA] */
493 /* Fly two bees close together */
494 for(N=0; N < 5000; N++){
495 for(i=0; i< 2; i++) {
496 v2 = Iterate(sp->p2+i, sp->ODE, sp->par2, sp->step2);
497 if(sp->yperiod > 0 && sp->p2[i].y > sp->yperiod)
498 sp->p2[i].y -= sp->yperiod;
500 if(fabs(sp->p2[i].x) > LOST_IN_SPACE ||
501 fabs(sp->p2[i].y) > LOST_IN_SPACE ||
502 fabs(sp->p2[i].z) > LOST_IN_SPACE) {
505 if(v2 > maxv2) maxv2 = v2; /* Track max v^2 */
508 /* find bounding box */
509 if ( sp->p2[0].x < min.x ) min.x = sp->p2[0].x;
510 else if ( sp->p2[0].x > max.x ) max.x = sp->p2[0].x;
511 if ( sp->p2[0].y < min.y ) min.y = sp->p2[0].y;
512 else if ( sp->p2[0].y > max.y ) max.y = sp->p2[0].y;
513 if ( sp->p2[0].z < min.z ) min.z = sp->p2[0].z;
514 else if ( sp->p2[0].z > max.z ) max.z = sp->p2[0].z;
516 /* Measure how much we have to pull the two bees to prevent
518 dl.x = sp->p2[1].x - sp->p2[0].x;
519 dl.y = sp->p2[1].y - sp->p2[0].y;
520 dl.z = sp->p2[1].z - sp->p2[0].z;
522 dl2 = dl.x*dl.x + dl.y*dl.y + dl.z*dl.z;
526 sp->p2[1].x = sp->p2[0].x + rs * dl.x;
527 sp->p2[1].y = sp->p2[0].y + rs * dl.y;
528 sp->p2[1].z = sp->p2[0].z + rs * dl.z;
529 lsum = lsum + log(df);
531 l = M_LOG2E / 2 * lsum / nl / sp->step2;
535 /* Anything that didn't explode has a finite attractor */
536 /* If Lyapunov is negative then it probably hit a fixed point or a
537 * limit cycle. Positive Lyapunov indicates a strange attractor. */
541 sp->size2 = max.x - min.x;
543 if(s > sp->size2) sp->size2 = s;
545 if(s > sp->size2) sp->size2 = s;
547 sp->mid2.x = (max.x + min.x) / 2;
548 sp->mid2.y = (max.y + min.y) / 2;
549 sp->mid2.z = (max.z + min.z) / 2;
551 if(sqrt(maxv2) > sp->size2 * 0.2) {
552 /* Flowing too fast, reduce step size. This
553 helps to eliminate high-speed limit cycles,
554 which can show +ve Lyapunov due to integration
561 /* Sets up initial conditions for a flow without all the extra baggage
562 that goes with init_flow */
564 restart_flow(ModeInfo * mi)
571 sp = &flows[MI_SCREEN(mi)];
574 /* Re-Initialize point positions, velocities, etc. */
575 for (b = 0; b < sp->beecount; b++) {
576 X(0, b) = Gauss_Rand(sp->range.x);
577 Y(0, b) = (sp->yperiod > 0)?
578 balance_rand(sp->range.y) : Gauss_Rand(sp->range.y);
579 Z(0, b) = Gauss_Rand(sp->range.z);
583 /* Returns true if line was behind a clip plane, or it clips the line */
584 /* nx,ny,nz is the normal to the plane. d is the distance from the origin */
585 /* s and e are the end points of the line to be clipped */
587 clip(double nx, double ny, double nz, double d, dvector *s, dvector *e)
593 front1 = (nx*s->x + ny*s->y + nz*s->z >= -d);
594 front2 = (nx*e->x + ny*e->y + nz*e->z >= -d);
595 if (!front1 && !front2) return 1;
596 if (front1 && front2) return 0;
601 /* Find t in line equation */
602 t = ( -d - nx*s->x - ny*s->y - nz*s->z) / ( nx*w.x + ny*w.y + nz*w.z);
604 p.x = s->x + w.x * t;
605 p.y = s->y + w.y * t;
606 p.z = s->z + w.z * t;
608 /* Move clipped point to the intersection */
622 init_flow (ModeInfo * mi)
628 if ((flows = (flowstruct *) calloc(MI_NUM_SCREENS(mi),
629 sizeof (flowstruct))) == NULL)
632 sp = &flows[MI_SCREEN(mi)];
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;
645 if(!rotatep && !ridep) rotatep = True; /* We need at least one viewpoint */
647 /* Start camera at Orbit or Bee */
653 sp->chasetime = 1; /* Go directly to target */
657 sp->step2 = INITIALSTEP;
659 /* Zero parameter set */
660 memset(sp->par2, 0, N_PARS * sizeof(dvector));
662 /* Set up standard examples */
663 switch (NRAND((periodicp) ? 5 : 3)) {
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 */
677 sp->par2[Z].z = - 2 + balance_rand(1*0); /* -c */
687 sp->par2[Z].x = -2 + balance_rand(1); /* a */
689 sp->par2[Y].y = 0.2 + balance_rand(0.1); /* b */
690 sp->par2[C].z = 0.2 + balance_rand(0.1); /* c */
692 sp->par2[Z].z = -5.7;
698 z' = 0.2 + z(x - 5.7)
700 name = "RosslerCone";
702 sp->par2[Z].x = -2; /* a */
704 sp->par2[Y].y = 0.2; /* b */
705 sp->par2[ZZ].y = -0.331 + balance_rand(0.01); /* c */
708 sp->par2[Z].z = -5.7;
714 z' = 0.7x + az(0.1 - x^2)
718 sp->par2[SINY].x = 0.35 + balance_rand(0.25); /* b */
719 sp->par2[C].y = 1.57; /* c */
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;
727 x' = -ax - z/2 - z^3/8 + b sin(y)
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 */
738 sp->yperiod = 2 * M_PI;
746 if(sp->yperiod > 0) {
748 /* periodic flows show either uniform distribution or a
749 snapshot on the 'time' axis */
750 sp->range.y = NRAND(2)? sp->yperiod : 0;
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. */
760 if(MI_IS_VERBOSE(mi))
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;
768 sp->step = sp->step2;
769 memcpy(sp->par, sp->par2, sizeof(sp->par2));
771 sp->count2 = 0; /* Reset search */
774 sp->beecount = MI_COUNT(mi);
776 sp->beecount = 1; /* The camera requires 1 or more */
777 } else if (sp->beecount < 0) { /* random variations */
778 sp->beecount = NRAND(-sp->beecount) + 1; /* Minimum 1 */
781 # ifdef HAVE_COCOA /* Don't second-guess Quartz's double-buffering */
785 if(dbufp) { /* Set up double buffer */
786 if (sp->buffer != None)
787 XFreePixmap(MI_DISPLAY(mi), sp->buffer);
788 sp->buffer = XCreatePixmap(MI_DISPLAY(mi), MI_WINDOW(mi),
789 MI_WIDTH(mi), MI_HEIGHT(mi), MI_DEPTH(mi));
791 sp->buffer = MI_WINDOW(mi);
793 /* no "NoExpose" events from XCopyArea wanted */
794 XSetGraphicsExposures(MI_DISPLAY(mi), MI_GC(mi), False);
796 /* Make sure we're using 'thin' lines */
797 XSetLineAttributes(MI_DISPLAY(mi), MI_GC(mi), 0, LineSolid, CapNotLast,
800 /* Clear the background (may be slow depending on user prefs). */
803 /* Allocate memory. */
804 if (sp->csegs == NULL) {
805 allocate(sp->csegs, XSegment,
806 (sp->beecount + BOX_L) * MI_NPIXELS(mi) * sp->taillen);
807 allocate(sp->cnsegs, int, MI_NPIXELS(mi));
808 allocate(sp->old_segs, XSegment, (sp->beecount + BOX_L) * sp->taillen);
809 allocate(sp->p, dvector, sp->beecount * sp->taillen);
812 /* Initialize point positions, velocities, etc. */
815 /* Set up camera tail */
816 X(1, 0) = sp->cam[1].x = 0;
817 Y(1, 0) = sp->cam[1].y = 0;
818 Z(1, 0) = sp->cam[1].z = 0;
822 draw_flow (ModeInfo * mi)
826 double M[3][3]; /* transformation matrix */
827 flowstruct *sp = NULL;
833 sp = &flows[MI_SCREEN(mi)];
834 if (sp->csegs == NULL)
837 #ifdef HAVE_COCOA /* Don't second-guess Quartz's double-buffering */
838 XClearWindow (MI_DISPLAY(mi), MI_WINDOW(mi));
841 /* multiplier for indexing segment arrays. Used in IX macro, etc. */
842 segindex = (sp->beecount + BOX_L) * sp->taillen;
845 if(sp->count2 == 0) { /* start new search */
846 sp->step2 = INITIALSTEP;
847 /* Pick random parameters. Actual range is irrelevant
848 since parameter scale determines flow speed but not
850 for(i=0; i< N_PARS; i++) {
851 sp->par2[i].x = Gauss_Rand(1.0);
852 sp->par2[i].y = Gauss_Rand(1.0);
853 sp->par2[i].z = Gauss_Rand(1.0);
856 if(!discover(mi)) { /* Flow exploded, reset. */
860 sp->count2 = 0; /* Attractor found, but it's not strange */
861 }else if(sp->count2 > 1000000) { /* This one will do */
862 sp->count2 = 0; /* Reset search */
863 if(MI_IS_VERBOSE(mi))
865 "flow: Lyapunov exponent: %g, step: %g, size: %g (unnamed)\n",
866 sp->lyap2, sp->step2, sp->size2);
867 /* Install new params */
868 sp->lyap = sp->lyap2;
869 sp->size = sp->size2;
871 sp->step = sp->step2;
872 memcpy(sp->par, sp->par2, sizeof(sp->par2));
874 /* If we're allowed to zoom out, do so now, so that we
875 get a look at the new attractor. */
876 if(sp->chaseto == BEE && rotatep) {
880 /* Reset initial conditions, so we don't get
881 misleading artifacts in the particle density. */
887 /* Reset segment buffers */
888 for (col = 0; col < MI_NPIXELS(mi); col++)
891 MI_IS_DRAWN(mi) = True;
893 /* Calculate circling POV [Chalky]*/
894 sp->circle[1] = sp->circle[0];
895 sp->circle[0].x = sp->size * 2 * sin(sp->count / 100.0) *
896 (-0.6 + 0.4 *cos(sp->count / 500.0)) + sp->mid.x;
897 sp->circle[0].y = sp->size * 2 * cos(sp->count / 100.0) *
898 (0.6 + 0.4 *cos(sp->count / 500.0)) + sp->mid.y;
899 sp->circle[0].z = sp->size * 2 * sin(sp->count / 421.0) + sp->mid.z;
901 /* Timed chase instead of Chalkie's Bistable oscillator [TDA] */
902 if(rotatep && ridep) {
903 if(sp->chaseto == BEE && NRAND(1000) == 0){
906 }else if(NRAND(4000) == 0){
912 /* Set up orientation matrix */
914 double x[3], p[3], x2=0, xp=0;
917 /* Chasetime is here to guarantee the camera makes it all the
918 way to the target in a finite number of steps. */
919 if(sp->chasetime > 1)
922 if(sp->chaseto == BEE){
923 /* Camera Head targets bee 0 */
924 sp->cam[0].x += (X(0, 0) - sp->cam[0].x)/sp->chasetime;
925 sp->cam[0].y += (Y(0, 0) - sp->cam[0].y)/sp->chasetime;
926 sp->cam[0].z += (Z(0, 0) - sp->cam[0].z)/sp->chasetime;
928 /* Camera Tail targets previous position of bee 0 */
929 sp->cam[1].x += (X(1, 0) - sp->cam[1].x)/sp->chasetime;
930 sp->cam[1].y += (Y(1, 0) - sp->cam[1].y)/sp->chasetime;
931 sp->cam[1].z += (Z(1, 0) - sp->cam[1].z)/sp->chasetime;
933 /* Camera Wing targets bee 1 */
934 sp->cam[2].x += (X(0, 1) - sp->cam[2].x)/sp->chasetime;
935 sp->cam[2].y += (Y(0, 1) - sp->cam[2].y)/sp->chasetime;
936 sp->cam[2].z += (Z(0, 1) - sp->cam[2].z)/sp->chasetime;
938 /* Camera Head targets Orbiter */
939 sp->cam[0].x += (sp->circle[0].x - sp->cam[0].x)/sp->chasetime;
940 sp->cam[0].y += (sp->circle[0].y - sp->cam[0].y)/sp->chasetime;
941 sp->cam[0].z += (sp->circle[0].z - sp->cam[0].z)/sp->chasetime;
943 /* Camera Tail targets diametrically opposite the middle
944 of the bounding box from the Orbiter */
946 (2*sp->circle[0].x - sp->mid.x - sp->cam[1].x)/sp->chasetime;
948 (2*sp->circle[0].y - sp->mid.y - sp->cam[1].y)/sp->chasetime;
950 (2*sp->circle[0].z - sp->mid.z - sp->cam[1].z)/sp->chasetime;
951 /* Camera Wing targets previous position of Orbiter */
952 sp->cam[2].x += (sp->circle[1].x - sp->cam[2].x)/sp->chasetime;
953 sp->cam[2].y += (sp->circle[1].y - sp->cam[2].y)/sp->chasetime;
954 sp->cam[2].z += (sp->circle[1].z - sp->cam[2].z)/sp->chasetime;
957 /* Viewpoint from Tail of camera */
958 sp->centre.x=sp->cam[1].x;
959 sp->centre.y=sp->cam[1].y;
960 sp->centre.z=sp->cam[1].z;
963 x[0] = sp->cam[0].x - sp->cam[1].x;
964 x[1] = sp->cam[0].y - sp->cam[1].y;
965 x[2] = sp->cam[0].z - sp->cam[1].z;
968 p[0] = sp->cam[2].x - sp->cam[1].x;
969 p[1] = sp->cam[2].y - sp->cam[1].y;
970 p[2] = sp->cam[2].z - sp->cam[1].z;
973 /* So long as X and P don't collide, these can be used to form
974 three mutually othogonal axes: X, (X x P) x X and X x P.
975 After being normalised to unit length, these form the
976 Orientation Matrix. */
979 x2+= x[i]*x[i]; /* X . X */
980 xp+= x[i]*p[i]; /* X . P */
981 M[0][i] = x[i]; /* X */
984 for(i=0; i<3; i++) /* (X x P) x X */
985 M[1][i] = x2*p[i] - xp*x[i]; /* == (X . X) P - (X . P) X */
987 M[2][0] = x[1]*p[2] - x[2]*p[1]; /* X x P */
988 M[2][1] = -x[0]*p[2] + x[2]*p[0];
989 M[2][2] = x[0]*p[1] - x[1]*p[0];
994 for(i=0; i<3; i++) A+=M[j][i]*M[j][i]; /* sum squares */
997 for(i=0; i<3; i++) M[j][i]/=A;
1000 if(sp->chaseto == BEE) {
1001 X(0, 1)=X(0, 0)+M[1][0]*sp->step; /* adjust neighbour */
1002 Y(0, 1)=Y(0, 0)+M[1][1]*sp->step;
1003 Z(0, 1)=Z(0, 0)+M[1][2]*sp->step;
1007 /* <=- Bounding Box -=> */
1009 for (b = 0; b < BOX_L; b++) {
1011 /* Chalky's clipping code, Only used for the box */
1012 /* clipping trails is slow and of little benefit. [TDA] */
1013 int p1 = lines[b][0];
1014 int p2 = lines[b][1];
1016 double x1=box[p1][0]* sp->size/2 + sp->mid.x - sp->centre.x;
1017 double y1=box[p1][1]* sp->size/2 + sp->mid.y - sp->centre.y;
1018 double z1=box[p1][2]* sp->size/2 + sp->mid.z - sp->centre.z;
1019 double x2=box[p2][0]* sp->size/2 + sp->mid.x - sp->centre.x;
1020 double y2=box[p2][1]* sp->size/2 + sp->mid.y - sp->centre.y;
1021 double z2=box[p2][2]* sp->size/2 + sp->mid.z - sp->centre.z;
1023 A1.x=M[0][0]*x1 + M[0][1]*y1 + M[0][2]*z1;
1024 A1.y=M[1][0]*x1 + M[1][1]*y1 + M[1][2]*z1;
1025 A1.z=M[2][0]*x1 + M[2][1]*y1 + M[2][2]*z1 + EYEHEIGHT * sp->size;
1026 A2.x=M[0][0]*x2 + M[0][1]*y2 + M[0][2]*z2;
1027 A2.y=M[1][0]*x2 + M[1][1]*y2 + M[1][2]*z2;
1028 A2.z=M[2][0]*x2 + M[2][1]*y2 + M[2][2]*z2 + EYEHEIGHT * sp->size;
1030 /* Clip in 3D before projecting down to 2D. A 2D clip
1031 after projection wouldn't be able to handle lines that
1033 if (clip(1, 0, 0,-1, &A1, &A2) || /* Screen */
1034 clip(1, 2, 0, 0, &A1, &A2) || /* Left */
1035 clip(1,-2, 0, 0, &A1, &A2) || /* Right */
1036 clip(1,0, 2.0*MI_WIDTH(mi)/MI_HEIGHT(mi), 0, &A1, &A2)||/*UP*/
1037 clip(1,0,-2.0*MI_WIDTH(mi)/MI_HEIGHT(mi), 0, &A1, &A2))/*Down*/
1040 /* Colour according to bee */
1041 col = b % (MI_NPIXELS(mi) - 1);
1043 sp->csegs[IX(col)].x1 = MI_WIDTH(mi)/2 + MI_WIDTH(mi) * A1.y/A1.x;
1044 sp->csegs[IX(col)].y1 = MI_HEIGHT(mi)/2 + MI_WIDTH(mi) * A1.z/A1.x;
1045 sp->csegs[IX(col)].x2 = MI_WIDTH(mi)/2 + MI_WIDTH(mi) * A2.y/A2.x;
1046 sp->csegs[IX(col)].y2 = MI_HEIGHT(mi)/2 + MI_WIDTH(mi) * A2.z/A2.x;
1052 for (b = 0; b < sp->beecount; b++) {
1053 if(fabs(X(0, b)) > LOST_IN_SPACE ||
1054 fabs(Y(0, b)) > LOST_IN_SPACE ||
1055 fabs(Z(0, b)) > LOST_IN_SPACE){
1056 if(sp->chaseto == BEE && b == 0){
1057 /* Lost camera bee. Need to replace it since
1058 rerunning init_flow could lose us a hard-won new
1059 attractor. Try moving it very close to a random
1060 other bee. This way we have a good chance of being
1061 close to the attractor and not forming a false
1062 artifact. If we've lost many bees this may need to
1064 /* Don't worry about camera wingbee. It stays close
1065 to the main camera bee no matter what happens. */
1066 int newb = 1 + NRAND(sp->beecount - 1);
1067 X(0, 0) = X(0, newb) + 0.001;
1068 Y(0, 0) = Y(0, newb);
1069 Z(0, 0) = Z(0, newb);
1070 if(MI_IS_VERBOSE(mi))
1072 "flow: resetting lost camera near bee %d\n",
1078 /* Age the tail. It's critical this be fast since
1079 beecount*taillen can be large. */
1080 memmove(B(1, b), B(0, b), (sp->taillen - 1) * sizeof(dvector));
1082 Iterate(B(0,b), sp->ODE, sp->par, sp->step);
1084 /* Don't show wingbee since he's not quite in the flow. */
1085 if(sp->chaseto == BEE && b == 1) continue;
1087 /* Colour according to bee */
1088 col = b % (MI_NPIXELS(mi) - 1);
1090 /* Fill the segment lists. */
1092 begin = 0; /* begin new trail */
1093 end = MIN(sp->taillen, sp->count); /* short trails at first */
1094 for(i=0; i < end; i++){
1095 double x = X(i,b)-sp->centre.x;
1096 double y = Y(i,b)*(sp->yperiod < 0? (sp->size/sp->yperiod) :1)
1098 double z = Z(i,b)-sp->centre.z;
1099 double XM=M[0][0]*x + M[0][1]*y + M[0][2]*z;
1100 double YM=M[1][0]*x + M[1][1]*y + M[1][2]*z;
1101 double ZM=M[2][0]*x + M[2][1]*y + M[2][2]*z + EYEHEIGHT * sp->size;
1104 swarm++; /* count the remaining bees */
1105 if(sp->yperiod > 0 && Y(i,b) > sp->yperiod){
1107 Y(i,b) -= sp->yperiod;
1108 /* hide tail to prevent streaks in Y. Streaks in X,Z
1109 are ok, they help to outline the Poincare'
1111 for(j = i; j < end; j++) Y(j,b) = Y(i,b);
1116 if(XM <= 0){ /* off screen - new trail */
1120 absx = MI_WIDTH(mi)/2 + MI_WIDTH(mi) * YM/XM;
1121 absy = MI_HEIGHT(mi)/2 + MI_WIDTH(mi) * ZM/XM;
1122 /* Performance bottleneck */
1123 if(absx <= 0 || absx >= MI_WIDTH(mi) ||
1124 absy <= 0 || absy >= MI_HEIGHT(mi)) {
1125 /* off screen - new trail */
1129 if(i > begin) { /* complete previous segment */
1130 sp->csegs[IX(col)].x2 = absx;
1131 sp->csegs[IX(col)].y2 = absy;
1135 if(i < end -1){ /* start new segment */
1136 sp->csegs[IX(col)].x1 = absx;
1137 sp->csegs[IX(col)].y1 = absy;
1143 XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_BLACK_PIXEL(mi));
1144 if (dbufp) { /* In Double Buffer case, prepare off-screen copy */
1145 /* For slow systems, this can be the single biggest bottleneck
1146 in the program. These systems may be better of not using
1147 the double buffer. */
1148 XFillRectangle(MI_DISPLAY(mi), sp->buffer, MI_GC(mi), 0, 0,
1149 MI_WIDTH(mi), MI_HEIGHT(mi));
1150 } else { /* Otherwise, erase previous segment list directly */
1151 XDrawSegments(MI_DISPLAY(mi), sp->buffer, MI_GC(mi),
1152 sp->old_segs, sp->nold_segs);
1156 if (MI_NPIXELS(mi) > 2){ /* colour */
1158 for (col = 0; col < MI_NPIXELS(mi) - 1; col++)
1159 if (sp->cnsegs[col] > 0) {
1160 if(sp->cnsegs[col] > mn) mn = sp->cnsegs[col];
1161 XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_PIXEL(mi, col+1));
1162 /* This is usually the biggest bottleneck on most
1163 systems. The maths load is insignificant compared
1165 XDrawSegments(MI_DISPLAY(mi), sp->buffer, MI_GC(mi),
1166 sp->csegs + col * segindex, sp->cnsegs[col]);
1168 } else { /* mono handled seperately since xlockmore uses '1' for
1170 XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_WHITE_PIXEL(mi));
1171 XDrawSegments(MI_DISPLAY(mi), sp->buffer, MI_GC(mi),
1172 sp->csegs, sp->cnsegs[0]);
1174 if (dbufp) { /* In Double Buffer case, this updates the screen */
1175 XCopyArea(MI_DISPLAY(mi), sp->buffer, MI_WINDOW(mi), MI_GC(mi), 0, 0,
1176 MI_WIDTH(mi), MI_HEIGHT(mi), 0, 0);
1177 } else { /* Otherwise, screen is already updated. Copy segments
1178 to erase-list to be erased directly next time. */
1180 for (col = 0; col < MI_NPIXELS(mi) - 1; col++) {
1181 memcpy(sp->old_segs + c, sp->csegs + col * segindex,
1182 sp->cnsegs[col] * sizeof(XSegment));
1183 c += sp->cnsegs[col];
1188 if(sp->count > 1 && swarm == 0) { /* all gone */
1189 if(MI_IS_VERBOSE(mi))
1190 fprintf(stdout, "flow: all gone at %d\n", sp->count);
1194 if(sp->count++ > MI_CYCLES(mi)){ /* Time's up. If we haven't
1195 found anything new by now we
1196 should pick a new standard
1203 reshape_flow(ModeInfo * mi, int width, int height)
1210 release_flow (ModeInfo * mi)
1212 if (flows != NULL) {
1215 for (screen = 0; screen < MI_NUM_SCREENS(mi); screen++)
1216 free_flow(&flows[screen]);
1218 flows = (flowstruct *) NULL;
1223 refresh_flow (ModeInfo * mi)
1225 if(!dbufp) MI_CLEARWINDOW(mi);
1229 flow_handle_event (ModeInfo *mi, XEvent *event)
1231 if (screenhack_event_helper (MI_DISPLAY(mi), MI_WINDOW(mi), event))
1240 XSCREENSAVER_MODULE ("Flow", flow)
1242 #endif /* MODE_flow */