#if 0
	cc -o chaos chaos.c -lm -lX11
	exit
#endif

/*
 * Chaos program.  Copyright 1991 Ken Shirriff
 * Send comments, etc. to shirriff@sprite.Berkeley.EDU
 *
 * This program draws the "chaotic attractor".
 * Usage: chaos [-r min max] [-i ilow ihigh] [-s step]\n");
 * -r gives the minimum and maximum values of the logistic equation parameter.
 * -i gives the iteration range to be plotted and used for the Lyapunov
 * exponent computation.  (The iterations below ilow are not used, so they
 * can be used to allow the function to settle.
 * -s gives the step size (i.e. how big the steps are in logistic equation
 * parameter space.  1=1 pixel.
 *
 * After drawing, left click prints the X position of the mouse in terms
 * of the logistic equation parameter.
 * Right click ends the program.
 *
 * Compile program with "sh chaos.c"
 *
 */


#include <X11/Xlib.h>
#include <stdio.h>
#include <math.h>

struct Point {
	double x,y;
};
typedef struct Point pos_t;

/*
 * SIZE is the size of the X window to use.
 */
#define SIZE 600

/*
 * The logistic equation.
 */
double logis(x,s)
double x,s;
{
    double xp;
    xp = s*x*(1-x);
    return xp;
}

main(argc,argv)
int argc;
char **argv;
{
    double s;	/* Our logistics equation driving parameter. */
    int i;
    double x;	/* Our iterated x value. */
    pos_t p;
    float lya;	/* Lyapunov exponent. */
    int N0=10, NT=100; /* Low and high iterations to plot. */
    float S0=0, S1=4.0; /* Low and high parameter values. */
    float div=1;	/* Plotting resolution in pixels. */
#define abs(x) ((x)>0?(x):-(x))

    if (argc < 1) {
usage:
	fprintf(stderr,"Usage: [-r min max] [-i ilow ihigh] [-s step]\n");
	exit(1);
    }
    while (1) {
	if (argc==1) {
	    break;
	} else if (strncmp("-h",argv[1],2)==0) {
	    goto usage;
	} else if (strcmp("-r",argv[1])==0 && argc>3) {
	    sscanf(argv[2],"%f", &S0);
	    sscanf(argv[3],"%f", &S1);
	    if (S0>=S1) {
		fprintf(stderr,"Invalid range\n");
		exit(-1);
	    }
	    argc -= 3;
	    argv += 3;
	} else if (strcmp("-i",argv[1])==0 && argc>3) {
	    sscanf(argv[2],"%d", &N0);
	    sscanf(argv[3],"%d", &NT);
	    argc -= 3;
	    argv += 3;
	} else if (strcmp("-s",argv[1])==0 && argc>2) {
	    sscanf(argv[2],"%f", &div);
	    argc -= 2;
	    argv += 2;
	} else if (strcmp("-i",argv[1])==0 && argc>3) {
	    sscanf(argv[2],"%d", &N0);
	    sscanf(argv[3],"%d", &NT);
	    argc -= 3;
	    argv += 3;
	} else {
	    goto usage;
	}
    }

    xinit();
    range(S0,-2.,S1,1.5);

    /*
     * Loop across the different 'C' values (in the variable s)
     */
    for (s=S0;s<S1;s+=(S1-S0)/SIZE*div) {
	x = .5;
	lya = 0;
	/*
	 * Do NT iterations of the logistics equation.
	 */
	for (i=0;i<NT;i++) {
	    x = logis(x,s);
	    if (i>N0) {
		/* 
		 * Plot the point.
		 */
		p.x = s;
		p.y = x;
		point(&p);
		/*
		 * Sum up the Lyapunov exponent.
		 */
		lya += log(abs(s-2*x*s))/(NT-N0);
	    }
	}
	/*
	 * Plot the Lyapunov exponent and 0 axis.
	 */
	p.x = s;
	p.y = lya/2-1; /* Scale lya to lower part of window. */
	point(&p);
	p.x = s;
	p.y = 0./2-1;
	point(&p);
    }

    /*
     * Flush the drawing buffer.
     */
    drawflush();
    printf("Done\n");
    click();
    endit();
}

/*
 * -------------------------------------------------------------
 * X drivers
 */


Display *dp;
Window w;
GC gc;

float Xxmin, Xxmax, Xymin, Xymax;

