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