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

   1 /* -*- Mode: C; tab-width: 4 -*- */
   2 /* flow --- flow of strange bees */
   3
   4 #if !defined( lint ) && !defined( SABER )
   5 static const char sccsid[] = "@(#)flow.c 4.10 98/04/24 xlockmore";
   6
   7 #endif
   8
   9 /*-
  10  * Copyright (c) 1996 by Tim Auckland <Tim.Auckland@Sun.COM>
  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  * 31-Nov-98: [TDA] Added Duffing  (what a strange day that was :) DAB)
  31  *   Duffing's forced oscillator has been added to the formula list and
  32  *   the parameters section has been updated to display it in Poincare'
  33  *   section.
  34  * 30-Nov-98: [TDA] Added travelling perspective option
  35  *   A more exciting point-of-view has been added to all autonomous flows.
  36  *   This views the flow as seen by a particle moving with the flow.  In the
  37  *   metaphor of the original code, I've attached a camera to one of the
  38  *   trained bees!
  39  * 30-Nov-98: [TDA] Much code cleanup.
  40  * 09-Apr-97: [TDA] Ported to xlockmore-4
  41  * 18-Jul-96: Adapted from swarm.c Copyright (c) 1991 by Patrick J. Naughton.
  42  * 31-Aug-90: Adapted from xswarm by Jeff Butterworth. (butterwo@ncsc.org)
  43  */
  44
  45 #ifdef STANDALONE
  46 # define PROGCLASS      "Flow"
  47 # define HACK_INIT      init_flow
  48 # define HACK_DRAW      draw_flow
  49 # define flow_opts      xlockmore_opts
  50 # define DEFAULTS       "*delay:                1000 \n" \
  51                                         "*count:                1024 \n" \
  52                                         "*cycles:               3000 \n" \
  53                                         "*ncolors:              200 \n"
  54 # define SMOOTH_COLORS
  55 # include "xlockmore.h"         /* in xscreensaver distribution */
  56 # include "erase.h"
  57
  58 #else /* STANDALONE */
  59 # include "xlock.h"             /* in xlockmore distribution */
  60 #endif /* STANDALONE */
  61
  62 ModeSpecOpt flow_opts = { 0, NULL, 0, NULL, NULL };
  63
  64 #ifdef USE_MODULES
  65 ModStruct   flow_description = {
  66         "flow", "init_flow", "draw_flow", "release_flow",
  67         "refresh_flow", "init_flow", NULL, &flow_opts,
  68         1000, 1024, 3000, 1, 64, 1.0, "",
  69         "Shows dynamic strange attractors", 0, NULL
  70 };
  71
  72 #endif
  73
  74 typedef struct {
  75         double      x;
  76         double      y;
  77         double      z;
  78 } dvector;
  79
  80 typedef struct {
  81         double      a, b, c;
  82 } Par;
  83
  84 /* Macros */
  85 #define X(t,b)  (sp->p[t][b].x)
  86 #define Y(t,b)  (sp->p[t][b].y)
  87 #define Z(t,b)  (sp->p[t][b].z)
  88 #define balance_rand(v) ((LRAND()/MAXRAND*(v))-((v)/2)) /* random around 0 */
  89 #define SCALE_X(A) (sp->width/2+sp->width/sp->size*(A))
  90 #define SCALE_Y(A) (sp->height/2+sp->height/sp->size*(A))
  91
  92 typedef struct {
  93         int         width;
  94         int         height;
  95         int         count;
  96         double      size;
  97
  98         int         beecount;   /* number of bees */
  99         XSegment   *csegs;          /* bee lines */
 100         int        *cnsegs;
 101         XSegment   *old_segs;   /* old bee lines */
 102         int         nold_segs;
 103         double      step;
 104         dvector     centre;             /* centre */
 105         struct {
 106                 double  depth;
 107                 double  height;
 108         }           view;
 109         dvector    *p[2];   /* bee positions x[time][bee#] */
 110         struct {
 111                 double  theta;
 112                 double  dtheta;
 113                 double  phi;
 114                 double  dphi;
 115         }           tumble;
 116         dvector  (*ODE) (Par par, double x, double y, double z);
 117         Par         par;
 118 } flowstruct;
 119
 120 static flowstruct *flows = NULL;
 121
 122 static dvector
 123 Lorentz(Par par, double x, double y, double z)
 124 {
 125         dvector d;
 126
 127         d.x = par.a * (y - x);
 128         d.y = x * (par.b - z) - y;
 129         d.z = x * y - par.c * z;
 130         return d;
 131 }
 132
 133 static dvector
 134 Rossler(Par par, double x, double y, double z)
 135 {
 136         dvector d;
 137
 138         d.x = -(y + par.a * z);
 139         d.y = x + y * par.b;
 140         d.z = par.c + z * (x - 5.7);
 141         return d;
 142 }
 143
 144 static dvector
 145 RosslerCone(Par par, double x, double y, double z)
 146 {
 147         dvector d;
 148
 149         d.x = -(y + par.a * z);
 150         d.y = x + y * par.b - z * z * par.c;
 151         d.z = 0.2 + z * (x - 5.7);
 152         return d;
 153 }
 154
 155 static dvector
 156 Birkhoff(Par par, double x, double y, double z)
 157 {
 158         dvector d;
 159
 160         d.x = -y + par.b * sin(z);
 161         d.y = 0.7 * x + par.a * y * (0.1 - x * x);
 162         d.z = par.c;
 163         return d;
 164 }
 165
 166 static dvector
 167 Duffing(Par par, double x, double y, double z)
 168 {
 169         dvector d;
 170
 171         d.x = -par.a * x - y/2 - y * y * y/8 + par.b * cos(z);
 172         d.y = 2*x;
 173         d.z = par.c;
 174         return d;
 175 }
 176
 177 void
 178 init_flow(ModeInfo * mi)
 179 {
 180         flowstruct *sp;
 181         int         b;
 182         double      beemult = 1;
 183         dvector     range;
 184         static int  allocated = 0;
 185
 186         if (flows == NULL) {
 187                 if ((flows = (flowstruct *) calloc(MI_NUM_SCREENS(mi),
 188                                                sizeof (flowstruct))) == NULL)
 189                         return;
 190         }
 191         sp = &flows[MI_SCREEN(mi)];
 192
 193         sp->count = 0;
 194
 195         sp->width = MI_WIDTH(mi);
 196         sp->height = MI_HEIGHT(mi);
 197
 198         sp->tumble.theta = balance_rand(M_PI);
 199         sp->tumble.phi = balance_rand(M_PI);
 200         sp->tumble.dtheta = 0.002;
 201         sp->tumble.dphi = 0.001;
 202         sp->view.height = 0;
 203         sp->view.depth = 0; /* no perspective view */
 204
 205         switch (NRAND(8)) {
 206         case 0:
 207                 sp->view.depth = 10;
 208                 sp->view.height = 0.2;
 209                 beemult = 3;
 210         case 1:
 211                 sp->ODE = Lorentz;
 212                 sp->step = 0.02;
 213                 sp->size = 60;
 214                 sp->centre.x = 0;
 215                 sp->centre.y = 0;
 216                 sp->centre.z = 24;
 217                 range.x = 5;
 218                 range.y = 5;
 219                 range.z = 1;
 220                 sp->par.a = 10 + balance_rand(5);
 221                 sp->par.b = 28 + balance_rand(5);
 222                 sp->par.c = 2 + balance_rand(1);
 223                 break;
 224         case 2:
 225                 sp->view.depth = 10;
 226                 sp->view.height = 0.1;
 227                 beemult = 4;
 228         case 3:
 229                 sp->ODE = Rossler;
 230                 sp->step = 0.05;
 231                 sp->size = 24;
 232                 sp->centre.x = 0;
 233                 sp->centre.y = 0;
 234                 sp->centre.z = 3;
 235                 range.x = 4;
 236                 range.y = 4;
 237                 range.z = 7;
 238                 sp->par.a = 2 + balance_rand(1);
 239                 sp->par.b = 0.2 + balance_rand(0.1);
 240                 sp->par.c = 0.2 + balance_rand(0.1);
 241                 break;
 242         case 4:
 243                 sp->view.depth = 10;
 244                 sp->view.height = 0.1;
 245                 beemult = 3;
 246         case 5:
 247                 sp->ODE = RosslerCone;
 248                 sp->step = 0.05;
 249                 sp->size = 24;
 250                 sp->centre.x = 0;
 251                 sp->centre.y = 0;
 252                 sp->centre.z = 3;
 253                 range.x = 4;
 254                 range.y = 4;
 255                 range.z = 4;
 256                 sp->par.a = 2;
 257                 sp->par.b = 0.2;
 258                 sp->par.c = 0.25 + balance_rand(0.09);
 259                 break;
 260         case 6:
 261                 sp->ODE = Birkhoff;
 262                 sp->step = 0.04;
 263                 sp->size = 2.6;
 264                 sp->centre.x = 0;
 265                 sp->centre.y = 0;
 266                 sp->centre.z = 0;
 267                 range.x = 3;
 268                 range.y = 4;
 269                 range.z = 0;
 270                 sp->par.a = 10 + balance_rand(5);
 271                 sp->par.b = 0.35 + balance_rand(0.25);
 272                 sp->par.c = 1.57;
 273                 sp->tumble.theta = 0;
 274                 sp->tumble.phi = 0;
 275                 sp->tumble.dtheta = 0;
 276                 sp->tumble.dphi = 0;
 277                 break;
 278         case 7:
 279         default:
 280                 sp->ODE = Duffing;
 281                 sp->step = 0.02;
 282                 sp->size = 30;
 283                 sp->centre.x = 0;
 284                 sp->centre.y = 0;
 285                 sp->centre.z = 0;
 286                 range.x = 20;
 287                 range.y = 20;
 288                 range.z = 0;
 289                 sp->par.a = 0.2 + balance_rand(0.1);
 290                 sp->par.b = 27.0 + balance_rand(3.0);
 291                 sp->par.c = 1.33;
 292                 sp->tumble.theta = 0;
 293                 sp->tumble.phi = 0;
 294                 sp->tumble.dtheta = -NRAND(2)*sp->par.c*sp->step;
 295                 sp->tumble.dphi = 0;
 296                 beemult = 0.5;
 297                 break;
 298         }
 299
 300         sp->beecount = beemult * MI_COUNT(mi);
 301         if (sp->beecount < 0)   /* random variations */
 302                 sp->beecount = NRAND(-sp->beecount) + 1; /* Minimum 1 */
 303
 304         /* Clear the background. */
 305         MI_CLEARWINDOW(mi);
 306
 307         if(!allocated || sp->beecount != allocated){ /* reallocate */
 308                 if (sp->csegs != NULL) {
 309                         (void) free((void *) sp->csegs);
 310                         sp->csegs = NULL;
 311                 }
 312                 if (sp->cnsegs != NULL) {
 313                         (void) free((void *) sp->cnsegs);
 314                         sp->cnsegs = NULL;
 315                 }
 316                 if (sp->old_segs != NULL) {
 317                         (void) free((void *) sp->old_segs);
 318                         sp->old_segs = NULL;
 319                 }
 320                 if (sp->p[0] != NULL) {
 321                         (void) free((void *) sp->p[0]);
 322                         sp->p[0] = NULL;
 323                 }
 324                 if (sp->p[1] != NULL) {
 325                         (void) free((void *) sp->p[1]);
 326                         sp->p[1] = NULL;
 327                 }
 328         }
 329
 330         /* Allocate memory. */
 331
 332         if (!sp->csegs) {
 333                 sp->csegs = (XSegment *) malloc(sizeof (XSegment) * sp->beecount
 334                                                 * MI_NPIXELS(mi));
 335                 sp->cnsegs = (int *) malloc(sizeof (int) * MI_NPIXELS(mi));
 336
 337                 sp->old_segs = (XSegment *) malloc(sizeof (XSegment) * sp->beecount);
 338                 sp->p[0] = (dvector *) malloc(sizeof (dvector) * sp->beecount);
 339                 sp->p[1] = (dvector *) malloc(sizeof (dvector) * sp->beecount);
 340         }
 341
 342         /* Initialize point positions, velocities, etc. */
 343
 344         for (b = 0; b < sp->beecount; b++) {
 345                 X(1, b) = X(0, b) = balance_rand(range.x);
 346                 Y(1, b) = Y(0, b) = balance_rand(range.y);
 347                 Z(1, b) = Z(0, b) = balance_rand(range.z);
 348         }
 349 }
 350
 351 void
 352 draw_flow(ModeInfo * mi)
 353 {
 354         Display    *display = MI_DISPLAY(mi);
 355         Window      window = MI_WINDOW(mi);
 356         GC          gc = MI_GC(mi);
 357         flowstruct *sp = &flows[MI_SCREEN(mi)];
 358         int         b, c, i;
 359         int         col, ix;
 360         double      M[3][3]; /* transformation matrix */
 361
 362         if(!sp->view.depth){ /* simple 3D tumble */
 363                 double      sint, cost, sinp, cosp;
 364                 sp->tumble.theta += sp->tumble.dtheta;
 365                 sp->tumble.phi += sp->tumble.dphi;
 366                 sint = sin(sp->tumble.theta);
 367                 cost = cos(sp->tumble.theta);
 368                 sinp = sin(sp->tumble.phi);
 369                 cosp = cos(sp->tumble.phi);
 370                 M[0][0]= cost; M[0][1]=-sint*cosp; M[0][2]= sint*sinp;
 371                 M[1][0]= sint; M[1][1]= cost*cosp; M[1][2]=-cost*sinp;
 372                 M[2][0]= 0;    M[2][1]= 0;         M[2][2]= 1;
 373         } else { /* initialize matrix */
 374                 M[0][0]= 0; M[0][1]= 0; M[0][2]= 0;
 375                 M[1][0]= 0; M[1][1]= 0; M[1][2]= 0;
 376                 M[2][0]= 0; M[2][1]= 0; M[2][2]= 0;
 377
 378         }
 379
 380         for (col = 0; col < MI_NPIXELS(mi); col++)
 381                 sp->cnsegs[col] = 0;
 382
 383         MI_IS_DRAWN(mi) = True;
 384
 385         /* <=- Bees -=> */
 386         for (b = 0; b < sp->beecount; b++) {
 387                 /* Age the arrays. */
 388                 X(1, b) = X(0, b);
 389                 Y(1, b) = Y(0, b);
 390                 Z(1, b) = Z(0, b);
 391
 392                 /* 2nd order Kunge Kutta */
 393                 {
 394                         dvector     k1, k2;
 395
 396                         k1 = sp->ODE(sp->par, X(1, b), Y(1, b), Z(1, b));
 397                         k1.x *= sp->step;
 398                         k1.y *= sp->step;
 399                         k1.z *= sp->step;
 400                         k2 = sp->ODE(sp->par, X(1, b) + k1.x, Y(1, b) + k1.y, Z(1, b) + k1.z);
 401                         k2.x *= sp->step;
 402                         k2.y *= sp->step;
 403                         k2.z *= sp->step;
 404                         X(0, b) = X(1, b) + (k1.x + k2.x) / 2.0;
 405                         Y(0, b) = Y(1, b) + (k1.y + k2.y) / 2.0;
 406                         Z(0, b) = Z(1, b) + (k1.z + k2.z) / 2.0;
 407                 }
 408
 409                 /* Colour according to bee */
 410                 col = b % (MI_NPIXELS(mi) - 1);
 411                 ix = col * sp->beecount + sp->cnsegs[col];
 412
 413                 /* Fill the segment lists. */
 414
 415                 if(sp->view.depth) /* perspective view has special points */
 416                         if(b==0){ /* point of view */
 417                                 sp->centre.x=X(0, b);
 418                                 sp->centre.y=Y(0, b);
 419                                 sp->centre.z=Z(0, b);
 420                         }else if(b==1){ /* neighbour: used to compute local axes */
 421                                 double x[3], p[3], x2=0, xp=0;
 422                                 int j;
 423
 424                                 /* forward */
 425                                 x[0] = X(0, 0) - X(1, 0);
 426                                 x[1] = Y(0, 0) - Y(1, 0);
 427                                 x[2] = Z(0, 0) - Z(1, 0);
 428
 429                                 /* neighbour */
 430                                 p[0] = X(0, 1) - X(1, 0);
 431                                 p[1] = Y(0, 1) - Y(1, 0);
 432                                 p[2] = Z(0, 1) - Z(1, 0);
 433
 434                                 for(i=0; i<3; i++){
 435                                         x2+= x[i]*x[i];    /* X . X */
 436                                         xp+= x[i]*p[i];    /* X . P */
 437                                         M[0][i] = x[i];    /* X */
 438                                 }
 439
 440                                 for(i=0; i<3; i++)               /* (X x P) x X */
 441                                         M[1][i] = x2*p[i] - xp*x[i]; /* == (X . X) P - (X . P) X */
 442
 443                                 M[2][0] =  x[1]*p[2] - x[2]*p[1]; /* X x P */
 444                                 M[2][1] = -x[0]*p[2] + x[2]*p[0];
 445                                 M[2][2] =  x[0]*p[1] - x[1]*p[0];
 446
 447                                 /* normalise axes */
 448                                 for(j=0; j<3; j++){
 449                                         double A=0;
 450                                         for(i=0; i<3; i++) A+=M[j][i]*M[j][i]; /* sum squares */
 451                                         A=sqrt(A);
 452                                         for(i=0; i<3; i++) M[j][i]/=A;
 453                                 }
 454
 455                                 X(0, 1)=X(0, 0)+M[1][0]; /* adjust neighbour */
 456                                 Y(0, 1)=Y(0, 0)+M[1][1];
 457                                 Z(0, 1)=Z(0, 0)+M[1][2];
 458
 459 #if 0  /* display local axes for testing */
 460                                 X(1, b)=X(0, 0);
 461                                 Y(1, b)=Y(0, 0);
 462                                 Z(1, b)=Z(0, 0);
 463                         }else if(b==2){
 464                                 X(0, b)=X(0, 0)+0.5*M[0][0];
 465                                 Y(0, b)=Y(0, 0)+0.5*M[0][1];
 466                                 Z(0, b)=Z(0, 0)+0.5*M[0][2];
 467                                 X(1, b)=X(0, 0);
 468                                 Y(1, b)=Y(0, 0);
 469                                 Z(1, b)=Z(0, 0);
 470                         }else if(b==3){
 471                                 X(0, b)=X(0, 0)+1.5*M[2][0];
 472                                 Y(0, b)=Y(0, 0)+1.5*M[2][1];
 473                                 Z(0, b)=Z(0, 0)+1.5*M[2][2];
 474                                 X(1, b)=X(0, 0);
 475                                 Y(1, b)=Y(0, 0);
 476                                 Z(1, b)=Z(0, 0);
 477 #endif
 478                         }
 479
 480                 for(i=0; i<2; i++){
 481                         double x=X(i,b)-sp->centre.x;
 482                         double y=Y(i,b)-sp->centre.y;
 483                         double z=Z(i,b)-sp->centre.z;
 484                         double X=M[0][0]*x + M[0][1]*y + M[0][2]*z;
 485                         double Y=M[1][0]*x + M[1][1]*y + M[1][2]*z;
 486                         double Z=M[2][0]*x + M[2][1]*y + M[2][2]*z+sp->view.height;
 487                         double absx, absy;
 488                         if(sp->view.depth){
 489                                 if(X <= 0) break;
 490                                 absx=SCALE_X(sp->view.depth*Y/X);
 491                                 absy=SCALE_Y(sp->view.depth*Z/X);
 492                                 if(absx < -sp->width || absx > 2*sp->width ||
 493                                    absy < -sp->height || absy > 2*sp->height)
 494                                         break;
 495                         }else{
 496                                 absx=SCALE_X(X);
 497                                 absy=SCALE_Y(Y);
 498                         }
 499                         if(i){
 500                                 sp->csegs[ix].x1 = (short) absx;
 501                                 sp->csegs[ix].y1 = (short) absy;
 502                         }else{
 503                                 sp->csegs[ix].x2 = (short) absx;
 504                                 sp->csegs[ix].y2 = (short) absy;
 505                         }
 506                 }
 507                 if(i == 2) /* both assigned */
 508                         sp->cnsegs[col]++;
 509     }
 510         if (sp->count) { /* erase */
 511                 XSetForeground(display, gc, MI_BLACK_PIXEL(mi));
 512                 XDrawSegments(display, window, gc, sp->old_segs, sp->nold_segs);
 513         }
 514
 515         if (MI_NPIXELS(mi) > 2){ /* render colour */
 516                 for (col = 0; col < MI_NPIXELS(mi); col++)
 517                         if (sp->cnsegs[col] > 0) {
 518                                 XSetForeground(display, gc, MI_PIXEL(mi, col));
 519                                 XDrawSegments(display, window, gc,
 520                                               sp->csegs + col * sp->beecount, sp->cnsegs[col]);
 521                         }
 522         } else {                /* render mono */
 523                 XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
 524                 XDrawSegments(display, window, gc,
 525                                           sp->csegs + col * sp->beecount, sp->cnsegs[col]);
 526         }
 527
 528         /* Copy to erase-list */
 529         for (col = 0, c = 0; col < MI_NPIXELS(mi); col++)
 530                 for (b = 0; b < sp->cnsegs[col]; b++)
 531                         sp->old_segs[c++] = (sp->csegs + col * sp->beecount)[b];
 532         sp->nold_segs = c;
 533
 534         if (++sp->count > MI_CYCLES(mi)) /* pick a new flow */
 535                 init_flow(mi);
 536 }
 537
 538 void
 539 release_flow(ModeInfo * mi)
 540 {
 541         if (flows != NULL) {
 542                 int         screen;
 543
 544                 for (screen = 0; screen < MI_NUM_SCREENS(mi); screen++) {
 545                         flowstruct *sp = &flows[screen];
 546
 547                         if (sp->csegs != NULL)
 548                                 (void) free((void *) sp->csegs);
 549                         if (sp->cnsegs != NULL)
 550                                 (void) free((void *) sp->cnsegs);
 551                         if (sp->old_segs != NULL)
 552                                 (void) free((void *) sp->old_segs);
 553                         if (sp->p[0] != NULL)
 554                                 (void) free((void *) sp->p[0]);
 555                         if (sp->p[1] != NULL)
 556                                 (void) free((void *) sp->p[1]);
 557                 }
 558                 (void) free((void *) flows);
 559                 flows = NULL;
 560         }
 561 }
 562
 563 void
 564 refresh_flow(ModeInfo * mi)
 565 {
 566         MI_CLEARWINDOW(mi);
 567 }