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

#if !defined( lint ) && !defined( SABER )
static const char sccsid[] = "@(#)flow.c 4.10 98/04/24 xlockmore";

#endif

/*-
 * Copyright (c) 1996 by Tim Auckland <Tim.Auckland@Sun.COM>
 *
 * 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:
 * 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:		1024 \n" \
					"*cycles:		3000 \n" \
					"*ncolors:		200 \n"
# define SMOOTH_COLORS
# include "xlockmore.h"		/* in xscreensaver distribution */
# include "erase.h"

#else /* STANDALONE */
# include "xlock.h"		/* in xlockmore distribution */
#endif /* STANDALONE */

ModeSpecOpt flow_opts = { 0, NULL, 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))

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;
	dvector     centre;		/* centre */
	struct {
		double  depth;
		double  height;
	}           view;
	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;
} 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_flow(ModeInfo * mi)
{
	flowstruct *sp;
	int         b;
	double      beemult = 1;
	dvector     range;
	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->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 */

	switch (NRAND(8)) {
	case 0:
		sp->view.depth = 10;
		sp->view.height = 0.2;
		beemult = 3;
	case 1:
		sp->ODE = Lorentz;
		sp->step = 0.02;
		sp->size = 60;
		sp->centre.x = 0;
		sp->centre.y = 0;
		sp->centre.z = 24;
		range.x = 5;
		range.y = 5;
		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 2:
		sp->view.depth = 10;
		sp->view.height = 0.1;
		beemult = 4;
	case 3:
		sp->ODE = Rossler;
		sp->step = 0.05;
		sp->size = 24;
		sp->centre.x = 0;
		sp->centre.y = 0;
		sp->centre.z = 3;
		range.x = 4;
		range.y = 4;
		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 4:
		sp->view.depth = 10;
		sp->view.height = 0.1;
		beemult = 3;
	case 5:
		sp->ODE = RosslerCone;
		sp->step = 0.05;
		sp->size = 24;
		sp->centre.x = 0;
		sp->centre.y = 0;
		sp->centre.z = 3;
		range.x = 4;
		range.y = 4;
		range.z = 4;
		sp->par.a = 2;
		sp->par.b = 0.2;
		sp->par.c = 0.25 + balance_rand(0.09);
		break;
	case 6:
		sp->ODE = Birkhoff;
		sp->step = 0.04;
		sp->size = 2.6;
		sp->centre.x = 0;
		sp->centre.y = 0;
		sp->centre.z = 0;
		range.x = 3;
		range.y = 4;
		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 7:
	default:
		sp->ODE = Duffing;
		sp->step = 0.02;
		sp->size = 30;
		sp->centre.x = 0;
		sp->centre.y = 0;
		sp->centre.z = 0;
		range.x = 20;
		range.y = 20;
		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->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;
		}
	}

	/* 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(range.x);
		Y(1, b) = Y(0, b) = balance_rand(range.y);
		Z(1, b) = Z(0, b) = balance_rand(range.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;
	double      M[3][3]; /* transformation matrix */

	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;

	/* <=- Bees -=> */
	for (b = 0; b < sp->beecount; b++) {
		/* Age the arrays. */
		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;

			k1 = sp->ODE(sp->par, X(1, b), Y(1, b), Z(1, b));
			k1.x *= sp->step;
			k1.y *= sp->step;
			k1.z *= sp->step;
			k2 = sp->ODE(sp->par, X(1, b) + k1.x, Y(1, b) + k1.y, Z(1, b) + k1.z);
			k2.x *= sp->step;
			k2.y *= sp->step;
			k2.z *= sp->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);
				sp->centre.y=Y(0, b);
				sp->centre.z=Z(0, b);
			}else if(b==1){ /* neighbour: used to compute local axes */
				double x[3], p[3], x2=0, xp=0;
				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];

#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
			}
		
		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]++;
    }
	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);
}

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