#define SCALEX(v) ((int)(((v)-(Xxmin))*SIZE/(Xxmax-(Xxmin))+.5))
#define ISCALEX(n) (Xxmin+(n)*(Xxmax-Xxmin)/SIZE)
#define SCALEY(v) (SIZE-((int)(((v)-(Xymin))*SIZE/(Xymax-(Xymin))+.5)))
#define ISCALEY(n) (Xymin+(SIZE-n)*(Xymax-Xymin)/SIZE)

/*
 * Initialize the X window.
 */

xinit()
{
    Window root;
    Screen *sc;
    int dscreen;
    XGCValues gcvals;
    XSetWindowAttributes watt;
    int depth;
    XEvent xevent;
    int done=0;

    if ((dp=XOpenDisplay(NULL))==NULL) {
	fprintf(stderr,"Could not open Display!0\n");
	exit(1);
    }
    XSynchronize(dp, True);
    dscreen = XDefaultScreen(dp);
    sc = ScreenOfDisplay(dp,dscreen);
    depth = DefaultDepth(dp,dscreen);
    root = DefaultRootWindow(dp);
    watt.background_pixel = WhitePixelOfScreen(sc);
    watt.bit_gravity = StaticGravity;
    w = XCreateWindow(dp,root,0,0,SIZE,SIZE,0,
       depth,InputOutput,CopyFromParent,CWBackPixel|CWBitGravity, &watt);
    XStoreName(dp,w,"Test");
    gcvals.background = BlackPixelOfScreen(sc);
    gcvals.foreground = BlackPixelOfScreen(sc);
    gc = XCreateGC(dp,w, GCForeground|GCBackground, &gcvals);
    XSetFunction(dp,gc,GXcopy);
    XSelectInput(dp,w,ExposureMask|EnterWindowMask|PointerMotionMask|
	ButtonPressMask|StructureNotifyMask|SubstructureNotifyMask);
    XMapWindow(dp,w);
    XSync(dp,False);
    XClearArea(dp,w,0,0,SIZE, SIZE, False);
    while (!done) {
	XNextEvent(dp, &xevent);
	switch (((XAnyEvent *)&xevent)->type) {
	    case ButtonPress:
		if ( ((XButtonEvent*)&xevent)->button == Button3) {
		    done = 1;
		}
		break;
	    case Expose:
		done = 1;
		break;
	    default:
		break;
	}
    }
}

/*
 * Set the plotting range.
 */
range(pxmin, pymin, pxmax, pymax)
float pxmin, pxmax, pymin, pymax;
{
    Xxmin = pxmin;
    Xxmax = pxmax;
    Xymin = pymin;
    Xymax = pymax;
}

/*
 * Plot the point.  Points are buffered and plotted in large bunches
 * to make it run faster.
 */
#define BUF 1000
XPoint buf[BUF];
int bufp=0;
point(p)
pos_t *p;
{
    if (bufp>=BUF) {
	XDrawPoints(dp,w,gc,buf,bufp,CoordModeOrigin);
	bufp=0;
    }
    buf[bufp].x = SCALEX(p->x);
    buf[bufp].y = SCALEY(p->y);
    bufp++;
}

/*
 * Flush the point drawing buffer.
 */
drawflush()
{
    if (bufp>0) {
	XDrawPoints(dp,w,gc,buf,bufp,CoordModeOrigin);
	bufp=0;
    }
}

/*
 * Wait for mouse click.
 */
click()
{
    XEvent xevent;
    int done=0;
    while (!done) {
	XNextEvent(dp, &xevent);
	switch (((XAnyEvent *)&xevent)->type) {
	    case ButtonPress:
		if ( ((XButtonEvent*)&xevent)->button == Button3) {
		    done = 1;
		} else if ( ((XButtonEvent*)&xevent)->button == Button1) {
		    /*
		     * If we get a click, print the mouse position.
		     */
		    int x,y;
		    float fy;
		    x =  ((XButtonEvent *)&xevent)->x;
		    y =  ((XButtonEvent *)&xevent)->y;
		    fy = ISCALEY(y);
		    /*
		     * Rescale the Lyapunov window coordinates.
		     */
		    if (fy<-.2) {
			printf("fy=%f\n", fy);
			fy = (fy+1)*2;
		    }
		    printf("(%f,%f)\n", ISCALEX(x), fy);
		}
		break;
	    default:
		break;
	}
    }
}

/*
 * Close the window.
 */
endit()
{
    XDestroyWindow(dp, w);
    XCloseDisplay(dp);
}