git.hungrycats.org Git - xscreensaver/blob - hacks/flow.c

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