Source Code 1994 March/Source_Code_CD-ROM_Walnut_Creek_March_1994.iso/compsrcs/x/volume21/lyap/part01

Text File | 1993-10-07 | 53.4 KB | 1,669 lines

Newsgroups: comp.sources.x
From: [e-mail] (Ronald Joe Record)
Subject: v21i004:  lyap - Calculate Lyapunov fractals, Part01/03
Message-ID: <csx-v21i004=[e-mail]>
X-Md4-Signature: 729ec9dae6c6a2e720a5fd49f36c7404
Sender: [e-mail] (Chris Olson)
Organization: Sterling Software
Date: Thu, 7 Oct 1993 17:46:45 GMT
Approved: [e-mail]

Submitted-by: [e-mail] (Ronald Joe Record)
Posting-number: Volume 21, Issue 4
Archive-name: lyap/part01
Environment: X11
Supersedes: lyapunov-xlib: Volume 17, Issue 48-49

The idea here is to calculate the Lyapunov exponent for a periodically
forced logistic map (later i added several other nonlinear maps of the unit
interval). In order to turn the 1-dimensional parameter space of the
logistic map into a 2-dimensional parameter space, select two parameter
values ('a' and 'b') then alternate the iterations of the logistic map using
first 'a' then 'b' as the parameter. This program accepts an argument to
specify a forcing function, so instead of just alternating 'a' and 'b', you
can use 'a' as the parameter for say 6 iterations, then 'b' for 6 iterations
and so on. An interesting forcing function to look at is abbabaab (the
Morse-Thue sequence, an aperiodic self-similar, self-generating sequence).
Anyway, step through all the values of 'a' and 'b' in the ranges you want,
calculating the Lyapunov exponent for each pair of values. The exponent
is calculated by iterating out a ways (specified by the variable "settle")
then on subsequent iterations calculating an average of the logarithm of
the absolute value of the derivative at that point. Points in parameter
space with a negative Lyapunov exponent are colored one way (using the
value of the exponent to index into a color map) while points with a
non-negative exponent are colored differently.

The algorithm was taken from the September 1991 Scientific American article
by A. K. Dewdney who gives credit to Mario Markus of the Max Planck Institute
for its creation. Additional information and ideas were gleaned from the
discussion on alt.fractals involving Stephen Hall, Ed Kubaitis, Dave Platt
and Baback Moghaddam. Assistance with colormaps and spinning color wheels
and X was gleaned from Hiram Clawson. Rubber banding code was adapted from
an existing Mandelbrot program written by Stacey Campbell.

#! /bin/sh
# This is a shell archive.  Remove anything before this line, then feed it
# into a shell via "sh file" or similar.  To overwrite existing files,
# type "sh file -c".
# Contents:  lyap lyap/lib lyap/lyap.c lyap/params lyap/testit
# Wrapped by chris@sparky on Thu Oct  7 12:38:39 1993
PATH=/bin:/usr/bin:/usr/ucb:/usr/local/bin:/usr/lbin ; export PATH
echo If this archive is complete, you will see the following message:
echo '          "shar: End of archive 1 (of 3)."'
if test ! -d 'lyap' ; then
    echo shar: Creating directory \"'lyap'\"
    mkdir 'lyap'
fi
if test ! -d 'lyap/lib' ; then
    echo shar: Creating directory \"'lyap/lib'\"
    mkdir 'lyap/lib'
fi
if test -f 'lyap/lyap.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'lyap/lyap.c'\"
else
  echo shar: Extracting \"'lyap/lyap.c'\" \(47187 characters\)
  sed "s/^X//" >'lyap/lyap.c' <<'END_OF_FILE'
