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