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

   1 /* -*- Mode: C; tab-width: 4 -*- */
   2 /* flow --- flow of strange bees */
   3
   4 #if 0
   5 static const char sccsid[] = "@(#)flow.c 4.10 98/04/24 xlockmore";
   6 #endif
   7
   8 /*-
   9  * Copyright (c) 1996 by Tim Auckland <Tim.Auckland@Sun.COM>
  10  * Portions added by Stephen Davies are Copyright (c) 2000 Stephen Davies
  11  *
  12  * Permission to use, copy, modify, and distribute this software and its
  13  * documentation for any purpose and without fee is hereby granted,
  14  * provided that the above copyright notice appear in all copies and that
  15  * both that copyright notice and this permission notice appear in
  16  * supporting documentation.
  17  *
  18  * This file is provided AS IS with no warranties of any kind.  The author
  19  * shall have no liability with respect to the infringement of copyrights,
  20  * trade secrets or any patents by this file or any part thereof.  In no
  21  * event will the author be liable for any lost revenue or profits or
  22  * other special, indirect and consequential damages.
  23  *
  24  * "flow" shows a variety of continuous phase-space flows around strange
  25  * attractors.  It includes the well-known Lorentz mask (the "Butterfly"
  26  * of chaos fame), two forms of Rossler's "Folded Band" and Poincare'
  27  * sections of the "Birkhoff Bagel" and Duffing's forced occilator.
  28  *
  29  * Revision History:
  30  * 21-Feb-00: Major hackage by Chalky (Stephen Davies, chalky@null.net)
  31  *            Forced perspective mode, added 3d box around attractor which
  32  *            involved coding 3d-planar-clipping against the view-frustrum
  33  *            thingy. Also made view alternate between piggybacking on a 'bee'
  34  *            to zooming around outside the attractor. Most bees slow down and
  35  *            stop, to make the structure of the attractor more obvious.
  36  * 31-Nov-98: [TDA] Added Duffing  (what a strange day that was :) DAB)
  37  *   Duffing's forced oscillator has been added to the formula list and
  38  *   the parameters section has been updated to display it in Poincare'
  39  *   section.
  40  * 30-Nov-98: [TDA] Added travelling perspective option
  41  *   A more exciting point-of-view has been added to all autonomous flows.
  42  *   This views the flow as seen by a particle moving with the flow.  In the
  43  *   metaphor of the original code, I've attached a camera to one of the
  44  *   trained bees!
  45  * 30-Nov-98: [TDA] Much code cleanup.
  46  * 09-Apr-97: [TDA] Ported to xlockmore-4
  47  * 18-Jul-96: Adapted from swarm.c Copyright (c) 1991 by Patrick J. Naughton.
  48  * 31-Aug-90: Adapted from xswarm by Jeff Butterworth. (butterwo@ncsc.org)
  49  */
  50
  51 #ifdef STANDALONE
  52 # define PROGCLASS      "Flow"
  53 # define HACK_INIT      init_flow
  54 # define HACK_DRAW      draw_flow
  55 # define flow_opts      xlockmore_opts
  56 # define DEFAULTS       "*delay:                1000 \n" \
  57                                         "*count:                500 \n" \
  58                                         "*cycles:               3000 \n" \
  59                                         "*ncolors:              200 \n" \
  60         "*rotate:         True \n" \
  61         "*ride:           True \n" \
  62         "*zoom:           True \n" \
  63         "*allow2d:        True \n" \
  64         "*box:            True \n" \
  65         "*slow:           True \n" \
  66         "*freeze:         True \n"
  67 # define SMOOTH_COLORS
  68 # include "xlockmore.h"         /* in xscreensaver distribution */
  69 # include "erase.h"
  70
  71 #else /* STANDALONE */
  72 # include "xlock.h"             /* in xlockmore distribution */
  73 #endif /* STANDALONE */
  74
  75 XrmOptionDescRec flow_options[];
  76 ModeSpecOpt flow_opts = { 7, flow_options, 0, NULL, NULL };
  77
  78 #ifdef USE_MODULES
  79 ModStruct   flow_description = {
  80         "flow", "init_flow", "draw_flow", "release_flow",
  81         "refresh_flow", "init_flow", NULL, &flow_opts,
  82         1000, 1024, 3000, 1, 64, 1.0, "",
  83         "Shows dynamic strange attractors", 0, NULL
  84 };
  85
  86 #endif
  87
  88 typedef struct {
  89         double      x;
  90         double      y;
  91         double      z;
  92 } dvector;
  93
  94 typedef struct {
  95         double      a, b, c;
  96 } Par;
  97
  98 /* Macros */
  99 #define X(t,b)  (sp->p[t][b].x)
 100 #define Y(t,b)  (sp->p[t][b].y)
 101 #define Z(t,b)  (sp->p[t][b].z)
 102 #define balance_rand(v) ((LRAND()/MAXRAND*(v))-((v)/2)) /* random around 0 */
 103 #define SCALE_X(A) (sp->width/2+sp->width/sp->size*(A))
 104 /*#define SCALE_Y(A) (sp->height/2+sp->height/sp->size*(A))*/
 105 #define SCALE_Y(A) (sp->height/2+sp->width/sp->size*(A))
 106
 107 /* Mode of operation. Rotate, ride and zoom are mutually exclusive */
 108 typedef enum {
 109         FLOW_ROTATE = 1, /* Rotate around attractor */
 110         FLOW_RIDE = 2,   /* Ride a trained bee */
 111         FLOW_ZOOM = 4,   /* Zoom in and out */
 112         FLOW_2D = 8,     /* Allow 2D attractors */
 113         FLOW_BOX = 16,    /* Compute a box around the attractor */
 114         FLOW_SLOW = 32,   /* Some bees are slower (and have antifreeze) */
 115         FLOW_FREEZE = 64  /* Freeze some of the bees in action */
 116 } FlowMode;
 117
 118 #define FLOW_DEFAULT (FLOW_ROTATE|FLOW_RIDE|FLOW_ZOOM|FLOW_2D|\
 119                 FLOW_BOX|FLOW_SLOW|FLOW_FREEZE)
 120
 121 typedef struct {
 122         int         width;
 123         int         height;
 124         int         count;
 125         double      size;
 126
 127         int         beecount;   /* number of bees */
 128         XSegment   *csegs;          /* bee lines */
 129         int        *cnsegs;
 130         XSegment   *old_segs;   /* old bee lines */
 131         int         nold_segs;
 132         double      step;
 133         double          slow;
 134         double          slow_view;
 135         dvector     centre;             /* centre */
 136         dvector         range;
 137         struct {
 138                 double  depth;
 139                 double  height;
 140         }           view;
 141         dvector         circle[2]; /* POV that circles around the scene */
 142         dvector    *p[2];   /* bee positions x[time][bee#] */
 143         struct {
 144                 double  theta;
 145                 double  dtheta;
 146                 double  phi;
 147                 double  dphi;
 148         }           tumble;
 149         dvector  (*ODE) (Par par, double x, double y, double z);
 150         Par         par;
 151         FlowMode                mode; /* Mode of operation */
 152 } flowstruct;
 153
 154 static flowstruct *flows = NULL;
 155
 156 static dvector
 157 Lorentz(Par par, double x, double y, double z)
 158 {
 159         dvector d;
 160
 161         d.x = par.a * (y - x);
 162         d.y = x * (par.b - z) - y;
 163         d.z = x * y - par.c * z;
 164         return d;
 165 }
 166
 167 static dvector
 168 Rossler(Par par, double x, double y, double z)
 169 {
 170         dvector d;
 171
 172         d.x = -(y + par.a * z);
 173         d.y = x + y * par.b;
 174         d.z = par.c + z * (x - 5.7);
 175         return d;
 176 }
 177
 178 static dvector
 179 RosslerCone(Par par, double x, double y, double z)
 180 {
 181         dvector d;
 182
 183         d.x = -(y + par.a * z);
 184         d.y = x + y * par.b - z * z * par.c;
 185         d.z = 0.2 + z * (x - 5.7);
 186         return d;
 187 }
 188
 189 static dvector
 190 Birkhoff(Par par, double x, double y, double z)
 191 {
 192         dvector d;
 193
 194         d.x = -y + par.b * sin(z);
 195         d.y = 0.7 * x + par.a * y * (0.1 - x * x);
 196         d.z = par.c;
 197         return d;
 198 }
 199
 200 static dvector
 201 Duffing(Par par, double x, double y, double z)
 202 {
 203         dvector d;
 204
 205         d.x = -par.a * x - y/2 - y * y * y/8 + par.b * cos(z);
 206         d.y = 2*x;
 207         d.z = par.c;
 208         return d;
 209 }
 210
 211 void init_clip(flowstruct *sp);
 212
 213 void
 214 init_flow(ModeInfo * mi)
 215 {
 216         flowstruct *sp;
 217         int         b;
 218         double      beemult = 1;
 219         static int  allocated = 0;
 220
 221         if (flows == NULL) {
 222                 if ((flows = (flowstruct *) calloc(MI_NUM_SCREENS(mi),
 223                                                sizeof (flowstruct))) == NULL)
 224                         return;
 225         }
 226         sp = &flows[MI_SCREEN(mi)];
 227
 228         sp->count = 0;
 229         sp->slow = 0.999;
 230         sp->slow_view = 0.90;
 231
 232         sp->width = MI_WIDTH(mi);
 233         sp->height = MI_HEIGHT(mi);
 234
 235         sp->tumble.theta = balance_rand(M_PI);
 236         sp->tumble.phi = balance_rand(M_PI);
 237         sp->tumble.dtheta = 0.002;
 238         sp->tumble.dphi = 0.001;
 239         sp->view.height = 0;
 240         sp->view.depth = 0; /* no perspective view */
 241         sp->mode = 0;
 242    if (get_boolean_resource ("rotate", "Boolean")) sp->mode |= FLOW_ROTATE;
 243    if (get_boolean_resource ("ride", "Boolean")) sp->mode |= FLOW_RIDE;
 244    if (get_boolean_resource ("zoom", "Boolean")) sp->mode |= FLOW_ZOOM;
 245    if (get_boolean_resource ("allow2d", "Boolean")) sp->mode |= FLOW_2D;
 246    if (get_boolean_resource ("slow", "Boolean")) sp->mode |= FLOW_SLOW;
 247    if (get_boolean_resource ("freeze", "Boolean")) sp->mode |= FLOW_FREEZE;
 248    if (get_boolean_resource ("box", "Boolean")) sp->mode |= FLOW_BOX;
 249
 250         b = (sp->mode & FLOW_2D) ? 5 : 3;
 251         b = NRAND(b);
 252
 253         /* If more than one of rotate, ride and zoom are set, choose one */
 254         if (b < 3) {
 255                 int num = 0, modes[3];
 256
 257                 if (sp->mode & FLOW_ROTATE) modes[num++] = FLOW_ROTATE;
 258                 if (sp->mode & FLOW_RIDE) modes[num++] = FLOW_RIDE;
 259                 if (sp->mode & FLOW_ZOOM) modes[num++] = FLOW_ZOOM;
 260
 261                 sp->mode &= ~(FLOW_ROTATE | FLOW_RIDE | FLOW_ZOOM);
 262
 263                 if (num) sp->mode |= modes[ NRAND(num) ];
 264                 else sp->mode |= FLOW_ZOOM;
 265         }
 266
 267         switch (b) {
 268         case 0:
 269                 sp->view.depth = 10;
 270                 sp->view.height = 0.2;
 271                 beemult = 3;
 272                 sp->ODE = Lorentz;
 273                 sp->step = 0.02;
 274                 sp->size = 60;
 275                 sp->centre.x = 0;
 276                 sp->centre.y = 0;
 277                 sp->centre.z = 24;
 278                 sp->range.x = 5;
 279                 sp->range.y = 5;
 280                 sp->range.z = 1;
 281                 sp->par.a = 10 + balance_rand(5);
 282                 sp->par.b = 28 + balance_rand(5);
 283                 sp->par.c = 2 + balance_rand(1);
 284                 break;
 285         case 1:
 286                 sp->view.depth = 10;
 287                 sp->view.height = 0.1;
 288                 beemult = 4;
 289                 sp->ODE = Rossler;
 290                 sp->step = 0.05;
 291                 sp->size = 24;
 292                 sp->centre.x = 0;
 293                 sp->centre.y = 0;
 294                 sp->centre.z = 3;
 295                 sp->range.x = 4;
 296                 sp->range.y = 4;
 297                 sp->range.z = 7;
 298                 sp->par.a = 2 + balance_rand(1);
 299                 sp->par.b = 0.2 + balance_rand(0.1);
 300                 sp->par.c = 0.2 + balance_rand(0.1);
 301                 break;
 302         case 2:
 303                 sp->view.depth = 10;
 304                 sp->view.height = 0.1;
 305                 beemult = 3;
 306                 sp->ODE = RosslerCone;
 307                 sp->step = 0.05;
 308                 sp->size = 24;
 309                 sp->centre.x = 0;
 310                 sp->centre.y = 0;
 311                 sp->centre.z = 3;
 312                 sp->range.x = 4;
 313                 sp->range.y = 4;
 314                 sp->range.z = 4;
 315                 sp->par.a = 2;
 316                 sp->par.b = 0.2;
 317                 sp->par.c = 0.25 + balance_rand(0.09);
 318                 break;
 319         case 3:
 320                 sp->ODE = Birkhoff;
 321                 sp->step = 0.04;
 322                 sp->size = 2.6;
 323                 sp->centre.x = 0;
 324                 sp->centre.y = 0;
 325                 sp->centre.z = 0;
 326                 sp->range.x = 3;
 327                 sp->range.y = 4;
 328                 sp->range.z = 0;
 329                 sp->par.a = 10 + balance_rand(5);
 330                 sp->par.b = 0.35 + balance_rand(0.25);
 331                 sp->par.c = 1.57;
 332                 sp->tumble.theta = 0;
 333                 sp->tumble.phi = 0;
 334                 sp->tumble.dtheta = 0;
 335                 sp->tumble.dphi = 0;
 336                 break;
 337         case 4:
 338         default:
 339                 sp->ODE = Duffing;
 340                 sp->step = 0.02;
 341                 sp->size = 30;
 342                 sp->centre.x = 0;
 343                 sp->centre.y = 0;
 344                 sp->centre.z = 0;
 345                 sp->range.x = 20;
 346                 sp->range.y = 20;
 347                 sp->range.z = 0;
 348                 sp->par.a = 0.2 + balance_rand(0.1);
 349                 sp->par.b = 27.0 + balance_rand(3.0);
 350                 sp->par.c = 1.33;
 351                 sp->tumble.theta = 0;
 352                 sp->tumble.phi = 0;
 353                 sp->tumble.dtheta = -NRAND(2)*sp->par.c*sp->step;
 354                 sp->tumble.dphi = 0;
 355                 beemult = 0.5;
 356                 break;
 357         }
 358
 359         sp->view.depth *= 4;
 360
 361         sp->beecount = beemult * MI_COUNT(mi);
 362         if (sp->beecount < 0)   /* random variations */
 363                 sp->beecount = NRAND(-sp->beecount) + 1; /* Minimum 1 */
 364
 365         /* Clear the background. */
 366         MI_CLEARWINDOW(mi);
 367
 368         if(!allocated || sp->beecount != allocated){ /* reallocate */
 369                 if (sp->csegs != NULL) {
 370                         (void) free((void *) sp->csegs);
 371                         sp->csegs = NULL;
 372                 }
 373                 if (sp->cnsegs != NULL) {
 374                         (void) free((void *) sp->cnsegs);
 375                         sp->cnsegs = NULL;
 376                 }
 377                 if (sp->old_segs != NULL) {
 378                         (void) free((void *) sp->old_segs);
 379                         sp->old_segs = NULL;
 380                 }
 381                 if (sp->p[0] != NULL) {
 382                         (void) free((void *) sp->p[0]);
 383                         sp->p[0] = NULL;
 384                 }
 385                 if (sp->p[1] != NULL) {
 386                         (void) free((void *) sp->p[1]);
 387                         sp->p[1] = NULL;
 388                 }
 389                 allocated = sp->beecount;
 390         }
 391
 392         /* Allocate memory. */
 393
 394         if (!sp->csegs) {
 395                 sp->csegs = (XSegment *) malloc(sizeof (XSegment) * sp->beecount
 396                                                 * MI_NPIXELS(mi));
 397                 sp->cnsegs = (int *) malloc(sizeof (int) * MI_NPIXELS(mi));
 398
 399                 sp->old_segs = (XSegment *) malloc(sizeof (XSegment) * sp->beecount);
 400                 sp->p[0] = (dvector *) malloc(sizeof (dvector) * sp->beecount);
 401                 sp->p[1] = (dvector *) malloc(sizeof (dvector) * sp->beecount);
 402         }
 403
 404         /* Initialize point positions, velocities, etc. */
 405
 406         for (b = 0; b < sp->beecount; b++) {
 407                 X(1, b) = X(0, b) = balance_rand(sp->range.x);
 408                 Y(1, b) = Y(0, b) = balance_rand(sp->range.y);
 409                 Z(1, b) = Z(0, b) = balance_rand(sp->range.z);
 410         }
 411
 412         init_clip(sp);
 413
 414 }
 415
 416 /* Clipping planes */
 417 #define PLANES 5
 418 static double plane_orig[][2][3] = {
 419         /* X goes into screen, Y goes right, Z goes down(up?) */
 420         /* {Normal}, {Point} */
 421         { {1.0, 0, 0}, {0.01, 0, 0} },
 422         { {1.0, 1.0, 0.0}, {0, 0, 0} },
 423         { {1.0,-1.0, 0.0}, {0, 0, 0} },
 424         { {1.0, 0.0, 1.0}, {0, 0, 0} },
 425         { {1.0, 0.0,-1.0}, {0, 0, 0} }
 426 };
 427 static double plane[PLANES][2][3];
 428 static double plane_d[PLANES];
 429
 430 #define BOX_P 32
 431 #define BOX_L 36
 432 #define MIN_BOX (3)
 433 #define MAX_BOX (MIN_BOX + BOX_L)
 434 /* Points that make up the box (normalized coordinates) */
 435 static double box_orig[][3] = {
 436         {1,1,1},   /* 0 */
 437         {1,1,-1},  /* 1 */
 438         {1,-1,-1}, /* 2 */
 439         {1,-1,1},  /* 3 */
 440         {-1,1,1},  /* 4 */
 441         {-1,1,-1}, /* 5 */
 442         {-1,-1,-1},/* 6 */
 443         {-1,-1,1}, /* 7 */
 444         {1, .8, .8},
 445         {1, .8,-.8},
 446         {1,-.8,-.8},
 447         {1,-.8, .8},
 448         { .8,1, .8},
 449         { .8,1,-.8},
 450         {-.8,1,-.8},
 451         {-.8,1, .8},
 452         { .8, .8,1},
 453         { .8,-.8,1},
 454         {-.8,-.8,1},
 455         {-.8, .8,1},
 456         {-1, .8, .8},
 457         {-1, .8,-.8},
 458         {-1,-.8,-.8},
 459         {-1,-.8, .8},
 460         { .8,-1, .8},
 461         { .8,-1,-.8},
 462         {-.8,-1,-.8},
 463         {-.8,-1, .8},
 464         { .8, .8,-1},
 465         { .8,-.8,-1},
 466         {-.8,-.8,-1},
 467         {-.8, .8,-1}
 468 };
 469
 470 /* Container for scaled box points */
 471 static double box[BOX_P][3];
 472
 473 /* Lines connecting the box dots */
 474 static double lines[][2] = {
 475         {0,1}, {1,2}, {2,3}, {3,0}, /* box */
 476         {4,5}, {5,6}, {6,7}, {7,4},
 477         {0,4}, {1,5}, {2,6}, {3,7},
 478         {4+4,5+4}, {5+4,6+4}, {6+4,7+4}, {7+4,4+4},
 479         {4+8,5+8}, {5+8,6+8}, {6+8,7+8}, {7+8,4+8},
 480         {4+12,5+12}, {5+12,6+12}, {6+12,7+12}, {7+12,4+12},
 481         {4+16,5+16}, {5+16,6+16}, {6+16,7+16}, {7+16,4+16},
 482         {4+20,5+20}, {5+20,6+20}, {6+20,7+20}, {7+20,4+20},
 483         {4+24,5+24}, {5+24,6+24}, {6+24,7+24}, {7+24,4+24},
 484 };
 485
 486 /* Boundaries of bees */
 487 double xmin, xmax;
 488 double ymin, ymax;
 489 double zmin, zmax;
 490
 491 void init_clip(flowstruct *sp)
 492 {
 493         int i;
 494
 495         /* Scale the planes to the screen. I had to invert the projection
 496          * algorithms so that when projected they would be right at the edge of the
 497          * screen. */
 498
 499     /* #### jwz: I'm not really sure what it means when sp->view.depth is 0
 500        in here -- what's the right thing to do? */
 501
 502         double width = (sp->view.depth
 503                     ? sp->size/sp->view.depth/2
 504                     : 1);
 505         double height = (sp->view.depth
 506                      ? (sp->size/sp->view.depth/2*
 507                         sp->view.height/sp->view.height)
 508                      : 1);
 509         for (i = 0; i < PLANES; i++) {
 510                 /* Copy orig planes into planes, expanding <-> clippings */
 511                 plane[i][0][0] = plane_orig[i][0][0];
 512                 plane[i][0][1] = plane_orig[i][0][1] / width;
 513                 plane[i][0][2] = plane_orig[i][0][2] / height;
 514                 plane[i][1][0] = plane_orig[i][1][0];
 515                 plane[i][1][1] = plane_orig[i][1][1];
 516                 plane[i][1][2] = plane_orig[i][1][2];
 517
 518                 /* Calculate the 'd' part of 'ax + by + cz = d' */
 519                 plane_d[i] = - plane[i][0][0] * plane[i][1][0];
 520                 plane_d[i] -= plane[i][0][1] * plane[i][1][1];
 521                 plane_d[i] -= plane[i][0][2] * plane[i][1][2];
 522         }
 523         xmin = X(0, i); xmax = X(0,i);
 524         ymin = Y(0, i); ymax = Y(0,i);
 525         zmin = Z(0, i); zmax = Z(0,i);
 526 }
 527
 528 /* Scale the box defined above to fit around all points */
 529 void create_box(flowstruct *sp)
 530 {
 531         int i = MAX_BOX;
 532         double xmid, ymid, zmid;
 533         double xsize, ysize, zsize;
 534         double size;
 535
 536         /* Count every 5th point for speed.. */
 537         for (; i < sp->beecount; i += 5) {
 538                 if ( X(0,i) < xmin ) xmin = X(0, i);
 539                 else if ( X(0,i) > xmax ) xmax = X(0, i);
 540                 if ( Y(0,i) < ymin ) ymin = Y(0, i);
 541                 else if ( Y(0,i) > ymax ) ymax = Y(0, i);
 542                 if ( Z(0,i) < zmin ) zmin = Z(0, i);
 543                 else if ( Z(0,i) > zmax ) zmax = Z(0, i);
 544         }
 545         xmid = (xmax+xmin)/2;
 546         ymid = (ymax+ymin)/2;
 547         zmid = (zmax+zmin)/2;
 548         xsize = xmax - xmin;
 549         ysize = ymax - ymin;
 550         zsize = zmax - zmin;
 551         size = xsize;
 552         if (ysize> size) size = ysize;
 553         if (zsize > size) size = zsize;
 554         size /= 2;
 555
 556         /* Scale box */
 557         for (i = 0; i < BOX_P; i++) {
 558                 box[i][0] = box_orig[i][0] * size + xmid;
 559                 box[i][1] = box_orig[i][1] * size + ymid;
 560                 box[i][2] = box_orig[i][2] * size + zmid;
 561         }
 562
 563 }
 564
 565 /* Returns true if point is infront of the plane (rather than behind) */
 566 int infront_of(double x, double y, double z, int i)
 567 {
 568         double sum = plane[i][0][0]*x + plane[i][0][1]*y + plane[i][0][2]*z + plane_d[i];
 569         return sum >= 0.0;
 570 }
 571
 572 /* Returns true if line was behind a clip plane, or clips the line */
 573 int clip(double *x1, double *y1, double *z1, double *x2, double *y2, double *z2)
 574 {
 575         int i;
 576         for (i = 0; i < PLANES; i++) {
 577                 double t;
 578                 double x, y, z; /* Intersection point */
 579                 double dx, dy, dz; /* line delta */
 580                 int front1, front2;
 581                 front1 = infront_of(*x1, *y1, *z1, i);
 582                 front2 = infront_of(*x2, *y2, *z2, i);
 583                 if (!front1 && !front2) return 1;
 584                 if (front1 && front2) continue;
 585
 586                 dx = *x2 - *x1;
 587                 dy = *y2 - *y1;
 588                 dz = *z2 - *z1;
 589
 590                 /* Find t in line equation */
 591                 t = ( plane_d[i] -
 592                                 plane[i][0][0]*(*x1) - plane[i][0][1]*(*y1) - plane[i][0][2]*(*z1) )
 593                                 /
 594                         ( plane[i][0][0]*dx + plane[i][0][1]*dy + plane[i][0][2]*dz );
 595
 596                 x = *x1 + dx * t;
 597                 y = *y1 + dy * t;
 598                 z = *z1 + dz * t;
 599                 /* Make point that was behind to be the intersect */
 600                 if (front2) {
 601                         *x1 = x;
 602                         *y1 = y;
 603                         *z1 = z;
 604                 } else {
 605                         *x2 = x;
 606                         *y2 = y;
 607                         *z2 = z;
 608                 }
 609         }
 610         return 0;
 611 }
 612
 613
 614 void
 615 draw_flow(ModeInfo * mi)
 616 {
 617         Display    *display = MI_DISPLAY(mi);
 618         Window      window = MI_WINDOW(mi);
 619         GC          gc = MI_GC(mi);
 620         flowstruct *sp = &flows[MI_SCREEN(mi)];
 621         int         b, c, i;
 622         int         col, ix;
 623         int                     new_view = 0;
 624         double      M[3][3]; /* transformation matrix */
 625         double          step_view = sp->step;
 626         double          step_bees = sp->step;
 627         double          step_slow = sp->step;
 628         double          pp, pc;
 629
 630         create_box(sp);
 631
 632         if(!sp->view.depth){ /* simple 3D tumble */
 633                 double      sint, cost, sinp, cosp;
 634                 sp->tumble.theta += sp->tumble.dtheta;
 635                 sp->tumble.phi += sp->tumble.dphi;
 636                 sint = sin(sp->tumble.theta);
 637                 cost = cos(sp->tumble.theta);
 638                 sinp = sin(sp->tumble.phi);
 639                 cosp = cos(sp->tumble.phi);
 640                 M[0][0]= cost; M[0][1]=-sint*cosp; M[0][2]= sint*sinp;
 641                 M[1][0]= sint; M[1][1]= cost*cosp; M[1][2]=-cost*sinp;
 642                 M[2][0]= 0;    M[2][1]= 0;         M[2][2]= 1;
 643         } else { /* initialize matrix */
 644                 M[0][0]= 0; M[0][1]= 0; M[0][2]= 0;
 645                 M[1][0]= 0; M[1][1]= 0; M[1][2]= 0;
 646                 M[2][0]= 0; M[2][1]= 0; M[2][2]= 0;
 647
 648         }
 649
 650         for (col = 0; col < MI_NPIXELS(mi); col++)
 651                 sp->cnsegs[col] = 0;
 652
 653         MI_IS_DRAWN(mi) = True;
 654
 655         /* Calculate circling POV */
 656         sp->circle[1] = sp->circle[0];
 657         sp->circle[0].x = sp->size * 2 * sin(sp->count / 40.0) * (0.6 + 0.4 *cos(sp->count / 100.0));
 658         sp->circle[0].y = sp->size * 2 * cos(sp->count / 40.0) * (0.6 + 0.4 *cos(sp->count / 100.0));
 659         sp->circle[0].z = sp->size * 2 * sin(sp->count / 421.0);
 660
 661         if (sp->mode & FLOW_ROTATE)
 662                 pp = 0;
 663         else if (sp->mode & FLOW_RIDE)
 664                 pp = 1;
 665         else /* ZOOM */
 666                 /* Bistable oscillator to switch between the trained bee and the circler */
 667                 pp = -sin(sin(sin(cos(sp->count / 150.0)*M_PI/2)*M_PI/2)*M_PI/2) *0.5 + 0.5;
 668         pc = 1 - pp;
 669
 670
 671         /* Slow down or speed up the bees / view: */
 672
 673         /* exponentially accelerate towards zero */
 674         sp->slow = sp->slow * 1.005 - 0.005;
 675         if (sp->slow < 0) sp->slow = 0;
 676
 677         sp->slow_view = sp->slow_view * 1.005 - 0.005;
 678         if (sp->slow_view < 0) sp->slow_view = 0;
 679
 680         /* View speeds up, slow bees slow to half speed, and other bees will
 681          * actually stop */
 682         step_view = step_view * (1.01 - sp->slow_view * sp->slow_view) * 0.2;
 683         step_slow = step_slow * (sp->slow + 0.5) / 2;
 684         if (sp->mode & FLOW_SLOW)
 685                 step_bees = step_bees * sp->slow;
 686         else
 687                 step_bees = step_slow;
 688
 689         /* <=- Bees -=> */
 690         for (b = 0; b < sp->beecount; b++) {
 691                 /* Calc if this bee is slow. Note normal bees are exempt from
 692                  * calculations once they slow to half speed, so that they remain as
 693                  * frozen lines rather than barely-visible points */
 694                 int slow = ((b & 0x7) == 0);
 695                 if ( !(sp->mode & FLOW_FREEZE) ) slow = 1;
 696                 /* Age the arrays. */
 697                 if (b < 2 || sp->slow > 0.5 || slow) {
 698                         X(1, b) = X(0, b);
 699                         Y(1, b) = Y(0, b);
 700                         Z(1, b) = Z(0, b);
 701
 702                         /* 2nd order Kunge Kutta */
 703                         {
 704                                 dvector     k1, k2;
 705                                 double          step;
 706
 707                                 if (b == 0 || b == 1) {
 708                                         step = step_view;
 709                                 } else if (slow) {
 710                                         step = step_slow;
 711                                 } else {
 712                                         step = step_bees;
 713                                 }
 714                                 k1 = sp->ODE(sp->par, X(1, b), Y(1, b), Z(1, b));
 715                                 k1.x *= step;
 716                                 k1.y *= step;
 717                                 k1.z *= step;
 718                                 k2 = sp->ODE(sp->par, X(1, b) + k1.x, Y(1, b) + k1.y, Z(1, b) + k1.z);
 719                                 k2.x *= step;
 720                                 k2.y *= step;
 721                                 k2.z *= step;
 722                                 X(0, b) = X(1, b) + (k1.x + k2.x) / 2.0;
 723                                 Y(0, b) = Y(1, b) + (k1.y + k2.y) / 2.0;
 724                                 Z(0, b) = Z(1, b) + (k1.z + k2.z) / 2.0;
 725                         }
 726                 }
 727
 728
 729                 /* Colour according to bee */
 730                 col = b % (MI_NPIXELS(mi) - 1);
 731                 ix = col * sp->beecount + sp->cnsegs[col];
 732
 733                 /* Fill the segment lists. */
 734
 735                 if(sp->view.depth) { /* perspective view has special points */
 736                         if(b==0){ /* point of view */
 737                                 sp->centre.x = X(0, b) * pp + sp->circle[0].x * pc;
 738                                 sp->centre.y = Y(0, b) * pp + sp->circle[0].y * pc;
 739                                 sp->centre.z = Z(0, b) * pp + sp->circle[0].z * pc;
 740                                 /*printf("center: (%3.3f,%3.3f,%3.3f)\n",sp->centre.x, sp->centre.y, sp->centre.z);*/
 741                         }else if(b==1){ /* neighbour: used to compute local axes */
 742                                 double x[3], p[3], x2=0, xp=0, C[3][3];
 743                                 int j;
 744
 745                                 /* forward */
 746                                 x[0] = X(0, 0) - X(1, 0);
 747                                 x[1] = Y(0, 0) - Y(1, 0);
 748                                 x[2] = Z(0, 0) - Z(1, 0);
 749
 750                                 /* neighbour */
 751                                 p[0] = X(0, 1) - X(1, 0);
 752                                 p[1] = Y(0, 1) - Y(1, 0);
 753                                 p[2] = Z(0, 1) - Z(1, 0);
 754
 755                                 for(i=0; i<3; i++){
 756                                         x2+= x[i]*x[i];    /* X . X */
 757                                         xp+= x[i]*p[i];    /* X . P */
 758                                         M[0][i] = x[i];    /* X */
 759                                 }
 760
 761                                 for(i=0; i<3; i++)               /* (X x P) x X */
 762                                         M[1][i] = x2*p[i] - xp*x[i]; /* == (X . X) P - (X . P) X */
 763
 764                                 M[2][0] =  x[1]*p[2] - x[2]*p[1]; /* X x P */
 765                                 M[2][1] = -x[0]*p[2] + x[2]*p[0];
 766                                 M[2][2] =  x[0]*p[1] - x[1]*p[0];
 767
 768                                 /* normalise axes */
 769                                 for(j=0; j<3; j++){
 770                                         double A=0;
 771                                         for(i=0; i<3; i++) A+=M[j][i]*M[j][i]; /* sum squares */
 772                                         A=sqrt(A);
 773                                         for(i=0; i<3; i++) M[j][i]/=A;
 774                                 }
 775
 776                                 X(0, 1)=X(0, 0)+M[1][0]; /* adjust neighbour */
 777                                 Y(0, 1)=Y(0, 0)+M[1][1];
 778                                 Z(0, 1)=Z(0, 0)+M[1][2];
 779
 780                                 /* Look at trained bee into C matrix */
 781                                 /* forward */
 782                                 x[0] = 0 - sp->circle[0].x;
 783                                 x[1] = 0 - sp->circle[0].y;
 784                                 x[2] = 0 - sp->circle[0].z;
 785
 786                                 /* neighbour */
 787                                 p[0] = sp->circle[0].x - sp->circle[1].x;
 788                                 p[1] = sp->circle[0].y - sp->circle[1].y;
 789                                 p[2] = sp->circle[0].z - sp->circle[1].z;
 790
 791                                 for(i=0; i<3; i++){
 792                                         x2+= x[i]*x[i];    /* X . X */
 793                                         xp+= x[i]*p[i];    /* X . P */
 794                                         C[0][i] = x[i];    /* X */
 795                                 }
 796
 797                                 for(i=0; i<3; i++)               /* (X x P) x X */
 798                                         C[1][i] = x2*p[i] - xp*x[i]; /* == (X . X) P - (X . P) X */
 799
 800                                 C[2][0] =  x[1]*p[2] - x[2]*p[1]; /* X x P */
 801                                 C[2][1] = -x[0]*p[2] + x[2]*p[0];
 802                                 C[2][2] =  x[0]*p[1] - x[1]*p[0];
 803
 804                                 /* normalise axes */
 805                                 for(j=0; j<3; j++){
 806                                         double A=0;
 807                                         for(i=0; i<3; i++) A+=C[j][i]*C[j][i]; /* sum squares */
 808                                         A=sqrt(A);
 809                     if (A != 0) /* #### is this right? */
 810                       for(i=0; i<3; i++) C[j][i]/=A;
 811                                 }
 812
 813                                 /* Interpolate between Center and Trained Bee matrices */
 814                                 /* This isn't very accurate and leads to weird transformations
 815                                  * (shearing, etc), but it works. Besides, sometimes they look
 816                                  * cool :) */
 817                                 pp = pp * pp; /* Don't follow bee's direction until very close */
 818                                 pc = 1 - pp;
 819                                 for (i = 0; i < 3; i++)
 820                                         for (j = 0; j < 3; j++)
 821                                                 M[i][j] = M[i][j] * pp + C[i][j] * pc;
 822
 823
 824 #if 0  /* display local axes for testing */
 825                                 X(1, b)=X(0, 0);
 826                                 Y(1, b)=Y(0, 0);
 827                                 Z(1, b)=Z(0, 0);
 828                         }else if(b==2){
 829                                 X(0, b)=X(0, 0)+0.5*M[0][0];
 830                                 Y(0, b)=Y(0, 0)+0.5*M[0][1];
 831                                 Z(0, b)=Z(0, 0)+0.5*M[0][2];
 832                                 X(1, b)=X(0, 0);
 833                                 Y(1, b)=Y(0, 0);
 834                                 Z(1, b)=Z(0, 0);
 835                         }else if(b==3){
 836                                 X(0, b)=X(0, 0)+1.5*M[2][0];
 837                                 Y(0, b)=Y(0, 0)+1.5*M[2][1];
 838                                 Z(0, b)=Z(0, 0)+1.5*M[2][2];
 839                                 X(1, b)=X(0, 0);
 840                                 Y(1, b)=Y(0, 0);
 841                                 Z(1, b)=Z(0, 0);
 842 #endif
 843                         /* Draw a box... */
 844                         }
 845                 }
 846                         if (b >= MIN_BOX && b < MAX_BOX) {
 847                                 int p1 = lines[b-MIN_BOX][0];
 848                                 int p2 = lines[b-MIN_BOX][1];
 849                                 X(0, b) = box[p1][0]; Y(0, b) = box[p1][1]; Z(0, b) = box[p1][2];
 850                                 X(1, b) = box[p2][0]; Y(1, b) = box[p2][1]; Z(1, b) = box[p2][2];
 851                         }
 852
 853
 854 #if 0   /* Original code */
 855                 for(i=0; i<2; i++){
 856                         double x=X(i,b)-sp->centre.x;
 857                         double y=Y(i,b)-sp->centre.y;
 858                         double z=Z(i,b)-sp->centre.z;
 859                         double X=M[0][0]*x + M[0][1]*y + M[0][2]*z;
 860                         double Y=M[1][0]*x + M[1][1]*y + M[1][2]*z;
 861                         double Z=M[2][0]*x + M[2][1]*y + M[2][2]*z+sp->view.height;
 862                         double absx, absy;
 863                         if(sp->view.depth){
 864                                 if(X <= 0) break;
 865                                 absx=SCALE_X(sp->view.depth*Y/X);
 866                                 absy=SCALE_Y(sp->view.depth*Z/X);
 867                                 if(absx < -sp->width || absx > 2*sp->width ||
 868                                    absy < -sp->height || absy > 2*sp->height)
 869                                         break;
 870                         }else{
 871                                 absx=SCALE_X(X);
 872                                 absy=SCALE_Y(Y);
 873                         }
 874                         if(i){
 875                                 sp->csegs[ix].x1 = (short) absx;
 876                                 sp->csegs[ix].y1 = (short) absy;
 877                         }else{
 878                                 sp->csegs[ix].x2 = (short) absx;
 879                                 sp->csegs[ix].y2 = (short) absy;
 880                         }
 881                 }
 882                 if(i == 2) /* both assigned */
 883                         sp->cnsegs[col]++;
 884 #else
 885                 /* Chalky's code w/ clipping */
 886                 if (b < ((sp->mode & FLOW_BOX) ? 2 : MAX_BOX))
 887                         continue;
 888                 do {
 889                         double x1=X(0,b)-sp->centre.x;
 890                         double y1=Y(0,b)-sp->centre.y;
 891                         double z1=Z(0,b)-sp->centre.z;
 892                         double X1=M[0][0]*x1 + M[0][1]*y1 + M[0][2]*z1;
 893                         double Y1=M[1][0]*x1 + M[1][1]*y1 + M[1][2]*z1;
 894                         double Z1=M[2][0]*x1 + M[2][1]*y1 + M[2][2]*z1+sp->view.height;
 895                         double absx1, absy1;
 896                         double x2=X(1,b)-sp->centre.x;
 897                         double y2=Y(1,b)-sp->centre.y;
 898                         double z2=Z(1,b)-sp->centre.z;
 899                         double X2=M[0][0]*x2 + M[0][1]*y2 + M[0][2]*z2;
 900                         double Y2=M[1][0]*x2 + M[1][1]*y2 + M[1][2]*z2;
 901                         double Z2=M[2][0]*x2 + M[2][1]*y2 + M[2][2]*z2+sp->view.height;
 902                         double absx2, absy2;
 903                         if(sp->view.depth){
 904                                 /* Need clipping if: is part of box, or close to viewer */
 905                                 if ( (b >= MIN_BOX && b < MAX_BOX) || X1 <= 0.1 || X2 < 0.1)
 906                                         if (clip(&X1, &Y1, &Z1, &X2, &Y2, &Z2))
 907                                                 break;
 908                                 if (X1 <= 0 || X2 <= 0) break;
 909                                 absx1=SCALE_X(sp->view.depth*Y1/X1);
 910                                 absy1=SCALE_Y(sp->view.depth*Z1/X1);
 911                                 if(absx1 < -sp->width || absx1 > 2*sp->width ||
 912                                    absy1 < -sp->height || absy1 > 2*sp->height)
 913                                         break;
 914                                 absx2=SCALE_X(sp->view.depth*Y2/X2);
 915                                 absy2=SCALE_Y(sp->view.depth*Z2/X2);
 916                                 if(absx2 < -sp->width || absx2 > 2*sp->width ||
 917                                    absy2 < -sp->height || absy2 > 2*sp->height)
 918                                         break;
 919                         }else{
 920                                 absx1=SCALE_X(X1);
 921                                 absy1=SCALE_Y(Y1);
 922                                 absx2=SCALE_X(X2);
 923                                 absy2=SCALE_Y(Y2);
 924                         }
 925
 926                         sp->csegs[ix].x1 = (short) absx1;
 927                         sp->csegs[ix].y1 = (short) absy1;
 928                         sp->csegs[ix].x2 = (short) absx2;
 929                         sp->csegs[ix].y2 = (short) absy2;
 930
 931                         sp->cnsegs[col]++;
 932                 } while (0);
 933 #endif
 934         }
 935
 936         if (sp->count) { /* erase */
 937                 XSetForeground(display, gc, MI_BLACK_PIXEL(mi));
 938                 XDrawSegments(display, window, gc, sp->old_segs, sp->nold_segs);
 939         }
 940
 941         if (MI_NPIXELS(mi) > 2){ /* render colour */
 942                 for (col = 0; col < MI_NPIXELS(mi); col++)
 943                         if (sp->cnsegs[col] > 0) {
 944                                 XSetForeground(display, gc, MI_PIXEL(mi, col));
 945                                 XDrawSegments(display, window, gc,
 946                                               sp->csegs + col * sp->beecount, sp->cnsegs[col]);
 947                         }
 948         } else {                /* render mono */
 949                 XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
 950                 XDrawSegments(display, window, gc,
 951                                           sp->csegs + col * sp->beecount, sp->cnsegs[col]);
 952         }
 953
 954         /* Copy to erase-list */
 955         for (col = 0, c = 0; col < MI_NPIXELS(mi); col++)
 956                 for (b = 0; b < sp->cnsegs[col]; b++)
 957                         sp->old_segs[c++] = (sp->csegs + col * sp->beecount)[b];
 958         sp->nold_segs = c;
 959
 960         if (++sp->count > MI_CYCLES(mi)) /* pick a new flow */
 961                 init_flow(mi);
 962
 963         if (sp->count % (MI_CYCLES(mi)/4) == 0) { /* pick a new view */
 964                 new_view = 0; /* change to 1 .. */
 965         }
 966
 967         if (X(0, 0) < xmin*2 || X(0, 0) > xmax*2) new_view = 1;
 968         if (Y(0, 0) < ymin*2 || Y(0, 0) > ymax*2) new_view = 1;
 969         if (Z(0, 0) < zmin*2 || Z(0, 0) > zmax*2) new_view = 1;
 970
 971         if (new_view) {
 972                 for (b = 0; b < 2; b++) {
 973                         X(1, b) = X(0, b) = balance_rand(sp->range.x*4);
 974                         Y(1, b) = Y(0, b) = balance_rand(sp->range.y*4);
 975                         Z(1, b) = Z(0, b) = balance_rand(sp->range.z*4);
 976                 }
 977                 sp->slow_view = 0.90;
 978         }
 979 }
 980
 981 void
 982 release_flow(ModeInfo * mi)
 983 {
 984         if (flows != NULL) {
 985                 int         screen;
 986
 987                 for (screen = 0; screen < MI_NUM_SCREENS(mi); screen++) {
 988                         flowstruct *sp = &flows[screen];
 989
 990                         if (sp->csegs != NULL)
 991                                 (void) free((void *) sp->csegs);
 992                         if (sp->cnsegs != NULL)
 993                                 (void) free((void *) sp->cnsegs);
 994                         if (sp->old_segs != NULL)
 995                                 (void) free((void *) sp->old_segs);
 996                         if (sp->p[0] != NULL)
 997                                 (void) free((void *) sp->p[0]);
 998                         if (sp->p[1] != NULL)
 999                                 (void) free((void *) sp->p[1]);
1000                 }
1001                 (void) free((void *) flows);
1002                 flows = NULL;
1003         }
1004 }
1005
1006 void
1007 refresh_flow(ModeInfo * mi)
1008 {
1009         MI_CLEARWINDOW(mi);
1010 }
1011
1012 XrmOptionDescRec flow_options[] =
1013 {
1014         {"-rotate",  ".rotate", XrmoptionSepArg, 0},
1015         {"-ride",  ".ride", XrmoptionSepArg, 0},
1016         {"-zoom",  ".zoom", XrmoptionSepArg, 0},
1017         {"-box",  ".box", XrmoptionSepArg, 0},
1018         {"-slow",  ".slow", XrmoptionSepArg, 0},
1019         {"-freeze",  ".freeze", XrmoptionSepArg, 0},
1020         {"-allow2d",  ".allow2d", XrmoptionSepArg, 0},
1021   { 0, 0, 0, 0 }
1022 };
1023
1024 /*
1025 char*        defaults[] =
1026 {
1027         "*rotate:         True",
1028         "*ride:           True",
1029         "*zoom:           True",
1030         "*allow2d:        True",
1031         "*box:            True",
1032         "*slow:           True",
1033         "*freeze:         True",
1034   0
1035 };
1036         */