X/*
X *    @(#) lyap.c 9.9 93/08/23 SCOINC
X */
X/*************************************************************************
X *                                                                       *
X *  Copyright (c) 1992, 1993 Ronald Joe Record                           *
X *                                                                       *
X *  All rights reserved. No part of this program or publication may be   *
X *  reproduced, transmitted, transcribed, stored in a retrieval system,  *
X *  or translated into any language or computer language, in any form or *
X *  by any means, electronic, mechanical, magnetic, optical, chemical,   *
X *  biological, or otherwise, without the prior written permission of:   *
X *                                                                       *
X *      Ronald Joe Record [phone]                                 *
X *      212 Owen St., Santa Cruz, California 95062 USA                   *
X *                                                                       *
X *************************************************************************/
X
X/* Lyap - calculate and display Lyapunov exponents */
X
X/* Written by Ronald Record ([e-mail]) 03 Sep 1991 */
X
X/* The idea here is to calculate the Lyapunov exponent for a periodically
X * forced logistic map (later i added several other nonlinear maps of the unit
X * interval). In order to turn the 1-dimensional parameter space of the
X * logistic map into a 2-dimensional parameter space, select two parameter
X * values ('a' and 'b') then alternate the iterations of the logistic map using
X * first 'a' then 'b' as the parameter. This program accepts an argument to
X * specify a forcing function, so instead of just alternating 'a' and 'b', you
X * can use 'a' as the parameter for say 6 iterations, then 'b' for 6 iterations
X * and so on. An interesting forcing function to look at is abbabaab (the
X * Morse-Thue sequence, an aperiodic self-similar, self-generating sequence).
X * Anyway, step through all the values of 'a' and 'b' in the ranges you want,
X * calculating the Lyapunov exponent for each pair of values. The exponent
X * is calculated by iterating out a ways (specified by the variable "settle")
X * then on subsequent iterations calculating an average of the logarithm of
X * the absolute value of the derivative at that point. Points in parameter
X * space with a negative Lyapunov exponent are colored one way (using the
X * value of the exponent to index into a color map) while points with a
X * non-negative exponent are colored differently.
X *
X * The algorithm was taken from the September 1991 Scientific American article
X * by A. K. Dewdney who gives credit to Mario Markus of the Max Planck Institute
X * for its creation. Additional information and ideas were gleaned from the
X * discussion on alt.fractals involving Stephen Hall, Ed Kubaitis, Dave Platt
X * and Baback Moghaddam. Assistance with colormaps and spinning color wheels
X * and X was gleaned from Hiram Clawson. Rubber banding code was adapted from
X * an existing Mandelbrot program written by Stacey Campbell.
X */
X
X#include "lyap.h"
X
Xstatic char *lyap_c_id = "@(#) lyap.c 9.9 93/08/23 SCOINC";
Xstatic char *version = LYAP_VERSION;
X
Xextern void BufferPoint(), InitBuffer(), FlushBuffer();
X
Xvoid
XCleanup() {
X    XCloseDisplay(dpy);
X}
X
Xmain(ac, av)
X    int ac;
X    char **av;
X{
X    int i, j;
X    XSizeHints hint;
X    extern void init_canvas(), init_data(), init_color(), parseargs();
X    extern void Clear(), restor_picture();
X
X    parseargs(ac, av);
X    dpy = XOpenDisplay("");
X    screen = DefaultScreen(dpy);
X    background = BlackPixel(dpy, screen);
X    if (width < 1)
X        width = XDisplayWidth(dpy, screen);
X    if (height < 1)
X        height = XDisplayHeight(dpy, screen);
X    setupmem();
X    init_data();
X    if (displayplanes > 1)
X        foreground = startcolor;
X    else
X        foreground = WhitePixel(dpy,XDefaultScreen(dpy));
X    if (xposition < 0)
X        xposition = Max(0,(XDisplayWidth(dpy, screen) - width)/2);
X    if (yposition < 0)
X        yposition = Max(0,(XDisplayHeight(dpy, screen) - height)/2);
X    hint.x = xposition;
X    hint.y = yposition;
X    hint.width = width;
X    hint.height = height;
X    hint.flags = PPosition | PSize;
X    /*
X     * Create the window to display the Lyapunov exponents
X     */
X    canvas = XCreateSimpleWindow(dpy, DefaultRootWindow(dpy),
X        hint.x, hint.y, hint.width, hint.height,
X        5, foreground, background);
X    XSetStandardProperties(dpy, canvas, "Lyap by Ron Record",
X        "Lyap", None, av, ac, &hint);
X    XChangeProperty( dpy, canvas, XA_WM_CLASS, XA_STRING, 8, PropModeReplace,
X                    "lyap", strlen("lyap"));
X    init_canvas();
X    XMapRaised(dpy, canvas);
X    /* Try to write into a new color map */
X    cmap = XCreateColormap(dpy, canvas, DefaultVisual(dpy, screen), AllocAll);
X    if (displayplanes > 1)
X        init_color(dpy, canvas, cmap, Colors, startcolor, mincolindex,
X                    numcolors, numwheels, "lyap", "Lyap", 0);
X    else
X        XQueryColors(dpy, DefaultColormap(dpy, DefaultScreen(dpy)),
X                Colors, numcolors);
X    pixmap = XCreatePixmap(dpy, DefaultRootWindow(dpy), width, height,
X                           DefaultDepth(dpy, screen));
X    rubber_data.band_cursor=XCreateFontCursor(dpy, XC_hand2);
X    CreateXorGC();
X    Clear();
X    if (restfile)
X        restor_picture();
X    if (demo)
X        for(;;) {
X            if (!run) {
X                DemoSpin(dpy, cmap, Colors, startcolor, maxcolor, delay, 2);
X                for (i=0; i<=MAXWHEELS; i++) {
X                    init_color(dpy, canvas, cmap, Colors, startcolor,
X                            mincolindex, numcolors, i, "lyap", "Lyap", 0);
X                    sleep(1);
X                }
X                if (!run) {
X                    Cleanup();
X                    exit(0);
X                }
X            }
X        }
X    else {
X        XSelectInput(dpy,canvas,KeyPressMask|ButtonPressMask|ButtonMotionMask|
X                ButtonReleaseMask|ExposureMask|StructureNotifyMask);
X        for(;;)
X            main_event();
X    }
X}
X
Xmain_event()
X{
X    int n;
X    XEvent event;
X
X    if (complyap() == TRUE)
X        run=0;
X    n = XEventsQueued(dpy, QueuedAfterFlush);
X    while (n--) {
X        XNextEvent(dpy, &event);
X        switch(event.type) {
X            case KeyPress: Getkey(&event); break;
X            case Expose: redisplay(canvas, &event); break;
X            case ConfigureNotify: resize(); break;
X            case ButtonPress: StartRubberBand(canvas, &rubber_data, &event);
X                 break;
X            case MotionNotify: TrackRubberBand(canvas, &rubber_data, &event);
X                 break;
X            case ButtonRelease: EndRubberBand(canvas, &rubber_data, &event);
X                 break;
X        }
X    }
X}
X
X/* complyap() is the guts of the program. This is where the Lyapunov exponent
X * is calculated. For each iteration (past some large number of iterations)
X * calculate the logarithm of the absolute value of the derivative at that
X * point. Then average them over some large number of iterations. Some small
X * speed up is achieved by utilizing the fact that log(a*b) = log(a) + log(b).
X */
Xcomplyap()
X{
X    register i, bindex, findex;
X    double total, prod, x, r;
X    extern void save_to_file(), store_to_file();
X
X    if (!run)
X        return TRUE;
X    if (a >= max_a)
X        if (sendpoint(lyapunov) == TRUE)
X            return FALSE;
X        else {
X            run=0;
X            FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points, 0, maxcolor);
X            if (savefile)
X                save_to_file();
X            else if (storefile)
X                store_to_file();
X            return TRUE;
X        }
X    if (b >= max_b) {
X        run=0;
X        FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points, 0, maxcolor);
X        if (savefile)
X            save_to_file();
X        else if (storefile)
X            store_to_file();
X        return TRUE;
X    }
X    prod = 1.0;
X    total = 0.0;
X    bindex = 0;
X    x = start_x;
X    r = (forcing[bindex]) ? b : a;
X#ifdef MAPS
X    findex = 0;
X    map = Maps[Forcing[findex]];
X#endif
X    for (i=0;i<settle;i++) {       /* Here's where we let the thing */
X        x = (*map)(x, r);       /* "settle down". There is usually */
X        if (++bindex >= maxindex) { /* some initial "noise" in the */
X            bindex = 0;       /* iterations. How can we optimize */
X            if (Rflag)        /* the value of settle ??? */
X                setforcing();
X        }
X        r = (forcing[bindex]) ? b : a;
X#ifdef MAPS
X        if (++findex >= funcmaxindex)
X            findex = 0;
X        map = Maps[Forcing[findex]];
X#endif
X    }
X#ifdef MAPS
X    deriv = Derivs[Forcing[findex]];
X#endif
X    if (useprod) {            /* using log(a*b) */
X        for (i=0;i<dwell;i++) {
X            x = (*map)(x, r);
X            prod *= ABS((*deriv)(x, r));
X            /* we need to prevent overflow and underflow */
X            if ((prod > 1.0e12) || (prod < 1.0e-12)) {
X                if (prod == 0.0)
X                    total += LN_MINDOUBLE;
X                else if (prod < 0.0)
X                    total += LN_MAXDOUBLE;
X                else
X                    total += log(prod);
X                prod = 1.0;
X            }
X            if (++bindex >= maxindex) {
X                bindex = 0;
X                if (Rflag)
X                    setforcing();
X            }
X            r = (forcing[bindex]) ? b : a;
X#ifdef MAPS
X            if (++findex >= funcmaxindex)
X                findex = 0;
X            map = Maps[Forcing[findex]];
X            deriv = Derivs[Forcing[findex]];
X#endif
X        }
X        if (prod == 0.0)
X            total += LN_MINDOUBLE;
X        else if (prod < 0.0)
X            total += LN_MAXDOUBLE;
X        else
X            total += log(prod);
X        lyapunov = (total * M_LOG2E) / (double)dwell;
X    }
X    else {                /* use log(a) + log(b) */
X        for (i=0;i<dwell;i++) {
X            x = (*map)(x, r);
X            total += log(ABS((*deriv)(x, r)));
X            if (++bindex >= maxindex) {
X                bindex = 0;
X                if (Rflag)
X                    setforcing();
X            }
X            r = (forcing[bindex]) ? b : a;
X#ifdef MAPS
X            if (++findex >= funcmaxindex)
X                findex = 0;
X            map = Maps[Forcing[findex]];
X            deriv = Derivs[Forcing[findex]];
X#endif
X        }
X        lyapunov = (total * M_LOG2E) / (double)dwell;
X    }
X    if (sendpoint(lyapunov) == TRUE)
X        return FALSE;
X    else {
X        run=0;
X        FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points, 0, maxcolor);
X        if (savefile)
X            save_to_file();
X        else if (storefile)
X            store_to_file();
X        return TRUE;
X    }
X}
X
Xdouble
Xlogistic(x, r)            /* the familiar logistic map */
Xdouble x, r;
X{
X    return(r * x * (1.0 - x));
X}
X
Xdouble
Xdlogistic(x, r)            /* the derivative of logistic map */
Xdouble x, r;
X{
X    return(r - (2.0 * r * x));
X}
X
Xdouble
Xcircle(x, r)            /* sin() hump or sorta like the circle map */
Xdouble x, r;
X{
X    extern double sin();
X
X    return(r * sin(M_PI * x));
X}
X
Xdouble
Xdcircle(x, r)            /* derivative of the "sin() hump" */
Xdouble x, r;
X{
X    extern double cos();
X
X    return(r * M_PI * cos(M_PI * x));
X}
X
Xdouble
Xleftlog(x, r)            /* left skewed logistic */
Xdouble x, r;
X{
X    double d;
X
X    d = 1.0 - x;
X    return(r * x * d * d);
X}
X
Xdouble
Xdleftlog(x, r)            /* derivative of the left skewed logistic */
Xdouble x, r;
X{
X    return(r * (1.0 - (4.0 * x) + (3.0 * x * x)));
X}
X
Xdouble
Xrightlog(x, r)            /* right skewed logistic */
Xdouble x, r;
X{
X    return(r * x * x * (1.0 - x));
X}
X
Xdouble
Xdrightlog(x, r)            /* derivative of the right skewed logistic */
Xdouble x, r;
X{
X    return(r * ((2.0 * x) - (3.0 * x * x)));
X}
X
Xdouble
Xdoublelog(x, r)            /* double logistic */
Xdouble x, r;
X{
X    double d;
X
X    d = 1.0 - x;
X    return(r * x * x * d * d);
X}
X
Xdouble
Xddoublelog(x, r)        /* derivative of the double logistic */
Xdouble x, r;
X{
X    double d;
X
X    d = x * x;
X    return(r * ((2.0 * x) - (6.0 * d) + (4.0 * x * d)));
X}
X
Xvoid
Xinit_data()
X{
X    static int i;
X
X    numcolors = XDisplayCells(dpy, XDefaultScreen(dpy));
X    displayplanes = DisplayPlanes(dpy, XDefaultScreen(dpy));
X    if (numcolors > maxcolor)
X        numcolors = maxcolor;
X    if (numcolors <= 16) {
X        maxcolor = 16; startcolor = 0;
X        color_offset = 0; mincolindex = 1;
X    }
X    numfreecols = numcolors - mincolindex;
X    lowrange = mincolindex - startcolor;
X    a_inc = a_range / (double)width;
X    b_inc = b_range / (double)height;
X    if (!restfile) {
X        point.x = point.y = 0;
X        a = min_a; b = min_b;
X    }
X    rubber_data.p_min = min_a;
X    rubber_data.q_min = min_b;
X    rubber_data.p_max = max_a;
X    rubber_data.q_max = max_b;
X    if (show)
X        show_defaults();
X    InitBuffer(&Points, maxcolor);
X    srand48(time(0));
X}
X
Xvoid
Xinit_canvas()
X{
X    static int i;
X
X    /*
X     * create default, writable, graphics contexts for the canvas.
X     */
X    for (i=0; i<maxcolor; i++) {
X        Data_GC[i] = XCreateGC(dpy, DefaultRootWindow(dpy),
X                (unsigned long) NULL, (XGCValues *) NULL);
X        /* set the background to black */
X        XSetBackground(dpy,Data_GC[i],BlackPixel(dpy,XDefaultScreen(dpy)));
X        /* set the foreground of the ith context to i */
X        XSetForeground(dpy, Data_GC[i], i);
X    }
X    if (displayplanes == 1) {
X        XSetForeground(dpy,Data_GC[0],BlackPixel(dpy,XDefaultScreen(dpy)));
X        XSetForeground(dpy,Data_GC[1],WhitePixel(dpy,XDefaultScreen(dpy)));
X    }
X}
X
Xvoid
Xparseargs(ac, av)
Xint ac;
Xchar **av;
X{
X    static int c;
X    static int i;
X    int bindex=0, findex;
X    char *ch;
X    extern int optind;
X    extern char *optarg;
X    extern double atof();
X    extern void check_params(), usage(), restor_params();
X
X    map = Maps[0];
X    deriv = Derivs[0];
X    maxexp=minlyap; minexp= -1.0 * minlyap;
X
X    while ((c=getopt(ac,av,
X        "Ldpuvc:i:m:C:W:H:I:M:N:O:R:S:a:b:D:F:f:o:w:h:s:x:y:"))!=EOF){
X        switch (c) {
X        case 'C':    mincolindex=atoi(optarg); break;
X        case 'D':    if (strcmp(optarg, "elay"))
X                        dwell=atoi(optarg);
X                     else
X                        delay = atoi(av[optind++]);
X                     break;
X#ifdef MAPS
X        case 'F':    funcmaxindex = strlen(optarg);
X                if (funcmaxindex > FUNCMAXINDEX)
X                    usage();
X                ch = optarg;
X                Force++;
X                for (findex=0;findex<funcmaxindex;findex++) {
X                    Forcing[findex] = (int)(*ch++ - '0');
X                    if (Forcing[findex] >= NUMMAPS)
X                        usage();
X                }
X                break;
X#endif
X        case 'H':    height=atoi(optarg); break;
X        case 'I':    start_x=atof(optarg); break;
X        case 'L':    useprod=0; break;
X        case 'M':    minlyap=ABS(atof(optarg));
X                maxexp=minlyap; minexp= -1.0 * minlyap; break;
X        case 'N':    maxcolor=ABS(atoi(optarg));
X                if ((maxcolor - startcolor) <= 0)
X                    startcolor = 0;
X                if ((maxcolor - mincolindex) <= 0) {
X                    mincolindex = 1;
X                    color_offset = 0;
X                }
X                break;
X        case 'O':    color_offset=atoi(optarg); break;
X        case 'R':    prob=atof(optarg); Rflag++; setforcing(); break;
X        case 'S':    settle=atoi(optarg); break;
X        case 'W':    width=atoi(optarg); break;
X        case 'a':    min_a=atof(optarg); aflag++; break;
X        case 'b':    min_b=atof(optarg); bflag++; break;
X        case 'c':    numwheels=atoi(optarg); break;
X        case 'd':    demo++; break;
X        case 'f':    maxindex = strlen(optarg);
X                if (maxindex > MAXINDEX)
X                    usage();
X                ch = optarg;
X                force++;
X                while (bindex < maxindex) {
X                    if (*ch == 'a')
X                        forcing[bindex++] = 0;
X                    else if (*ch == 'b')
X                        forcing[bindex++] = 1;
X                    else
X                        usage();
X                    ch++;
X                }
X                break;
X        case 'h':    b_range=atof(optarg); hflag++; break;
X        case 'i':    restfile++; inname=optarg; break;
X        case 'm':    mapindex=atoi(optarg);
X                if ((mapindex >= NUMMAPS) || (mapindex < 0))
X                    usage();
X                map = Maps[mapindex];
X                deriv = Derivs[mapindex];
X                if (!aflag)
X                    min_a = amins[mapindex];
X                if (!wflag)
X                    a_range = aranges[mapindex];
X                if (!bflag)
X                    min_b = bmins[mapindex];
X                if (!hflag)
X                    b_range = branges[mapindex];
X                if (!Force)
X                    for (i=0;i<FUNCMAXINDEX;i++)
X                        Forcing[i] = mapindex;
X                break;
X        case 'o':    savefile++; outname=optarg; break;
X        case 'p':    negative--; break;
X        case 's':    storefile++; savname=optarg; break;
X        case 'u':    usage(); break;
X        case 'v':    show=1; break;
X        case 'w':    a_range=atof(optarg); wflag++; break;
X        case 'x':    xposition=atoi(optarg); break;
X        case 'y':    yposition=atoi(optarg); break;
X        case '?':    usage(); break;
X        }
X    }
X    if (restfile)
X        restor_params();
X    else {
X        max_a = min_a + a_range;
X        max_b = min_b + b_range;
X    }
X    a_minimums[0] = min_a; b_minimums[0] = min_b;
X    a_maximums[0] = max_a; b_maximums[0] = max_b;
X    if (Force)
X        if (maxindex == funcmaxindex)
X            for (findex=0;findex<funcmaxindex;findex++)
X                check_params(Forcing[findex],forcing[findex]);
X        else
X            fprintf(stderr, "Warning! Unable to check parameters\n");
X    else
X        check_params(mapindex,2);
X}
X
Xvoid
Xcheck_params(mapnum, parnum)
Xint mapnum;
Xint parnum;
X{
X
X    if (parnum != 1) {
X        if ((max_a > pmaxs[mapnum]) || (min_a < pmins[mapnum])) {
X            fprintf(stderr, "Warning! Parameter 'a' out of range.\n");
X            fprintf(stderr, "You have requested a range of (%lf,%lf).\n",
X                    min_a,max_a);
X            fprintf(stderr, "Valid range is (%lf,%lf).\n",
X                    pmins[mapnum],pmaxs[mapnum]);
X        }
X    }
X    if (parnum != 0) {
X        if ((max_b > pmaxs[mapnum]) || (min_b < pmins[mapnum])) {
X            fprintf(stderr, "Warning! Parameter 'b' out of range.\n");
X            fprintf(stderr, "You have requested a range of (%lf,%lf).\n",
X                    min_b,max_b);
X            fprintf(stderr, "Valid range is (%lf,%lf).\n",
X                    pmins[mapnum],pmaxs[mapnum]);
X        }
X    }
X}
X
Xvoid
Xusage()
X{
X    fprintf(stderr,"lyap [-Ls][-W#][-H#][-a#][-b#][-w#][-h#][-x xstart]\n");
X    fprintf(stderr,"\t[-M#][-S#][-D#][-f string][-r#][-O#][-C#][-c#][-m#]\n");
X#ifdef MAPS
X    fprintf(stderr,"\t[-F string]\n");
X#endif
X    fprintf(stderr,"\tWhere: -C# specifies the minimum color index\n");
X    fprintf(stderr,"\t       -r# specifies the maxzimum rgb value\n");
X    fprintf(stderr,"\t       -u displays this message\n");
X    fprintf(stderr,"\t       -a# specifies the minimum horizontal parameter\n");
X    fprintf(stderr,"\t       -b# specifies the minimum vertical parameter\n");
X    fprintf(stderr,"\t       -w# specifies the horizontal parameter range\n");
X    fprintf(stderr,"\t       -h# specifies the vertical parameter range\n");
X    fprintf(stderr,"\t       -D# specifies the dwell\n");
X    fprintf(stderr,"\t       -S# specifies the settle\n");
X    fprintf(stderr,"\t       -H# specifies the initial window height\n");
X    fprintf(stderr,"\t       -W# specifies the initial window width\n");
X    fprintf(stderr,"\t       -O# specifies the color offset\n");
X    fprintf(stderr,"\t       -c# specifies the desired color wheel\n");
X    fprintf(stderr,"\t       -m# specifies the desired map (0-4)\n");
X    fprintf(stderr,"\t       -f aabbb specifies a forcing function of 00111\n");
X#ifdef MAPS
X    fprintf(stderr,"\t       -F 00111 specifies the function forcing function\n");
X#endif
X    fprintf(stderr,"\t       -L indicates use log(x)+log(y) rather than log(xy)\n");
X    fprintf(stderr,"\tDuring display :\n");
X    fprintf(stderr,"\t     Use the mouse to zoom in on an area\n");
X    fprintf(stderr,"\t     e or E recalculates color indices\n");
X    fprintf(stderr,"\t     f or F saves exponents to a file\n");
X    fprintf(stderr,"\t     KJmn increase/decrease minimum negative exponent\n");
X    fprintf(stderr,"\t     r or R redraws\n");
X    fprintf(stderr,"\t     s or S spins the colorwheel\n");
X    fprintf(stderr,"\t     w or W changes the color wheel\n");
X    fprintf(stderr,"\t     x or X clears the window\n");
X    fprintf(stderr,"\t     q or Q exits\n");
X    exit(1);
X}
X
XCycle_frames()
X{
X    static int i;
X    extern void redraw();
X
X    for (i=0;i<=maxframe;i++)
X        redraw(exponents[i], expind[i], 1);
X}
X
XGetkey(event)
XXKeyEvent *event;
X{
X    unsigned char key;
X    static int i, running, spindir=0, spinning=0;
X    extern void init_color(), recalc(), print_values(), print_help();
X    extern void go_init(), go_back(), go_down(), Clear();
X    extern void save_to_file(), Redraw(), redraw(), store_to_file();
X
X    if (XLookupString(event, (char *)&key, sizeof(key), (KeySym *)0,
X            (XComposeStatus *) 0) > 0)
X        switch (key) {
X            case '(': delay -= 25; if (delay < 0) delay = 0; break;
X            case ')': delay += 25; break;
X            case '<': dwell /= 2; if (dwell < 1) dwell = 1; break;
X            case '>': dwell *= 2; break;
X            case '[': settle /= 2; if (settle < 1) settle = 0; break;
X            case ']': settle *= 2; if (settle < 1) settle = 1; break;
X            case 'd': go_down(); break;
X            case 'D': FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points,
X                                0, maxcolor);
X                    break;
X            case 'e':
X            case 'E': FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points,
X                                0, maxcolor);
X                  dorecalc = (!dorecalc);
X                  if (dorecalc)
X                      recalc();
X                  else {
X                      maxexp = minlyap; minexp = -1.0 * minlyap;
X                  }
X                  redraw(exponents[frame], expind[frame], 1);
X                  break;
X            case 'f': FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points,
X                        0, maxcolor);
X                      store_to_file(); break;
X            case 'F': FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points,
X                        0, maxcolor);
X                      save_to_file(); break;
X            case 'i':
X                  if (stripe_interval > 0) {
X                      stripe_interval--;
X                      if (displayplanes > 1) {
X                          running = run; run = 0;
X                          init_color(dpy,canvas,cmap,Colors,startcolor,
X                               mincolindex,numcolors,numwheels,"lyap","Lyap",0);
X                          run = running;
X                      }
X                  }
X                  break;
X            case 'I': stripe_interval++;
X                  if (displayplanes > 1) {
X                      running = run; run = 0;
X                      init_color(dpy,canvas,cmap,Colors,startcolor,
X                            mincolindex,numcolors,numwheels,"lyap","Lyap",0);
X                      run = running;
X                  }
X                  break;
X            case 'K': if (minlyap > 0.05)
X                          minlyap -= 0.05;
X                   break;
X            case 'J': minlyap += 0.05;
X                   break;
X            case 'm': mapindex++;
X                  if (mapindex >= NUMMAPS)
X                      mapindex=0;
X                  map = Maps[mapindex];
X                  deriv = Derivs[mapindex];
X                  if (!aflag)
X                      min_a = amins[mapindex];
X                  if (!wflag)
X                      a_range = aranges[mapindex];
X                  if (!bflag)
X                      min_b = bmins[mapindex];
X                  if (!hflag)
X                      b_range = branges[mapindex];
X                  if (!Force)
X                      for (i=0;i<FUNCMAXINDEX;i++)
X                           Forcing[i] = mapindex;
X                  max_a = min_a + a_range;
X                  max_b = min_b + b_range;
X                  a_minimums[0] = min_a; b_minimums[0] = min_b;
X                  a_maximums[0] = max_a; b_maximums[0] = max_b;
X                  a_inc = a_range / (double)width;
X                  b_inc = b_range / (double)height;
X                  point.x = 0;
X                  point.y = 0;
X                  a = rubber_data.p_min = min_a;
X                  b = rubber_data.q_min = min_b;
X                  rubber_data.p_max = max_a;
X                  rubber_data.q_max = max_b;
X                  Clear();
X                  run = 1;
X                  break;
X            case 'M': if (minlyap > 0.005)
X                      minlyap -= 0.005;
X                  break;
X            case 'N': minlyap += 0.005;
X                  break;
X            case 'p':
X            case 'P': negative = (!negative);
X                  FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points,
X                                0, maxcolor);
X                  redraw(exponents[frame], expind[frame], 1);
X                  break;
X            case 'r': FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points,
X                            0, maxcolor);
X                  redraw(exponents[frame],expind[frame],1);
X                  break;
X            case 'R': FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points,
X                            0, maxcolor);
X                  Redraw();
X                  break;
X            case 'S': spinning=0; break;
X            case 's': spinning=1; spindir=(!spindir);
X                  Spin(dpy,cmap,Colors,startcolor,numcolors,delay,spindir);
X                  break;
X            case 'u': go_back(); break;
X            case 'U': go_init(); break;
X            case 'v':
X            case 'V': print_values(); break;
X            case '\027': /* (ctrl-W) read color palette from $HOME/.lyapmap */
X                  numwheels = 0;
X                  if (displayplanes > 1) {
X                      running = run; run = 0;
X                      init_color(dpy,canvas,cmap,Colors,startcolor,
X                            mincolindex,numcolors,numwheels,"lyap","Lyap",0);
X                      run = running;
X                  }
X                  break;
X            case 'W': if (numwheels < MAXWHEELS)
X                          numwheels++;
X                      else
X                          numwheels = 0;
X                      if (displayplanes > 1) {
X                          running = run; run = 0;
X                            init_color(dpy,canvas,cmap,Colors,startcolor,
X                               mincolindex,numcolors,numwheels,"lyap","Lyap",0);
X                          run = running;
X                      }
X                      break;
X            case 'w': if (numwheels > 0)
X                          numwheels--;
X                      else
X                          numwheels = MAXWHEELS;
X                      if (displayplanes > 1) {
X                          running = run; run = 0;
X                            init_color(dpy,canvas,cmap,Colors,startcolor,
X                               mincolindex,numcolors,numwheels,"lyap","Lyap",0);
X                          run = running;
X                      }
X                      break;
X            case 'x': Clear(); break;
X            case 'X': Destroy_frame(); break;
X            case 'z': running = run; run = 0;
X                  Cycle_frames();
X                  run = running; redraw(exponents[frame], expind[frame], 1);
X                  break;
X            case 'Z': running = run; run = 0;
X                  while (!XPending(dpy)) Cycle_frames();
X                  run = running; redraw(exponents[frame], expind[frame], 1);
X                  break;
X            case 'q':
X            case 'Q': Cleanup(); exit(0); break;
X            case '?':
X            case 'h':
X            case 'H': print_help(); break;
X            default:  break;
X        }
X        if (spinning)
X            Spin(dpy,cmap,Colors,startcolor,numcolors,delay,spindir);
X}
X
X/* Here's where we index into a color map. After the Lyapunov exponent is
X * calculated, it is used to determine what color to use for that point.
X * I suppose there are a lot of ways to do this. I used the following :
X * if it's non-negative then there's a reserved area at the lower range
X * of the color map that i index into. The ratio of some "minimum exponent
X * value" and the calculated value is used as a ratio of how high to index
X * into this reserved range. Usually these colors are dark red (see init_color).
X * If the exponent is negative, the same ratio (expo/minlyap) is used to index
X * into the remaining portion of the colormap (which is usually some light
X * shades of color or a rainbow wheel). The coloring scheme can actually make
X * a great deal of difference in the quality of the picture. Different colormaps
X * bring out different details of the dynamics while different indexing
X * algorithms also greatly effect what details are seen. Play around with this.
X */
Xsendpoint(expo)
Xdouble expo;
X{
X    static int index;
X    static double tmpexpo;
X    extern double exp();
X
X    tmpexpo = (negative) ? expo : -1.0 * expo;
X    if (tmpexpo > 0) {
X        if (displayplanes >1) {
X            index = (int)(tmpexpo*lowrange/maxexp);
X            index = (index % lowrange) + startcolor;
X        }
X        else
X            index = 0;
X    }
X    else {
X        if (displayplanes >1) {
X            index = (int)(tmpexpo*numfreecols/minexp);
X            index = (index % numfreecols) + mincolindex;
X        }
X        else
X            index = 1;
X    }
X    BufferPoint(dpy, canvas, pixmap, Data_GC, &Points, index, point.x, point.y);
X    point.x++;
X    a += a_inc;
X    if (save)
X        exponents[frame][expind[frame]++] = expo;
X    if (point.x >= width) {
X        point.y++;
X        point.x = 0;
X        if (save) {
X            b += b_inc;
X            a = min_a;
X        }
X        if (point.y >= height)
X            return FALSE;
X        else
X            return TRUE;
X    }
X    return TRUE;
X}
X
Xvoid
Xredisplay (w, event)
XWindow          w;
XXExposeEvent    *event;
X{
X    /*
X     * Extract the exposed area from the event and copy
X     * from the saved pixmap to the window.
X     */
X    XCopyArea(dpy, pixmap, canvas, Data_GC[0],
X           event->x, event->y, event->width, event->height,
X           event->x, event->y);
X}
X
Xvoid
Xresize()
X{
X    Window r;
X    int n, x, y;
X    unsigned int bw, d, new_w, new_h;
X    extern void Clear(), Redraw();
X
X    XGetGeometry(dpy,canvas,&r,&x,&y,&new_w,&new_h,&bw,&d);
X    if ((new_w == width) && (new_h == height))
X        return;
X    freemem();
X    width = new_w; height = new_h;
X    XClearWindow(dpy, canvas);
X    if (pixmap)
X        XFreePixmap(dpy, pixmap);
X    pixmap = XCreatePixmap(dpy, DefaultRootWindow(dpy),
X            width, height, DefaultDepth(dpy, screen));
X    a_inc = a_range / (double)width;
X    b_inc = b_range / (double)height;
X    point.x = 0;
X    point.y = 0;
X    run = 1;
X    a = rubber_data.p_min = min_a;
X    b = rubber_data.q_min = min_b;
X    rubber_data.p_max = max_a;
X    rubber_data.q_max = max_b;
X    setupmem();
X    for (n=0;n<MAXFRAMES;n++)
X        if ((n <= maxframe) && (n != frame))
X            resized[n] = 1;
X    InitBuffer(&Points, maxcolor);
X    Clear();
X    Redraw();
X}
X
Xvoid
Xredraw(exparray, index, cont)
Xdouble *exparray;
Xint index, cont;
X{
X    static int i;
X    static int x_sav, y_sav;
X
X    x_sav = point.x;
X    y_sav = point.y;
X
X    point.x = 0;
X    point.y = 0;
X
X    save=0;
X    for (i=0;i<index;i++)
X        sendpoint(exparray[i]);
X    save=1;
X
X    if (cont) {
X        point.x = x_sav;
X        point.y = y_sav;
X    }
X    else {
X        a = point.x * a_inc + min_a;
X        b = point.y * b_inc + min_b;
X    }
X    FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points, 0, maxcolor);
X}
X
Xvoid
XRedraw()
X{
X    FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points, 0, maxcolor);
X    point.x = 0;
X    point.y = 0;
X    run = 1;
X    a = min_a;
X    b = min_b;
X    expind[frame] = 0;
X    resized[frame] = 0;
X}
X
X/* Store colormap indices in file so we can read them in later */
Xvoid
Xstore_to_file()
X{
X    FILE *outfile;
X    unsigned char c;
X    XImage *ximage;
X    static int i,j;
X
X    outfile = fopen(savname,"w");
X    if(!outfile) {
X        perror(savname);
X        Cleanup();
X        exit(-1);
X    }
X    fprintf(outfile,"# width=%d height=%d \n",width,height);
X    fprintf(outfile,"# settle=%d dwell=%d start_x=%lf \n",settle,dwell,start_x);
X    fprintf(outfile,"# min_a=%lf a_rng=%lf max_a=%lf \n",min_a,a_range,max_a);
X    fprintf(outfile,"# min_b=%lf b_rng=%lf max_b=%lf \n",min_b,b_range,max_b);
X    fprintf(outfile,"# run=%d point.x=%d point.y=%d \n",run,point.x,point.y);
X    fprintf(outfile,"# negative=%d mincolindex=%d startcolor=%d \n",
X            negative,mincolindex,startcolor);
X    fprintf(outfile,"# minexp=%lf maxexp=%lf a=%lf b=%lf \n",
X            minexp,maxexp,a,b);
X    fprintf(outfile,"# Force=%d force=%d useprod=%d \n",Force,force,useprod);
X    fprintf(outfile,"# mapindex=%d maxindex=%d \n", mapindex, maxindex);
X    fprintf(outfile,"# numwheels=%d Rflag=%d color_offset=%d \n",
X            numwheels,Rflag,color_offset);
X    fprintf(outfile,"# maxcolor=%d start_x=%lf minlyap=%lf prob=%lf\n",
X            maxcolor,start_x,minlyap,prob);
X    fprintf(outfile,"# maxindex=%d periodic forcing=", maxindex);
X        for (i=0;i<maxindex;i++) {
X            fprintf(outfile," %d",forcing[i]);
X        }
X    fprintf(outfile," \n");
X    fprintf(outfile,"# funcmaxindex=%d function forcing=", funcmaxindex);
X    for (i=0;i<funcmaxindex;i++)
X         fprintf(outfile," %d",Forcing[i]);
X    fprintf(outfile," \n");
X    fprintf(outfile,"# numcolors=%d\n",numcolors);
X
X    ximage=XGetImage(dpy, pixmap, 0, 0,
X            (unsigned int)width, (unsigned int)height, AllPlanes, XYPixmap);
X
X    for (j=0;j<=point.y;j++)
X        for (i=0;i<width;i++) {
X            c = (unsigned char)XGetPixel(ximage,i,j);
X            fwrite((char *)&c,sizeof c,1,outfile);
X        }
X    fclose(outfile);
X}
X
X/* Store color pics in PPM format and monochrome in PGM */
Xvoid
Xsave_to_file()
X{
X    FILE *outfile;
X    unsigned char c;
X    XImage *ximage;
X    static int i,j;
X    struct Colormap {
X        unsigned char red;
X        unsigned char green;
X        unsigned char blue;
X    };
X    struct Colormap *colormap=NULL;
X
X    if (colormap)
X        free(colormap);
X    if ((colormap=
X        (struct Colormap *)malloc(sizeof(struct Colormap)*maxcolor))
X            == NULL) {
X        fprintf(stderr,"Error malloc'ing colormap array\n");
X        Cleanup();
X        exit(-1);
X    }
X    outfile = fopen(outname,"w");
X    if(!outfile) {
X        perror(outname);
X        Cleanup();
X        exit(-1);
X    }
X
X    ximage=XGetImage(dpy, pixmap, 0, 0, width, height, AllPlanes, XYPixmap);
X
X    if (displayplanes > 1) {
X        for (i=0;i<maxcolor;i++) {
X            colormap[i].red=(unsigned char)(Colors[i].red >> 8);
X            colormap[i].green=(unsigned char)(Colors[i].green >> 8);
X            colormap[i].blue =(unsigned char)(Colors[i].blue >> 8);
X        }
X        fprintf(outfile,"P%d %d %d\n",6,width,height);
X    }
X    else
X        fprintf(outfile,"P%d %d %d\n",5,width,height);
X    fprintf(outfile,"# settle=%d  dwell=%d start_x=%lf\n",settle,dwell,start_x);
X    fprintf(outfile,"# min_a=%lf  a_rng=%lf  max_a=%lf\n",min_a,a_range,max_a);
X    fprintf(outfile,"# min_b=%lf  b_rng=%lf  max_b=%lf\n",min_b,b_range,max_b);
X    if (Rflag)
X        fprintf(outfile,"# pseudo-random forcing\n");
X    else if (force) {
X        fprintf(outfile,"# periodic forcing=");
X        for (i=0;i<maxindex;i++) {
X            fprintf(outfile,"%d",forcing[i]);
X        }
X        fprintf(outfile,"\n");
X    }
X    else
X         fprintf(outfile,"# periodic forcing=01\n");
X    if (Force) {
X         fprintf(outfile,"# function forcing=");
X         for (i=0;i<funcmaxindex;i++)
X             fprintf(outfile,"%d",Forcing[i]);
X         fprintf(outfile,"\n");
X    }
X    fprintf(outfile,"%d\n",numcolors-1);
X
X    for (j=0;j<height;j++)
X        for (i=0;i<width;i++) {
X            c = (unsigned char)XGetPixel(ximage,i,j);
X            if (displayplanes > 1)
X                fwrite((char *)&colormap[c],sizeof colormap[0],1,outfile);
X            else
X                fwrite((char *)&c,sizeof c,1,outfile);
X        }
X    fclose(outfile);
X}
X
X/* Read saved file with parameters */
Xvoid
Xrestor_params()
X{
X    unsigned char s[16];
X    static int i;
X
X    infile = fopen(inname,"r");
X    if(!infile) {
X        perror(inname);
X        Cleanup();
X        exit(-1);
X    }
X
X    fscanf(infile,"%1s%6s%d%7s%d",s,s,&width,s,&height);
X    fscanf(infile,"%1s%7s%d%6s%d%8s%lf",s,s,&settle,s,&dwell,s,&start_x);
X    fscanf(infile,"%1s%6s%lf%6s%lf%6s%lf", s,s,&min_a,s,&a_range,s,&max_a);
X    fscanf(infile,"%1s%6s%lf%6s%lf%6s%lf", s,s,&min_b,s,&b_range,s,&max_b);
X    fscanf(infile,"%1s%4s%d%8s%d%8s%d",s,s,&run,s,&point.x,s,&point.y);
X    fscanf(infile,"%1s%9s%d%12s%d%11s%d",
X            s,s,&negative,s,&mincolindex,s,&startcolor);
X    fscanf(infile,"%1s%7s%lf%7s%lf%2s%lf%2s%lf",
X            s,s,&minexp,s,&maxexp,s,&a,s,&b);
X    fscanf(infile,"%1s%6s%d%6s%d%8s%d",s,s,&Force,s,&force,s,&useprod);
X    fscanf(infile,"%1s%9s%d%9s%d",
X            s,s,&mapindex,s,&maxindex);
X    fscanf(infile,"%1s%10s%d%6s%d%13s%d",
X            s,s,&numwheels,s,&Rflag,s,&color_offset);
X    fscanf(infile,"%1s%9s%d%8s%lf%8s%lf%5s%lf",
X            s,s,&maxcolor,s,&start_x,s,&minlyap,s,&prob);
X    fscanf(infile,"%1s%9s%d%s%s", s,s,&maxindex,s,s);
X    for(i=0;i<maxindex;i++) {
X        fscanf(infile, "%d", &forcing[i]);
X    }
X    fscanf(infile,"%1s%13s%d%s%s", s,s,&funcmaxindex,s,s);
X    for(i=0;i<funcmaxindex;i++)
X        fscanf(infile, "%d", &Forcing[i]);
X    fscanf(infile, "%1s%10s%d", s,s,&numcolors);
X    fread(s, sizeof s[0], 1, infile);
X    map = Maps[mapindex];
X    deriv = Derivs[mapindex];
X}
X
Xvoid
Xrestor_picture() {
X    static int i, j, k;
X    unsigned char c;
X
X    for (j=0;j<point.y;j++)
X        for (i=0;i<width;i++) {
X            k = fread(&c, sizeof c, 1, infile);
X            if (k) {
X                BufferPoint(dpy,canvas,pixmap,Data_GC,&Points,(int)c,i,j);
X                if (k < mincolindex)
X                    exponents[frame][expind[frame]++] =
X                        -(double)k/(double)mincolindex;
X                else
X                    exponents[frame][expind[frame]++] =
X                        (double)k/(double)numfreecols;
X            }
X            else
X                break;
X        }
X    for (j=0;j<point.x;j++) {
X        k = fread(&c, sizeof c, 1, infile);
X        if (k) {
X            BufferPoint(dpy,canvas,pixmap,Data_GC,&Points,(int)c,j,point.y);
X            if (k < mincolindex)
X                exponents[frame][expind[frame]++] =
X                    -(double)k/(double)mincolindex;
X            else
X                exponents[frame][expind[frame]++] =
X                    (double)k/(double)numfreecols;
X        }
X        else
X            break;
X    }
X    fclose(infile);
X    FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points, 0, maxcolor);
X}
X
Xvoid
Xrecalc()
X{
X    static int i, index, x, y;
X
X    minexp = maxexp = 0.0;
X    x = y = 0;
X    for (i=0;i<expind[frame];i++) {
X        if (exponents[frame][i] < minexp)
X            minexp = exponents[frame][i];
X        if (exponents[frame][i] > maxexp)
X            maxexp = exponents[frame][i];
X    }
X}
X
Xvoid
XClear()
X{
X    XFillRectangle(dpy, pixmap, Data_GC[0], 0, 0, width, height);
X    XCopyArea(dpy, pixmap, canvas, Data_GC[0],
X                0, 0, width, height, 0, 0);
X    InitBuffer(&Points, maxcolor);
X}
X
Xvoid
Xshow_defaults()
X{
X
X    printf("Width=%d  Height=%d  numcolors=%d  settle=%d  dwell=%d\n",
X        width,height,numcolors,settle,dwell);
X    printf("min_a=%lf  a_range=%lf  max_a=%lf\n", min_a,a_range,max_a);
X    printf("min_b=%lf  b_range=%lf  max_b=%lf\n", min_b,b_range,max_b);
X    printf("minlyap=%lf  minexp=%lf  maxexp=%lf\n", minlyap,minexp,maxexp);
X    exit(0);
X}
X
Xvoid
XCreateXorGC()
X{
X    XGCValues values;
X
X    values.foreground = foreground;
X    values.line_style = LineSolid;
X    values.function = GXxor;
X    RubberGC = XCreateGC(dpy, DefaultRootWindow(dpy),
X          GCForeground | GCBackground | GCFunction | GCLineStyle, &values);
X}
X
Xvoid
XStartRubberBand(w, data, event)
XWindow w;
Ximage_data_t *data;
XXEvent *event;
X{
X    XPoint corners[5];
X    extern void SetupCorners();
X
X    nostart = 0;
X    data->rubber_band.last_x = data->rubber_band.start_x = event->xbutton.x;
X    data->rubber_band.last_y = data->rubber_band.start_y = event->xbutton.y;
X    SetupCorners(corners, data);
X    XDrawLines(dpy, canvas, RubberGC,
X        corners, sizeof(corners) / sizeof(corners[0]), CoordModeOrigin);
X}
X
Xvoid
XSetupCorners(corners, data)
XXPoint *corners;
Ximage_data_t *data;
X{
X    corners[0].x = data->rubber_band.start_x;
X    corners[0].y = data->rubber_band.start_y;
X    corners[1].x = data->rubber_band.start_x;
X    corners[1].y = data->rubber_band.last_y;
X    corners[2].x = data->rubber_band.last_x;
X    corners[2].y = data->rubber_band.last_y;
X    corners[3].x = data->rubber_band.last_x;
X    corners[3].y = data->rubber_band.start_y;
X    corners[4] = corners[0];
X}
X
Xvoid
XTrackRubberBand(w, data, event)
XWindow w;
Ximage_data_t *data;
XXEvent *event;
X{
X    XPoint corners[5];
X    int xdiff, ydiff, diff;
X    extern void SetupCorners();
X
X    if (nostart)
X        return;
X    SetupCorners(corners, data);
X    XDrawLines(dpy, canvas, RubberGC, corners,
X               sizeof(corners) / sizeof(corners[0]), CoordModeOrigin);
X    ydiff = event->xbutton.y - data->rubber_band.start_y;
X    xdiff = event->xbutton.x - data->rubber_band.start_x;
X    data->rubber_band.last_x = data->rubber_band.start_x + xdiff;
X    data->rubber_band.last_y = data->rubber_band.start_y + ydiff;
X    if (data->rubber_band.last_y < data->rubber_band.start_y ||
X        data->rubber_band.last_x < data->rubber_band.start_x) {
X            data->rubber_band.last_y = data->rubber_band.start_y;
X            data->rubber_band.last_x = data->rubber_band.start_x;
X    }
X    SetupCorners(corners, data);
X    XDrawLines(dpy, canvas, RubberGC, corners,
X               sizeof(corners) / sizeof(corners[0]), CoordModeOrigin);
X}
X
Xvoid
XEndRubberBand(w, data, event)
XWindow w;
Ximage_data_t *data;
XXEvent *event;
X{
X    XPoint corners[5];
X    XPoint top, bot;
X    double delta, diff;
X    extern void set_new_params(), SetupCorners();
X
X    nostart = 1;
X    SetupCorners(corners, data);
X    XDrawLines(dpy, canvas, RubberGC,
X        corners, sizeof(corners) / sizeof(corners[0]), CoordModeOrigin);
X    if (data->rubber_band.start_x >= data->rubber_band.last_x ||
X        data->rubber_band.start_y >= data->rubber_band.last_y)
X        return;
X    top.x = data->rubber_band.start_x;
X    bot.x = data->rubber_band.last_x;
X    top.y = data->rubber_band.start_y;
X    bot.y = data->rubber_band.last_y;
X    diff = data->q_max - data->q_min;
X    delta = (double)top.y / (double)height;
X    data->q_min += diff * delta;
X    delta = (double)(height - bot.y) / (double)height;
X    data->q_max -= diff * delta;
X    diff = data->p_max - data->p_min;
X    delta = (double)top.x / (double)width;
X    data->p_min += diff * delta;
X    delta = (double)(width - bot.x) / (double)width;
X    data->p_max -= diff * delta;
X    fflush(stdout);
X    set_new_params(w, data);
X}
X
Xvoid
Xset_new_params(w, data)
XWindow w;
Ximage_data_t *data;
X{
X    extern void Clear();
X
X    frame = (maxframe + 1) % MAXFRAMES;
X    if (frame > maxframe)
X        maxframe = frame;
X    a_range = data->p_max - data->p_min;
X    b_range = data->q_max - data->q_min;
X    a_minimums[frame] = min_a = data->p_min;
X    b_minimums[frame] = min_b = data->q_min;
X    a_inc = a_range / (double)width;
X    b_inc = b_range / (double)height;
X    point.x = 0;
X    point.y = 0;
X    run = 1;
X    a = min_a;
X    b = min_b;
X    a_maximums[frame] = max_a = data->p_max;
X    b_maximums[frame] = max_b = data->q_max;
X    expind[frame] = 0;
X    Clear();
X}
X
Xvoid
Xgo_down()
X{
X    static int i;
X
X    frame++;
X    if (frame > maxframe)
X        frame = 0;
X    jumpwin();
X}
X
Xvoid
Xgo_back()
X{
X    static int i;
X
X    frame--;
X    if (frame < 0)
X        frame = maxframe;
X    jumpwin();
X}
X
Xjumpwin()
X{
X    rubber_data.p_min = min_a = a_minimums[frame];
X    rubber_data.q_min = min_b = b_minimums[frame];
X    rubber_data.p_max = max_a = a_maximums[frame];
X    rubber_data.q_max = max_b = b_maximums[frame];
X    a_range = max_a - min_a;
X    b_range = max_b - min_b;
X    a_inc = a_range / (double)width;
X    b_inc = b_range / (double)height;
X    point.x = 0;
X    point.y = 0;
X    a = min_a;
X    b = min_b;
X    Clear();
X    if (resized[frame])
X        Redraw();
X    else
X        redraw(exponents[frame], expind[frame], 0);
X}
X
Xvoid
Xgo_init()
X{
X    static int i;
X
X    frame = 0;
X    jumpwin();
X}
X
XDestroy_frame()
X{
X    static int i;
X
X    for (i=frame; i<maxframe; i++) {
X        exponents[frame] = exponents[frame+1];
X        expind[frame] = expind[frame+1];
X        a_minimums[frame] = a_minimums[frame+1];
X        b_minimums[frame] = b_minimums[frame+1];
X        a_maximums[frame] = a_maximums[frame+1];
X        b_maximums[frame] = b_maximums[frame+1];
X    }
X    maxframe--;
X    go_back();
X}
X
Xvoid
Xprint_help()
X{
X    printf("During run-time, interactive control can be exerted via : \n");
X    printf("Mouse buttons allow rubber-banding of a zoom box\n");
X    printf("< halves the 'dwell', > doubles the 'dwell'\n");
X    printf("[ halves the 'settle', ] doubles the 'settle'\n");
X    printf("D flushes the drawing buffer\n");
X    printf("d descends a ladder of windows created via the u or U keys\n");
X    printf("e or E recalculates color indices\n");
X    printf("f saves program state to a file\n");
X    printf("F saves window to a PPM file\n");
X    printf("h or H or ? displays this message\n");
X    printf("i decrements, I increments the stripe interval\n");
X    printf("KJMN increase/decrease minimum negative exponent\n");
X    printf("m increments the map index, changing maps\n");
X    printf("p or P reverses the colormap for negative/positive exponents\n");
X    printf("r redraws without recalculating\n");
X    printf("R redraws, recalculating with new dwell and settle values\n");
X    printf("s or S spins the colorwheel\n");
X    printf("u pops back up to the last zoom\n");
X    printf("U pops back up to the first picture\n");
X    printf("v or V displays the values of various settings\n");
X    printf("w decrements, W increments the color wheel index\n");
X    printf("x or X clears the window\n");
X    printf("q or Q exits\n");
X}
X
Xvoid
Xprint_values()
X{
X    static int i;
X
X    printf("\nminlyap=%lf minexp=%lf maxexp=%lf\n",minlyap,minexp,maxexp);
X    printf("width=%d height=%d\n",width,height);
X    printf("settle=%d  dwell=%d start_x=%lf\n",settle,dwell, start_x);
X    printf("min_a=%lf  a_rng=%lf  max_a=%lf\n",min_a,a_range,max_a);
X    printf("min_b=%lf  b_rng=%lf  max_b=%lf\n",min_b,b_range,max_b);
X    if (Rflag)
X        printf("pseudo-random forcing\n");
X    else if (force) {
X        printf("periodic forcing=");
X        for (i=0;i<maxindex;i++)
X            printf("%d",forcing[i]);
X        printf("\n");
X    }
X    else
X        printf("periodic forcing=01\n");
X    if (Force) {
X        printf("function forcing=");
X        for (i=0;i<funcmaxindex;i++)
X            printf("%d",Forcing[i]);
X        printf("\n");
X    }
X    printf("numcolors=%d\n",numcolors);
X}
X
Xfreemem()
X{
X    static int i;
X    for (i=0;i<MAXFRAMES;i++)
X        free(exponents[i]);
X}
X
Xsetupmem()
X{
X    static int i;
X
X    for (i=0;i<MAXFRAMES;i++) {
X        if((exponents[i]=
X            (double *)malloc(sizeof(double)*width*height))==NULL){
X                fprintf(stderr,"Error malloc'ing exponent array.\n");
X                exit(-1);
X        }
X    }
X}
X
Xsetforcing()
X{
X    static int i;
X    extern double drand48();
X
X    for (i=0;i<MAXINDEX;i++)
X        forcing[i] = (drand48() > prob) ? 0 : 1;
X}
END_OF_FILE
  if test 47187 -ne `wc -c <'lyap/lyap.c'`; then
    echo shar: \"'lyap/lyap.c'\" unpacked with wrong size!
  fi
  # end of 'lyap/lyap.c'
