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 static const char sccsid[] = "@(#)flow.c        5.00 2000/11/01 xlockmore";
   6 #endif
   7
   8 /*-
   9  * Copyright (c) 1996 by Tim Auckland <tda10.geo@yahoo.com>
  10  * Incorporating some code from Stephen Davies Copyright (c) 2000
  11  *
  12  * Search code based on techniques described in "Strange Attractors:
  13  * Creating Patterns in Chaos" by Julien C. Sprott
  14  *
  15  * Permission to use, copy, modify, and distribute this software and its
  16  * documentation for any purpose and without fee is hereby granted,
  17  * provided that the above copyright notice appear in all copies and that
  18  * both that copyright notice and this permission notice appear in
  19  * supporting documentation.
  20  *
  21  * This file is provided AS IS with no warranties of any kind.  The author
  22  * shall have no liability with respect to the infringement of copyrights,
  23  * trade secrets or any patents by this file or any part thereof.  In no
  24  * event will the author be liable for any lost revenue or profits or
  25  * other special, indirect and consequential damages.
  26  *
  27  * "flow" shows a variety of continuous phase-space flows around strange
  28  * attractors.  It includes the well-known Lorentz mask (the "Butterfly"
  29  * of chaos fame), two forms of Rossler's "Folded Band" and Poincare'
  30  * sections of the "Birkhoff Bagel" and Duffing's forced occilator.  "flow"
  31  * can now discover new attractors.
  32  *
  33  * Revision History:
  34  *
  35  * 29-Oct-2004: [TDA] Discover Attractors unknown to science.
  36  *   Replace 2D rendering of Periodic Attractors with a 3D
  37  *   'interrupted' rendering.  Replace "-/+allow2d" with "-/+periodic"
  38  *   Replace all ODE formulae with completely generic forms.
  39  *   Add '-search' option to perform background high-speed discovery
  40  *   for completely new attractors without impacting rendering
  41  *   performance.
  42  *   Use gaussian distribution for initial point positions and for
  43  *   parameter search.
  44  *   Add "+dbuf" option to allow Double-Buffering to be turned off on
  45  *   slow X servers.
  46  *   Remove redundant '-zoom' option.  Now automatically zooms if both
  47  *   rotation and riding are permitted.
  48  *   Replace dynamic bounding box with static one pre-calculated
  49  *   during discovery phase.
  50  *   Simplify and fix bounding box clipping code.  Should now be safe
  51  *   to run without double buffer on all XFree86 servers if desired.
  52  * 12-Oct-2004: [TDA] Merge Xscreensaver and Xlockmore branches
  53  *   Added Chalky's orbital camera, but made the zooming work by
  54  *   flying the camera rather than interpolating the view transforms.
  55  *   Added Chalky's Bounding Box, but time-averaged the boundaries to
  56  *   let the lost bees escape.
  57  *   Added Chalky's 'view-frustrum' clipping, but only applying it to
  58  *   the Bounding Box.  Trails make clipping less useful.
  59  *   Added Chalky's "-slow" and "-freeze" options for compatibility,
  60  *   but haven't implemented the features, since the results are ugly
  61  *   and make no mathematical contribution.
  62  *   Added Double-Buffering as a work-around for a persistent XFree86
  63  *   bug that left debris on the screen.
  64  * 21-Mar-2003: [TDA] Trails added (XLockmore branch)
  65  * 01-Nov-2000: [TDA] Allocation checks (XLockmore branch)
  66  * 21-Feb-2000: [Chalky] Major hackage (Stephen Davies, chalky@null.net)
  67  *   (Xscreensaver branch)
  68  *   Forced perspective mode, added 3d box around attractor which
  69  *   involved coding 3d-planar-clipping against the view-frustrum
  70  *   thingy. Also made view alternate between piggybacking on a 'bee'
  71  *   to zooming around outside the attractor. Most bees slow down and
  72  *   stop, to make the structure of the attractor more obvious.
  73 * 28-Jan-1999: [TDA] Catch 'lost' bees in flow.c and disable them.
  74  *   (XLockmore branch)
  75  *   I chose to disable them rather than reinitialise them because
  76  *   reinitialising can produce fake attractors.
  77  *   This has allowed me to relax some of the parameters and initial
  78  *   conditions slightly to catch some of the more extreme cases.  As a
  79  *   result you may see some bees fly away at the start - these are the ones
  80  *   that 'missed' the attractor.  If the bee with the camera should fly
  81  *   away the mode will restart  :-)
  82  * 31-Nov-1998: [TDA] Added Duffing  (what a strange day that was :) DAB)
  83  *   Duffing's forced oscillator has been added to the formula list and
  84  *   the parameters section has been updated to display it in Poincare'
  85  *   section.
  86  * 30-Nov-1998: [TDA] Added travelling perspective option
  87  *   A more exciting point-of-view has been added to all autonomous flows.
  88  *   This views the flow as seen by a particle moving with the flow.  In the
  89  *   metaphor of the original code, I've attached a camera to one of the
  90  *   trained bees!
  91  * 30-Nov-1998: [TDA] Much code cleanup.
  92  * 09-Apr-1997: [TDA] Ported to xlockmore-4
  93  * 18-Jul-1996: Adapted from swarm.c Copyright (c) 1991 by Patrick J. Naughton.
  94  * 31-Aug-1990: Adapted from xswarm by Jeff Butterworth. (butterwo@ncsc.org).
  95  */
  96
  97 #ifdef STANDALONE
  98 # define MODE_flow
  99 # define DEFAULTS   "*delay:       10000 \n" \
 100                                         "*count:       3000  \n" \
 101                                         "*size:        -10   \n" \
 102                                         "*cycles:      10000 \n" \
 103                                         "*ncolors:     200   \n"
 104
 105 # define flow_handle_event 0
 106 # include "xlockmore.h"         /* in xscreensaver distribution */
 107 #else /* STANDALONE */
 108 # include "xlock.h"             /* in xlockmore distribution */
 109 #endif /* STANDALONE */
 110
 111 #ifdef MODE_flow
 112
 113 #define DEF_ROTATE   "TRUE"
 114 #define DEF_RIDE     "TRUE"
 115 #define DEF_BOX      "TRUE"
 116 #define DEF_PERIODIC "TRUE"
 117 #define DEF_SEARCH   "TRUE"
 118 #define DEF_DBUF     "TRUE"
 119
 120 static Bool rotatep;
 121 static Bool ridep;
 122 static Bool boxp;
 123 static Bool periodicp;
 124 static Bool searchp;
 125 static Bool dbufp;
 126
 127 static XrmOptionDescRec opts[] = {
 128         {"-rotate",   ".flow.rotate",   XrmoptionNoArg, "on"},
 129         {"+rotate",   ".flow.rotate",   XrmoptionNoArg, "off"},
 130         {"-ride",     ".flow.ride",     XrmoptionNoArg, "on"},
 131         {"+ride",     ".flow.ride",     XrmoptionNoArg, "off"},
 132         {"-box",      ".flow.box",      XrmoptionNoArg, "on"},
 133         {"+box",      ".flow.box",      XrmoptionNoArg, "off"},
 134         {"-periodic", ".flow.periodic", XrmoptionNoArg, "on"},
 135         {"+periodic", ".flow.periodic", XrmoptionNoArg, "off"},
 136         {"-search",   ".flow.search",   XrmoptionNoArg, "on"},
 137         {"+search",   ".flow.search",   XrmoptionNoArg, "off"},
 138         {"-dbuf",     ".flow.dbuf",     XrmoptionNoArg, "on"},
 139         {"+dbuf",     ".flow.dbuf",     XrmoptionNoArg, "off"},
 140 };
 141
 142 static argtype vars[] = {
 143     {&rotatep,   "rotate",   "Rotate",   DEF_ROTATE,   t_Bool},
 144     {&ridep,     "ride",     "Ride",     DEF_RIDE,     t_Bool},
 145     {&boxp,      "box",      "Box",      DEF_BOX,      t_Bool},
 146     {&periodicp, "periodic", "Periodic", DEF_PERIODIC, t_Bool},
 147     {&searchp,   "search",   "Search",   DEF_SEARCH,   t_Bool},
 148     {&dbufp,     "dbuf",     "Dbuf",     DEF_DBUF,     t_Bool},
 149 };
 150
 151 static OptionStruct desc[] = {
 152     {"-/+rotate",   "turn on/off rotating around attractor."},
 153     {"-/+ride",     "turn on/off ride in the flow."},
 154     {"-/+box",      "turn on/off bounding box."},
 155     {"-/+periodic", "turn on/off periodic attractors."},
 156     {"-/+search",   "turn on/off search for new attractors."},
 157     {"-/+dbuf",     "turn on/off double buffering."},
 158 };
 159
 160 ENTRYPOINT ModeSpecOpt flow_opts =
 161 {sizeof opts / sizeof opts[0], opts,
 162  sizeof vars / sizeof vars[0], vars, desc};
 163
 164 #ifdef USE_MODULES
 165 ModStruct   flow_description = {
 166         "flow", "init_flow", "draw_flow", "release_flow",
 167         "refresh_flow", "init_flow", NULL, &flow_opts,
 168         1000, 1024, 10000, -10, 200, 1.0, "",
 169         "Shows dynamic strange attractors", 0, NULL
 170 };
 171
 172 #endif
 173
 174 typedef struct { double x, y, z; } dvector;
 175
 176 #define N_PARS 20 /* Enough for Full Cubic or Periodic Cubic */
 177 typedef dvector Par[N_PARS];
 178 enum { /* Name the parameter indices to make it easier to write
 179                   standard examples */
 180         C,
 181         X,XX,XXX,XXY,XXZ,XY,XYY,XYZ,XZ,XZZ,
 182         Y,YY,YYY,YYZ,YZ,YZZ,
 183         Z,ZZ,ZZZ,
 184         SINY = XY /* OK to overlap in this case */
 185 };
 186
 187 /* Camera target [TDA] */
 188 typedef enum {
 189         ORBIT = 0,
 190         BEE = 1
 191 } Chaseto;
 192
 193 /* Macros */
 194 #define IX(C) ((C) * segindex + sp->cnsegs[(C)])
 195 #define B(t,b)  (sp->p + (t) + (b) * sp->taillen)
 196 #define X(t,b)  (B((t),(b))->x)
 197 #define Y(t,b)  (B((t),(b))->y)
 198 #define Z(t,b)  (B((t),(b))->z)
 199 #define balance_rand(v) ((LRAND()/MAXRAND*(v))-((v)/2)) /* random around 0 */
 200 #define LOST_IN_SPACE 2000.0
 201 #define INITIALSTEP 0.04
 202 #define EYEHEIGHT   0.005
 203 #define MINTRAIL 2
 204 #define BOX_L 36
 205
 206 /* Points that make up the box (normalized coordinates) */
 207 static const double box[][3] = {
 208         {1,1,1},   /* 0 */
 209         {1,1,-1},  /* 1 */
 210         {1,-1,-1}, /* 2 */
 211         {1,-1,1},  /* 3 */
 212         {-1,1,1},  /* 4 */
 213         {-1,1,-1}, /* 5 */
 214         {-1,-1,-1},/* 6 */
 215         {-1,-1,1}, /* 7 */
 216         {1, .8, .8},
 217         {1, .8,-.8},
 218         {1,-.8,-.8},
 219         {1,-.8, .8},
 220         { .8,1, .8},
 221         { .8,1,-.8},
 222         {-.8,1,-.8},
 223         {-.8,1, .8},
 224         { .8, .8,1},
 225         { .8,-.8,1},
 226         {-.8,-.8,1},
 227         {-.8, .8,1},
 228         {-1, .8, .8},
 229         {-1, .8,-.8},
 230         {-1,-.8,-.8},
 231         {-1,-.8, .8},
 232         { .8,-1, .8},
 233         { .8,-1,-.8},
 234         {-.8,-1,-.8},
 235         {-.8,-1, .8},
 236         { .8, .8,-1},
 237         { .8,-.8,-1},
 238         {-.8,-.8,-1},
 239         {-.8, .8,-1}
 240 };
 241
 242 /* Lines connecting the box dots */
 243 static const double lines[][2] = {
 244         {0,1}, {1,2}, {2,3}, {3,0}, /* box */
 245         {4,5}, {5,6}, {6,7}, {7,4},
 246         {0,4}, {1,5}, {2,6}, {3,7},
 247         {4+4,5+4}, {5+4,6+4}, {6+4,7+4}, {7+4,4+4},
 248         {4+8,5+8}, {5+8,6+8}, {6+8,7+8}, {7+8,4+8},
 249         {4+12,5+12}, {5+12,6+12}, {6+12,7+12}, {7+12,4+12},
 250         {4+16,5+16}, {5+16,6+16}, {6+16,7+16}, {7+16,4+16},
 251         {4+20,5+20}, {5+20,6+20}, {6+20,7+20}, {7+20,4+20},
 252         {4+24,5+24}, {5+24,6+24}, {6+24,7+24}, {7+24,4+24},
 253 };
 254
 255 typedef struct {
 256         /* Variables used in rendering */
 257         dvector     cam[3]; /* camera flight path */
 258         int         chasetime;
 259         Chaseto         chaseto;
 260         Pixmap      buffer; /* Double Buffer */
 261         dvector         circle[2]; /* POV that circles around the scene */
 262         dvector     centre;             /* centre */
 263         int         beecount;   /* number of bees */
 264         XSegment   *csegs;          /* bee lines */
 265         int        *cnsegs;
 266         XSegment   *old_segs;   /* old bee lines */
 267         int         nold_segs;
 268         int         taillen;
 269
 270         /* Variables common to iterators */
 271         dvector  (*ODE) (Par par, double x, double y, double z);
 272         dvector  range; /* Initial conditions */
 273         double   yperiod; /* ODE's where Y is periodic. */
 274
 275         /* Variables used in iterating main flow */
 276         Par         par;
 277         dvector    *p;   /* bee positions x[time][bee#] */
 278         int         count;
 279         double      lyap;
 280         double      size;
 281         dvector     mid; /* Effective bounding box */
 282         double      step;
 283
 284         /* second set of variables, used for parallel search */
 285         Par         par2;
 286         dvector     p2[2];
 287         int         count2;
 288         double      lyap2;
 289         double      size2;
 290         dvector     mid2;
 291         double      step2;
 292
 293 } flowstruct;
 294
 295 static flowstruct *flows = (flowstruct *) NULL;
 296
 297 /*
 298  * Private functions
 299  */
 300
 301
 302 /* ODE functions */
 303
 304 /* Generic 3D Cubic Polynomial.  Includes all the Quadratics (Lorentz,
 305    Rossler) and much more! */
 306
 307 /* I considered offering a seperate 'Quadratic' option, since Cubic is
 308    clearly overkill for the standard examples, but the performance
 309    difference is too small to measure.  The compute time is entirely
 310    dominated by the XDrawSegments calls anyway. [TDA] */
 311 static dvector
 312 Cubic(Par a, double x, double y, double z)
 313 {
 314         dvector d;
 315         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 +
 316                 a[XXZ].x*x*x*z + a[XY].x*x*y + a[XYY].x*x*y*y + a[XYZ].x*x*y*z +
 317                 a[XZ].x*x*z + a[XZZ].x*x*z*z + a[Y].x*y + a[YY].x*y*y +
 318                 a[YYY].x*y*y*y + a[YYZ].x*y*y*z + a[YZ].x*y*z + a[YZZ].x*y*z*z +
 319                 a[Z].x*z + a[ZZ].x*z*z + a[ZZZ].x*z*z*z;
 320
 321         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 +
 322                 a[XXZ].y*x*x*z + a[XY].y*x*y + a[XYY].y*x*y*y + a[XYZ].y*x*y*z +
 323                 a[XZ].y*x*z + a[XZZ].y*x*z*z + a[Y].y*y + a[YY].y*y*y +
 324                 a[YYY].y*y*y*y + a[YYZ].y*y*y*z + a[YZ].y*y*z + a[YZZ].y*y*z*z +
 325                 a[Z].y*z + a[ZZ].y*z*z + a[ZZZ].y*z*z*z;
 326
 327         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 +
 328                 a[XXZ].z*x*x*z + a[XY].z*x*y + a[XYY].z*x*y*y + a[XYZ].z*x*y*z +
 329                 a[XZ].z*x*z + a[XZZ].z*x*z*z + a[Y].z*y + a[YY].z*y*y +
 330                 a[YYY].z*y*y*y + a[YYZ].z*y*y*z + a[YZ].z*y*z + a[YZZ].z*y*z*z +
 331                 a[Z].z*z + a[ZZ].z*z*z + a[ZZZ].z*z*z*z;
 332
 333         return d;
 334 }
 335
 336 /* 3D Cubic in (x,z) with periodic sinusoidal forcing term in x.  y is
 337    the independent periodic (time) axis.  This includes Birkhoff's
 338    Bagel and Duffing's Attractor */
 339 static dvector
 340 Periodic(Par a, double x, double y, double z)
 341 {
 342         dvector d;
 343
 344         d.x = a[C].x + a[X].x*x + a[XX].x*x*x + a[XXX].x*x*x*x +
 345                 a[XXZ].x*x*x*z + a[XZ].x*x*z + a[XZZ].x*x*z*z + a[Z].x*z +
 346                 a[ZZ].x*z*z + a[ZZZ].x*z*z*z + a[SINY].x*sin(y);
 347
 348         d.y = a[C].y;
 349
 350         d.z = a[C].z + a[X].z*x + a[XX].z*x*x + a[XXX].z*x*x*x +
 351                 a[XXZ].z*x*x*z + a[XZ].z*x*z + a[XZZ].z*x*z*z + a[Z].z*z +
 352                 a[ZZ].z*z*z + a[ZZZ].z*z*z*z;
 353
 354         return d;
 355 }
 356
 357 /* Numerical integration of the ODE using 2nd order Runge Kutta.
 358    Returns length^2 of the update, so that we can detect if the step
 359    size needs reducing. */
 360 static double
 361 Iterate(dvector *p, dvector(*ODE)(Par par, double x, double y, double z),
 362                 Par par, double step)
 363 {
 364         dvector     k1, k2, k3;
 365
 366         k1 = ODE(par, p->x, p->y, p->z);
 367         k1.x *= step;
 368         k1.y *= step;
 369         k1.z *= step;
 370         k2 = ODE(par, p->x + k1.x, p->y + k1.y, p->z + k1.z);
 371         k2.x *= step;
 372         k2.y *= step;
 373         k2.z *= step;
 374         k3.x = (k1.x + k2.x) / 2.0;
 375         k3.y = (k1.y + k2.y) / 2.0;
 376         k3.z = (k1.z + k2.z) / 2.0;
 377
 378         p->x += k3.x;
 379         p->y += k3.y;
 380         p->z += k3.z;
 381
 382         return k3.x*k3.x + k3.y*k3.y + k3.z*k3.z;
 383 }
 384
 385 /* Memory functions */
 386
 387 #define deallocate(p,t) if (p!=NULL) {free(p); p=(t*)NULL; }
 388 #define allocate(p,t,s) if ((p=(t*)malloc(sizeof(t)*s))==NULL)\
 389 {free_flow(sp);return;}
 390
 391 static void
 392 free_flow(flowstruct *sp)
 393 {
 394         deallocate(sp->csegs, XSegment);
 395         deallocate(sp->cnsegs, int);
 396         deallocate(sp->old_segs, XSegment);
 397         deallocate(sp->p, dvector);
 398 }
 399
 400 /* Generate Gaussian random number: mean 0, "amplitude" A (actually
 401    A is 3*standard deviation). */
 402
 403 /* Note this generates a pair of gaussian variables, so it saves one
 404    to give out next time it's called */
 405 static double
 406 Gauss_Rand(double A)
 407 {
 408         static double d;
 409         static Bool ready = 0;
 410         if(ready) {
 411                 ready = 0;
 412                 return A/3 * d;
 413         } else {
 414                 double x, y, w;
 415                 do {
 416                         x = 2.0 * (double)LRAND() / MAXRAND - 1.0;
 417                         y = 2.0 * (double)LRAND() / MAXRAND - 1.0;
 418                         w = x*x + y*y;
 419                 } while(w >= 1.0);
 420
 421                 w = sqrt((-2 * log(w))/w);
 422                 ready = 1;
 423                 d =          x * w;
 424                 return A/3 * y * w;
 425         }
 426 }
 427
 428 /* Attempt to discover new atractors by sending a pair of bees on a
 429    fast trip through the new flow and computing their Lyapunov
 430    exponent.  Returns False if the bees fly away.
 431
 432    If the bees stay bounded, the new bounds and the Lyapunov exponent
 433    are stored in sp and the function returns True.
 434
 435    Repeat invocations continue the flow and improve the accuracy of
 436    the bounds and the Lyapunov exponent.  Set sp->count2 to zero to
 437    start a new test.
 438
 439    Acts on alternate variable set, so that it can be run in parallel
 440    with the main flow */
 441
 442 static Bool
 443 discover(ModeInfo * mi)
 444 {
 445         flowstruct *sp;
 446         double l = 0;
 447         dvector dl;
 448         dvector max, min;
 449         double dl2, df, rs, lsum = 0, s, maxv2 = 0, v2;
 450
 451         int N, i, nl = 0;
 452
 453         if (flows == NULL)
 454                 return 0;
 455         sp = &flows[MI_SCREEN(mi)];
 456
 457         if(sp->count2 == 0) {
 458                 /* initial conditions */
 459                 sp->p2[0].x = Gauss_Rand(sp->range.x);
 460                 sp->p2[0].y = (sp->yperiod > 0)?
 461                         balance_rand(sp->range.y) : Gauss_Rand(sp->range.y);
 462                 sp->p2[0].z = Gauss_Rand(sp->range.z);
 463
 464                 /* 1000 steps to find an attractor */
 465                 /* Most cases explode out here */
 466                 for(N=0; N < 1000; N++){
 467                         Iterate(sp->p2, sp->ODE, sp->par2, sp->step2);
 468                         if(sp->yperiod > 0 && sp->p2[0].y > sp->yperiod)
 469                                 sp->p2[0].y -= sp->yperiod;
 470                         if(fabs(sp->p2[0].x) > LOST_IN_SPACE ||
 471                            fabs(sp->p2[0].y) > LOST_IN_SPACE ||
 472                            fabs(sp->p2[0].z) > LOST_IN_SPACE) {
 473                                 return 0;
 474                         }
 475                         sp->count2++;
 476                 }
 477                 /* Small perturbation */
 478                 sp->p2[1].x = sp->p2[0].x + 0.000001;
 479                 sp->p2[1].y = sp->p2[0].y;
 480                 sp->p2[1].z = sp->p2[0].z;
 481         }
 482
 483         /* Reset bounding box */
 484         max.x = min.x = sp->p2[0].x;
 485         max.y = min.y = sp->p2[0].y;
 486         max.z = min.z = sp->p2[0].z;
 487
 488         /* Compute Lyapunov Exponent */
 489
 490         /* (Technically, we're only estimating the largest Lyapunov
 491            Exponent, but that's all we need to know to determine if we
 492            have a strange attractor.) [TDA] */
 493
 494         /* Fly two bees close together */
 495         for(N=0; N < 5000; N++){
 496                 for(i=0; i< 2; i++) {
 497                         v2 = Iterate(sp->p2+i, sp->ODE, sp->par2, sp->step2);
 498                         if(sp->yperiod > 0 && sp->p2[i].y > sp->yperiod)
 499                                 sp->p2[i].y -= sp->yperiod;
 500
 501                         if(fabs(sp->p2[i].x) > LOST_IN_SPACE ||
 502                            fabs(sp->p2[i].y) > LOST_IN_SPACE ||
 503                            fabs(sp->p2[i].z) > LOST_IN_SPACE) {
 504                                 return 0;
 505                         }
 506                         if(v2 > maxv2) maxv2 = v2; /* Track max v^2 */
 507                 }
 508
 509                 /* find bounding box */
 510                 if ( sp->p2[0].x < min.x )      min.x = sp->p2[0].x;
 511                 else if ( sp->p2[0].x > max.x ) max.x = sp->p2[0].x;
 512                 if ( sp->p2[0].y < min.y )      min.y = sp->p2[0].y;
 513                 else if ( sp->p2[0].y > max.y ) max.y = sp->p2[0].y;
 514                 if ( sp->p2[0].z < min.z )      min.z = sp->p2[0].z;
 515                 else if ( sp->p2[0].z > max.z ) max.z = sp->p2[0].z;
 516
 517                 /* Measure how much we have to pull the two bees to prevent
 518                    them diverging. */
 519                 dl.x = sp->p2[1].x - sp->p2[0].x;
 520                 dl.y = sp->p2[1].y - sp->p2[0].y;
 521                 dl.z = sp->p2[1].z - sp->p2[0].z;
 522
 523                 dl2 = dl.x*dl.x + dl.y*dl.y + dl.z*dl.z;
 524                 if(dl2 > 0) {
 525                         df = 1e12 * dl2;
 526                         rs = 1/sqrt(df);
 527                         sp->p2[1].x = sp->p2[0].x + rs * dl.x;
 528                         sp->p2[1].y = sp->p2[0].y + rs * dl.y;
 529                         sp->p2[1].z = sp->p2[0].z + rs * dl.z;
 530                         lsum = lsum + log(df);
 531                         nl = nl + 1;
 532                         l = M_LOG2E / 2 * lsum / nl / sp->step2;
 533                 }
 534                 sp->count2++;
 535         }
 536         /* Anything that didn't explode has a finite attractor */
 537         /* If Lyapunov is negative then it probably hit a fixed point or a
 538      * limit cycle.  Positive Lyapunov indicates a strange attractor. */
 539
 540         sp->lyap2 = l;
 541
 542         sp->size2 = max.x - min.x;
 543         s = max.y - min.y;
 544         if(s > sp->size2) sp->size2 = s;
 545         s = max.z - min.z;
 546         if(s > sp->size2) sp->size2 = s;
 547
 548         sp->mid2.x = (max.x + min.x) / 2;
 549         sp->mid2.y = (max.y + min.y) / 2;
 550         sp->mid2.z = (max.z + min.z) / 2;
 551
 552         if(sqrt(maxv2) > sp->size2 * 0.2) {
 553                 /* Flowing too fast, reduce step size.  This
 554                    helps to eliminate high-speed limit cycles,
 555                    which can show +ve Lyapunov due to integration
 556                    inaccuracy. */
 557                 sp->step2 /= 2;
 558         }
 559         return 1;
 560 }
 561
 562 /* Sets up initial conditions for a flow without all the extra baggage
 563    that goes with init_flow */
 564 static void
 565 restart_flow(ModeInfo * mi)
 566 {
 567         flowstruct *sp;
 568         int         b;
 569
 570         if (flows == NULL)
 571                 return;
 572         sp = &flows[MI_SCREEN(mi)];
 573         sp->count = 0;
 574
 575         /* Re-Initialize point positions, velocities, etc. */
 576         for (b = 0; b < sp->beecount; b++) {
 577                 X(0, b) = Gauss_Rand(sp->range.x);
 578                 Y(0, b) = (sp->yperiod > 0)?
 579                         balance_rand(sp->range.y) : Gauss_Rand(sp->range.y);
 580                 Z(0, b) = Gauss_Rand(sp->range.z);
 581         }
 582 }
 583
 584 /* Returns true if line was behind a clip plane, or it clips the line */
 585 /* nx,ny,nz is the normal to the plane.   d is the distance from the origin */
 586 /* s and e are the end points of the line to be clipped */
 587 static int
 588 clip(double nx, double ny, double nz, double d, dvector *s, dvector *e)
 589 {
 590         int front1, front2;
 591         dvector w, p;
 592         double t;
 593
 594         front1 = (nx*s->x + ny*s->y + nz*s->z >= -d);
 595         front2 = (nx*e->x + ny*e->y + nz*e->z >= -d);
 596         if (!front1 && !front2) return 1;
 597         if (front1 && front2) return 0;
 598         w.x = e->x - s->x;
 599         w.y = e->y - s->y;
 600         w.z = e->z - s->z;
 601
 602         /* Find t in line equation */
 603         t = ( -d - nx*s->x - ny*s->y - nz*s->z) / ( nx*w.x + ny*w.y + nz*w.z);
 604
 605         p.x = s->x + w.x * t;
 606         p.y = s->y + w.y * t;
 607         p.z = s->z + w.z * t;
 608
 609         /* Move clipped point to the intersection */
 610         if (front2) {
 611                 *s = p;
 612         } else {
 613                 *e = p;
 614         }
 615         return 0;
 616 }
 617
 618 /*
 619  * Public functions
 620  */
 621
 622 ENTRYPOINT void
 623 init_flow (ModeInfo * mi)
 624 {
 625         flowstruct *sp;
 626         char       *name;
 627
 628         if (flows == NULL) {
 629                 if ((flows = (flowstruct *) calloc(MI_NUM_SCREENS(mi),
 630                                                                                    sizeof (flowstruct))) == NULL)
 631                         return;
 632         }
 633         sp = &flows[MI_SCREEN(mi)];
 634
 635         sp->count2 = 0;
 636
 637         sp->taillen = MI_SIZE(mi);
 638         if (sp->taillen < -MINTRAIL) {
 639                 /* Change by sqrt so it seems more variable */
 640                 sp->taillen = NRAND((int)sqrt((double) (-sp->taillen - MINTRAIL + 1)));
 641                 sp->taillen = sp->taillen * sp->taillen + MINTRAIL;
 642         } else if (sp->taillen < MINTRAIL) {
 643                 sp->taillen = MINTRAIL;
 644         }
 645
 646         if(!rotatep && !ridep) rotatep = True; /* We need at least one viewpoint */
 647
 648         /* Start camera at Orbit or Bee */
 649         if(rotatep) {
 650                 sp->chaseto = ORBIT;
 651         } else {
 652                 sp->chaseto = BEE;
 653         }
 654         sp->chasetime = 1; /* Go directly to target */
 655
 656         sp->lyap = 0;
 657         sp->yperiod = 0;
 658         sp->step2 = INITIALSTEP;
 659
 660         /* Zero parameter set */
 661         memset(sp->par2, 0, N_PARS * sizeof(dvector));
 662
 663         /* Set up standard examples */
 664         switch (NRAND((periodicp) ? 5 : 3)) {
 665         case 0:
 666                 /*
 667                   x' = a(y - x)
 668                   y' = x(b - z) - y
 669                   z' = xy - cz
 670                  */
 671                 name = "Lorentz";
 672                 sp->par2[Y].x = 10 + balance_rand(5*0); /* a */
 673                 sp->par2[X].x = - sp->par2[Y].x;        /* -a */
 674                 sp->par2[X].y = 28 + balance_rand(5*0); /* b */
 675                 sp->par2[XZ].y = -1;
 676                 sp->par2[Y].y = -1;
 677                 sp->par2[XY].z = 1;
 678                 sp->par2[Z].z = - 2 + balance_rand(1*0); /* -c */
 679                 break;
 680         case 1:
 681                 /*
 682                   x' = -(y + az)
 683                   y' = x + by
 684                   z' = c + z(x - 5.7)
 685                  */
 686                 name = "Rossler";
 687                 sp->par2[Y].x = -1;
 688                 sp->par2[Z].x = -2 + balance_rand(1); /* a */
 689                 sp->par2[X].y = 1;
 690                 sp->par2[Y].y = 0.2 + balance_rand(0.1); /* b */
 691                 sp->par2[C].z = 0.2 + balance_rand(0.1); /* c */
 692                 sp->par2[XZ].z = 1;
 693                 sp->par2[Z].z = -5.7;
 694                 break;
 695         case 2:
 696                 /*
 697                   x' = -(y + az)
 698                   y' = x + by - cz^2
 699                   z' = 0.2 + z(x - 5.7)
 700                  */
 701                 name = "RosslerCone";
 702                 sp->par2[Y].x = -1;
 703                 sp->par2[Z].x = -2; /* a */
 704                 sp->par2[X].y = 1;
 705                 sp->par2[Y].y = 0.2; /* b */
 706                 sp->par2[ZZ].y = -0.331 + balance_rand(0.01); /* c */
 707                 sp->par2[C].z = 0.2;
 708                 sp->par2[XZ].z = 1;
 709                 sp->par2[Z].z = -5.7;
 710                 break;
 711         case 3:
 712                 /*
 713                   x' = -z + b sin(y)
 714                   y' = c
 715                   z' = 0.7x + az(0.1 - x^2)
 716                  */
 717                 name = "Birkhoff";
 718                 sp->par2[Z].x = -1;
 719                 sp->par2[SINY].x = 0.35 + balance_rand(0.25); /* b */
 720                 sp->par2[C].y = 1.57; /* c */
 721                 sp->par2[X].z = 0.7;
 722                 sp->par2[Z].z = 1 + balance_rand(0.5); /* a/10 */
 723                 sp->par2[XXZ].z = -10 * sp->par2[Z].z; /* -a */
 724                 sp->yperiod = 2 * M_PI;
 725                 break;
 726         default:
 727                 /*
 728                   x' = -ax - z/2 - z^3/8 + b sin(y)
 729                   y' = c
 730                   z' = 2x
 731                  */
 732                 name = "Duffing";
 733                 sp->par2[X].x = -0.2 + balance_rand(0.1); /* a */
 734                 sp->par2[Z].x = -0.5;
 735                 sp->par2[ZZZ].x = -0.125;
 736                 sp->par2[SINY].x = 27.0 + balance_rand(3.0); /* b */
 737                 sp->par2[C].y = 1.33; /* c */
 738                 sp->par2[X].z = 2;
 739                 sp->yperiod = 2 * M_PI;
 740                 break;
 741
 742         }
 743
 744         sp->range.x = 5;
 745         sp->range.z = 5;
 746
 747         if(sp->yperiod > 0) {
 748                 sp->ODE = Periodic;
 749                 /* periodic flows show either uniform distribution or a
 750            snapshot on the 'time' axis */
 751                 sp->range.y = NRAND(2)? sp->yperiod : 0;
 752         } else {
 753                 sp->range.y = 5;
 754                 sp->ODE = Cubic;
 755         }
 756
 757         /* Run discoverer to set up bounding box, etc.  Lyapunov will
 758            probably be innaccurate, since we're only running it once, but
 759            we're using known strange attractors so it should be ok. */
 760         discover(mi);
 761         if(MI_IS_VERBOSE(mi))
 762                 fprintf(stdout,
 763                                 "flow: Lyapunov exponent: %g, step: %g, size: %g (%s)\n",
 764                                 sp->lyap2, sp->step2, sp->size2, name);
 765         /* Install new params */
 766         sp->lyap = sp->lyap2;
 767         sp->size = sp->size2;
 768         sp->mid = sp->mid2;
 769         sp->step = sp->step2;
 770         memcpy(sp->par, sp->par2, sizeof(sp->par2));
 771
 772         sp->count2 = 0; /* Reset search */
 773
 774         free_flow(sp);
 775         sp->beecount = MI_COUNT(mi);
 776         if (sp->beecount < 0) { /* random variations */
 777                 sp->beecount = NRAND(-sp->beecount) + 1; /* Minimum 1 */
 778         }
 779
 780 # ifdef HAVE_COCOA      /* Don't second-guess Quartz's double-buffering */
 781   dbufp = False;
 782 # endif
 783
 784         if(dbufp) { /* Set up double buffer */
 785                 if (sp->buffer != None)
 786                         XFreePixmap(MI_DISPLAY(mi), sp->buffer);
 787                 sp->buffer = XCreatePixmap(MI_DISPLAY(mi), MI_WINDOW(mi),
 788                                                                  MI_WIDTH(mi), MI_HEIGHT(mi), MI_DEPTH(mi));
 789         } else {
 790                 sp->buffer = MI_WINDOW(mi);
 791         }
 792         /* no "NoExpose" events from XCopyArea wanted */
 793         XSetGraphicsExposures(MI_DISPLAY(mi), MI_GC(mi), False);
 794
 795         /* Make sure we're using 'thin' lines */
 796         XSetLineAttributes(MI_DISPLAY(mi), MI_GC(mi), 0, LineSolid, CapNotLast,
 797                                            JoinMiter);
 798
 799         /* Clear the background (may be slow depending on user prefs). */
 800         MI_CLEARWINDOW(mi);
 801
 802         /* Allocate memory. */
 803         if (sp->csegs == NULL) {
 804                 allocate(sp->csegs, XSegment,
 805                                  (sp->beecount + BOX_L) * MI_NPIXELS(mi) * sp->taillen);
 806                 allocate(sp->cnsegs, int, MI_NPIXELS(mi));
 807                 allocate(sp->old_segs, XSegment, sp->beecount * sp->taillen);
 808                 allocate(sp->p, dvector, sp->beecount * sp->taillen);
 809         }
 810
 811         /* Initialize point positions, velocities, etc. */
 812         restart_flow(mi);
 813
 814         /* Set up camera tail */
 815         X(1, 0) = sp->cam[1].x = 0;
 816         Y(1, 0) = sp->cam[1].y = 0;
 817         Z(1, 0) = sp->cam[1].z = 0;
 818 }
 819
 820 ENTRYPOINT void
 821 draw_flow (ModeInfo * mi)
 822 {
 823         int         b, i;
 824         int         col, begin, end;
 825         double      M[3][3]; /* transformation matrix */
 826         flowstruct *sp = NULL;
 827         int         swarm = 0;
 828         int         segindex;
 829
 830         if (flows == NULL)
 831                 return;
 832         sp = &flows[MI_SCREEN(mi)];
 833         if (sp->csegs == NULL)
 834                 return;
 835
 836 #ifdef HAVE_COCOA       /* Don't second-guess Quartz's double-buffering */
 837     XClearWindow (MI_DISPLAY(mi), MI_WINDOW(mi));
 838 #endif
 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 ENTRYPOINT void
1202 reshape_flow(ModeInfo * mi, int width, int height)
1203 {
1204   init_flow (mi);
1205 }
1206
1207
1208 ENTRYPOINT void
1209 release_flow (ModeInfo * mi)
1210 {
1211         if (flows != NULL) {
1212                 int         screen;
1213
1214                 for (screen = 0; screen < MI_NUM_SCREENS(mi); screen++)
1215                         free_flow(&flows[screen]);
1216                 free(flows);
1217                 flows = (flowstruct *) NULL;
1218         }
1219 }
1220
1221 ENTRYPOINT void
1222 refresh_flow (ModeInfo * mi)
1223 {
1224         if(!dbufp) MI_CLEARWINDOW(mi);
1225 }
1226
1227 XSCREENSAVER_MODULE ("Flow", flow)
1228
1229 #endif /* MODE_flow */