X-Git-Url: http://git.hungrycats.org/cgi-bin/gitweb.cgi?p=xscreensaver;a=blobdiff_plain;f=hacks%2Fflow.c;h=df9aadc366efc82d5866e3e75c18a3d0b2e849fc;hp=607994d0387b1494e8de74aee346fc5906640367;hb=4361b69d3178d7fc98d0388f9a223af6c2651aba;hpb=6cee540bdbb571485cd5e519f89f389faebd0495 diff --git a/hacks/flow.c b/hacks/flow.c index 607994d0..df9aadc3 100644 --- a/hacks/flow.c +++ b/hacks/flow.c @@ -1,14 +1,16 @@ -/* -*- Mode: C; tab-width: 4 -*- */ +/* -*- Mode: C; tab-width: 4; c-basic-offset: 4 -*- */ /* flow --- flow of strange bees */ -#if !defined( lint ) && !defined( SABER ) -static const char sccsid[] = "@(#)flow.c 4.10 98/04/24 xlockmore"; - +#if 0 +static const char sccsid[] = "@(#)flow.c 5.00 2000/11/01 xlockmore"; #endif /*- - * Copyright (c) 1996 by Tim Auckland - * Portions added by Stephen Davies are Copyright (c) 2000 Stephen Davies + * Copyright (c) 1996 by Tim Auckland + * Incorporating some code from Stephen Davies Copyright (c) 2000 + * + * Search code based on techniques described in "Strange Attractors: + * Creating Patterns in Chaos" by Julien C. Sprott * * Permission to use, copy, modify, and distribute this software and its * documentation for any purpose and without fee is hereby granted, @@ -25,1013 +27,1200 @@ static const char sccsid[] = "@(#)flow.c 4.10 98/04/24 xlockmore"; * "flow" shows a variety of continuous phase-space flows around strange * attractors. It includes the well-known Lorentz mask (the "Butterfly" * of chaos fame), two forms of Rossler's "Folded Band" and Poincare' - * sections of the "Birkhoff Bagel" and Duffing's forced occilator. + * sections of the "Birkhoff Bagel" and Duffing's forced occilator. "flow" + * can now discover new attractors. * * Revision History: - * 21-Feb-00: Major hackage by Chalky (Stephen Davies, chalky@null.net) - * Forced perspective mode, added 3d box around attractor which - * involved coding 3d-planar-clipping against the view-frustrum - * thingy. Also made view alternate between piggybacking on a 'bee' - * to zooming around outside the attractor. Most bees slow down and - * stop, to make the structure of the attractor more obvious. - * 31-Nov-98: [TDA] Added Duffing (what a strange day that was :) DAB) + * + * 29-Oct-2004: [TDA] Discover Attractors unknown to science. + * Replace 2D rendering of Periodic Attractors with a 3D + * 'interrupted' rendering. Replace "-/+allow2d" with "-/+periodic" + * Replace all ODE formulae with completely generic forms. + * Add '-search' option to perform background high-speed discovery + * for completely new attractors without impacting rendering + * performance. + * Use gaussian distribution for initial point positions and for + * parameter search. + * Add "+dbuf" option to allow Double-Buffering to be turned off on + * slow X servers. + * Remove redundant '-zoom' option. Now automatically zooms if both + * rotation and riding are permitted. + * Replace dynamic bounding box with static one pre-calculated + * during discovery phase. + * Simplify and fix bounding box clipping code. Should now be safe + * to run without double buffer on all XFree86 servers if desired. + * 12-Oct-2004: [TDA] Merge Xscreensaver and Xlockmore branches + * Added Chalky's orbital camera, but made the zooming work by + * flying the camera rather than interpolating the view transforms. + * Added Chalky's Bounding Box, but time-averaged the boundaries to + * let the lost bees escape. + * Added Chalky's 'view-frustrum' clipping, but only applying it to + * the Bounding Box. Trails make clipping less useful. + * Added Chalky's "-slow" and "-freeze" options for compatibility, + * but haven't implemented the features, since the results are ugly + * and make no mathematical contribution. + * Added Double-Buffering as a work-around for a persistent XFree86 + * bug that left debris on the screen. + * 21-Mar-2003: [TDA] Trails added (XLockmore branch) + * 01-Nov-2000: [TDA] Allocation checks (XLockmore branch) + * 21-Feb-2000: [Chalky] Major hackage (Stephen Davies, chalky@null.net) + * (Xscreensaver branch) + * Forced perspective mode, added 3d box around attractor which + * involved coding 3d-planar-clipping against the view-frustrum + * thingy. Also made view alternate between piggybacking on a 'bee' + * to zooming around outside the attractor. Most bees slow down and + * stop, to make the structure of the attractor more obvious. +* 28-Jan-1999: [TDA] Catch 'lost' bees in flow.c and disable them. + * (XLockmore branch) + * I chose to disable them rather than reinitialise them because + * reinitialising can produce fake attractors. + * This has allowed me to relax some of the parameters and initial + * conditions slightly to catch some of the more extreme cases. As a + * result you may see some bees fly away at the start - these are the ones + * that 'missed' the attractor. If the bee with the camera should fly + * away the mode will restart :-) + * 31-Nov-1998: [TDA] Added Duffing (what a strange day that was :) DAB) * Duffing's forced oscillator has been added to the formula list and * the parameters section has been updated to display it in Poincare' * section. - * 30-Nov-98: [TDA] Added travelling perspective option + * 30-Nov-1998: [TDA] Added travelling perspective option * A more exciting point-of-view has been added to all autonomous flows. * This views the flow as seen by a particle moving with the flow. In the * metaphor of the original code, I've attached a camera to one of the * trained bees! - * 30-Nov-98: [TDA] Much code cleanup. - * 09-Apr-97: [TDA] Ported to xlockmore-4 - * 18-Jul-96: Adapted from swarm.c Copyright (c) 1991 by Patrick J. Naughton. - * 31-Aug-90: Adapted from xswarm by Jeff Butterworth. (butterwo@ncsc.org) + * 30-Nov-1998: [TDA] Much code cleanup. + * 09-Apr-1997: [TDA] Ported to xlockmore-4 + * 18-Jul-1996: Adapted from swarm.c Copyright (c) 1991 by Patrick J. Naughton. + * 31-Aug-1990: Adapted from xswarm by Jeff Butterworth. (butterwo@ncsc.org). */ #ifdef STANDALONE -# define PROGCLASS "Flow" -# define HACK_INIT init_flow -# define HACK_DRAW draw_flow -# define flow_opts xlockmore_opts -# define DEFAULTS "*delay: 1000 \n" \ - "*count: 500 \n" \ - "*cycles: 3000 \n" \ - "*ncolors: 200 \n" \ - "*rotate: True \n" \ - "*ride: True \n" \ - "*zoom: True \n" \ - "*allow2d: True \n" \ - "*box: True \n" \ - "*slow: True \n" \ - "*freeze: True \n" -# define SMOOTH_COLORS +# define MODE_flow +# define DEFAULTS "*delay: 10000 \n" \ + "*count: 3000 \n" \ + "*size: -10 \n" \ + "*cycles: 10000 \n" \ + "*ncolors: 200 \n" + +# define release_flow 0 # include "xlockmore.h" /* in xscreensaver distribution */ -# include "erase.h" - #else /* STANDALONE */ # include "xlock.h" /* in xlockmore distribution */ #endif /* STANDALONE */ -XrmOptionDescRec flow_options[]; -ModeSpecOpt flow_opts = { 7, flow_options, 0, NULL, NULL }; +#ifdef MODE_flow + +#define DEF_ROTATE "TRUE" +#define DEF_RIDE "TRUE" +#define DEF_BOX "TRUE" +#define DEF_PERIODIC "TRUE" +#define DEF_SEARCH "TRUE" +#define DEF_DBUF "TRUE" + +static Bool rotatep; +static Bool ridep; +static Bool boxp; +static Bool periodicp; +static Bool searchp; +static Bool dbufp; + +static XrmOptionDescRec opts[] = { + {"-rotate", ".flow.rotate", XrmoptionNoArg, "on"}, + {"+rotate", ".flow.rotate", XrmoptionNoArg, "off"}, + {"-ride", ".flow.ride", XrmoptionNoArg, "on"}, + {"+ride", ".flow.ride", XrmoptionNoArg, "off"}, + {"-box", ".flow.box", XrmoptionNoArg, "on"}, + {"+box", ".flow.box", XrmoptionNoArg, "off"}, + {"-periodic", ".flow.periodic", XrmoptionNoArg, "on"}, + {"+periodic", ".flow.periodic", XrmoptionNoArg, "off"}, + {"-search", ".flow.search", XrmoptionNoArg, "on"}, + {"+search", ".flow.search", XrmoptionNoArg, "off"}, + {"-dbuf", ".flow.dbuf", XrmoptionNoArg, "on"}, + {"+dbuf", ".flow.dbuf", XrmoptionNoArg, "off"}, +}; + +static argtype vars[] = { + {&rotatep, "rotate", "Rotate", DEF_ROTATE, t_Bool}, + {&ridep, "ride", "Ride", DEF_RIDE, t_Bool}, + {&boxp, "box", "Box", DEF_BOX, t_Bool}, + {&periodicp, "periodic", "Periodic", DEF_PERIODIC, t_Bool}, + {&searchp, "search", "Search", DEF_SEARCH, t_Bool}, + {&dbufp, "dbuf", "Dbuf", DEF_DBUF, t_Bool}, +}; + +static OptionStruct desc[] = { + {"-/+rotate", "turn on/off rotating around attractor."}, + {"-/+ride", "turn on/off ride in the flow."}, + {"-/+box", "turn on/off bounding box."}, + {"-/+periodic", "turn on/off periodic attractors."}, + {"-/+search", "turn on/off search for new attractors."}, + {"-/+dbuf", "turn on/off double buffering."}, +}; + +ENTRYPOINT ModeSpecOpt flow_opts = +{sizeof opts / sizeof opts[0], opts, + sizeof vars / sizeof vars[0], vars, desc}; #ifdef USE_MODULES ModStruct flow_description = { "flow", "init_flow", "draw_flow", "release_flow", "refresh_flow", "init_flow", NULL, &flow_opts, - 1000, 1024, 3000, 1, 64, 1.0, "", + 1000, 1024, 10000, -10, 200, 1.0, "", "Shows dynamic strange attractors", 0, NULL }; #endif -typedef struct { - double x; - double y; - double z; -} dvector; +typedef struct { double x, y, z; } dvector; + +#define N_PARS 20 /* Enough for Full Cubic or Periodic Cubic */ +typedef dvector Par[N_PARS]; +enum { /* Name the parameter indices to make it easier to write + standard examples */ + C, + X,XX,XXX,XXY,XXZ,XY,XYY,XYZ,XZ,XZZ, + Y,YY,YYY,YYZ,YZ,YZZ, + Z,ZZ,ZZZ, + SINY = XY /* OK to overlap in this case */ +}; -typedef struct { - double a, b, c; -} Par; +/* Camera target [TDA] */ +typedef enum { + ORBIT = 0, + BEE = 1 +} Chaseto; /* Macros */ -#define X(t,b) (sp->p[t][b].x) -#define Y(t,b) (sp->p[t][b].y) -#define Z(t,b) (sp->p[t][b].z) +#define IX(C) ((C) * segindex + sp->cnsegs[(C)]) +#define B(t,b) (sp->p + (t) + (b) * sp->taillen) +#define X(t,b) (B((t),(b))->x) +#define Y(t,b) (B((t),(b))->y) +#define Z(t,b) (B((t),(b))->z) #define balance_rand(v) ((LRAND()/MAXRAND*(v))-((v)/2)) /* random around 0 */ -#define SCALE_X(A) (sp->width/2+sp->width/sp->size*(A)) -/*#define SCALE_Y(A) (sp->height/2+sp->height/sp->size*(A))*/ -#define SCALE_Y(A) (sp->height/2+sp->width/sp->size*(A)) +#define LOST_IN_SPACE 2000.0 +#define INITIALSTEP 0.04 +#define EYEHEIGHT 0.005 +#define MINTRAIL 2 +#define BOX_L 36 -/* Mode of operation. Rotate, ride and zoom are mutually exclusive */ -typedef enum { - FLOW_ROTATE = 1, /* Rotate around attractor */ - FLOW_RIDE = 2, /* Ride a trained bee */ - FLOW_ZOOM = 4, /* Zoom in and out */ - FLOW_2D = 8, /* Allow 2D attractors */ - FLOW_BOX = 16, /* Compute a box around the attractor */ - FLOW_SLOW = 32, /* Some bees are slower (and have antifreeze) */ - FLOW_FREEZE = 64 /* Freeze some of the bees in action */ -} FlowMode; - -#define FLOW_DEFAULT (FLOW_ROTATE|FLOW_RIDE|FLOW_ZOOM|FLOW_2D|\ - FLOW_BOX|FLOW_SLOW|FLOW_FREEZE) +/* Points that make up the box (normalized coordinates) */ +static const double box[][3] = { + {1,1,1}, /* 0 */ + {1,1,-1}, /* 1 */ + {1,-1,-1}, /* 2 */ + {1,-1,1}, /* 3 */ + {-1,1,1}, /* 4 */ + {-1,1,-1}, /* 5 */ + {-1,-1,-1},/* 6 */ + {-1,-1,1}, /* 7 */ + {1, .8, .8}, + {1, .8,-.8}, + {1,-.8,-.8}, + {1,-.8, .8}, + { .8,1, .8}, + { .8,1,-.8}, + {-.8,1,-.8}, + {-.8,1, .8}, + { .8, .8,1}, + { .8,-.8,1}, + {-.8,-.8,1}, + {-.8, .8,1}, + {-1, .8, .8}, + {-1, .8,-.8}, + {-1,-.8,-.8}, + {-1,-.8, .8}, + { .8,-1, .8}, + { .8,-1,-.8}, + {-.8,-1,-.8}, + {-.8,-1, .8}, + { .8, .8,-1}, + { .8,-.8,-1}, + {-.8,-.8,-1}, + {-.8, .8,-1} +}; + +/* Lines connecting the box dots */ +static const double lines[][2] = { + {0,1}, {1,2}, {2,3}, {3,0}, /* box */ + {4,5}, {5,6}, {6,7}, {7,4}, + {0,4}, {1,5}, {2,6}, {3,7}, + {4+4,5+4}, {5+4,6+4}, {6+4,7+4}, {7+4,4+4}, + {4+8,5+8}, {5+8,6+8}, {6+8,7+8}, {7+8,4+8}, + {4+12,5+12}, {5+12,6+12}, {6+12,7+12}, {7+12,4+12}, + {4+16,5+16}, {5+16,6+16}, {6+16,7+16}, {7+16,4+16}, + {4+20,5+20}, {5+20,6+20}, {6+20,7+20}, {7+20,4+20}, + {4+24,5+24}, {5+24,6+24}, {6+24,7+24}, {7+24,4+24}, +}; typedef struct { - int width; - int height; - int count; - double size; - + /* Variables used in rendering */ + dvector cam[3]; /* camera flight path */ + int chasetime; + Chaseto chaseto; + Pixmap buffer; /* Double Buffer */ + dvector circle[2]; /* POV that circles around the scene */ + dvector centre; /* centre */ int beecount; /* number of bees */ XSegment *csegs; /* bee lines */ int *cnsegs; XSegment *old_segs; /* old bee lines */ int nold_segs; - double step; - double slow; - double slow_view; - dvector centre; /* centre */ - dvector range; - struct { - double depth; - double height; - } view; - dvector circle[2]; /* POV that circles around the scene */ - dvector *p[2]; /* bee positions x[time][bee#] */ - struct { - double theta; - double dtheta; - double phi; - double dphi; - } tumble; + int taillen; + + /* Variables common to iterators */ dvector (*ODE) (Par par, double x, double y, double z); + dvector range; /* Initial conditions */ + double yperiod; /* ODE's where Y is periodic. */ + + /* Variables used in iterating main flow */ Par par; - FlowMode mode; /* Mode of operation */ + dvector *p; /* bee positions x[time][bee#] */ + int count; + double lyap; + double size; + dvector mid; /* Effective bounding box */ + double step; + + /* second set of variables, used for parallel search */ + Par par2; + dvector p2[2]; + int count2; + double lyap2; + double size2; + dvector mid2; + double step2; + } flowstruct; -static flowstruct *flows = NULL; +static flowstruct *flows = (flowstruct *) NULL; + +/* + * Private functions + */ + + +/* ODE functions */ +/* Generic 3D Cubic Polynomial. Includes all the Quadratics (Lorentz, + Rossler) and much more! */ + +/* I considered offering a seperate 'Quadratic' option, since Cubic is + clearly overkill for the standard examples, but the performance + difference is too small to measure. The compute time is entirely + dominated by the XDrawSegments calls anyway. [TDA] */ static dvector -Lorentz(Par par, double x, double y, double z) +Cubic(Par a, double x, double y, double z) { dvector d; + 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 + + a[XXZ].x*x*x*z + a[XY].x*x*y + a[XYY].x*x*y*y + a[XYZ].x*x*y*z + + a[XZ].x*x*z + a[XZZ].x*x*z*z + a[Y].x*y + a[YY].x*y*y + + a[YYY].x*y*y*y + a[YYZ].x*y*y*z + a[YZ].x*y*z + a[YZZ].x*y*z*z + + a[Z].x*z + a[ZZ].x*z*z + a[ZZZ].x*z*z*z; + + 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 + + a[XXZ].y*x*x*z + a[XY].y*x*y + a[XYY].y*x*y*y + a[XYZ].y*x*y*z + + a[XZ].y*x*z + a[XZZ].y*x*z*z + a[Y].y*y + a[YY].y*y*y + + a[YYY].y*y*y*y + a[YYZ].y*y*y*z + a[YZ].y*y*z + a[YZZ].y*y*z*z + + a[Z].y*z + a[ZZ].y*z*z + a[ZZZ].y*z*z*z; + + 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 + + a[XXZ].z*x*x*z + a[XY].z*x*y + a[XYY].z*x*y*y + a[XYZ].z*x*y*z + + a[XZ].z*x*z + a[XZZ].z*x*z*z + a[Y].z*y + a[YY].z*y*y + + a[YYY].z*y*y*y + a[YYZ].z*y*y*z + a[YZ].z*y*z + a[YZZ].z*y*z*z + + a[Z].z*z + a[ZZ].z*z*z + a[ZZZ].z*z*z*z; - d.x = par.a * (y - x); - d.y = x * (par.b - z) - y; - d.z = x * y - par.c * z; return d; } +/* 3D Cubic in (x,z) with periodic sinusoidal forcing term in x. y is + the independent periodic (time) axis. This includes Birkhoff's + Bagel and Duffing's Attractor */ static dvector -Rossler(Par par, double x, double y, double z) +Periodic(Par a, double x, double y, double z) { dvector d; - d.x = -(y + par.a * z); - d.y = x + y * par.b; - d.z = par.c + z * (x - 5.7); + d.x = a[C].x + a[X].x*x + a[XX].x*x*x + a[XXX].x*x*x*x + + a[XXZ].x*x*x*z + a[XZ].x*x*z + a[XZZ].x*x*z*z + a[Z].x*z + + a[ZZ].x*z*z + a[ZZZ].x*z*z*z + a[SINY].x*sin(y); + + d.y = a[C].y; + + d.z = a[C].z + a[X].z*x + a[XX].z*x*x + a[XXX].z*x*x*x + + a[XXZ].z*x*x*z + a[XZ].z*x*z + a[XZZ].z*x*z*z + a[Z].z*z + + a[ZZ].z*z*z + a[ZZZ].z*z*z*z; + return d; } -static dvector -RosslerCone(Par par, double x, double y, double z) +/* Numerical integration of the ODE using 2nd order Runge Kutta. + Returns length^2 of the update, so that we can detect if the step + size needs reducing. */ +static double +Iterate(dvector *p, dvector(*ODE)(Par par, double x, double y, double z), + Par par, double step) +{ + dvector k1, k2, k3; + + k1 = ODE(par, p->x, p->y, p->z); + k1.x *= step; + k1.y *= step; + k1.z *= step; + k2 = ODE(par, p->x + k1.x, p->y + k1.y, p->z + k1.z); + k2.x *= step; + k2.y *= step; + k2.z *= step; + k3.x = (k1.x + k2.x) / 2.0; + k3.y = (k1.y + k2.y) / 2.0; + k3.z = (k1.z + k2.z) / 2.0; + + p->x += k3.x; + p->y += k3.y; + p->z += k3.z; + + return k3.x*k3.x + k3.y*k3.y + k3.z*k3.z; +} + +/* Memory functions */ + +#define deallocate(p,t) if (p!=NULL) {free(p); p=(t*)NULL; } +#define allocate(p,t,s) if ((p=(t*)malloc(sizeof(t)*s))==NULL)\ +{free_flow(mi);return;} + +static void +free_flow(ModeInfo * mi) { - dvector d; + flowstruct *sp = &flows[MI_SCREEN(mi)]; + deallocate(sp->csegs, XSegment); + deallocate(sp->cnsegs, int); + deallocate(sp->old_segs, XSegment); + deallocate(sp->p, dvector); +} - d.x = -(y + par.a * z); - d.y = x + y * par.b - z * z * par.c; - d.z = 0.2 + z * (x - 5.7); - return d; +/* Generate Gaussian random number: mean 0, "amplitude" A (actually + A is 3*standard deviation). */ + +/* Note this generates a pair of gaussian variables, so it saves one + to give out next time it's called */ +static double +Gauss_Rand(double A) +{ + static double d; + static Bool ready = 0; + if(ready) { + ready = 0; + return A/3 * d; + } else { + double x, y, w; + do { + x = 2.0 * (double)LRAND() / MAXRAND - 1.0; + y = 2.0 * (double)LRAND() / MAXRAND - 1.0; + w = x*x + y*y; + } while(w >= 1.0); + + w = sqrt((-2 * log(w))/w); + ready = 1; + d = x * w; + return A/3 * y * w; + } } -static dvector -Birkhoff(Par par, double x, double y, double z) +/* Attempt to discover new atractors by sending a pair of bees on a + fast trip through the new flow and computing their Lyapunov + exponent. Returns False if the bees fly away. + + If the bees stay bounded, the new bounds and the Lyapunov exponent + are stored in sp and the function returns True. + + Repeat invocations continue the flow and improve the accuracy of + the bounds and the Lyapunov exponent. Set sp->count2 to zero to + start a new test. + + Acts on alternate variable set, so that it can be run in parallel + with the main flow */ + +static Bool +discover(ModeInfo * mi) { - dvector d; + flowstruct *sp; + double l = 0; + dvector dl; + dvector max, min; + double dl2, df, rs, lsum = 0, s, maxv2 = 0, v2; - d.x = -y + par.b * sin(z); - d.y = 0.7 * x + par.a * y * (0.1 - x * x); - d.z = par.c; - return d; + int N, i, nl = 0; + + if (flows == NULL) + return 0; + sp = &flows[MI_SCREEN(mi)]; + + if(sp->count2 == 0) { + /* initial conditions */ + sp->p2[0].x = Gauss_Rand(sp->range.x); + sp->p2[0].y = (sp->yperiod > 0)? + balance_rand(sp->range.y) : Gauss_Rand(sp->range.y); + sp->p2[0].z = Gauss_Rand(sp->range.z); + + /* 1000 steps to find an attractor */ + /* Most cases explode out here */ + for(N=0; N < 1000; N++){ + Iterate(sp->p2, sp->ODE, sp->par2, sp->step2); + if(sp->yperiod > 0 && sp->p2[0].y > sp->yperiod) + sp->p2[0].y -= sp->yperiod; + if(fabs(sp->p2[0].x) > LOST_IN_SPACE || + fabs(sp->p2[0].y) > LOST_IN_SPACE || + fabs(sp->p2[0].z) > LOST_IN_SPACE) { + return 0; + } + sp->count2++; + } + /* Small perturbation */ + sp->p2[1].x = sp->p2[0].x + 0.000001; + sp->p2[1].y = sp->p2[0].y; + sp->p2[1].z = sp->p2[0].z; + } + + /* Reset bounding box */ + max.x = min.x = sp->p2[0].x; + max.y = min.y = sp->p2[0].y; + max.z = min.z = sp->p2[0].z; + + /* Compute Lyapunov Exponent */ + + /* (Technically, we're only estimating the largest Lyapunov + Exponent, but that's all we need to know to determine if we + have a strange attractor.) [TDA] */ + + /* Fly two bees close together */ + for(N=0; N < 5000; N++){ + for(i=0; i< 2; i++) { + v2 = Iterate(sp->p2+i, sp->ODE, sp->par2, sp->step2); + if(sp->yperiod > 0 && sp->p2[i].y > sp->yperiod) + sp->p2[i].y -= sp->yperiod; + + if(fabs(sp->p2[i].x) > LOST_IN_SPACE || + fabs(sp->p2[i].y) > LOST_IN_SPACE || + fabs(sp->p2[i].z) > LOST_IN_SPACE) { + return 0; + } + if(v2 > maxv2) maxv2 = v2; /* Track max v^2 */ + } + + /* find bounding box */ + if ( sp->p2[0].x < min.x ) min.x = sp->p2[0].x; + else if ( sp->p2[0].x > max.x ) max.x = sp->p2[0].x; + if ( sp->p2[0].y < min.y ) min.y = sp->p2[0].y; + else if ( sp->p2[0].y > max.y ) max.y = sp->p2[0].y; + if ( sp->p2[0].z < min.z ) min.z = sp->p2[0].z; + else if ( sp->p2[0].z > max.z ) max.z = sp->p2[0].z; + + /* Measure how much we have to pull the two bees to prevent + them diverging. */ + dl.x = sp->p2[1].x - sp->p2[0].x; + dl.y = sp->p2[1].y - sp->p2[0].y; + dl.z = sp->p2[1].z - sp->p2[0].z; + + dl2 = dl.x*dl.x + dl.y*dl.y + dl.z*dl.z; + if(dl2 > 0) { + df = 1e12 * dl2; + rs = 1/sqrt(df); + sp->p2[1].x = sp->p2[0].x + rs * dl.x; + sp->p2[1].y = sp->p2[0].y + rs * dl.y; + sp->p2[1].z = sp->p2[0].z + rs * dl.z; + lsum = lsum + log(df); + nl = nl + 1; + l = M_LOG2E / 2 * lsum / nl / sp->step2; + } + sp->count2++; + } + /* Anything that didn't explode has a finite attractor */ + /* If Lyapunov is negative then it probably hit a fixed point or a + * limit cycle. Positive Lyapunov indicates a strange attractor. */ + + sp->lyap2 = l; + + sp->size2 = max.x - min.x; + s = max.y - min.y; + if(s > sp->size2) sp->size2 = s; + s = max.z - min.z; + if(s > sp->size2) sp->size2 = s; + + sp->mid2.x = (max.x + min.x) / 2; + sp->mid2.y = (max.y + min.y) / 2; + sp->mid2.z = (max.z + min.z) / 2; + + if(sqrt(maxv2) > sp->size2 * 0.2) { + /* Flowing too fast, reduce step size. This + helps to eliminate high-speed limit cycles, + which can show +ve Lyapunov due to integration + inaccuracy. */ + sp->step2 /= 2; + } + return 1; } -static dvector -Duffing(Par par, double x, double y, double z) +/* Sets up initial conditions for a flow without all the extra baggage + that goes with init_flow */ +static void +restart_flow(ModeInfo * mi) { - dvector d; + flowstruct *sp; + int b; - d.x = -par.a * x - y/2 - y * y * y/8 + par.b * cos(z); - d.y = 2*x; - d.z = par.c; - return d; + if (flows == NULL) + return; + sp = &flows[MI_SCREEN(mi)]; + sp->count = 0; + + /* Re-Initialize point positions, velocities, etc. */ + for (b = 0; b < sp->beecount; b++) { + X(0, b) = Gauss_Rand(sp->range.x); + Y(0, b) = (sp->yperiod > 0)? + balance_rand(sp->range.y) : Gauss_Rand(sp->range.y); + Z(0, b) = Gauss_Rand(sp->range.z); + } } -void init_clip(flowstruct *sp); +/* Returns true if line was behind a clip plane, or it clips the line */ +/* nx,ny,nz is the normal to the plane. d is the distance from the origin */ +/* s and e are the end points of the line to be clipped */ +static int +clip(double nx, double ny, double nz, double d, dvector *s, dvector *e) +{ + int front1, front2; + dvector w, p; + double t; + + front1 = (nx*s->x + ny*s->y + nz*s->z >= -d); + front2 = (nx*e->x + ny*e->y + nz*e->z >= -d); + if (!front1 && !front2) return 1; + if (front1 && front2) return 0; + w.x = e->x - s->x; + w.y = e->y - s->y; + w.z = e->z - s->z; + + /* Find t in line equation */ + t = ( -d - nx*s->x - ny*s->y - nz*s->z) / ( nx*w.x + ny*w.y + nz*w.z); + + p.x = s->x + w.x * t; + p.y = s->y + w.y * t; + p.z = s->z + w.z * t; + + /* Move clipped point to the intersection */ + if (front2) { + *s = p; + } else { + *e = p; + } + return 0; +} -void -init_flow(ModeInfo * mi) +/* + * Public functions + */ + +ENTRYPOINT void +init_flow (ModeInfo * mi) { flowstruct *sp; - int b; - double beemult = 1; - static int allocated = 0; + char *name; + + MI_INIT (mi, flows, free_flow); + sp = &flows[MI_SCREEN(mi)]; + + sp->count2 = 0; - if (flows == NULL) { - if ((flows = (flowstruct *) calloc(MI_NUM_SCREENS(mi), - sizeof (flowstruct))) == NULL) - return; + sp->taillen = MI_SIZE(mi); + if (sp->taillen < -MINTRAIL) { + /* Change by sqrt so it seems more variable */ + sp->taillen = NRAND((int)sqrt((double) (-sp->taillen - MINTRAIL + 1))); + sp->taillen = sp->taillen * sp->taillen + MINTRAIL; + } else if (sp->taillen < MINTRAIL) { + sp->taillen = MINTRAIL; } - sp = &flows[MI_SCREEN(mi)]; - sp->count = 0; - sp->slow = 0.999; - sp->slow_view = 0.90; - - sp->width = MI_WIDTH(mi); - sp->height = MI_HEIGHT(mi); - - sp->tumble.theta = balance_rand(M_PI); - sp->tumble.phi = balance_rand(M_PI); - sp->tumble.dtheta = 0.002; - sp->tumble.dphi = 0.001; - sp->view.height = 0; - sp->view.depth = 0; /* no perspective view */ - sp->mode = 0; - if (get_boolean_resource ("rotate", "Boolean")) sp->mode |= FLOW_ROTATE; - if (get_boolean_resource ("ride", "Boolean")) sp->mode |= FLOW_RIDE; - if (get_boolean_resource ("zoom", "Boolean")) sp->mode |= FLOW_ZOOM; - if (get_boolean_resource ("allow2d", "Boolean")) sp->mode |= FLOW_2D; - if (get_boolean_resource ("slow", "Boolean")) sp->mode |= FLOW_SLOW; - if (get_boolean_resource ("freeze", "Boolean")) sp->mode |= FLOW_FREEZE; - if (get_boolean_resource ("box", "Boolean")) sp->mode |= FLOW_BOX; - - b = (sp->mode & FLOW_2D) ? 5 : 3; - b = NRAND(b); - - /* If more than one of rotate, ride and zoom are set, choose one */ - if (b < 3) { - int num = 0, modes[3]; - - if (sp->mode & FLOW_ROTATE) modes[num++] = FLOW_ROTATE; - if (sp->mode & FLOW_RIDE) modes[num++] = FLOW_RIDE; - if (sp->mode & FLOW_ZOOM) modes[num++] = FLOW_ZOOM; - - sp->mode &= ~(FLOW_ROTATE | FLOW_RIDE | FLOW_ZOOM); - - if (num) sp->mode |= modes[ NRAND(num) ]; - else sp->mode |= FLOW_ZOOM; + if(!rotatep && !ridep) rotatep = True; /* We need at least one viewpoint */ + + /* Start camera at Orbit or Bee */ + if(rotatep) { + sp->chaseto = ORBIT; + } else { + sp->chaseto = BEE; } - - switch (b) { + sp->chasetime = 1; /* Go directly to target */ + + sp->lyap = 0; + sp->yperiod = 0; + sp->step2 = INITIALSTEP; + + /* Zero parameter set */ + memset(sp->par2, 0, N_PARS * sizeof(dvector)); + + /* Set up standard examples */ + switch (NRAND((periodicp) ? 5 : 3)) { case 0: - sp->view.depth = 10; - sp->view.height = 0.2; - beemult = 3; - sp->ODE = Lorentz; - sp->step = 0.02; - sp->size = 60; - sp->centre.x = 0; - sp->centre.y = 0; - sp->centre.z = 24; - sp->range.x = 5; - sp->range.y = 5; - sp->range.z = 1; - sp->par.a = 10 + balance_rand(5); - sp->par.b = 28 + balance_rand(5); - sp->par.c = 2 + balance_rand(1); + /* + x' = a(y - x) + y' = x(b - z) - y + z' = xy - cz + */ + name = "Lorentz"; + sp->par2[Y].x = 10 + balance_rand(5*0); /* a */ + sp->par2[X].x = - sp->par2[Y].x; /* -a */ + sp->par2[X].y = 28 + balance_rand(5*0); /* b */ + sp->par2[XZ].y = -1; + sp->par2[Y].y = -1; + sp->par2[XY].z = 1; + sp->par2[Z].z = - 2 + balance_rand(1*0); /* -c */ break; case 1: - sp->view.depth = 10; - sp->view.height = 0.1; - beemult = 4; - sp->ODE = Rossler; - sp->step = 0.05; - sp->size = 24; - sp->centre.x = 0; - sp->centre.y = 0; - sp->centre.z = 3; - sp->range.x = 4; - sp->range.y = 4; - sp->range.z = 7; - sp->par.a = 2 + balance_rand(1); - sp->par.b = 0.2 + balance_rand(0.1); - sp->par.c = 0.2 + balance_rand(0.1); + /* + x' = -(y + az) + y' = x + by + z' = c + z(x - 5.7) + */ + name = "Rossler"; + sp->par2[Y].x = -1; + sp->par2[Z].x = -2 + balance_rand(1); /* a */ + sp->par2[X].y = 1; + sp->par2[Y].y = 0.2 + balance_rand(0.1); /* b */ + sp->par2[C].z = 0.2 + balance_rand(0.1); /* c */ + sp->par2[XZ].z = 1; + sp->par2[Z].z = -5.7; break; - case 2: - sp->view.depth = 10; - sp->view.height = 0.1; - beemult = 3; - sp->ODE = RosslerCone; - sp->step = 0.05; - sp->size = 24; - sp->centre.x = 0; - sp->centre.y = 0; - sp->centre.z = 3; - sp->range.x = 4; - sp->range.y = 4; - sp->range.z = 4; - sp->par.a = 2; - sp->par.b = 0.2; - sp->par.c = 0.25 + balance_rand(0.09); + case 2: + /* + x' = -(y + az) + y' = x + by - cz^2 + z' = 0.2 + z(x - 5.7) + */ + name = "RosslerCone"; + sp->par2[Y].x = -1; + sp->par2[Z].x = -2; /* a */ + sp->par2[X].y = 1; + sp->par2[Y].y = 0.2; /* b */ + sp->par2[ZZ].y = -0.331 + balance_rand(0.01); /* c */ + sp->par2[C].z = 0.2; + sp->par2[XZ].z = 1; + sp->par2[Z].z = -5.7; break; case 3: - sp->ODE = Birkhoff; - sp->step = 0.04; - sp->size = 2.6; - sp->centre.x = 0; - sp->centre.y = 0; - sp->centre.z = 0; - sp->range.x = 3; - sp->range.y = 4; - sp->range.z = 0; - sp->par.a = 10 + balance_rand(5); - sp->par.b = 0.35 + balance_rand(0.25); - sp->par.c = 1.57; - sp->tumble.theta = 0; - sp->tumble.phi = 0; - sp->tumble.dtheta = 0; - sp->tumble.dphi = 0; + /* + x' = -z + b sin(y) + y' = c + z' = 0.7x + az(0.1 - x^2) + */ + name = "Birkhoff"; + sp->par2[Z].x = -1; + sp->par2[SINY].x = 0.35 + balance_rand(0.25); /* b */ + sp->par2[C].y = 1.57; /* c */ + sp->par2[X].z = 0.7; + sp->par2[Z].z = 1 + balance_rand(0.5); /* a/10 */ + sp->par2[XXZ].z = -10 * sp->par2[Z].z; /* -a */ + sp->yperiod = 2 * M_PI; break; - case 4: default: - sp->ODE = Duffing; - sp->step = 0.02; - sp->size = 30; - sp->centre.x = 0; - sp->centre.y = 0; - sp->centre.z = 0; - sp->range.x = 20; - sp->range.y = 20; - sp->range.z = 0; - sp->par.a = 0.2 + balance_rand(0.1); - sp->par.b = 27.0 + balance_rand(3.0); - sp->par.c = 1.33; - sp->tumble.theta = 0; - sp->tumble.phi = 0; - sp->tumble.dtheta = -NRAND(2)*sp->par.c*sp->step; - sp->tumble.dphi = 0; - beemult = 0.5; + /* + x' = -ax - z/2 - z^3/8 + b sin(y) + y' = c + z' = 2x + */ + name = "Duffing"; + sp->par2[X].x = -0.2 + balance_rand(0.1); /* a */ + sp->par2[Z].x = -0.5; + sp->par2[ZZZ].x = -0.125; + sp->par2[SINY].x = 27.0 + balance_rand(3.0); /* b */ + sp->par2[C].y = 1.33; /* c */ + sp->par2[X].z = 2; + sp->yperiod = 2 * M_PI; break; + } - sp->view.depth *= 4; + sp->range.x = 5; + sp->range.z = 5; - sp->beecount = beemult * MI_COUNT(mi); - if (sp->beecount < 0) /* random variations */ - sp->beecount = NRAND(-sp->beecount) + 1; /* Minimum 1 */ + if(sp->yperiod > 0) { + sp->ODE = Periodic; + /* periodic flows show either uniform distribution or a + snapshot on the 'time' axis */ + sp->range.y = NRAND(2)? sp->yperiod : 0; + } else { + sp->range.y = 5; + sp->ODE = Cubic; + } - /* Clear the background. */ - MI_CLEARWINDOW(mi); + /* Run discoverer to set up bounding box, etc. Lyapunov will + probably be innaccurate, since we're only running it once, but + we're using known strange attractors so it should be ok. */ + discover(mi); + if(MI_IS_VERBOSE(mi)) + fprintf(stdout, + "flow: Lyapunov exponent: %g, step: %g, size: %g (%s)\n", + sp->lyap2, sp->step2, sp->size2, name); + /* Install new params */ + sp->lyap = sp->lyap2; + sp->size = sp->size2; + sp->mid = sp->mid2; + sp->step = sp->step2; + memcpy(sp->par, sp->par2, sizeof(sp->par2)); + + sp->count2 = 0; /* Reset search */ + + sp->beecount = MI_COUNT(mi); + if (!sp->beecount) { + sp->beecount = 1; /* The camera requires 1 or more */ + } else if (sp->beecount < 0) { /* random variations */ + sp->beecount = NRAND(-sp->beecount) + 1; /* Minimum 1 */ + } - if(!allocated || sp->beecount != allocated){ /* reallocate */ - if (sp->csegs != NULL) { - (void) free((void *) sp->csegs); - sp->csegs = NULL; - } - if (sp->cnsegs != NULL) { - (void) free((void *) sp->cnsegs); - sp->cnsegs = NULL; - } - if (sp->old_segs != NULL) { - (void) free((void *) sp->old_segs); - sp->old_segs = NULL; - } - if (sp->p[0] != NULL) { - (void) free((void *) sp->p[0]); - sp->p[0] = NULL; - } - if (sp->p[1] != NULL) { - (void) free((void *) sp->p[1]); - sp->p[1] = NULL; - } - allocated = sp->beecount; +# ifdef HAVE_JWXYZ /* Don't second-guess Quartz's double-buffering */ + dbufp = False; +# endif + + if(dbufp) { /* Set up double buffer */ + if (sp->buffer != None) + XFreePixmap(MI_DISPLAY(mi), sp->buffer); + sp->buffer = XCreatePixmap(MI_DISPLAY(mi), MI_WINDOW(mi), + MI_WIDTH(mi), MI_HEIGHT(mi), MI_DEPTH(mi)); + } else { + sp->buffer = MI_WINDOW(mi); } + /* no "NoExpose" events from XCopyArea wanted */ + XSetGraphicsExposures(MI_DISPLAY(mi), MI_GC(mi), False); - /* Allocate memory. */ + /* Make sure we're using 'thin' lines */ + XSetLineAttributes(MI_DISPLAY(mi), MI_GC(mi), 0, LineSolid, CapNotLast, + JoinMiter); - if (!sp->csegs) { - sp->csegs = (XSegment *) malloc(sizeof (XSegment) * sp->beecount - * MI_NPIXELS(mi)); - sp->cnsegs = (int *) malloc(sizeof (int) * MI_NPIXELS(mi)); + /* Clear the background (may be slow depending on user prefs). */ + MI_CLEARWINDOW(mi); - sp->old_segs = (XSegment *) malloc(sizeof (XSegment) * sp->beecount); - sp->p[0] = (dvector *) malloc(sizeof (dvector) * sp->beecount); - sp->p[1] = (dvector *) malloc(sizeof (dvector) * sp->beecount); + /* Allocate memory. */ + if (sp->csegs == NULL) { + allocate(sp->csegs, XSegment, + (sp->beecount + BOX_L) * MI_NPIXELS(mi) * sp->taillen); + allocate(sp->cnsegs, int, MI_NPIXELS(mi)); + allocate(sp->old_segs, XSegment, (sp->beecount + BOX_L) * sp->taillen); + allocate(sp->p, dvector, sp->beecount * sp->taillen); } /* Initialize point positions, velocities, etc. */ + restart_flow(mi); - for (b = 0; b < sp->beecount; b++) { - X(1, b) = X(0, b) = balance_rand(sp->range.x); - Y(1, b) = Y(0, b) = balance_rand(sp->range.y); - Z(1, b) = Z(0, b) = balance_rand(sp->range.z); - } - - init_clip(sp); - + /* Set up camera tail */ + X(1, 0) = sp->cam[1].x = 0; + Y(1, 0) = sp->cam[1].y = 0; + Z(1, 0) = sp->cam[1].z = 0; } -/* Clipping planes */ -#define PLANES 5 -static double plane_orig[][2][3] = { - /* X goes into screen, Y goes right, Z goes down(up?) */ - /* {Normal}, {Point} */ - { {1.0, 0, 0}, {0.01, 0, 0} }, - { {1.0, 1.0, 0.0}, {0, 0, 0} }, - { {1.0,-1.0, 0.0}, {0, 0, 0} }, - { {1.0, 0.0, 1.0}, {0, 0, 0} }, - { {1.0, 0.0,-1.0}, {0, 0, 0} } -}; -static double plane[PLANES][2][3]; -static double plane_d[PLANES]; +ENTRYPOINT void +draw_flow (ModeInfo * mi) +{ + int b, i; + int col, begin, end; + double M[3][3]; /* transformation matrix */ + flowstruct *sp = NULL; + int swarm = 0; + int segindex; -#define BOX_P 32 -#define BOX_L 36 -#define MIN_BOX (3) -#define MAX_BOX (MIN_BOX + BOX_L) -/* Points that make up the box (normalized coordinates) */ -static double box_orig[][3] = { - {1,1,1}, /* 0 */ - {1,1,-1}, /* 1 */ - {1,-1,-1}, /* 2 */ - {1,-1,1}, /* 3 */ - {-1,1,1}, /* 4 */ - {-1,1,-1}, /* 5 */ - {-1,-1,-1},/* 6 */ - {-1,-1,1}, /* 7 */ - {1, .8, .8}, - {1, .8,-.8}, - {1,-.8,-.8}, - {1,-.8, .8}, - { .8,1, .8}, - { .8,1,-.8}, - {-.8,1,-.8}, - {-.8,1, .8}, - { .8, .8,1}, - { .8,-.8,1}, - {-.8,-.8,1}, - {-.8, .8,1}, - {-1, .8, .8}, - {-1, .8,-.8}, - {-1,-.8,-.8}, - {-1,-.8, .8}, - { .8,-1, .8}, - { .8,-1,-.8}, - {-.8,-1,-.8}, - {-.8,-1, .8}, - { .8, .8,-1}, - { .8,-.8,-1}, - {-.8,-.8,-1}, - {-.8, .8,-1} -}; + if (flows == NULL) + return; + sp = &flows[MI_SCREEN(mi)]; + if (sp->csegs == NULL) + return; -/* Container for scaled box points */ -static double box[BOX_P][3]; +#ifdef HAVE_JWXYZ /* Don't second-guess Quartz's double-buffering */ + XClearWindow (MI_DISPLAY(mi), MI_WINDOW(mi)); +#endif -/* Lines connecting the box dots */ -static double lines[][2] = { - {0,1}, {1,2}, {2,3}, {3,0}, /* box */ - {4,5}, {5,6}, {6,7}, {7,4}, - {0,4}, {1,5}, {2,6}, {3,7}, - {4+4,5+4}, {5+4,6+4}, {6+4,7+4}, {7+4,4+4}, - {4+8,5+8}, {5+8,6+8}, {6+8,7+8}, {7+8,4+8}, - {4+12,5+12}, {5+12,6+12}, {6+12,7+12}, {7+12,4+12}, - {4+16,5+16}, {5+16,6+16}, {6+16,7+16}, {7+16,4+16}, - {4+20,5+20}, {5+20,6+20}, {6+20,7+20}, {7+20,4+20}, - {4+24,5+24}, {5+24,6+24}, {6+24,7+24}, {7+24,4+24}, -}; + /* multiplier for indexing segment arrays. Used in IX macro, etc. */ + segindex = (sp->beecount + BOX_L) * sp->taillen; + + if(searchp){ + if(sp->count2 == 0) { /* start new search */ + sp->step2 = INITIALSTEP; + /* Pick random parameters. Actual range is irrelevant + since parameter scale determines flow speed but not + structure. */ + for(i=0; i< N_PARS; i++) { + sp->par2[i].x = Gauss_Rand(1.0); + sp->par2[i].y = Gauss_Rand(1.0); + sp->par2[i].z = Gauss_Rand(1.0); + } + } + if(!discover(mi)) { /* Flow exploded, reset. */ + sp->count2 = 0; + } else { + if(sp->lyap2 < 0) { + sp->count2 = 0; /* Attractor found, but it's not strange */ + }else if(sp->count2 > 1000000) { /* This one will do */ + sp->count2 = 0; /* Reset search */ + if(MI_IS_VERBOSE(mi)) + fprintf(stdout, + "flow: Lyapunov exponent: %g, step: %g, size: %g (unnamed)\n", + sp->lyap2, sp->step2, sp->size2); + /* Install new params */ + sp->lyap = sp->lyap2; + sp->size = sp->size2; + sp->mid = sp->mid2; + sp->step = sp->step2; + memcpy(sp->par, sp->par2, sizeof(sp->par2)); + + /* If we're allowed to zoom out, do so now, so that we + get a look at the new attractor. */ + if(sp->chaseto == BEE && rotatep) { + sp->chaseto = ORBIT; + sp->chasetime = 100; + } + /* Reset initial conditions, so we don't get + misleading artifacts in the particle density. */ + restart_flow(mi); + } + } + } -/* Boundaries of bees */ -double xmin, xmax; -double ymin, ymax; -double zmin, zmax; + /* Reset segment buffers */ + for (col = 0; col < MI_NPIXELS(mi); col++) + sp->cnsegs[col] = 0; -void init_clip(flowstruct *sp) -{ - int i; - - /* Scale the planes to the screen. I had to invert the projection - * algorithms so that when projected they would be right at the edge of the - * screen. */ - - /* #### jwz: I'm not really sure what it means when sp->view.depth is 0 - in here -- what's the right thing to do? */ - - double width = (sp->view.depth - ? sp->size/sp->view.depth/2 - : 1); - double height = (sp->view.depth - ? (sp->size/sp->view.depth/2* - sp->view.height/sp->view.height) - : 1); - for (i = 0; i < PLANES; i++) { - /* Copy orig planes into planes, expanding <-> clippings */ - plane[i][0][0] = plane_orig[i][0][0]; - plane[i][0][1] = plane_orig[i][0][1] / width; - plane[i][0][2] = plane_orig[i][0][2] / height; - plane[i][1][0] = plane_orig[i][1][0]; - plane[i][1][1] = plane_orig[i][1][1]; - plane[i][1][2] = plane_orig[i][1][2]; - - /* Calculate the 'd' part of 'ax + by + cz = d' */ - plane_d[i] = - plane[i][0][0] * plane[i][1][0]; - plane_d[i] -= plane[i][0][1] * plane[i][1][1]; - plane_d[i] -= plane[i][0][2] * plane[i][1][2]; - } - xmin = X(0, i); xmax = X(0,i); - ymin = Y(0, i); ymax = Y(0,i); - zmin = Z(0, i); zmax = Z(0,i); -} + MI_IS_DRAWN(mi) = True; -/* Scale the box defined above to fit around all points */ -void create_box(flowstruct *sp) -{ - int i = MAX_BOX; - double xmid, ymid, zmid; - double xsize, ysize, zsize; - double size; - - /* Count every 5th point for speed.. */ - for (; i < sp->beecount; i += 5) { - if ( X(0,i) < xmin ) xmin = X(0, i); - else if ( X(0,i) > xmax ) xmax = X(0, i); - if ( Y(0,i) < ymin ) ymin = Y(0, i); - else if ( Y(0,i) > ymax ) ymax = Y(0, i); - if ( Z(0,i) < zmin ) zmin = Z(0, i); - else if ( Z(0,i) > zmax ) zmax = Z(0, i); - } - xmid = (xmax+xmin)/2; - ymid = (ymax+ymin)/2; - zmid = (zmax+zmin)/2; - xsize = xmax - xmin; - ysize = ymax - ymin; - zsize = zmax - zmin; - size = xsize; - if (ysize> size) size = ysize; - if (zsize > size) size = zsize; - size /= 2; - - /* Scale box */ - for (i = 0; i < BOX_P; i++) { - box[i][0] = box_orig[i][0] * size + xmid; - box[i][1] = box_orig[i][1] * size + ymid; - box[i][2] = box_orig[i][2] * size + zmid; + /* Calculate circling POV [Chalky]*/ + sp->circle[1] = sp->circle[0]; + sp->circle[0].x = sp->size * 2 * sin(sp->count / 100.0) * + (-0.6 + 0.4 *cos(sp->count / 500.0)) + sp->mid.x; + sp->circle[0].y = sp->size * 2 * cos(sp->count / 100.0) * + (0.6 + 0.4 *cos(sp->count / 500.0)) + sp->mid.y; + sp->circle[0].z = sp->size * 2 * sin(sp->count / 421.0) + sp->mid.z; + + /* Timed chase instead of Chalkie's Bistable oscillator [TDA] */ + if(rotatep && ridep) { + if(sp->chaseto == BEE && NRAND(1000) == 0){ + sp->chaseto = ORBIT; + sp->chasetime = 100; + }else if(NRAND(4000) == 0){ + sp->chaseto = BEE; + sp->chasetime = 100; + } } -} - -/* Returns true if point is infront of the plane (rather than behind) */ -int infront_of(double x, double y, double z, int i) -{ - double sum = plane[i][0][0]*x + plane[i][0][1]*y + plane[i][0][2]*z + plane_d[i]; - return sum >= 0.0; -} - -/* Returns true if line was behind a clip plane, or clips the line */ -int clip(double *x1, double *y1, double *z1, double *x2, double *y2, double *z2) -{ - int i; - for (i = 0; i < PLANES; i++) { - double t; - double x, y, z; /* Intersection point */ - double dx, dy, dz; /* line delta */ - int front1, front2; - front1 = infront_of(*x1, *y1, *z1, i); - front2 = infront_of(*x2, *y2, *z2, i); - if (!front1 && !front2) return 1; - if (front1 && front2) continue; - - dx = *x2 - *x1; - dy = *y2 - *y1; - dz = *z2 - *z1; - - /* Find t in line equation */ - t = ( plane_d[i] - - plane[i][0][0]*(*x1) - plane[i][0][1]*(*y1) - plane[i][0][2]*(*z1) ) - / - ( plane[i][0][0]*dx + plane[i][0][1]*dy + plane[i][0][2]*dz ); - - x = *x1 + dx * t; - y = *y1 + dy * t; - z = *z1 + dz * t; - /* Make point that was behind to be the intersect */ - if (front2) { - *x1 = x; - *y1 = y; - *z1 = z; + /* Set up orientation matrix */ + { + double x[3], p[3], x2=0, xp=0; + int j; + + /* Chasetime is here to guarantee the camera makes it all the + way to the target in a finite number of steps. */ + if(sp->chasetime > 1) + sp->chasetime--; + + if(sp->chaseto == BEE){ + /* Camera Head targets bee 0 */ + sp->cam[0].x += (X(0, 0) - sp->cam[0].x)/sp->chasetime; + sp->cam[0].y += (Y(0, 0) - sp->cam[0].y)/sp->chasetime; + sp->cam[0].z += (Z(0, 0) - sp->cam[0].z)/sp->chasetime; + + /* Camera Tail targets previous position of bee 0 */ + sp->cam[1].x += (X(1, 0) - sp->cam[1].x)/sp->chasetime; + sp->cam[1].y += (Y(1, 0) - sp->cam[1].y)/sp->chasetime; + sp->cam[1].z += (Z(1, 0) - sp->cam[1].z)/sp->chasetime; + + /* Camera Wing targets bee 1 */ + sp->cam[2].x += (X(0, 1) - sp->cam[2].x)/sp->chasetime; + sp->cam[2].y += (Y(0, 1) - sp->cam[2].y)/sp->chasetime; + sp->cam[2].z += (Z(0, 1) - sp->cam[2].z)/sp->chasetime; } else { - *x2 = x; - *y2 = y; - *z2 = z; + /* Camera Head targets Orbiter */ + sp->cam[0].x += (sp->circle[0].x - sp->cam[0].x)/sp->chasetime; + sp->cam[0].y += (sp->circle[0].y - sp->cam[0].y)/sp->chasetime; + sp->cam[0].z += (sp->circle[0].z - sp->cam[0].z)/sp->chasetime; + + /* Camera Tail targets diametrically opposite the middle + of the bounding box from the Orbiter */ + sp->cam[1].x += + (2*sp->circle[0].x - sp->mid.x - sp->cam[1].x)/sp->chasetime; + sp->cam[1].y += + (2*sp->circle[0].y - sp->mid.y - sp->cam[1].y)/sp->chasetime; + sp->cam[1].z += + (2*sp->circle[0].z - sp->mid.z - sp->cam[1].z)/sp->chasetime; + /* Camera Wing targets previous position of Orbiter */ + sp->cam[2].x += (sp->circle[1].x - sp->cam[2].x)/sp->chasetime; + sp->cam[2].y += (sp->circle[1].y - sp->cam[2].y)/sp->chasetime; + sp->cam[2].z += (sp->circle[1].z - sp->cam[2].z)/sp->chasetime; } - } - return 0; -} + + /* Viewpoint from Tail of camera */ + sp->centre.x=sp->cam[1].x; + sp->centre.y=sp->cam[1].y; + sp->centre.z=sp->cam[1].z; + + /* forward vector */ + x[0] = sp->cam[0].x - sp->cam[1].x; + x[1] = sp->cam[0].y - sp->cam[1].y; + x[2] = sp->cam[0].z - sp->cam[1].z; + + /* side */ + p[0] = sp->cam[2].x - sp->cam[1].x; + p[1] = sp->cam[2].y - sp->cam[1].y; + p[2] = sp->cam[2].z - sp->cam[1].z; -void -draw_flow(ModeInfo * mi) -{ - Display *display = MI_DISPLAY(mi); - Window window = MI_WINDOW(mi); - GC gc = MI_GC(mi); - flowstruct *sp = &flows[MI_SCREEN(mi)]; - int b, c, i; - int col, ix; - int new_view = 0; - double M[3][3]; /* transformation matrix */ - double step_view = sp->step; - double step_bees = sp->step; - double step_slow = sp->step; - double pp, pc; - - create_box(sp); - - if(!sp->view.depth){ /* simple 3D tumble */ - double sint, cost, sinp, cosp; - sp->tumble.theta += sp->tumble.dtheta; - sp->tumble.phi += sp->tumble.dphi; - sint = sin(sp->tumble.theta); - cost = cos(sp->tumble.theta); - sinp = sin(sp->tumble.phi); - cosp = cos(sp->tumble.phi); - M[0][0]= cost; M[0][1]=-sint*cosp; M[0][2]= sint*sinp; - M[1][0]= sint; M[1][1]= cost*cosp; M[1][2]=-cost*sinp; - M[2][0]= 0; M[2][1]= 0; M[2][2]= 1; - } else { /* initialize matrix */ - M[0][0]= 0; M[0][1]= 0; M[0][2]= 0; - M[1][0]= 0; M[1][1]= 0; M[1][2]= 0; - M[2][0]= 0; M[2][1]= 0; M[2][2]= 0; + /* So long as X and P don't collide, these can be used to form + three mutually othogonal axes: X, (X x P) x X and X x P. + After being normalised to unit length, these form the + Orientation Matrix. */ + + for(i=0; i<3; i++){ + x2+= x[i]*x[i]; /* X . X */ + xp+= x[i]*p[i]; /* X . P */ + M[0][i] = x[i]; /* X */ + } + + for(i=0; i<3; i++) /* (X x P) x X */ + M[1][i] = x2*p[i] - xp*x[i]; /* == (X . X) P - (X . P) X */ + + M[2][0] = x[1]*p[2] - x[2]*p[1]; /* X x P */ + M[2][1] = -x[0]*p[2] + x[2]*p[0]; + M[2][2] = x[0]*p[1] - x[1]*p[0]; + + /* normalise axes */ + for(j=0; j<3; j++){ + double A=0; + for(i=0; i<3; i++) A+=M[j][i]*M[j][i]; /* sum squares */ + A=sqrt(A); + if(A>0) + for(i=0; i<3; i++) M[j][i]/=A; + } + if(sp->chaseto == BEE) { + X(0, 1)=X(0, 0)+M[1][0]*sp->step; /* adjust neighbour */ + Y(0, 1)=Y(0, 0)+M[1][1]*sp->step; + Z(0, 1)=Z(0, 0)+M[1][2]*sp->step; + } } - for (col = 0; col < MI_NPIXELS(mi); col++) - sp->cnsegs[col] = 0; - - MI_IS_DRAWN(mi) = True; - - /* Calculate circling POV */ - sp->circle[1] = sp->circle[0]; - sp->circle[0].x = sp->size * 2 * sin(sp->count / 40.0) * (0.6 + 0.4 *cos(sp->count / 100.0)); - sp->circle[0].y = sp->size * 2 * cos(sp->count / 40.0) * (0.6 + 0.4 *cos(sp->count / 100.0)); - sp->circle[0].z = sp->size * 2 * sin(sp->count / 421.0); - - if (sp->mode & FLOW_ROTATE) - pp = 0; - else if (sp->mode & FLOW_RIDE) - pp = 1; - else /* ZOOM */ - /* Bistable oscillator to switch between the trained bee and the circler */ - pp = -sin(sin(sin(cos(sp->count / 150.0)*M_PI/2)*M_PI/2)*M_PI/2) *0.5 + 0.5; - pc = 1 - pp; - - - /* Slow down or speed up the bees / view: */ - - /* exponentially accelerate towards zero */ - sp->slow = sp->slow * 1.005 - 0.005; - if (sp->slow < 0) sp->slow = 0; - - sp->slow_view = sp->slow_view * 1.005 - 0.005; - if (sp->slow_view < 0) sp->slow_view = 0; - - /* View speeds up, slow bees slow to half speed, and other bees will - * actually stop */ - step_view = step_view * (1.01 - sp->slow_view * sp->slow_view) * 0.2; - step_slow = step_slow * (sp->slow + 0.5) / 2; - if (sp->mode & FLOW_SLOW) - step_bees = step_bees * sp->slow; - else - step_bees = step_slow; + /* <=- Bounding Box -=> */ + if(boxp) { + for (b = 0; b < BOX_L; b++) { + + /* Chalky's clipping code, Only used for the box */ + /* clipping trails is slow and of little benefit. [TDA] */ + int p1 = lines[b][0]; + int p2 = lines[b][1]; + dvector A1, A2; + double x1=box[p1][0]* sp->size/2 + sp->mid.x - sp->centre.x; + double y1=box[p1][1]* sp->size/2 + sp->mid.y - sp->centre.y; + double z1=box[p1][2]* sp->size/2 + sp->mid.z - sp->centre.z; + double x2=box[p2][0]* sp->size/2 + sp->mid.x - sp->centre.x; + double y2=box[p2][1]* sp->size/2 + sp->mid.y - sp->centre.y; + double z2=box[p2][2]* sp->size/2 + sp->mid.z - sp->centre.z; + + A1.x=M[0][0]*x1 + M[0][1]*y1 + M[0][2]*z1; + A1.y=M[1][0]*x1 + M[1][1]*y1 + M[1][2]*z1; + A1.z=M[2][0]*x1 + M[2][1]*y1 + M[2][2]*z1 + EYEHEIGHT * sp->size; + A2.x=M[0][0]*x2 + M[0][1]*y2 + M[0][2]*z2; + A2.y=M[1][0]*x2 + M[1][1]*y2 + M[1][2]*z2; + A2.z=M[2][0]*x2 + M[2][1]*y2 + M[2][2]*z2 + EYEHEIGHT * sp->size; + + /* Clip in 3D before projecting down to 2D. A 2D clip + after projection wouldn't be able to handle lines that + cross x=0 */ + if (clip(1, 0, 0,-1, &A1, &A2) || /* Screen */ + clip(1, 2, 0, 0, &A1, &A2) || /* Left */ + clip(1,-2, 0, 0, &A1, &A2) || /* Right */ + clip(1,0, 2.0*MI_WIDTH(mi)/MI_HEIGHT(mi), 0, &A1, &A2)||/*UP*/ + clip(1,0,-2.0*MI_WIDTH(mi)/MI_HEIGHT(mi), 0, &A1, &A2))/*Down*/ + continue; + + /* Colour according to bee */ + col = b % (MI_NPIXELS(mi) - 1); + + sp->csegs[IX(col)].x1 = MI_WIDTH(mi)/2 + MI_WIDTH(mi) * A1.y/A1.x; + sp->csegs[IX(col)].y1 = MI_HEIGHT(mi)/2 + MI_WIDTH(mi) * A1.z/A1.x; + sp->csegs[IX(col)].x2 = MI_WIDTH(mi)/2 + MI_WIDTH(mi) * A2.y/A2.x; + sp->csegs[IX(col)].y2 = MI_HEIGHT(mi)/2 + MI_WIDTH(mi) * A2.z/A2.x; + sp->cnsegs[col]++; + } + } /* <=- Bees -=> */ for (b = 0; b < sp->beecount; b++) { - /* Calc if this bee is slow. Note normal bees are exempt from - * calculations once they slow to half speed, so that they remain as - * frozen lines rather than barely-visible points */ - int slow = ((b & 0x7) == 0); - if ( !(sp->mode & FLOW_FREEZE) ) slow = 1; - /* Age the arrays. */ - if (b < 2 || sp->slow > 0.5 || slow) { - X(1, b) = X(0, b); - Y(1, b) = Y(0, b); - Z(1, b) = Z(0, b); - - /* 2nd order Kunge Kutta */ - { - dvector k1, k2; - double step; - - if (b == 0 || b == 1) { - step = step_view; - } else if (slow) { - step = step_slow; - } else { - step = step_bees; - } - k1 = sp->ODE(sp->par, X(1, b), Y(1, b), Z(1, b)); - k1.x *= step; - k1.y *= step; - k1.z *= step; - k2 = sp->ODE(sp->par, X(1, b) + k1.x, Y(1, b) + k1.y, Z(1, b) + k1.z); - k2.x *= step; - k2.y *= step; - k2.z *= step; - X(0, b) = X(1, b) + (k1.x + k2.x) / 2.0; - Y(0, b) = Y(1, b) + (k1.y + k2.y) / 2.0; - Z(0, b) = Z(1, b) + (k1.z + k2.z) / 2.0; + if(fabs(X(0, b)) > LOST_IN_SPACE || + fabs(Y(0, b)) > LOST_IN_SPACE || + fabs(Z(0, b)) > LOST_IN_SPACE){ + if(sp->chaseto == BEE && b == 0){ + /* Lost camera bee. Need to replace it since + rerunning init_flow could lose us a hard-won new + attractor. Try moving it very close to a random + other bee. This way we have a good chance of being + close to the attractor and not forming a false + artifact. If we've lost many bees this may need to + be repeated. */ + /* Don't worry about camera wingbee. It stays close + to the main camera bee no matter what happens. */ + int newb = 1 + NRAND(sp->beecount - 1); + X(0, 0) = X(0, newb) + 0.001; + Y(0, 0) = Y(0, newb); + Z(0, 0) = Z(0, newb); + if(MI_IS_VERBOSE(mi)) + fprintf(stdout, + "flow: resetting lost camera near bee %d\n", + newb); } + continue; } + /* Age the tail. It's critical this be fast since + beecount*taillen can be large. */ + memmove(B(1, b), B(0, b), (sp->taillen - 1) * sizeof(dvector)); + + Iterate(B(0,b), sp->ODE, sp->par, sp->step); + + /* Don't show wingbee since he's not quite in the flow. */ + if(sp->chaseto == BEE && b == 1) continue; /* Colour according to bee */ col = b % (MI_NPIXELS(mi) - 1); - ix = col * sp->beecount + sp->cnsegs[col]; - + /* Fill the segment lists. */ - - if(sp->view.depth) { /* perspective view has special points */ - if(b==0){ /* point of view */ - sp->centre.x = X(0, b) * pp + sp->circle[0].x * pc; - sp->centre.y = Y(0, b) * pp + sp->circle[0].y * pc; - sp->centre.z = Z(0, b) * pp + sp->circle[0].z * pc; - /*printf("center: (%3.3f,%3.3f,%3.3f)\n",sp->centre.x, sp->centre.y, sp->centre.z);*/ - }else if(b==1){ /* neighbour: used to compute local axes */ - double x[3], p[3], x2=0, xp=0, C[3][3]; + + begin = 0; /* begin new trail */ + end = MIN(sp->taillen, sp->count); /* short trails at first */ + for(i=0; i < end; i++){ + double x = X(i,b)-sp->centre.x; + double y = Y(i,b)*(sp->yperiod < 0? (sp->size/sp->yperiod) :1) + -sp->centre.y; + double z = Z(i,b)-sp->centre.z; + double XM=M[0][0]*x + M[0][1]*y + M[0][2]*z; + double YM=M[1][0]*x + M[1][1]*y + M[1][2]*z; + double ZM=M[2][0]*x + M[2][1]*y + M[2][2]*z + EYEHEIGHT * sp->size; + short absx, absy; + + swarm++; /* count the remaining bees */ + if(sp->yperiod > 0 && Y(i,b) > sp->yperiod){ int j; - - /* forward */ - x[0] = X(0, 0) - X(1, 0); - x[1] = Y(0, 0) - Y(1, 0); - x[2] = Z(0, 0) - Z(1, 0); - - /* neighbour */ - p[0] = X(0, 1) - X(1, 0); - p[1] = Y(0, 1) - Y(1, 0); - p[2] = Z(0, 1) - Z(1, 0); - - for(i=0; i<3; i++){ - x2+= x[i]*x[i]; /* X . X */ - xp+= x[i]*p[i]; /* X . P */ - M[0][i] = x[i]; /* X */ - } - - for(i=0; i<3; i++) /* (X x P) x X */ - M[1][i] = x2*p[i] - xp*x[i]; /* == (X . X) P - (X . P) X */ - - M[2][0] = x[1]*p[2] - x[2]*p[1]; /* X x P */ - M[2][1] = -x[0]*p[2] + x[2]*p[0]; - M[2][2] = x[0]*p[1] - x[1]*p[0]; - - /* normalise axes */ - for(j=0; j<3; j++){ - double A=0; - for(i=0; i<3; i++) A+=M[j][i]*M[j][i]; /* sum squares */ - A=sqrt(A); - for(i=0; i<3; i++) M[j][i]/=A; - } - - X(0, 1)=X(0, 0)+M[1][0]; /* adjust neighbour */ - Y(0, 1)=Y(0, 0)+M[1][1]; - Z(0, 1)=Z(0, 0)+M[1][2]; - - /* Look at trained bee into C matrix */ - /* forward */ - x[0] = 0 - sp->circle[0].x; - x[1] = 0 - sp->circle[0].y; - x[2] = 0 - sp->circle[0].z; + Y(i,b) -= sp->yperiod; + /* hide tail to prevent streaks in Y. Streaks in X,Z + are ok, they help to outline the Poincare' + slice. */ + for(j = i; j < end; j++) Y(j,b) = Y(i,b); + /*begin = i + 1;*/ + break; + } - /* neighbour */ - p[0] = sp->circle[0].x - sp->circle[1].x; - p[1] = sp->circle[0].y - sp->circle[1].y; - p[2] = sp->circle[0].z - sp->circle[1].z; - - for(i=0; i<3; i++){ - x2+= x[i]*x[i]; /* X . X */ - xp+= x[i]*p[i]; /* X . P */ - C[0][i] = x[i]; /* X */ - } - - for(i=0; i<3; i++) /* (X x P) x X */ - C[1][i] = x2*p[i] - xp*x[i]; /* == (X . X) P - (X . P) X */ - - C[2][0] = x[1]*p[2] - x[2]*p[1]; /* X x P */ - C[2][1] = -x[0]*p[2] + x[2]*p[0]; - C[2][2] = x[0]*p[1] - x[1]*p[0]; - - /* normalise axes */ - for(j=0; j<3; j++){ - double A=0; - for(i=0; i<3; i++) A+=C[j][i]*C[j][i]; /* sum squares */ - A=sqrt(A); - if (A != 0) /* #### is this right? */ - for(i=0; i<3; i++) C[j][i]/=A; - } - - /* Interpolate between Center and Trained Bee matrices */ - /* This isn't very accurate and leads to weird transformations - * (shearing, etc), but it works. Besides, sometimes they look - * cool :) */ - pp = pp * pp; /* Don't follow bee's direction until very close */ - pc = 1 - pp; - for (i = 0; i < 3; i++) - for (j = 0; j < 3; j++) - M[i][j] = M[i][j] * pp + C[i][j] * pc; - - -#if 0 /* display local axes for testing */ - X(1, b)=X(0, 0); - Y(1, b)=Y(0, 0); - Z(1, b)=Z(0, 0); - }else if(b==2){ - X(0, b)=X(0, 0)+0.5*M[0][0]; - Y(0, b)=Y(0, 0)+0.5*M[0][1]; - Z(0, b)=Z(0, 0)+0.5*M[0][2]; - X(1, b)=X(0, 0); - Y(1, b)=Y(0, 0); - Z(1, b)=Z(0, 0); - }else if(b==3){ - X(0, b)=X(0, 0)+1.5*M[2][0]; - Y(0, b)=Y(0, 0)+1.5*M[2][1]; - Z(0, b)=Z(0, 0)+1.5*M[2][2]; - X(1, b)=X(0, 0); - Y(1, b)=Y(0, 0); - Z(1, b)=Z(0, 0); -#endif - /* Draw a box... */ + if(XM <= 0){ /* off screen - new trail */ + begin = i + 1; + continue; } - } - if (b >= MIN_BOX && b < MAX_BOX) { - int p1 = lines[b-MIN_BOX][0]; - int p2 = lines[b-MIN_BOX][1]; - X(0, b) = box[p1][0]; Y(0, b) = box[p1][1]; Z(0, b) = box[p1][2]; - X(1, b) = box[p2][0]; Y(1, b) = box[p2][1]; Z(1, b) = box[p2][2]; + absx = MI_WIDTH(mi)/2 + MI_WIDTH(mi) * YM/XM; + absy = MI_HEIGHT(mi)/2 + MI_WIDTH(mi) * ZM/XM; + /* Performance bottleneck */ + if(absx <= 0 || absx >= MI_WIDTH(mi) || + absy <= 0 || absy >= MI_HEIGHT(mi)) { + /* off screen - new trail */ + begin = i + 1; + continue; } - - -#if 0 /* Original code */ - for(i=0; i<2; i++){ - double x=X(i,b)-sp->centre.x; - double y=Y(i,b)-sp->centre.y; - double z=Z(i,b)-sp->centre.z; - double X=M[0][0]*x + M[0][1]*y + M[0][2]*z; - double Y=M[1][0]*x + M[1][1]*y + M[1][2]*z; - double Z=M[2][0]*x + M[2][1]*y + M[2][2]*z+sp->view.height; - double absx, absy; - if(sp->view.depth){ - if(X <= 0) break; - absx=SCALE_X(sp->view.depth*Y/X); - absy=SCALE_Y(sp->view.depth*Z/X); - if(absx < -sp->width || absx > 2*sp->width || - absy < -sp->height || absy > 2*sp->height) - break; - }else{ - absx=SCALE_X(X); - absy=SCALE_Y(Y); + if(i > begin) { /* complete previous segment */ + sp->csegs[IX(col)].x2 = absx; + sp->csegs[IX(col)].y2 = absy; + sp->cnsegs[col]++; } - if(i){ - sp->csegs[ix].x1 = (short) absx; - sp->csegs[ix].y1 = (short) absy; - }else{ - sp->csegs[ix].x2 = (short) absx; - sp->csegs[ix].y2 = (short) absy; + + if(i < end -1){ /* start new segment */ + sp->csegs[IX(col)].x1 = absx; + sp->csegs[IX(col)].y1 = absy; } } - if(i == 2) /* both assigned */ - sp->cnsegs[col]++; -#else - /* Chalky's code w/ clipping */ - if (b < ((sp->mode & FLOW_BOX) ? 2 : MAX_BOX)) - continue; - do { - double x1=X(0,b)-sp->centre.x; - double y1=Y(0,b)-sp->centre.y; - double z1=Z(0,b)-sp->centre.z; - double X1=M[0][0]*x1 + M[0][1]*y1 + M[0][2]*z1; - double Y1=M[1][0]*x1 + M[1][1]*y1 + M[1][2]*z1; - double Z1=M[2][0]*x1 + M[2][1]*y1 + M[2][2]*z1+sp->view.height; - double absx1, absy1; - double x2=X(1,b)-sp->centre.x; - double y2=Y(1,b)-sp->centre.y; - double z2=Z(1,b)-sp->centre.z; - double X2=M[0][0]*x2 + M[0][1]*y2 + M[0][2]*z2; - double Y2=M[1][0]*x2 + M[1][1]*y2 + M[1][2]*z2; - double Z2=M[2][0]*x2 + M[2][1]*y2 + M[2][2]*z2+sp->view.height; - double absx2, absy2; - if(sp->view.depth){ - /* Need clipping if: is part of box, or close to viewer */ - if ( (b >= MIN_BOX && b < MAX_BOX) || X1 <= 0.1 || X2 < 0.1) - if (clip(&X1, &Y1, &Z1, &X2, &Y2, &Z2)) - break; - if (X1 <= 0 || X2 <= 0) break; - absx1=SCALE_X(sp->view.depth*Y1/X1); - absy1=SCALE_Y(sp->view.depth*Z1/X1); - if(absx1 < -sp->width || absx1 > 2*sp->width || - absy1 < -sp->height || absy1 > 2*sp->height) - break; - absx2=SCALE_X(sp->view.depth*Y2/X2); - absy2=SCALE_Y(sp->view.depth*Z2/X2); - if(absx2 < -sp->width || absx2 > 2*sp->width || - absy2 < -sp->height || absy2 > 2*sp->height) - break; - }else{ - absx1=SCALE_X(X1); - absy1=SCALE_Y(Y1); - absx2=SCALE_X(X2); - absy2=SCALE_Y(Y2); - } - - sp->csegs[ix].x1 = (short) absx1; - sp->csegs[ix].y1 = (short) absy1; - sp->csegs[ix].x2 = (short) absx2; - sp->csegs[ix].y2 = (short) absy2; - - sp->cnsegs[col]++; - } while (0); -#endif } - if (sp->count) { /* erase */ - XSetForeground(display, gc, MI_BLACK_PIXEL(mi)); - XDrawSegments(display, window, gc, sp->old_segs, sp->nold_segs); + /* Erase */ + XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_BLACK_PIXEL(mi)); + if (dbufp) { /* In Double Buffer case, prepare off-screen copy */ + /* For slow systems, this can be the single biggest bottleneck + in the program. These systems may be better of not using + the double buffer. */ + XFillRectangle(MI_DISPLAY(mi), sp->buffer, MI_GC(mi), 0, 0, + MI_WIDTH(mi), MI_HEIGHT(mi)); + } else { /* Otherwise, erase previous segment list directly */ + XDrawSegments(MI_DISPLAY(mi), sp->buffer, MI_GC(mi), + sp->old_segs, sp->nold_segs); } - if (MI_NPIXELS(mi) > 2){ /* render colour */ - for (col = 0; col < MI_NPIXELS(mi); col++) + /* Render */ + if (MI_NPIXELS(mi) > 2){ /* colour */ + int mn = 0; + for (col = 0; col < MI_NPIXELS(mi) - 1; col++) if (sp->cnsegs[col] > 0) { - XSetForeground(display, gc, MI_PIXEL(mi, col)); - XDrawSegments(display, window, gc, - sp->csegs + col * sp->beecount, sp->cnsegs[col]); + if(sp->cnsegs[col] > mn) mn = sp->cnsegs[col]; + XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_PIXEL(mi, col+1)); + /* This is usually the biggest bottleneck on most + systems. The maths load is insignificant compared + to this. */ + XDrawSegments(MI_DISPLAY(mi), sp->buffer, MI_GC(mi), + sp->csegs + col * segindex, sp->cnsegs[col]); } - } else { /* render mono */ - XSetForeground(display, gc, MI_WHITE_PIXEL(mi)); - XDrawSegments(display, window, gc, - sp->csegs + col * sp->beecount, sp->cnsegs[col]); + } else { /* mono handled seperately since xlockmore uses '1' for + mono white! */ + XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_WHITE_PIXEL(mi)); + XDrawSegments(MI_DISPLAY(mi), sp->buffer, MI_GC(mi), + sp->csegs, sp->cnsegs[0]); + } + if (dbufp) { /* In Double Buffer case, this updates the screen */ + XCopyArea(MI_DISPLAY(mi), sp->buffer, MI_WINDOW(mi), MI_GC(mi), 0, 0, + MI_WIDTH(mi), MI_HEIGHT(mi), 0, 0); + } else { /* Otherwise, screen is already updated. Copy segments + to erase-list to be erased directly next time. */ + int c = 0; + for (col = 0; col < MI_NPIXELS(mi) - 1; col++) { + memcpy(sp->old_segs + c, sp->csegs + col * segindex, + sp->cnsegs[col] * sizeof(XSegment)); + c += sp->cnsegs[col]; + } + sp->nold_segs = c; } - /* Copy to erase-list */ - for (col = 0, c = 0; col < MI_NPIXELS(mi); col++) - for (b = 0; b < sp->cnsegs[col]; b++) - sp->old_segs[c++] = (sp->csegs + col * sp->beecount)[b]; - sp->nold_segs = c; - - if (++sp->count > MI_CYCLES(mi)) /* pick a new flow */ + if(sp->count > 1 && swarm == 0) { /* all gone */ + if(MI_IS_VERBOSE(mi)) + fprintf(stdout, "flow: all gone at %d\n", sp->count); init_flow(mi); - - if (sp->count % (MI_CYCLES(mi)/4) == 0) { /* pick a new view */ - new_view = 0; /* change to 1 .. */ } - if (X(0, 0) < xmin*2 || X(0, 0) > xmax*2) new_view = 1; - if (Y(0, 0) < ymin*2 || Y(0, 0) > ymax*2) new_view = 1; - if (Z(0, 0) < zmin*2 || Z(0, 0) > zmax*2) new_view = 1; - - if (new_view) { - for (b = 0; b < 2; b++) { - X(1, b) = X(0, b) = balance_rand(sp->range.x*4); - Y(1, b) = Y(0, b) = balance_rand(sp->range.y*4); - Z(1, b) = Z(0, b) = balance_rand(sp->range.z*4); - } - sp->slow_view = 0.90; + if(sp->count++ > MI_CYCLES(mi)){ /* Time's up. If we haven't + found anything new by now we + should pick a new standard + flow */ + init_flow(mi); } } -void -release_flow(ModeInfo * mi) +ENTRYPOINT void +reshape_flow(ModeInfo * mi, int width, int height) { - if (flows != NULL) { - int screen; - - for (screen = 0; screen < MI_NUM_SCREENS(mi); screen++) { - flowstruct *sp = &flows[screen]; - - if (sp->csegs != NULL) - (void) free((void *) sp->csegs); - if (sp->cnsegs != NULL) - (void) free((void *) sp->cnsegs); - if (sp->old_segs != NULL) - (void) free((void *) sp->old_segs); - if (sp->p[0] != NULL) - (void) free((void *) sp->p[0]); - if (sp->p[1] != NULL) - (void) free((void *) sp->p[1]); - } - (void) free((void *) flows); - flows = NULL; - } + init_flow (mi); } -void -refresh_flow(ModeInfo * mi) + +ENTRYPOINT void +refresh_flow (ModeInfo * mi) { - MI_CLEARWINDOW(mi); + if(!dbufp) MI_CLEARWINDOW(mi); } -XrmOptionDescRec flow_options[] = +ENTRYPOINT Bool +flow_handle_event (ModeInfo *mi, XEvent *event) { - {"-rotate", ".rotate", XrmoptionSepArg, 0}, - {"-ride", ".ride", XrmoptionSepArg, 0}, - {"-zoom", ".zoom", XrmoptionSepArg, 0}, - {"-box", ".box", XrmoptionSepArg, 0}, - {"-slow", ".slow", XrmoptionSepArg, 0}, - {"-freeze", ".freeze", XrmoptionSepArg, 0}, - {"-allow2d", ".allow2d", XrmoptionSepArg, 0}, - { 0, 0, 0, 0 } -}; + if (screenhack_event_helper (MI_DISPLAY(mi), MI_WINDOW(mi), event)) + { + init_flow (mi); + return True; + } + return False; +} -/* -char* defaults[] = -{ - "*rotate: True", - "*ride: True", - "*zoom: True", - "*allow2d: True", - "*box: True", - "*slow: True", - "*freeze: True", - 0 -}; - */ + +XSCREENSAVER_MODULE ("Flow", flow) + +#endif /* MODE_flow */