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