source: trunk/third/xscreensaver/hacks/flow.c @ 20148

/* -*- Mode: C; tab-width: 4 -*- */
/* flow --- flow of strange bees */

#if 0
static const char sccsid[] = "@(#)flow.c 4.10 98/04/24 xlockmore";
#endif

/*-
 * Copyright (c) 1996 by Tim Auckland <Tim.Auckland@Sun.COM>
 * Portions added by Stephen Davies are Copyright (c) 2000 Stephen Davies
 *
 * Permission to use, copy, modify, and distribute this software and its
 * documentation for any purpose and without fee is hereby granted,
 * provided that the above copyright notice appear in all copies and that
 * both that copyright notice and this permission notice appear in
 * supporting documentation.
 *
 * This file is provided AS IS with no warranties of any kind.  The author
 * shall have no liability with respect to the infringement of copyrights,
 * trade secrets or any patents by this file or any part thereof.  In no
 * event will the author be liable for any lost revenue or profits or
 * other special, indirect and consequential damages.
 *
 * "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.
 *
 * 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)
 *   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
 *   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)
 */

#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
# 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 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, "",
        "Shows dynamic strange attractors", 0, NULL
};

#endif

typedef struct {
        double      x;
        double      y;
        double      z;
} dvector;

typedef struct {
        double      a, b, c;
} Par;

/* 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 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))

/* 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)

typedef struct {
        int         width;
        int         height;
        int         count;
        double      size;
        
        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;
        dvector  (*ODE) (Par par, double x, double y, double z);
        Par         par;
        FlowMode                mode; /* Mode of operation */
} flowstruct;

static flowstruct *flows = NULL;

static dvector
Lorentz(Par par, double x, double y, double z)
{
        dvector d;

        d.x = par.a * (y - x);
        d.y = x * (par.b - z) - y;
        d.z = x * y - par.c * z;
        return d;
}

static dvector
Rossler(Par par, 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);
        return d;
}

static dvector
RosslerCone(Par par, double x, double y, double z)
{
        dvector d;

        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;
}

static dvector
Birkhoff(Par par, double x, double y, double z)
{
        dvector d;

        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;
}

static dvector
Duffing(Par par, double x, double y, double z)
{
        dvector d;

        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;
}

void init_clip(flowstruct *sp);

void
init_flow(ModeInfo * mi)
{
        flowstruct *sp;
        int         b;
        double      beemult = 1;
        static int  allocated = 0;

        if (flows == NULL) {
                if ((flows = (flowstruct *) calloc(MI_NUM_SCREENS(mi),
                                               sizeof (flowstruct))) == NULL)
                        return;
        }
        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;
        }
        
        switch (b) {
        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);
                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);
                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);
                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;
                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;
                break;
        }

        sp->view.depth *= 4;

        sp->beecount = beemult * MI_COUNT(mi);
        if (sp->beecount < 0)   /* random variations */ 
                sp->beecount = NRAND(-sp->beecount) + 1; /* Minimum 1 */

        /* Clear the background. */
        MI_CLEARWINDOW(mi);

        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;
        }

        /* Allocate memory. */

        if (!sp->csegs) {
                sp->csegs = (XSegment *) malloc(sizeof (XSegment) * sp->beecount
                                                * MI_NPIXELS(mi));
                sp->cnsegs = (int *) malloc(sizeof (int) * MI_NPIXELS(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);
        }

        /* Initialize point positions, velocities, etc. */

        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);

}

/* 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];

#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}
};

/* Container for scaled box points */
static double box[BOX_P][3];

/* 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},
};
        
/* Boundaries of bees */
double xmin, xmax;
double ymin, ymax;
double zmin, zmax;

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);
}

/* 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;
        }

}

/* 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;
                } else {
                        *x2 = x;
                        *y2 = y;
                        *z2 = z;
                }
        }
        return 0;
}       


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;

        }

        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;

        /* <=- 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;
                        }
                }


                /* 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];
                                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;
                        
                                /* 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 (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];
                        }
                
                
#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){
                                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 == 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);
        }

        if (MI_NPIXELS(mi) > 2){ /* render colour */
                for (col = 0; col < MI_NPIXELS(mi); 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]);
                        }
        } else {                /* render mono */
                XSetForeground(display, gc, MI_WHITE_PIXEL(mi));
                XDrawSegments(display, window, gc,
                                          sp->csegs + col * sp->beecount, sp->cnsegs[col]);
        }

        /* 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 */
                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;
        }
}

void
release_flow(ModeInfo * mi)
{
        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;
        }
}

void
refresh_flow(ModeInfo * mi)
{
        MI_CLEARWINDOW(mi);
}

XrmOptionDescRec flow_options[] =
{
        {"-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 }
};

/*
char*        defaults[] =
{
        "*rotate:         True",
        "*ride:           True",
        "*zoom:           True",
        "*allow2d:        True",
        "*box:            True",
        "*slow:           True",
        "*freeze:         True",
  0
};
        */