1 /* -*- Mode: C; tab-width: 4; c-basic-offset: 4 -*- */
2 /* flow --- flow of strange bees */
5 #if !defined( lint ) && !defined( SABER )
6 static const char sccsid[] = "@(#)flow.c 5.00 2000/11/01 xlockmore";
11 * Copyright (c) 1996 by Tim Auckland <tda10.geo@yahoo.com>
12 * Incorporating some code from Stephen Davies Copyright (c) 2000
14 * Search code based on techniques described in "Strange Attractors:
15 * Creating Patterns in Chaos" by Julien C. Sprott
17 * Permission to use, copy, modify, and distribute this software and its
18 * documentation for any purpose and without fee is hereby granted,
19 * provided that the above copyright notice appear in all copies and that
20 * both that copyright notice and this permission notice appear in
21 * supporting documentation.
23 * This file is provided AS IS with no warranties of any kind. The author
24 * shall have no liability with respect to the infringement of copyrights,
25 * trade secrets or any patents by this file or any part thereof. In no
26 * event will the author be liable for any lost revenue or profits or
27 * other special, indirect and consequential damages.
29 * "flow" shows a variety of continuous phase-space flows around strange
30 * attractors. It includes the well-known Lorentz mask (the "Butterfly"
31 * of chaos fame), two forms of Rossler's "Folded Band" and Poincare'
32 * sections of the "Birkhoff Bagel" and Duffing's forced occilator. "flow"
33 * can now discover new attractors.
37 * 29-Oct-2004: [TDA] Discover Attractors unknown to science.
38 * Replace 2D rendering of Periodic Attractors with a 3D
39 * 'interrupted' rendering. Replace "-/+allow2d" with "-/+periodic"
40 * Replace all ODE formulae with completely generic forms.
41 * Add '-search' option to perform background high-speed discovery
42 * for completely new attractors without impacting rendering
44 * Use gaussian distribution for initial point positions and for
46 * Add "+dbuf" option to allow Double-Buffering to be turned off on
48 * Remove redundant '-zoom' option. Now automatically zooms if both
49 * rotation and riding are permitted.
50 * Replace dynamic bounding box with static one pre-calculated
51 * during discovery phase.
52 * Simplify and fix bounding box clipping code. Should now be safe
53 * to run without double buffer on all XFree86 servers if desired.
54 * 12-Oct-2004: [TDA] Merge Xscreensaver and Xlockmore branches
55 * Added Chalky's orbital camera, but made the zooming work by
56 * flying the camera rather than interpolating the view transforms.
57 * Added Chalky's Bounding Box, but time-averaged the boundaries to
58 * let the lost bees escape.
59 * Added Chalky's 'view-frustrum' clipping, but only applying it to
60 * the Bounding Box. Trails make clipping less useful.
61 * Added Chalky's "-slow" and "-freeze" options for compatibility,
62 * but haven't implemented the features, since the results are ugly
63 * and make no mathematical contribution.
64 * Added Double-Buffering as a work-around for a persistent XFree86
65 * bug that left debris on the screen.
66 * 21-Mar-2003: [TDA] Trails added (XLockmore branch)
67 * 01-Nov-2000: [TDA] Allocation checks (XLockmore branch)
68 * 21-Feb-2000: [Chalky] Major hackage (Stephen Davies, chalky@null.net)
69 * (Xscreensaver branch)
70 * Forced perspective mode, added 3d box around attractor which
71 * involved coding 3d-planar-clipping against the view-frustrum
72 * thingy. Also made view alternate between piggybacking on a 'bee'
73 * to zooming around outside the attractor. Most bees slow down and
74 * stop, to make the structure of the attractor more obvious.
75 * 28-Jan-1999: [TDA] Catch 'lost' bees in flow.c and disable them.
77 * I chose to disable them rather than reinitialise them because
78 * reinitialising can produce fake attractors.
79 * This has allowed me to relax some of the parameters and initial
80 * conditions slightly to catch some of the more extreme cases. As a
81 * result you may see some bees fly away at the start - these are the ones
82 * that 'missed' the attractor. If the bee with the camera should fly
83 * away the mode will restart :-)
84 * 31-Nov-1998: [TDA] Added Duffing (what a strange day that was :) DAB)
85 * Duffing's forced oscillator has been added to the formula list and
86 * the parameters section has been updated to display it in Poincare'
88 * 30-Nov-1998: [TDA] Added travelling perspective option
89 * A more exciting point-of-view has been added to all autonomous flows.
90 * This views the flow as seen by a particle moving with the flow. In the
91 * metaphor of the original code, I've attached a camera to one of the
93 * 30-Nov-1998: [TDA] Much code cleanup.
94 * 09-Apr-1997: [TDA] Ported to xlockmore-4
95 * 18-Jul-1996: Adapted from swarm.c Copyright (c) 1991 by Patrick J. Naughton.
96 * 31-Aug-1990: Adapted from xswarm by Jeff Butterworth. (butterwo@ncsc.org).
101 # define DEFAULTS "*delay: 10000 \n" \
104 "*cycles: 10000 \n" \
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 */
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"
125 static Bool periodicp;
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"},
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},
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."},
162 ENTRYPOINT ModeSpecOpt flow_opts =
163 {sizeof opts / sizeof opts[0], opts,
164 sizeof vars / sizeof vars[0], vars, desc};
167 ModStruct flow_description = {
168 "flow", "init_flow", "draw_flow", "release_flow",
169 "refresh_flow", "init_flow", NULL, &flow_opts,
170 1000, 1024, 10000, -10, 200, 1.0, "",
171 "Shows dynamic strange attractors", 0, NULL
176 typedef struct { double x, y, z; } dvector;
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
183 X,XX,XXX,XXY,XXZ,XY,XYY,XYZ,XZ,XZZ,
186 SINY = XY /* OK to overlap in this case */
189 /* Camera target [TDA] */
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
208 /* Points that make up the box (normalized coordinates) */
209 static const double box[][3] = {
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},
258 /* Variables used in rendering */
259 dvector cam[3]; /* camera flight path */
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 */
268 XSegment *old_segs; /* old bee lines */
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. */
277 /* Variables used in iterating main flow */
279 dvector *p; /* bee positions x[time][bee#] */
283 dvector mid; /* Effective bounding box */
286 /* second set of variables, used for parallel search */
297 static flowstruct *flows = (flowstruct *) NULL;
306 /* Generic 3D Cubic Polynomial. Includes all the Quadratics (Lorentz,
307 Rossler) and much more! */
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] */
314 Cubic(Par a, double x, double y, double z)
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;
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;
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;
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 */
342 Periodic(Par a, double x, double y, double z)
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);
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;
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. */
363 Iterate(dvector *p, dvector(*ODE)(Par par, double x, double y, double z),
364 Par par, double step)
368 k1 = ODE(par, p->x, p->y, p->z);
372 k2 = ODE(par, p->x + k1.x, p->y + k1.y, p->z + k1.z);
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;
384 return k3.x*k3.x + k3.y*k3.y + k3.z*k3.z;
387 /* Memory functions */
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(sp);return;}
394 free_flow(flowstruct *sp)
396 deallocate(sp->csegs, XSegment);
397 deallocate(sp->cnsegs, int);
398 deallocate(sp->old_segs, XSegment);
399 deallocate(sp->p, dvector);
402 /* Generate Gaussian random number: mean 0, "amplitude" A (actually
403 A is 3*standard deviation). */
405 /* Note this generates a pair of gaussian variables, so it saves one
406 to give out next time it's called */
411 static Bool ready = 0;
418 x = 2.0 * (double)LRAND() / MAXRAND - 1.0;
419 y = 2.0 * (double)LRAND() / MAXRAND - 1.0;
423 w = sqrt((-2 * log(w))/w);
430 /* Attempt to discover new atractors by sending a pair of bees on a
431 fast trip through the new flow and computing their Lyapunov
432 exponent. Returns False if the bees fly away.
434 If the bees stay bounded, the new bounds and the Lyapunov exponent
435 are stored in sp and the function returns True.
437 Repeat invocations continue the flow and improve the accuracy of
438 the bounds and the Lyapunov exponent. Set sp->count2 to zero to
441 Acts on alternate variable set, so that it can be run in parallel
442 with the main flow */
445 discover(ModeInfo * mi)
451 double dl2, df, rs, lsum = 0, s, maxv2 = 0, v2;
457 sp = &flows[MI_SCREEN(mi)];
459 if(sp->count2 == 0) {
460 /* initial conditions */
461 sp->p2[0].x = Gauss_Rand(sp->range.x);
462 sp->p2[0].y = (sp->yperiod > 0)?
463 balance_rand(sp->range.y) : Gauss_Rand(sp->range.y);
464 sp->p2[0].z = Gauss_Rand(sp->range.z);
466 /* 1000 steps to find an attractor */
467 /* Most cases explode out here */
468 for(N=0; N < 1000; N++){
469 Iterate(sp->p2, sp->ODE, sp->par2, sp->step2);
470 if(sp->yperiod > 0 && sp->p2[0].y > sp->yperiod)
471 sp->p2[0].y -= sp->yperiod;
472 if(fabs(sp->p2[0].x) > LOST_IN_SPACE ||
473 fabs(sp->p2[0].y) > LOST_IN_SPACE ||
474 fabs(sp->p2[0].z) > LOST_IN_SPACE) {
479 /* Small perturbation */
480 sp->p2[1].x = sp->p2[0].x + 0.000001;
481 sp->p2[1].y = sp->p2[0].y;
482 sp->p2[1].z = sp->p2[0].z;
485 /* Reset bounding box */
486 max.x = min.x = sp->p2[0].x;
487 max.y = min.y = sp->p2[0].y;
488 max.z = min.z = sp->p2[0].z;
490 /* Compute Lyapunov Exponent */
492 /* (Technically, we're only estimating the largest Lyapunov
493 Exponent, but that's all we need to know to determine if we
494 have a strange attractor.) [TDA] */
496 /* Fly two bees close together */
497 for(N=0; N < 5000; N++){
498 for(i=0; i< 2; i++) {
499 v2 = Iterate(sp->p2+i, sp->ODE, sp->par2, sp->step2);
500 if(sp->yperiod > 0 && sp->p2[i].y > sp->yperiod)
501 sp->p2[i].y -= sp->yperiod;
503 if(fabs(sp->p2[i].x) > LOST_IN_SPACE ||
504 fabs(sp->p2[i].y) > LOST_IN_SPACE ||
505 fabs(sp->p2[i].z) > LOST_IN_SPACE) {
508 if(v2 > maxv2) maxv2 = v2; /* Track max v^2 */
511 /* find bounding box */
512 if ( sp->p2[0].x < min.x ) min.x = sp->p2[0].x;
513 else if ( sp->p2[0].x > max.x ) max.x = sp->p2[0].x;
514 if ( sp->p2[0].y < min.y ) min.y = sp->p2[0].y;
515 else if ( sp->p2[0].y > max.y ) max.y = sp->p2[0].y;
516 if ( sp->p2[0].z < min.z ) min.z = sp->p2[0].z;
517 else if ( sp->p2[0].z > max.z ) max.z = sp->p2[0].z;
519 /* Measure how much we have to pull the two bees to prevent
521 dl.x = sp->p2[1].x - sp->p2[0].x;
522 dl.y = sp->p2[1].y - sp->p2[0].y;
523 dl.z = sp->p2[1].z - sp->p2[0].z;
525 dl2 = dl.x*dl.x + dl.y*dl.y + dl.z*dl.z;
529 sp->p2[1].x = sp->p2[0].x + rs * dl.x;
530 sp->p2[1].y = sp->p2[0].y + rs * dl.y;
531 sp->p2[1].z = sp->p2[0].z + rs * dl.z;
532 lsum = lsum + log(df);
534 l = M_LOG2E / 2 * lsum / nl / sp->step2;
538 /* Anything that didn't explode has a finite attractor */
539 /* If Lyapunov is negative then it probably hit a fixed point or a
540 * limit cycle. Positive Lyapunov indicates a strange attractor. */
544 sp->size2 = max.x - min.x;
546 if(s > sp->size2) sp->size2 = s;
548 if(s > sp->size2) sp->size2 = s;
550 sp->mid2.x = (max.x + min.x) / 2;
551 sp->mid2.y = (max.y + min.y) / 2;
552 sp->mid2.z = (max.z + min.z) / 2;
554 if(sqrt(maxv2) > sp->size2 * 0.2) {
555 /* Flowing too fast, reduce step size. This
556 helps to eliminate high-speed limit cycles,
557 which can show +ve Lyapunov due to integration
564 /* Sets up initial conditions for a flow without all the extra baggage
565 that goes with init_flow */
567 restart_flow(ModeInfo * mi)
574 sp = &flows[MI_SCREEN(mi)];
577 /* Re-Initialize point positions, velocities, etc. */
578 for (b = 0; b < sp->beecount; b++) {
579 X(0, b) = Gauss_Rand(sp->range.x);
580 Y(0, b) = (sp->yperiod > 0)?
581 balance_rand(sp->range.y) : Gauss_Rand(sp->range.y);
582 Z(0, b) = Gauss_Rand(sp->range.z);
586 /* Returns true if line was behind a clip plane, or it clips the line */
587 /* nx,ny,nz is the normal to the plane. d is the distance from the origin */
588 /* s and e are the end points of the line to be clipped */
590 clip(double nx, double ny, double nz, double d, dvector *s, dvector *e)
596 front1 = (nx*s->x + ny*s->y + nz*s->z >= -d);
597 front2 = (nx*e->x + ny*e->y + nz*e->z >= -d);
598 if (!front1 && !front2) return 1;
599 if (front1 && front2) return 0;
604 /* Find t in line equation */
605 t = ( -d - nx*s->x - ny*s->y - nz*s->z) / ( nx*w.x + ny*w.y + nz*w.z);
607 p.x = s->x + w.x * t;
608 p.y = s->y + w.y * t;
609 p.z = s->z + w.z * t;
611 /* Move clipped point to the intersection */
625 init_flow (ModeInfo * mi)
631 if ((flows = (flowstruct *) calloc(MI_NUM_SCREENS(mi),
632 sizeof (flowstruct))) == NULL)
635 sp = &flows[MI_SCREEN(mi)];
639 sp->taillen = MI_SIZE(mi);
640 if (sp->taillen < -MINTRAIL) {
641 /* Change by sqrt so it seems more variable */
642 sp->taillen = NRAND((int)sqrt((double) (-sp->taillen - MINTRAIL + 1)));
643 sp->taillen = sp->taillen * sp->taillen + MINTRAIL;
644 } else if (sp->taillen < MINTRAIL) {
645 sp->taillen = MINTRAIL;
648 if(!rotatep && !ridep) rotatep = True; /* We need at least one viewpoint */
650 /* Start camera at Orbit or Bee */
656 sp->chasetime = 1; /* Go directly to target */
660 sp->step2 = INITIALSTEP;
662 /* Zero parameter set */
663 memset(sp->par2, 0, N_PARS * sizeof(dvector));
665 /* Set up standard examples */
666 switch (NRAND((periodicp) ? 5 : 3)) {
674 sp->par2[Y].x = 10 + balance_rand(5*0); /* a */
675 sp->par2[X].x = - sp->par2[Y].x; /* -a */
676 sp->par2[X].y = 28 + balance_rand(5*0); /* b */
680 sp->par2[Z].z = - 2 + balance_rand(1*0); /* -c */
690 sp->par2[Z].x = -2 + balance_rand(1); /* a */
692 sp->par2[Y].y = 0.2 + balance_rand(0.1); /* b */
693 sp->par2[C].z = 0.2 + balance_rand(0.1); /* c */
695 sp->par2[Z].z = -5.7;
701 z' = 0.2 + z(x - 5.7)
703 name = "RosslerCone";
705 sp->par2[Z].x = -2; /* a */
707 sp->par2[Y].y = 0.2; /* b */
708 sp->par2[ZZ].y = -0.331 + balance_rand(0.01); /* c */
711 sp->par2[Z].z = -5.7;
717 z' = 0.7x + az(0.1 - x^2)
721 sp->par2[SINY].x = 0.35 + balance_rand(0.25); /* b */
722 sp->par2[C].y = 1.57; /* c */
724 sp->par2[Z].z = 1 + balance_rand(0.5); /* a/10 */
725 sp->par2[XXZ].z = -10 * sp->par2[Z].z; /* -a */
726 sp->yperiod = 2 * M_PI;
730 x' = -ax - z/2 - z^3/8 + b sin(y)
735 sp->par2[X].x = -0.2 + balance_rand(0.1); /* a */
736 sp->par2[Z].x = -0.5;
737 sp->par2[ZZZ].x = -0.125;
738 sp->par2[SINY].x = 27.0 + balance_rand(3.0); /* b */
739 sp->par2[C].y = 1.33; /* c */
741 sp->yperiod = 2 * M_PI;
749 if(sp->yperiod > 0) {
751 /* periodic flows show either uniform distribution or a
752 snapshot on the 'time' axis */
753 sp->range.y = NRAND(2)? sp->yperiod : 0;
759 /* Run discoverer to set up bounding box, etc. Lyapunov will
760 probably be innaccurate, since we're only running it once, but
761 we're using known strange attractors so it should be ok. */
763 if(MI_IS_VERBOSE(mi))
765 "flow: Lyapunov exponent: %g, step: %g, size: %g (%s)\n",
766 sp->lyap2, sp->step2, sp->size2, name);
767 /* Install new params */
768 sp->lyap = sp->lyap2;
769 sp->size = sp->size2;
771 sp->step = sp->step2;
772 memcpy(sp->par, sp->par2, sizeof(sp->par2));
774 sp->count2 = 0; /* Reset search */
777 sp->beecount = MI_COUNT(mi);
778 if (sp->beecount < 0) { /* random variations */
779 sp->beecount = NRAND(-sp->beecount) + 1; /* Minimum 1 */
782 # ifdef HAVE_COCOA /* Don't second-guess Quartz's double-buffering */
786 if(dbufp) { /* Set up double buffer */
787 if (sp->buffer != None)
788 XFreePixmap(MI_DISPLAY(mi), sp->buffer);
789 sp->buffer = XCreatePixmap(MI_DISPLAY(mi), MI_WINDOW(mi),
790 MI_WIDTH(mi), MI_HEIGHT(mi), MI_DEPTH(mi));
792 sp->buffer = MI_WINDOW(mi);
794 /* no "NoExpose" events from XCopyArea wanted */
795 XSetGraphicsExposures(MI_DISPLAY(mi), MI_GC(mi), False);
797 /* Make sure we're using 'thin' lines */
798 XSetLineAttributes(MI_DISPLAY(mi), MI_GC(mi), 0, LineSolid, CapNotLast,
801 /* Clear the background (may be slow depending on user prefs). */
804 /* Allocate memory. */
805 if (sp->csegs == NULL) {
806 allocate(sp->csegs, XSegment,
807 (sp->beecount + BOX_L) * MI_NPIXELS(mi) * sp->taillen);
808 allocate(sp->cnsegs, int, MI_NPIXELS(mi));
809 allocate(sp->old_segs, XSegment, sp->beecount * sp->taillen);
810 allocate(sp->p, dvector, sp->beecount * sp->taillen);
813 /* Initialize point positions, velocities, etc. */
816 /* Set up camera tail */
817 X(1, 0) = sp->cam[1].x = 0;
818 Y(1, 0) = sp->cam[1].y = 0;
819 Z(1, 0) = sp->cam[1].z = 0;
823 draw_flow (ModeInfo * mi)
827 double M[3][3]; /* transformation matrix */
828 flowstruct *sp = NULL;
834 sp = &flows[MI_SCREEN(mi)];
835 if (sp->csegs == NULL)
838 #ifdef HAVE_COCOA /* Don't second-guess Quartz's double-buffering */
839 XClearWindow (MI_DISPLAY(mi), MI_WINDOW(mi));
842 /* multiplier for indexing segment arrays. Used in IX macro, etc. */
843 segindex = (sp->beecount + BOX_L) * sp->taillen;
846 if(sp->count2 == 0) { /* start new search */
847 sp->step2 = INITIALSTEP;
848 /* Pick random parameters. Actual range is irrelevant
849 since parameter scale determines flow speed but not
851 for(i=0; i< N_PARS; i++) {
852 sp->par2[i].x = Gauss_Rand(1.0);
853 sp->par2[i].y = Gauss_Rand(1.0);
854 sp->par2[i].z = Gauss_Rand(1.0);
857 if(!discover(mi)) { /* Flow exploded, reset. */
861 sp->count2 = 0; /* Attractor found, but it's not strange */
862 }else if(sp->count2 > 1000000) { /* This one will do */
863 sp->count2 = 0; /* Reset search */
864 if(MI_IS_VERBOSE(mi))
866 "flow: Lyapunov exponent: %g, step: %g, size: %g (unnamed)\n",
867 sp->lyap2, sp->step2, sp->size2);
868 /* Install new params */
869 sp->lyap = sp->lyap2;
870 sp->size = sp->size2;
872 sp->step = sp->step2;
873 memcpy(sp->par, sp->par2, sizeof(sp->par2));
875 /* If we're allowed to zoom out, do so now, so that we
876 get a look at the new attractor. */
877 if(sp->chaseto == BEE && rotatep) {
881 /* Reset initial conditions, so we don't get
882 misleading artifacts in the particle density. */
888 /* Reset segment buffers */
889 for (col = 0; col < MI_NPIXELS(mi); col++)
892 MI_IS_DRAWN(mi) = True;
894 /* Calculate circling POV [Chalky]*/
895 sp->circle[1] = sp->circle[0];
896 sp->circle[0].x = sp->size * 2 * sin(sp->count / 100.0) *
897 (-0.6 + 0.4 *cos(sp->count / 500.0)) + sp->mid.x;
898 sp->circle[0].y = sp->size * 2 * cos(sp->count / 100.0) *
899 (0.6 + 0.4 *cos(sp->count / 500.0)) + sp->mid.y;
900 sp->circle[0].z = sp->size * 2 * sin(sp->count / 421.0) + sp->mid.z;
902 /* Timed chase instead of Chalkie's Bistable oscillator [TDA] */
903 if(rotatep && ridep) {
904 if(sp->chaseto == BEE && NRAND(1000) == 0){
907 }else if(NRAND(4000) == 0){
913 /* Set up orientation matrix */
915 double x[3], p[3], x2=0, xp=0;
918 /* Chasetime is here to guarantee the camera makes it all the
919 way to the target in a finite number of steps. */
920 if(sp->chasetime > 1)
923 if(sp->chaseto == BEE){
924 /* Camera Head targets bee 0 */
925 sp->cam[0].x += (X(0, 0) - sp->cam[0].x)/sp->chasetime;
926 sp->cam[0].y += (Y(0, 0) - sp->cam[0].y)/sp->chasetime;
927 sp->cam[0].z += (Z(0, 0) - sp->cam[0].z)/sp->chasetime;
929 /* Camera Tail targets previous position of bee 0 */
930 sp->cam[1].x += (X(1, 0) - sp->cam[1].x)/sp->chasetime;
931 sp->cam[1].y += (Y(1, 0) - sp->cam[1].y)/sp->chasetime;
932 sp->cam[1].z += (Z(1, 0) - sp->cam[1].z)/sp->chasetime;
934 /* Camera Wing targets bee 1 */
935 sp->cam[2].x += (X(0, 1) - sp->cam[2].x)/sp->chasetime;
936 sp->cam[2].y += (Y(0, 1) - sp->cam[2].y)/sp->chasetime;
937 sp->cam[2].z += (Z(0, 1) - sp->cam[2].z)/sp->chasetime;
939 /* Camera Head targets Orbiter */
940 sp->cam[0].x += (sp->circle[0].x - sp->cam[0].x)/sp->chasetime;
941 sp->cam[0].y += (sp->circle[0].y - sp->cam[0].y)/sp->chasetime;
942 sp->cam[0].z += (sp->circle[0].z - sp->cam[0].z)/sp->chasetime;
944 /* Camera Tail targets diametrically opposite the middle
945 of the bounding box from the Orbiter */
947 (2*sp->circle[0].x - sp->mid.x - sp->cam[1].x)/sp->chasetime;
949 (2*sp->circle[0].y - sp->mid.y - sp->cam[1].y)/sp->chasetime;
951 (2*sp->circle[0].z - sp->mid.z - sp->cam[1].z)/sp->chasetime;
952 /* Camera Wing targets previous position of Orbiter */
953 sp->cam[2].x += (sp->circle[1].x - sp->cam[2].x)/sp->chasetime;
954 sp->cam[2].y += (sp->circle[1].y - sp->cam[2].y)/sp->chasetime;
955 sp->cam[2].z += (sp->circle[1].z - sp->cam[2].z)/sp->chasetime;
958 /* Viewpoint from Tail of camera */
959 sp->centre.x=sp->cam[1].x;
960 sp->centre.y=sp->cam[1].y;
961 sp->centre.z=sp->cam[1].z;
964 x[0] = sp->cam[0].x - sp->cam[1].x;
965 x[1] = sp->cam[0].y - sp->cam[1].y;
966 x[2] = sp->cam[0].z - sp->cam[1].z;
969 p[0] = sp->cam[2].x - sp->cam[1].x;
970 p[1] = sp->cam[2].y - sp->cam[1].y;
971 p[2] = sp->cam[2].z - sp->cam[1].z;
974 /* So long as X and P don't collide, these can be used to form
975 three mutually othogonal axes: X, (X x P) x X and X x P.
976 After being normalised to unit length, these form the
977 Orientation Matrix. */
980 x2+= x[i]*x[i]; /* X . X */
981 xp+= x[i]*p[i]; /* X . P */
982 M[0][i] = x[i]; /* X */
985 for(i=0; i<3; i++) /* (X x P) x X */
986 M[1][i] = x2*p[i] - xp*x[i]; /* == (X . X) P - (X . P) X */
988 M[2][0] = x[1]*p[2] - x[2]*p[1]; /* X x P */
989 M[2][1] = -x[0]*p[2] + x[2]*p[0];
990 M[2][2] = x[0]*p[1] - x[1]*p[0];
995 for(i=0; i<3; i++) A+=M[j][i]*M[j][i]; /* sum squares */
998 for(i=0; i<3; i++) M[j][i]/=A;
1001 if(sp->chaseto == BEE) {
1002 X(0, 1)=X(0, 0)+M[1][0]*sp->step; /* adjust neighbour */
1003 Y(0, 1)=Y(0, 0)+M[1][1]*sp->step;
1004 Z(0, 1)=Z(0, 0)+M[1][2]*sp->step;
1008 /* <=- Bounding Box -=> */
1010 for (b = 0; b < BOX_L; b++) {
1012 /* Chalky's clipping code, Only used for the box */
1013 /* clipping trails is slow and of little benefit. [TDA] */
1014 int p1 = lines[b][0];
1015 int p2 = lines[b][1];
1017 double x1=box[p1][0]* sp->size/2 + sp->mid.x - sp->centre.x;
1018 double y1=box[p1][1]* sp->size/2 + sp->mid.y - sp->centre.y;
1019 double z1=box[p1][2]* sp->size/2 + sp->mid.z - sp->centre.z;
1020 double x2=box[p2][0]* sp->size/2 + sp->mid.x - sp->centre.x;
1021 double y2=box[p2][1]* sp->size/2 + sp->mid.y - sp->centre.y;
1022 double z2=box[p2][2]* sp->size/2 + sp->mid.z - sp->centre.z;
1024 A1.x=M[0][0]*x1 + M[0][1]*y1 + M[0][2]*z1;
1025 A1.y=M[1][0]*x1 + M[1][1]*y1 + M[1][2]*z1;
1026 A1.z=M[2][0]*x1 + M[2][1]*y1 + M[2][2]*z1 + EYEHEIGHT * sp->size;
1027 A2.x=M[0][0]*x2 + M[0][1]*y2 + M[0][2]*z2;
1028 A2.y=M[1][0]*x2 + M[1][1]*y2 + M[1][2]*z2;
1029 A2.z=M[2][0]*x2 + M[2][1]*y2 + M[2][2]*z2 + EYEHEIGHT * sp->size;
1031 /* Clip in 3D before projecting down to 2D. A 2D clip
1032 after projection wouldn't be able to handle lines that
1034 if (clip(1, 0, 0,-1, &A1, &A2) || /* Screen */
1035 clip(1, 2, 0, 0, &A1, &A2) || /* Left */
1036 clip(1,-2, 0, 0, &A1, &A2) || /* Right */
1037 clip(1,0, 2.0*MI_WIDTH(mi)/MI_HEIGHT(mi), 0, &A1, &A2)||/*UP*/
1038 clip(1,0,-2.0*MI_WIDTH(mi)/MI_HEIGHT(mi), 0, &A1, &A2))/*Down*/
1041 /* Colour according to bee */
1042 col = b % (MI_NPIXELS(mi) - 1);
1044 sp->csegs[IX(col)].x1 = MI_WIDTH(mi)/2 + MI_WIDTH(mi) * A1.y/A1.x;
1045 sp->csegs[IX(col)].y1 = MI_HEIGHT(mi)/2 + MI_WIDTH(mi) * A1.z/A1.x;
1046 sp->csegs[IX(col)].x2 = MI_WIDTH(mi)/2 + MI_WIDTH(mi) * A2.y/A2.x;
1047 sp->csegs[IX(col)].y2 = MI_HEIGHT(mi)/2 + MI_WIDTH(mi) * A2.z/A2.x;
1053 for (b = 0; b < sp->beecount; b++) {
1054 if(fabs(X(0, b)) > LOST_IN_SPACE ||
1055 fabs(Y(0, b)) > LOST_IN_SPACE ||
1056 fabs(Z(0, b)) > LOST_IN_SPACE){
1057 if(sp->chaseto == BEE && b == 0){
1058 /* Lost camera bee. Need to replace it since
1059 rerunning init_flow could lose us a hard-won new
1060 attractor. Try moving it very close to a random
1061 other bee. This way we have a good chance of being
1062 close to the attractor and not forming a false
1063 artifact. If we've lost many bees this may need to
1065 /* Don't worry about camera wingbee. It stays close
1066 to the main camera bee no matter what happens. */
1067 int newb = 1 + NRAND(sp->beecount - 1);
1068 X(0, 0) = X(0, newb) + 0.001;
1069 Y(0, 0) = Y(0, newb);
1070 Z(0, 0) = Z(0, newb);
1071 if(MI_IS_VERBOSE(mi))
1073 "flow: resetting lost camera near bee %d\n",
1079 /* Age the tail. It's critical this be fast since
1080 beecount*taillen can be large. */
1081 memmove(B(1, b), B(0, b), (sp->taillen - 1) * sizeof(dvector));
1083 Iterate(B(0,b), sp->ODE, sp->par, sp->step);
1085 /* Don't show wingbee since he's not quite in the flow. */
1086 if(sp->chaseto == BEE && b == 1) continue;
1088 /* Colour according to bee */
1089 col = b % (MI_NPIXELS(mi) - 1);
1091 /* Fill the segment lists. */
1093 begin = 0; /* begin new trail */
1094 end = MIN(sp->taillen, sp->count); /* short trails at first */
1095 for(i=0; i < end; i++){
1096 double x = X(i,b)-sp->centre.x;
1097 double y = Y(i,b)*(sp->yperiod < 0? (sp->size/sp->yperiod) :1)
1099 double z = Z(i,b)-sp->centre.z;
1100 double XM=M[0][0]*x + M[0][1]*y + M[0][2]*z;
1101 double YM=M[1][0]*x + M[1][1]*y + M[1][2]*z;
1102 double ZM=M[2][0]*x + M[2][1]*y + M[2][2]*z + EYEHEIGHT * sp->size;
1105 swarm++; /* count the remaining bees */
1106 if(sp->yperiod > 0 && Y(i,b) > sp->yperiod){
1108 Y(i,b) -= sp->yperiod;
1109 /* hide tail to prevent streaks in Y. Streaks in X,Z
1110 are ok, they help to outline the Poincare'
1112 for(j = i; j < end; j++) Y(j,b) = Y(i,b);
1117 if(XM <= 0){ /* off screen - new trail */
1121 absx = MI_WIDTH(mi)/2 + MI_WIDTH(mi) * YM/XM;
1122 absy = MI_HEIGHT(mi)/2 + MI_WIDTH(mi) * ZM/XM;
1123 /* Performance bottleneck */
1124 if(absx <= 0 || absx >= MI_WIDTH(mi) ||
1125 absy <= 0 || absy >= MI_HEIGHT(mi)) {
1126 /* off screen - new trail */
1130 if(i > begin) { /* complete previous segment */
1131 sp->csegs[IX(col)].x2 = absx;
1132 sp->csegs[IX(col)].y2 = absy;
1136 if(i < end -1){ /* start new segment */
1137 sp->csegs[IX(col)].x1 = absx;
1138 sp->csegs[IX(col)].y1 = absy;
1144 XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_BLACK_PIXEL(mi));
1145 if (dbufp) { /* In Double Buffer case, prepare off-screen copy */
1146 /* For slow systems, this can be the single biggest bottleneck
1147 in the program. These systems may be better of not using
1148 the double buffer. */
1149 XFillRectangle(MI_DISPLAY(mi), sp->buffer, MI_GC(mi), 0, 0,
1150 MI_WIDTH(mi), MI_HEIGHT(mi));
1151 } else { /* Otherwise, erase previous segment list directly */
1152 XDrawSegments(MI_DISPLAY(mi), sp->buffer, MI_GC(mi),
1153 sp->old_segs, sp->nold_segs);
1157 if (MI_NPIXELS(mi) > 2){ /* colour */
1159 for (col = 0; col < MI_NPIXELS(mi) - 1; col++)
1160 if (sp->cnsegs[col] > 0) {
1161 if(sp->cnsegs[col] > mn) mn = sp->cnsegs[col];
1162 XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_PIXEL(mi, col+1));
1163 /* This is usually the biggest bottleneck on most
1164 systems. The maths load is insignificant compared
1166 XDrawSegments(MI_DISPLAY(mi), sp->buffer, MI_GC(mi),
1167 sp->csegs + col * segindex, sp->cnsegs[col]);
1169 } else { /* mono handled seperately since xlockmore uses '1' for
1171 XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_WHITE_PIXEL(mi));
1172 XDrawSegments(MI_DISPLAY(mi), sp->buffer, MI_GC(mi),
1173 sp->csegs, sp->cnsegs[0]);
1175 if (dbufp) { /* In Double Buffer case, this updates the screen */
1176 XCopyArea(MI_DISPLAY(mi), sp->buffer, MI_WINDOW(mi), MI_GC(mi), 0, 0,
1177 MI_WIDTH(mi), MI_HEIGHT(mi), 0, 0);
1178 } else { /* Otherwise, screen is already updated. Copy segments
1179 to erase-list to be erased directly next time. */
1181 for (col = 0; col < MI_NPIXELS(mi) - 1; col++) {
1182 memcpy(sp->old_segs + c, sp->csegs + col * segindex,
1183 sp->cnsegs[col] * sizeof(XSegment));
1184 c += sp->cnsegs[col];
1189 if(sp->count > 1 && swarm == 0) { /* all gone */
1190 if(MI_IS_VERBOSE(mi))
1191 fprintf(stdout, "flow: all gone at %d\n", sp->count);
1195 if(sp->count++ > MI_CYCLES(mi)){ /* Time's up. If we haven't
1196 found anything new by now we
1197 should pick a new standard
1204 reshape_flow(ModeInfo * mi, int width, int height)
1211 release_flow (ModeInfo * mi)
1213 if (flows != NULL) {
1216 for (screen = 0; screen < MI_NUM_SCREENS(mi); screen++)
1217 free_flow(&flows[screen]);
1219 flows = (flowstruct *) NULL;
1224 refresh_flow (ModeInfo * mi)
1226 if(!dbufp) MI_CLEARWINDOW(mi);
1229 XSCREENSAVER_MODULE ("Flow", flow)
1231 #endif /* MODE_flow */