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

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