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

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