fi
if test ! -d 'lyap/params' ; then
    echo shar: Creating directory \"'lyap/params'\"
    mkdir 'lyap/params'
fi
if test -f 'lyap/testit' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'lyap/testit'\"
else
  echo shar: Extracting \"'lyap/testit'\" \(1389 characters\)
  sed "s/^X//" >'lyap/testit' <<'END_OF_FILE'
X#! /bin/sh
X# Test Harness for Lyapunov Space Exploring Program
X
XS=4
XD=8
XW=128
XH=120
XV=$1
X
X[ "$V" = "demo" ] && {
X    S=40
X    D=80
X    W=350
X    H=300
X}
X
X[ "$V" = "vista" ] &&
X    echo "vista:hits 1.3" >> hits
X
X# testing function
Xtest_case()
X{
X
X    [ "$V" = "vista" ] && echo partition "$1" >> hits
X
X    if [ "$V" = "demo" ] ; then
X        $1
X    else
X        $1 > /dev/null 2>&1
X    fi
X    if [ $? = 0 ] ; then
X    [ "$V" = "demo" ] || echo $1 : PASS
X    else
X    [ "$V" = "demo" ] || echo $1 : FAIL
X    fi
X}
X
X[ "$V" = "demo" ] || {
X    echo ""
X    echo "    Running Suite of Test Cases for Lyap"
X    echo ""
X}
X
X[ "$V" = "demo" ] || {
X    echo "Usage message should report as test FAIL"
X    test_case "./lyap -u -d -S $S -D $D -W $W -H $H"
X    echo ""
X}
X
Xtest_case "./lyap -d -S $S -D $D -W $W -H $H"
Xtest_case "./lyap -d -f aaaaaabbbbbb -S $S -D $D -W $W -H $H"
Xtest_case "./lyap -d -f aaaaaabbbbbb -F 000000111111 -S $S -D $D -W $W -H $H"
Xtest_case "./lyap -d -f abbabaab -m 1 -S $S -D $D -W $W -H $H"
Xtest_case "./lyap -d -m 2 -S $S -D $D -W $W -H $H"
Xtest_case "./lyap -d -m 3 -S $S -D $D -W $W -H $H"
Xtest_case "./lyap -d -m 4 -S $S -D $D -W $W -H $H"
X
X[ "$V" = "demo" ] || {
X    echo ""
X    echo "Map index too high should report as test FAIL"
X    test_case "./lyap -d -m 5 -S $S -D $D -W $W -H $H"
X    echo ""
X    echo "    Test Cases Complete "
X    echo ""
X}
X
X[ "$V" = "vista" ] &&
X    echo "  Run ViSTA Reports MetGen, AudGen & CovGen"
Xecho ""
END_OF_FILE
  if test 1389 -ne `wc -c <'lyap/testit'`; then
    echo shar: \"'lyap/testit'\" unpacked with wrong size!
  fi
  chmod +x 'lyap/testit'
  # end of 'lyap/testit'
fi
echo shar: End of archive 1 \(of 3\).
cp /dev/null ark1isdone
MISSING=""
for I in 1 2 3 ; do
    if test ! -f ark${I}isdone ; then
    MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 3 archives.
    rm -f ark[1-9]isdone
else
    echo You still must unpack the following archives:
    echo "        " ${MISSING}
fi
exit 0
exit 0 # Just in case...
--
  // [e-mail]           | Send comp.sources.x submissions to:
\X/  Amiga - The only way to fly! |    [e-mail]
 "It's intuitively obvious to the |
  most casual observer..."        | GCS d+/-- p+ c++ l+ m+ s++/+ g+ w+ t+ r+ x+