git.hungrycats.org Git - xscreensaver/blob - hacks/xlyap.c

   1 /* Lyap - calculate and display Lyapunov exponents */
   2
   3 /* Written by Ron Record (rr@sco) 03 Sep 1991 */
   4
   5 /* The idea here is to calculate the Lyapunov exponent for a periodically
   6  * forced logistic map (later i added several other nonlinear maps of the unit
   7  * interval). In order to turn the 1-dimensional parameter space of the
   8  * logistic map into a 2-dimensional parameter space, select two parameter
   9  * values ('a' and 'b') then alternate the iterations of the logistic map using
  10  * first 'a' then 'b' as the parameter. This program accepts an argument to
  11  * specify a forcing function, so instead of just alternating 'a' and 'b', you
  12  * can use 'a' as the parameter for say 6 iterations, then 'b' for 6 iterations
  13  * and so on. An interesting forcing function to look at is abbabaab (the
  14  * Morse-Thue sequence, an aperiodic self-similar, self-generating sequence).
  15  * Anyway, step through all the values of 'a' and 'b' in the ranges you want,
  16  * calculating the Lyapunov exponent for each pair of values. The exponent
  17  * is calculated by iterating out a ways (specified by the variable "settle")
  18  * then on subsequent iterations calculating an average of the logarithm of
  19  * the absolute value of the derivative at that point. Points in parameter
  20  * space with a negative Lyapunov exponent are colored one way (using the
  21  * value of the exponent to index into a color map) while points with a
  22  * non-negative exponent are colored differently.
  23  *
  24  * The algorithm was taken from the September 1991 Scientific American article
  25  * by A. K. Dewdney who gives credit to Mario Markus of the Max Planck Institute
  26  * for its creation. Additional information and ideas were gleaned from the
  27  * discussion on alt.fractals involving Stephen Hall, Ed Kubaitis, Dave Platt
  28  * and Baback Moghaddam. Assistance with colormaps and spinning color wheels
  29  * and X was gleaned from Hiram Clawson. Rubber banding code was adapted from
  30  * an existing Mandelbrot program written by Stacey Campbell.
  31  */
  32
  33 #define LYAP_PATCHLEVEL 4
  34 #define LYAP_VERSION "#(@) lyap 2.3 2/20/92"
  35
  36 #include <assert.h>
  37 #include <math.h>
  38
  39 #include "screenhack.h"
  40 #include "yarandom.h"
  41 #include "hsv.h"
  42 #include "vroot.h"
  43
  44 #include <X11/cursorfont.h>
  45 #include <X11/Xutil.h>
  46
  47 char *progclass = "XLyap";
  48
  49 char *defaults [] = {
  50   ".background:         black",
  51   ".foreground:         white",
  52   "*randomize:          false",
  53   "*builtin:            -1",
  54   "*minColor:           1",
  55   "*maxColor:           256",
  56   "*dwell:              50",
  57   "*useLog:             false",
  58   "*colorExponent:      1.0",
  59   "*colorOffset:        0",
  60   "*randomForce:        ",              /* 0.5 */
  61   "*settle:             50",
  62   "*minA:               2.0",
  63   "*minB:               2.0",
  64   "*wheels:             7",
  65   "*function:           10101010",
  66   "*forcingFunction:    abbabaab",
  67   "*bRange:             ",              /* 2.0 */
  68   "*startX:             0.65",
  69   "*mapIndex:           ",              /* 0 */
  70   "*outputFile:         ",
  71   "*beNegative:         false",
  72   "*rgbMax:             65000",
  73   "*spinLength:         256",
  74   "*show:               false",
  75   "*aRange:             ",              /* 2.0 */
  76   0
  77 };
  78
  79 XrmOptionDescRec options [] = {
  80   { "-randomize", ".randomize", XrmoptionNoArg, "true" },
  81   { "-builtin",   ".builtin",   XrmoptionSepArg, 0 },
  82   { "-C", ".minColor",          XrmoptionSepArg, 0 },   /* n */
  83   { "-D", ".dwell",             XrmoptionSepArg, 0 },   /* n */
  84   { "-L", ".useLog",            XrmoptionNoArg, "true" },
  85   { "-M", ".colorExponent",     XrmoptionSepArg, 0 },   /* r */
  86   { "-O", ".colorOffset",       XrmoptionSepArg, 0 },   /* n */
  87   { "-R", ".randomForce",       XrmoptionSepArg, 0 },   /* p */
  88   { "-S", ".settle",            XrmoptionSepArg, 0 },   /* n */
  89   { "-a", ".minA",              XrmoptionSepArg, 0 },   /* r */
  90   { "-b", ".minB",              XrmoptionSepArg, 0 },   /* n */
  91   { "-c", ".wheels",            XrmoptionSepArg, 0 },   /* n */
  92   { "-F", ".function",          XrmoptionSepArg, 0 },   /* 10101010 */
  93   { "-f", ".forcingFunction",   XrmoptionSepArg, 0 },   /* abbabaab */
  94   { "-h", ".bRange",            XrmoptionSepArg, 0 },   /* r */
  95   { "-i", ".startX",            XrmoptionSepArg, 0 },   /* r */
  96   { "-m", ".mapIndex",          XrmoptionSepArg, 0 },   /* n */
  97   { "-o", ".outputFile",        XrmoptionSepArg, 0 },   /* filename */
  98   { "-p", ".beNegative",        XrmoptionNoArg, "true" },
  99   { "-r", ".rgbMax",            XrmoptionSepArg, 0 },   /* n */
 100   { "-s", ".spinLength",        XrmoptionSepArg, 0 },   /* n */
 101   { "-v", ".show",              XrmoptionNoArg, "true" },
 102   { "-w", ".aRange",            XrmoptionSepArg, 0 },   /* r */
 103   { 0, 0, 0, 0 }
 104 };
 105
 106
 107 #define ABS(a)  (((a)<0) ? (0-(a)) : (a) )
 108 #define Min(x,y) ((x < y)?x:y)
 109 #define Max(x,y) ((x > y)?x:y)
 110
 111 #ifdef SIXTEEN_COLORS
 112 #define MAXPOINTS  128
 113 #ifdef BIGMEM
 114 #define MAXFRAMES 4
 115 #else
 116 #define MAXFRAMES 2
 117 #endif
 118 #define MAXCOLOR 16
 119 static int maxcolor=16, startcolor=0, color_offset=0, mincolindex=1;
 120 static int dwell=50, settle=25;
 121 static int width=128, height=128, xposition=128, yposition=128;
 122 #else
 123 #define MAXPOINTS  256
 124 #ifdef BIGMEM
 125 #define MAXFRAMES 8
 126 #else
 127 #define MAXFRAMES 2
 128 #endif
 129 #define MAXCOLOR 256
 130 static int maxcolor=256, startcolor=17, color_offset=96, mincolindex=33;
 131 static int dwell=100, settle=50;
 132 static int width=256, height=256;
 133 #endif
 134
 135 #ifndef TRUE
 136 #define TRUE 1
 137 #define FALSE 0
 138 #endif
 139
 140 static int screen;
 141 static Display* dpy;
 142 static Visual *visual;
 143
 144 static unsigned long foreground, background;
 145
 146 static Window canvas;
 147
 148 typedef struct {
 149         int x, y;
 150 } xy_t;
 151
 152 typedef struct {
 153         int start_x, start_y;
 154         int last_x, last_y;
 155         } rubber_band_data_t;
 156
 157 typedef struct {
 158         Cursor band_cursor;
 159         double p_min, p_max, q_min, q_max;
 160         rubber_band_data_t rubber_band;
 161         } image_data_t;
 162
 163 typedef struct points_t {
 164         XPoint data[MAXCOLOR][MAXPOINTS];
 165         int npoints[MAXCOLOR];
 166         } points_t;
 167
 168 static points_t Points;
 169 static image_data_t rubber_data;
 170
 171 #ifndef TRUE
 172 #define TRUE 1
 173 #define FALSE 0
 174 #endif
 175
 176 static GC Data_GC[MAXCOLOR], RubberGC;
 177
 178 #define MAXINDEX 64
 179 #define FUNCMAXINDEX 16
 180 #define MAXWHEELS 7
 181 #define NUMMAPS 5
 182
 183 typedef double (*PFD)(double,double);
 184
 185 static double logistic(double,double), circle(double,double), leftlog(double,double), rightlog(double,double), doublelog(double,double);
 186 static double dlogistic(double,double), dcircle(double,double), dleftlog(double,double), drightlog(double,double), ddoublelog(double,double);
 187 static PFD map, deriv;
 188 static PFD Maps[NUMMAPS] = { logistic, circle, leftlog, rightlog, doublelog };
 189 static PFD Derivs[NUMMAPS] = { dlogistic, dcircle, dleftlog, drightlog, ddoublelog };
 190
 191 static int aflag=0, bflag=0, wflag=0, hflag=0, Rflag=0;
 192 static double pmins[NUMMAPS] = { 2.0, 0.0, 0.0, 0.0, 0.0 };
 193 static double pmaxs[NUMMAPS] = { 4.0, 1.0, 6.75, 6.75, 16.0 };
 194 static double amins[NUMMAPS] = { 2.0, 0.0, 0.0, 0.0, 0.0 };
 195 static double aranges[NUMMAPS] = { 2.0, 1.0, 6.75, 6.75, 16.0 };
 196 static double bmins[NUMMAPS] = { 2.0, 0.0, 0.0, 0.0, 0.0 };
 197 static double branges[NUMMAPS] = { 2.0, 1.0, 6.75, 6.75, 16.0 };
 198
 199 static int   forcing[MAXINDEX] = { 0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,
 200                         0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,
 201                         0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1 };
 202 static int   Forcing[FUNCMAXINDEX] = { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 };
 203
 204 static int   maxindex = MAXINDEX;
 205 static int   funcmaxindex = FUNCMAXINDEX;
 206 static double   min_a=2.0, min_b=2.0, a_range=2.0, b_range=2.0, minlyap=1.0;
 207 static double  max_a=4.0, max_b=4.0;
 208 static double  start_x=0.65, lyapunov, a_inc, b_inc, a, b;
 209 static int      numcolors=16, numfreecols, displayplanes, lowrange;
 210 static xy_t     point;
 211 static Pixmap  pixmap;
 212 static Colormap cmap;
 213 static XColor   Colors[MAXCOLOR];
 214 static double  *exponents[MAXFRAMES];
 215 static double  a_minimums[MAXFRAMES], b_minimums[MAXFRAMES];
 216 static double  a_maximums[MAXFRAMES], b_maximums[MAXFRAMES];
 217 static double  minexp, maxexp, prob=0.5;
 218 static int     expind[MAXFRAMES]={0}, resized[MAXFRAMES]={0};
 219 static int      numwheels=MAXWHEELS, force=0, Force=0, negative=1;
 220 static int     rgb_max=65000, nostart=1, stripe_interval=7;
 221 static int      save=1, show=0, useprod=1, spinlength=256, savefile=0;
 222 static int      maxframe=0, frame=0, dorecalc=0, mapindex=0, run=1;
 223 static char     *outname="lyap.out";
 224
 225
 226 const char * const version = LYAP_VERSION;
 227
 228 static void resize(void);
 229 static void redisplay(Window w, XExposeEvent *event);
 230 static void Spin(Window w);
 231 static void show_defaults(void);
 232 static void StartRubberBand(Window w, image_data_t *data, XEvent *event);
 233 static void TrackRubberBand(Window w, image_data_t *data, XEvent *event);
 234 static void EndRubberBand(Window w, image_data_t *data, XEvent *event);
 235 static void CreateXorGC(void);
 236 static void InitBuffer(void);
 237 static void BufferPoint(Display *display, Window window, int color,
 238                         int x, int y);
 239 static void FlushBuffer(void);
 240 static void init_canvas(void);
 241 static void init_data(void);
 242 static void init_color(void);
 243 static void parseargs(void);
 244 static void Clear(void);
 245 static void setupmem(void);
 246 static void main_event(void);
 247 static int complyap(void);
 248 static void Getkey(XKeyEvent *event);
 249 static int sendpoint(double expo);
 250 static void save_to_file(void);
 251 static void setforcing(void);
 252 static void check_params(int mapnum, int parnum);
 253 static void usage(void);
 254 static void Destroy_frame(void);
 255 static void freemem(void);
 256 static void Redraw(void);
 257 static void redraw(double *exparray, int index, int cont);
 258 static void recalc(void);
 259 static void SetupCorners(XPoint *corners, image_data_t *data);
 260 static void set_new_params(Window w, image_data_t *data);
 261 static void go_down(void);
 262 static void go_back(void);
 263 static void go_init(void);
 264 static void jumpwin(void);
 265 static void print_help(void);
 266 static void print_values(void);
 267
 268
 269 void
 270 screenhack (Display *d, Window window)
 271 {
 272   XWindowAttributes xgwa;
 273   int builtin = -1;
 274   dpy = d;
 275   XGetWindowAttributes (dpy, window, &xgwa);
 276   width = xgwa.width;
 277   height = xgwa.height;
 278   visual = xgwa.visual;
 279   cmap = xgwa.colormap;
 280
 281   parseargs();
 282
 283   if (get_boolean_resource("randomize", "Boolean"))
 284     builtin = random() % 22;
 285   else {
 286     char *s = get_string_resource("builtin", "Integer");
 287     if (s && *s)
 288       builtin = atoi(s);
 289     if (s) free (s);
 290   }
 291
 292   if (builtin >= 0)
 293     {
 294       char *ff = 0;
 295       switch (builtin) {
 296       case 0:
 297         min_a = 3.75; aflag++;
 298         min_b = 3.299999; bflag++;
 299         a_range = 0.05; wflag++;
 300         b_range = 0.05; hflag++;
 301         dwell = 200;
 302         settle = 100;
 303         ff = "abaabbaaabbb";
 304         break;
 305
 306       case 1:
 307         min_a = 3.8; aflag++;
 308         min_b = 3.2; bflag++;
 309         b_range = .05; hflag++;
 310         a_range = .05; wflag++;
 311         ff = "bbbbbaaaaa";
 312         break;
 313
 314       case 2:
 315         min_a =  3.4; aflag++;
 316         min_b =  3.04; bflag++;
 317         a_range =  .5; wflag++;
 318         b_range =  .5; hflag++;
 319         ff = "abbbbbbbbb";
 320         settle = 500;
 321         dwell = 1000;
 322         break;
 323
 324       case 3:
 325         min_a = 3.5; aflag++;
 326         min_b = 3.0; bflag++;
 327         a_range = 0.2; wflag++;
 328         b_range = 0.2; hflag++;
 329         dwell = 600;
 330         settle = 300;
 331         ff = "aaabbbab";
 332         break;
 333
 334       case 4:
 335         min_a = 3.55667; aflag++;
 336         min_b = 3.2; bflag++;
 337         b_range = .05; hflag++;
 338         a_range = .05; wflag++;
 339         ff = "bbbbbaaaaa";
 340         break;
 341
 342       case 5:
 343         min_a = 3.79; aflag++;
 344         min_b = 3.22; bflag++;
 345         b_range = .02999; hflag++;
 346         a_range = .02999; wflag++;
 347         ff = "bbbbbaaaaa";
 348         break;
 349
 350       case 6:
 351         min_a = 3.7999; aflag++;
 352         min_b = 3.299999; bflag++;
 353         a_range = 0.2; wflag++;
 354         b_range = 0.2; hflag++;
 355         dwell = 300;
 356         settle = 150;
 357         ff = "abaabbaaabbb";
 358         break;
 359
 360       case 7:
 361         min_a = 3.89; aflag++;
 362         min_b = 3.22; bflag++;
 363         b_range = .028; hflag++;
 364         a_range = .02999; wflag++;
 365         ff = "bbbbbaaaaa";
 366         settle = 600;
 367         dwell = 1000;
 368         break;
 369
 370       case 8:
 371         min_a = 3.2; aflag++;
 372         min_b = 3.7; bflag++;
 373         a_range = 0.05; wflag++;
 374         b_range = .005; hflag++;
 375         ff = "abbbbaa";
 376         break;
 377
 378       case 9:
 379         ff = "aaaaaabbbbbb";
 380         mapindex = 1;
 381         dwell =  400;
 382         settle =  200;
 383         minlyap = maxexp = ABS(-0.85);
 384         minexp = -1.0 * minlyap;
 385         break;
 386
 387       case 10:
 388         ff = "aaaaaabbbbbb";
 389         mapindex = 1;
 390         dwell =  400;
 391         settle = 200;
 392         minlyap = maxexp = ABS(-0.85);
 393         minexp = -1.0 * minlyap;
 394         break;
 395
 396       case 11:
 397         mapindex = 1;
 398         dwell =  400;
 399         settle = 200;
 400         minlyap = maxexp = ABS(-0.85);
 401         minexp = -1.0 * minlyap;
 402         break;
 403
 404       case 12:
 405         ff = "abbb";
 406         mapindex = 1;
 407         dwell =  400;
 408         settle = 200;
 409         minlyap = maxexp = ABS(-0.85);
 410         minexp = -1.0 * minlyap;
 411         break;
 412
 413       case 13:
 414         ff = "abbabaab";
 415         mapindex = 1;
 416         dwell =  400;
 417         settle = 200;
 418         minlyap = maxexp = ABS(-0.85);
 419         minexp = -1.0 * minlyap;
 420         break;
 421
 422       case 14:
 423         ff = "abbabaab";
 424         dwell =  800;
 425         settle = 200;
 426         minlyap = maxexp = ABS(-0.85);
 427         minexp = -1.0 * minlyap;
 428         /* ####  -x 0.05 */
 429         min_a = 3.91; aflag++;
 430         a_range =  0.0899999999; wflag++;
 431         min_b =  3.28; bflag++;
 432         b_range =  0.35; hflag++;
 433         break;
 434
 435       case 15:
 436         ff = "aaaaaabbbbbb";
 437         dwell =  400;
 438         settle = 200;
 439         minlyap = maxexp = ABS(-0.85);
 440         minexp = -1.0 * minlyap;
 441         break;
 442
 443       case 16:
 444         dwell =  400;
 445         settle = 200;
 446         minlyap = maxexp = ABS(-0.85);
 447         minexp = -1.0 * minlyap;
 448         break;
 449
 450       case 17:
 451         ff = "abbb";
 452         dwell =  400;
 453         settle = 200;
 454         minlyap = maxexp = ABS(-0.85);
 455         minexp = -1.0 * minlyap;
 456         break;
 457
 458       case 18:
 459         ff = "abbabaab";
 460         dwell =  400;
 461         settle = 200;
 462         minlyap = maxexp = ABS(-0.85);
 463         minexp = -1.0 * minlyap;
 464         break;
 465
 466       case 19:
 467         mapindex = 2;
 468         ff = "aaaaaabbbbbb";
 469         dwell =  400;
 470         settle = 200;
 471         minlyap = maxexp = ABS(-0.85);
 472         minexp = -1.0 * minlyap;
 473         break;
 474
 475       case 20:
 476         mapindex = 2;
 477         dwell =  400;
 478         settle = 200;
 479         minlyap = maxexp = ABS(-0.85);
 480         minexp = -1.0 * minlyap;
 481         break;
 482
 483       case 21:
 484         mapindex = 2;
 485         ff = "abbb";
 486         dwell =  400;
 487         settle = 200;
 488         minlyap = maxexp = ABS(-0.85);
 489         minexp = -1.0 * minlyap;
 490         break;
 491
 492       case 22:
 493         mapindex = 2;
 494         ff = "abbabaab";
 495         dwell =  400;
 496         settle = 200;
 497         minlyap = maxexp = ABS(-0.85);
 498         minexp = -1.0 * minlyap;
 499         break;
 500       }
 501
 502       if (ff) {
 503         char *ch;
 504         int bindex = 0;
 505         maxindex = strlen(ff);
 506         if (maxindex > MAXINDEX)
 507           usage();
 508         ch = ff;
 509         force++;
 510         while (bindex < maxindex) {
 511           if (*ch == 'a')
 512             forcing[bindex++] = 0;
 513           else if (*ch == 'b')
 514             forcing[bindex++] = 1;
 515           else
 516             usage();
 517           ch++;
 518         }
 519       }
 520     }
 521
 522   screen = DefaultScreen(dpy);
 523   background = BlackPixel(dpy, screen);
 524   setupmem();
 525   init_data();
 526   if (displayplanes > 1)
 527     foreground = startcolor;
 528   else
 529     foreground = WhitePixel(dpy, screen);
 530
 531   /*
 532   * Create the window to display the Lyapunov exponents
 533   */
 534   canvas = window;
 535   init_canvas();
 536
 537   if (window != DefaultRootWindow(dpy))
 538     XSelectInput(dpy,canvas,KeyPressMask|ButtonPressMask|ButtonMotionMask|
 539                  ButtonReleaseMask|ExposureMask|StructureNotifyMask);
 540   if (displayplanes > 1) {
 541     init_color();
 542   } else {
 543     XQueryColors(dpy, DefaultColormap(dpy, DefaultScreen(dpy)),
 544         Colors, numcolors);
 545   }
 546   pixmap = XCreatePixmap(dpy, DefaultRootWindow(dpy),
 547                          width, height, DefaultDepth(dpy, screen));
 548   rubber_data.band_cursor = XCreateFontCursor(dpy, XC_hand2);
 549   CreateXorGC();
 550   Clear();
 551   for(;;)
 552       main_event();
 553 }
 554
 555 static void
 556 main_event(void)
 557 {
 558   int n;
 559   XEvent event;
 560
 561   if (complyap() == TRUE)
 562       run=0;
 563   n = XEventsQueued(dpy, QueuedAfterFlush);
 564   while (n--) {
 565           XNextEvent(dpy, &event);
 566             switch(event.type)
 567             {
 568             case KeyPress:
 569     Getkey(&event.xkey);
 570     break;
 571             case Expose:
 572     redisplay(canvas, &event.xexpose);
 573           break;
 574             case ConfigureNotify:
 575     resize();
 576           break;
 577             case ButtonPress:
 578     StartRubberBand(canvas, &rubber_data, &event);
 579           break;
 580             case MotionNotify:
 581     TrackRubberBand(canvas, &rubber_data, &event);
 582           break;
 583             case ButtonRelease:
 584     EndRubberBand(canvas, &rubber_data, &event);
 585           break;
 586             default:
 587     screenhack_handle_event (dpy, &event);
 588           break;
 589             }
 590         }
 591 }
 592
 593 /* complyap() is the guts of the program. This is where the Lyapunov exponent
 594  * is calculated. For each iteration (past some large number of iterations)
 595  * calculate the logarithm of the absolute value of the derivative at that
 596  * point. Then average them over some large number of iterations. Some small
 597  * speed up is achieved by utilizing the fact that log(a*b) = log(a) + log(b).
 598  */
 599 static int
 600 complyap(void)
 601 {
 602   int i, bindex;
 603   double total, prod, x, dx, r;
 604
 605   if (!run)
 606     return TRUE;
 607   a += a_inc;
 608   if (a >= max_a)
 609     if (sendpoint(lyapunov) == TRUE)
 610       return FALSE;
 611     else {
 612       FlushBuffer();
 613       if (savefile)
 614         save_to_file();
 615       return TRUE;
 616     }
 617   if (b >= max_b) {
 618     FlushBuffer();
 619     if (savefile)
 620       save_to_file();
 621     return TRUE;
 622   }
 623   prod = 1.0;
 624   total = 0.0;
 625   bindex = 0;
 626   x = start_x;
 627   r = (forcing[bindex]) ? b : a;
 628 #ifdef MAPS
 629   findex = 0;
 630   map = Maps[Forcing[findex]];
 631 #endif
 632   for (i=0;i<settle;i++) {     /* Here's where we let the thing */
 633     x = (*map)(x, r);    /* "settle down". There is usually */
 634     if (++bindex >= maxindex) { /* some initial "noise" in the */
 635       bindex = 0;    /* iterations. How can we optimize */
 636       if (Rflag)      /* the value of settle ??? */
 637           setforcing();
 638     }
 639     r = (forcing[bindex]) ? b : a;
 640 #ifdef MAPS
 641     if (++findex >= funcmaxindex)
 642       findex = 0;
 643     map = Maps[Forcing[findex]];
 644 #endif
 645   }
 646 #ifdef MAPS
 647   deriv = Derivs[Forcing[findex]];
 648 #endif
 649   if (useprod) {      /* using log(a*b) */
 650     for (i=0;i<dwell;i++) {
 651       x = (*map)(x, r);
 652       dx = (*deriv)(x, r); /* ABS is a macro, so don't be fancy */
 653       dx = ABS(dx);
 654       if (dx == 0.0) /* log(0) is nasty so break out. */
 655       {
 656         i++;
 657         break;
 658       }
 659       prod *= dx;
 660       /* we need to prevent overflow and underflow */
 661       if ((prod > 1.0e12) || (prod < 1.0e-12)) {
 662         total += log(prod);
 663         prod = 1.0;
 664       }
 665       if (++bindex >= maxindex) {
 666         bindex = 0;
 667         if (Rflag)
 668           setforcing();
 669       }
 670       r = (forcing[bindex]) ? b : a;
 671 #ifdef MAPS
 672       if (++findex >= funcmaxindex)
 673         findex = 0;
 674       map = Maps[Forcing[findex]];
 675       deriv = Derivs[Forcing[findex]];
 676 #endif
 677     }
 678     total += log(prod);
 679     lyapunov = (total * M_LOG2E) / (double)i;
 680   }
 681   else {        /* use log(a) + log(b) */
 682     for (i=0;i<dwell;i++) {
 683       x = (*map)(x, r);
 684       dx = (*deriv)(x, r); /* ABS is a macro, so don't be fancy */
 685       dx = ABS(dx);
 686       if (x == 0.0)  /* log(0) check */
 687       {
 688         i++;
 689         break;
 690       }
 691       total += log(dx);
 692       if (++bindex >= maxindex) {
 693         bindex = 0;
 694         if (Rflag)
 695           setforcing();
 696       }
 697       r = (forcing[bindex]) ? b : a;
 698 #ifdef MAPS
 699       if (++findex >= funcmaxindex)
 700         findex = 0;
 701       map = Maps[Forcing[findex]];
 702       deriv = Derivs[Forcing[findex]];
 703 #endif
 704     }
 705     lyapunov = (total * M_LOG2E) / (double)i;
 706   }
 707
 708   if (sendpoint(lyapunov) == TRUE)
 709     return FALSE;
 710   else {
 711     FlushBuffer();
 712     if (savefile)
 713       save_to_file();
 714     return TRUE;
 715   }
 716 }
 717
 718 static double
 719 logistic(double x, double r)        /* the familiar logistic map */
 720 {
 721   return(r * x * (1.0 - x));
 722 }
 723
 724 static double
 725 dlogistic(double x, double r)       /* the derivative of logistic map */
 726 {
 727   return(r - (2.0 * r * x));
 728 }
 729
 730 static double
 731 circle(double x, double r)        /* sin() hump or sorta like the circle map */
 732 {
 733   return(r * sin(M_PI * x));
 734 }
 735
 736 static double
 737 dcircle(double x, double r)       /* derivative of the "sin() hump" */
 738 {
 739   return(r * M_PI * cos(M_PI * x));
 740 }
 741
 742 static double
 743 leftlog(double x, double r)       /* left skewed logistic */
 744 {
 745   double d;
 746
 747   d = 1.0 - x;
 748   return(r * x * d * d);
 749 }
 750
 751 static double
 752 dleftlog(double x, double r)    /* derivative of the left skewed logistic */
 753 {
 754   return(r * (1.0 - (4.0 * x) + (3.0 * x * x)));
 755 }
 756
 757 static double
 758 rightlog(double x, double r)    /* right skewed logistic */
 759 {
 760   return(r * x * x * (1.0 - x));
 761 }
 762
 763 static double
 764 drightlog(double x, double r)    /* derivative of the right skewed logistic */
 765 {
 766   return(r * ((2.0 * x) - (3.0 * x * x)));
 767 }
 768
 769 static double
 770 doublelog(double x, double r)    /* double logistic */
 771 {
 772   double d;
 773
 774   d = 1.0 - x;
 775   return(r * x * x * d * d);
 776 }
 777
 778 static double
 779 ddoublelog(double x, double r)   /* derivative of the double logistic */
 780 {
 781   double d;
 782
 783   d = x * x;
 784   return(r * ((2.0 * x) - (6.0 * d) + (4.0 * x * d)));
 785 }
 786
 787 static void
 788 init_data(void)
 789 {
 790   numcolors = XDisplayCells(dpy, XDefaultScreen(dpy));
 791   displayplanes = DisplayPlanes(dpy, XDefaultScreen(dpy));
 792   if (numcolors > maxcolor)
 793     numcolors = maxcolor;
 794   numfreecols = numcolors - mincolindex;
 795   lowrange = mincolindex - startcolor;
 796   a_inc = a_range / (double)width;
 797   b_inc = b_range / (double)height;
 798   point.x = -1;
 799   point.y = 0;
 800   a = rubber_data.p_min = min_a;
 801   b = rubber_data.q_min = min_b;
 802   rubber_data.p_max = max_a;
 803   rubber_data.q_max = max_b;
 804   if (show)
 805     show_defaults();
 806   InitBuffer();
 807 }
 808
 809 static void
 810 init_canvas(void)
 811 {
 812   static int i;
 813
 814   /*
 815   * create default, writable, graphics contexts for the canvas.
 816   */
 817         for (i=0; i<maxcolor; i++) {
 818             Data_GC[i] = XCreateGC(dpy, DefaultRootWindow(dpy),
 819                 (unsigned long) NULL, (XGCValues *) NULL);
 820             /* set the background to black */
 821             XSetBackground(dpy,Data_GC[i],BlackPixel(dpy,XDefaultScreen(dpy)));
 822             /* set the foreground of the ith context to i */
 823             XSetForeground(dpy, Data_GC[i], i);
 824         }
 825         if (displayplanes == 1) {
 826             XSetForeground(dpy,Data_GC[0],BlackPixel(dpy,XDefaultScreen(dpy)));
 827             XSetForeground(dpy,Data_GC[1],WhitePixel(dpy,XDefaultScreen(dpy)));
 828         }
 829 }
 830
 831 #if 0
 832 static void
 833 hls2rgb(int hue_light_sat[3],
 834         int rgb[3])             /*      Each in range [0..65535]        */
 835 {
 836   unsigned short r, g, b;
 837   hsv_to_rgb((int) (hue_light_sat[0] / 10),             /* 0-3600 -> 0-360 */
 838              (int) ((hue_light_sat[2]/1000.0) * 64435), /* 0-1000 -> 0-65535 */
 839              (int) ((hue_light_sat[1]/1000.0) * 64435), /* 0-1000 -> 0-65535 */
 840              &r, &g, &b);
 841   rgb[0] = r;
 842   rgb[1] = g;
 843   rgb[2] = b;
 844 }
 845 #endif /* 0 */
 846
 847
 848 static void
 849 init_color(void)
 850 {
 851 #if 1
 852
 853   int i;
 854   XColor colors[256];
 855   int ncolors = maxcolor;
 856   Bool writable = False;
 857   make_smooth_colormap(dpy, visual, cmap,
 858                         colors, &ncolors, True, &writable, True);
 859
 860   for (i = 0; i < maxcolor; i++)
 861     XSetForeground(dpy, Data_GC[i],
 862                    colors[((int) ((i / ((float)maxcolor)) * ncolors))].pixel);
 863
 864 #else
 865   static int i, j, colgap, leg, step;
 866   static Visual *visual;
 867   Colormap def_cmap;
 868   int hls[3], rgb[3];
 869
 870   def_cmap = DefaultColormap(dpy, DefaultScreen(dpy));
 871   for (i=0; i<numcolors; i++) {
 872     Colors[i].pixel = i;
 873     Colors[i].flags = DoRed|DoGreen|DoBlue;
 874   }
 875
 876   /* Try to write into a new color map */
 877   visual = DefaultVisual(dpy, DefaultScreen(dpy));
 878   cmap = XCreateColormap(dpy, canvas, visual, AllocAll);
 879   XQueryColors(dpy, def_cmap, Colors, numcolors);
 880   if (mincolindex)
 881     colgap = rgb_max / mincolindex;
 882   else
 883     colgap = rgb_max;
 884   hls[0] = 50;  /* Hue in low range */
 885   hls[2] = 1000;  /* Fully saturated */
 886   for (i=startcolor; i<lowrange + startcolor; i++) {
 887     hls[1] = 1000L * (i-startcolor) / lowrange;
 888     hls2rgb(hls, rgb);
 889     Colors[i].red = rgb[0];
 890     Colors[i].green = rgb[1];
 891     Colors[i].blue = rgb[2];
 892   }
 893   colgap = rgb_max / numcolors;
 894   if (numwheels == 0)
 895     XQueryColors(dpy, def_cmap, Colors, numcolors);
 896   else if (numwheels == 1) {
 897     colgap = 2*rgb_max/(numcolors - color_offset);
 898     for (i=mincolindex; i<(numcolors/2); i++) {
 899       Colors[i].blue = 0;
 900       Colors[i].green=((i+color_offset)*colgap);
 901       Colors[i].red=((i+color_offset)*colgap);
 902     }
 903     for (i=(numcolors/2); i<(numcolors); i++) {
 904       Colors[i].blue = 0;
 905       Colors[i].green=(((numcolors-i)+color_offset)*colgap);
 906       Colors[i].red=(((numcolors-i)+color_offset)*colgap);
 907     }
 908   }
 909   else if (numwheels == 2) {
 910           hls[0] = 800; /* Hue in mid range */
 911           hls[2] = 1000;  /* Fully saturated */
 912           for (i=startcolor; i<lowrange + startcolor; i++) {
 913       hls[1] = 1000L * (i-startcolor) / lowrange;
 914       hls2rgb(hls, rgb);
 915       Colors[i].red = rgb[0];
 916       Colors[i].green = rgb[1];
 917       Colors[i].blue = rgb[2];
 918           }
 919     for (i=mincolindex; i<(numcolors/2); i++) {
 920       Colors[i].blue = rgb_max;
 921       Colors[i].green = 0;
 922       Colors[i].red=(i*2*rgb_max/numcolors);
 923     }
 924     for (i=(numcolors/2); i<numcolors; i++) {
 925       Colors[i].blue = rgb_max;
 926       Colors[i].green = 0;
 927       Colors[i].red=((numcolors - i)*2*rgb_max/numcolors);
 928     }
 929   }
 930   else if (numwheels == 3) {
 931           hls[0] = 800; /* Hue in mid range */
 932           hls[2] = 1000;  /* Fully saturated */
 933           for (i=startcolor; i<lowrange + startcolor; i++) {
 934       hls[1] = 1000L * (i-startcolor) / lowrange;
 935       hls2rgb(hls, rgb);
 936       Colors[i].red = rgb[0];
 937       Colors[i].green = rgb[1];
 938       Colors[i].blue = rgb[2];
 939           }
 940     colgap = 4*rgb_max/numcolors;
 941     for (i=mincolindex; i<(numcolors/4); i++) {
 942       Colors[i].blue = rgb_max;
 943       Colors[i].green = 0;
 944       Colors[i].red=(i*colgap);
 945     }
 946     for (i=(numcolors/4); i<(numcolors/2); i++) {
 947       Colors[i].red = rgb_max;
 948       Colors[i].green = 0;
 949       Colors[i].blue=((numcolors/2) - i) * colgap;
 950     }
 951     for (i=(numcolors/2); i<(0.75*numcolors); i++) {
 952       Colors[i].red = rgb_max;
 953       Colors[i].blue=(i * colgap);
 954       Colors[i].green = 0;
 955     }
 956     for (i=(0.75*numcolors); i<numcolors; i++) {
 957       Colors[i].blue = rgb_max;
 958       Colors[i].green = 0;
 959       Colors[i].red=(numcolors-i)*colgap;
 960     }
 961   }
 962   else if (numwheels == 4) {
 963           hls[0] = 800; /* Hue in mid range */
 964           hls[2] = 1000;  /* Fully saturated */
 965           for (i=startcolor; i<lowrange + startcolor; i++) {
 966       hls[1] = 1000L * (i-startcolor) / lowrange;
 967       hls2rgb(hls, rgb);
 968       Colors[i].red = rgb[0];
 969       Colors[i].green = rgb[1];
 970       Colors[i].blue = rgb[2];
 971           }
 972     colgap = numwheels * rgb_max / numcolors;
 973     for (i=mincolindex; i<(numcolors/numwheels); i++) {
 974       Colors[i].blue = rgb_max;
 975       Colors[i].green = 0;
 976       Colors[i].red=(i*colgap);
 977     }
 978     for (i=(numcolors/numwheels); i<(2*numcolors/numwheels); i++) {
 979       Colors[i].red = rgb_max;
 980       Colors[i].green = 0;
 981       Colors[i].blue=((2*numcolors/numwheels) - i) * colgap;
 982     }
 983     for (i=(2*numcolors/numwheels); i<numcolors; i++) {
 984       Colors[i].red = rgb_max;
 985       Colors[i].green=(i - (2*numcolors/numwheels)) * colgap;
 986       Colors[i].blue = 0;
 987     }
 988   }
 989   else if (numwheels == 5) {
 990     hls[1] = 700; /* Lightness in midrange */
 991     hls[2] = 1000;  /* Fully saturated */
 992     for (i=mincolindex; i<numcolors; i++) {
 993       hls[0] = 3600L * i / numcolors;
 994       hls2rgb(hls, rgb);
 995       Colors[i].red = rgb[0];
 996       Colors[i].green = rgb[1];
 997       Colors[i].blue = rgb[2];
 998     }
 999     for (i=mincolindex; i<numcolors; i+=stripe_interval) {
1000       hls[0] = 3600L * i / numcolors;
1001       hls2rgb(hls, rgb);
1002       Colors[i].red = rgb[0] / 2;
1003       Colors[i].green = rgb[1] / 2;
1004       Colors[i].blue = rgb[2] / 2;
1005     }
1006   }
1007   else if (numwheels == 6) {
1008       hls[0] = 800; /* Hue in mid range */
1009       hls[2] = 1000;  /* Fully saturated */
1010       for (i=startcolor; i<lowrange + startcolor; i++) {
1011     hls[1] = 1000L * (i-startcolor) / lowrange;
1012     hls2rgb(hls, rgb);
1013     Colors[i].red = rgb[0];
1014     Colors[i].green = rgb[1];
1015     Colors[i].blue = rgb[2];
1016       }
1017       step = numfreecols / 3;
1018       leg = step+mincolindex;
1019       for (i = mincolindex; i < leg; ++i)
1020       {
1021     Colors[i].pixel = i;
1022     Colors[i].red = fabs(65535 - (double)i / step * 65535.0);
1023     Colors[i].blue = (double)i / step * 65535.0;
1024     Colors[i].green = 0;
1025     Colors[i].flags = DoRed | DoGreen | DoBlue;
1026       }
1027       for (j = 0, i = leg, leg += step; i < leg; ++i, ++j)
1028       {
1029     Colors[i].pixel = i;
1030     Colors[i].red = (double)j / step * 65535.0;
1031     Colors[i].blue = 65535;
1032     Colors[i].green = Colors[i].red;
1033     Colors[i].flags = DoRed | DoGreen | DoBlue;
1034       }
1035       for (j = 0, i = leg, leg += step; i < leg; ++i, ++j)
1036       {
1037     Colors[i].pixel = i;
1038     Colors[i].red = 65535;
1039     Colors[i].blue = fabs(65535 - (double)j / step * 65535.0);
1040     Colors[i].green = Colors[i].blue;
1041     Colors[i].flags = DoRed | DoGreen | DoBlue;
1042       }
1043   }
1044   else if (numwheels == MAXWHEELS) {  /* rainbow palette */
1045     hls[1] = 500; /* Lightness in midrange */
1046     hls[2] = 1000;  /* Fully saturated */
1047     for (i=mincolindex; i<numcolors; i++) {
1048       hls[0] = 3600L * i / numcolors;
1049       hls2rgb(hls, rgb);
1050       Colors[i].red = rgb[0];
1051       Colors[i].green = rgb[1];
1052       Colors[i].blue = rgb[2];
1053     }
1054   }
1055   XStoreColors(dpy, cmap, Colors, numcolors);
1056
1057   XSetWindowColormap(dpy, canvas, cmap);
1058 #endif
1059 }
1060
1061 static void
1062 parseargs()
1063 {
1064   static int i;
1065   int bindex=0, findex;
1066   char *s, *ch;
1067
1068   map = Maps[0];
1069   deriv = Derivs[0];
1070   maxexp=minlyap; minexp= -1.0 * minlyap;
1071
1072   mincolindex = get_integer_resource("minColor", "Integer");
1073   dwell = get_integer_resource("dwell", "Integer");
1074 #ifdef MAPS
1075   {
1076     char *optarg = get_string_resource("function", "String");
1077     funcmaxindex = strlen(optarg);
1078     if (funcmaxindex > FUNCMAXINDEX)
1079       usage();
1080     ch = optarg;
1081     Force++;
1082     for (findex=0;findex<funcmaxindex;findex++) {
1083       Forcing[findex] = (int)(*ch++ - '0');;
1084       if (Forcing[findex] >= NUMMAPS)
1085         usage();
1086     }
1087   }
1088 #endif
1089   if (get_boolean_resource("useLog", "Boolean"))
1090     useprod=0;
1091
1092   minlyap=ABS(get_float_resource("colorExponent", "Float"));
1093   maxexp=minlyap;
1094   minexp= -1.0 * minlyap;
1095
1096   color_offset = get_integer_resource("colorOffset", "Integer");
1097
1098   maxcolor=ABS(get_integer_resource("maxColor", "Integer"));
1099   if ((maxcolor - startcolor) <= 0)
1100     startcolor = 0;
1101   if ((maxcolor - mincolindex) <= 0) {
1102     mincolindex = 1;
1103     color_offset = 0;
1104   }
1105
1106   s = get_string_resource("randomForce", "Float");
1107   if (s && *s) {
1108     prob=atof(s); Rflag++; setforcing();
1109   }
1110
1111   settle = get_integer_resource("settle", "Integer");
1112
1113   s = get_string_resource("minA", "Float");
1114   if (s && *s) {
1115     min_a = atof(s);
1116     aflag++;
1117   }
1118
1119   s = get_string_resource("minB", "Float");
1120   if (s && *s) {
1121     min_b=atof(s); bflag++;
1122   }
1123
1124   numwheels = get_integer_resource("wheels", "Integer");
1125
1126   s = get_string_resource("forcingFunction", "String");
1127   if (s && *s) {
1128     maxindex = strlen(s);
1129     if (maxindex > MAXINDEX)
1130       usage();
1131     ch = s;
1132     force++;
1133     while (bindex < maxindex) {
1134       if (*ch == 'a')
1135         forcing[bindex++] = 0;
1136       else if (*ch == 'b')
1137         forcing[bindex++] = 1;
1138       else
1139         usage();
1140       ch++;
1141     }
1142   }
1143
1144   s = get_string_resource("bRange", "Float");
1145   if (s && *s) {
1146     b_range = atof(s);
1147     hflag++;
1148   }
1149
1150   start_x = get_float_resource("startX", "Float");
1151
1152   s = get_string_resource("mapIndex", "Integer");
1153   if (s && *s) {
1154     mapindex=atoi(s);
1155     if ((mapindex >= NUMMAPS) || (mapindex < 0))
1156       usage();
1157     map = Maps[mapindex];
1158     deriv = Derivs[mapindex];
1159     if (!aflag)
1160       min_a = amins[mapindex];
1161     if (!wflag)
1162       a_range = aranges[mapindex];
1163     if (!bflag)
1164       min_b = bmins[mapindex];
1165     if (!hflag)
1166       b_range = branges[mapindex];
1167     if (!Force)
1168       for (i=0;i<FUNCMAXINDEX;i++)
1169         Forcing[i] = mapindex;
1170   }
1171
1172   outname = get_string_resource("outputFile", "Integer");
1173
1174   if (get_boolean_resource("beNegative", "Boolean"))
1175     negative--;
1176
1177   rgb_max = get_integer_resource("rgbMax", "Integer");
1178   spinlength = get_integer_resource("spinLength", "Integer");
1179   show = get_boolean_resource("show", "Boolean");
1180
1181   s = get_string_resource("aRange", "Float");
1182   if (s && *s) {
1183     a_range = atof(s); wflag++;
1184   }
1185
1186   max_a = min_a + a_range;
1187   max_b = min_b + b_range;
1188
1189   a_minimums[0] = min_a; b_minimums[0] = min_b;
1190   a_maximums[0] = max_a; b_maximums[0] = max_b;
1191
1192   if (Force)
1193     if (maxindex == funcmaxindex)
1194       for (findex=0;findex<funcmaxindex;findex++)
1195         check_params(Forcing[findex],forcing[findex]);
1196     else
1197       fprintf(stderr, "Warning! Unable to check parameters\n");
1198   else
1199     check_params(mapindex,2);
1200 }
1201
1202 static void
1203 check_params(int mapnum, int parnum)
1204 {
1205
1206   if (parnum != 1) {
1207       if ((max_a > pmaxs[mapnum]) || (min_a < pmins[mapnum])) {
1208     fprintf(stderr, "Warning! Parameter 'a' out of range.\n");
1209     fprintf(stderr, "You have requested a range of (%f,%f).\n",
1210       min_a,max_a);
1211     fprintf(stderr, "Valid range is (%f,%f).\n",
1212       pmins[mapnum],pmaxs[mapnum]);
1213       }
1214   }
1215   if (parnum != 0) {
1216       if ((max_b > pmaxs[mapnum]) || (min_b < pmins[mapnum])) {
1217     fprintf(stderr, "Warning! Parameter 'b' out of range.\n");
1218     fprintf(stderr, "You have requested a range of (%f,%f).\n",
1219       min_b,max_b);
1220     fprintf(stderr, "Valid range is (%f,%f).\n",
1221       pmins[mapnum],pmaxs[mapnum]);
1222       }
1223   }
1224 }
1225
1226 static void
1227 usage(void)
1228 {
1229     fprintf(stderr,"lyap [-BLs][-W#][-H#][-a#][-b#][-w#][-h#][-x xstart]\n");
1230     fprintf(stderr,"\t[-M#][-S#][-D#][-f string][-r#][-O#][-C#][-c#][-m#]\n");
1231 #ifdef MAPS
1232     fprintf(stderr,"\t[-F string]\n");
1233 #endif
1234     fprintf(stderr,"\tWhere: -C# specifies the minimum color index\n");
1235     fprintf(stderr,"\t       -r# specifies the maxzimum rgb value\n");
1236     fprintf(stderr,"\t       -u displays this message\n");
1237     fprintf(stderr,"\t       -a# specifies the minimum horizontal parameter\n");
1238     fprintf(stderr,"\t       -b# specifies the minimum vertical parameter\n");
1239     fprintf(stderr,"\t       -w# specifies the horizontal parameter range\n");
1240     fprintf(stderr,"\t       -h# specifies the vertical parameter range\n");
1241     fprintf(stderr,"\t       -D# specifies the dwell\n");
1242     fprintf(stderr,"\t       -S# specifies the settle\n");
1243     fprintf(stderr,"\t       -H# specifies the initial window height\n");
1244     fprintf(stderr,"\t       -W# specifies the initial window width\n");
1245     fprintf(stderr,"\t       -O# specifies the color offset\n");
1246     fprintf(stderr,"\t       -c# specifies the desired color wheel\n");
1247     fprintf(stderr,"\t       -m# specifies the desired map (0-4)\n");
1248     fprintf(stderr,"\t       -f aabbb specifies a forcing function of 00111\n");
1249 #ifdef MAPS
1250     fprintf(stderr,"\t       -F 00111 specifies the function forcing function\n");
1251 #endif
1252     fprintf(stderr,"\t       -L indicates use log(x)+log(y) rather than log(xy)\n");
1253     fprintf(stderr,"\tDuring display :\n");
1254     fprintf(stderr,"\t     Use the mouse to zoom in on an area\n");
1255     fprintf(stderr,"\t     e or E recalculates color indices\n");
1256     fprintf(stderr,"\t     f or F saves exponents to a file\n");
1257     fprintf(stderr,"\t     KJmn increase/decrease minimum negative exponent\n");
1258     fprintf(stderr,"\t     r or R redraws\n");
1259     fprintf(stderr,"\t     s or S spins the colorwheel\n");
1260     fprintf(stderr,"\t     w or W changes the color wheel\n");
1261     fprintf(stderr,"\t     x or X clears the window\n");
1262     fprintf(stderr,"\t     q or Q exits\n");
1263     exit(1);
1264 }
1265
1266 static void
1267 Cycle_frames(void)
1268 {
1269   static int i;
1270   for (i=0;i<=maxframe;i++)
1271     redraw(exponents[i], expind[i], 1);
1272 }
1273
1274 static void
1275 Spin(Window w)
1276 {
1277   static int i, j;
1278   long tmpxcolor;
1279
1280   if (displayplanes > 1) {
1281     for (j=0;j<spinlength;j++) {
1282       tmpxcolor = Colors[mincolindex].pixel;
1283       for (i=mincolindex;i<numcolors-1;i++)
1284         Colors[i].pixel = Colors[i+1].pixel;
1285       Colors[numcolors-1].pixel = tmpxcolor;
1286       XStoreColors(dpy, cmap, Colors, numcolors);
1287     }
1288     for (j=0;j<spinlength;j++) {
1289       tmpxcolor = Colors[numcolors-1].pixel;
1290       for (i=numcolors-1;i>mincolindex;i--)
1291         Colors[i].pixel = Colors[i-1].pixel;
1292       Colors[mincolindex].pixel = tmpxcolor;
1293       XStoreColors(dpy, cmap, Colors, numcolors);
1294     }
1295   }
1296 }
1297
1298 static void
1299 Getkey(XKeyEvent *event)
1300 {
1301   unsigned char key;
1302   static int i;
1303   if (XLookupString(event, (char *)&key, sizeof(key), (KeySym *)0,
1304             (XComposeStatus *) 0) > 0)
1305     switch (key) {
1306   case '<': dwell /= 2; if (dwell < 1) dwell = 1; break;
1307   case '>': dwell *= 2; break;
1308   case '[': settle /= 2; if (settle < 1) settle = 1; break;
1309   case ']': settle *= 2; break;
1310   case 'd': go_down(); break;
1311   case 'D': FlushBuffer(); break;
1312   case 'e':
1313   case 'E': FlushBuffer();
1314       dorecalc = (!dorecalc);
1315       if (dorecalc)
1316       recalc();
1317       else {
1318       maxexp = minlyap; minexp = -1.0 * minlyap;
1319       }
1320       redraw(exponents[frame], expind[frame], 1);
1321       break;
1322   case 'f':
1323   case 'F': save_to_file(); break;
1324   case 'i': if (stripe_interval > 0) {
1325       stripe_interval--;
1326         if (displayplanes > 1) {
1327             init_color();
1328         }
1329       }
1330       break;
1331   case 'I': stripe_interval++;
1332       if (displayplanes > 1) {
1333         init_color();
1334       }
1335       break;
1336   case 'K': if (minlyap > 0.05)
1337       minlyap -= 0.05;
1338        break;
1339   case 'J': minlyap += 0.05;
1340        break;
1341   case 'm': mapindex++;
1342                   if (mapindex >= NUMMAPS)
1343                         mapindex=0;
1344                   map = Maps[mapindex];
1345                   deriv = Derivs[mapindex];
1346       if (!aflag)
1347                         min_a = amins[mapindex];
1348                   if (!wflag)
1349                         a_range = aranges[mapindex];
1350                   if (!bflag)
1351                         min_b = bmins[mapindex];
1352                   if (!hflag)
1353                         b_range = branges[mapindex];
1354                   if (!Force)
1355                         for (i=0;i<FUNCMAXINDEX;i++)
1356                              Forcing[i] = mapindex;
1357             max_a = min_a + a_range;
1358             max_b = min_b + b_range;
1359             a_minimums[0] = min_a; b_minimums[0] = min_b;
1360             a_maximums[0] = max_a; b_maximums[0] = max_b;
1361             a_inc = a_range / (double)width;
1362             b_inc = b_range / (double)height;
1363             point.x = -1;
1364             point.y = 0;
1365             a = rubber_data.p_min = min_a;
1366             b = rubber_data.q_min = min_b;
1367             rubber_data.p_max = max_a;
1368             rubber_data.q_max = max_b;
1369                   Clear();
1370                   break;
1371   case 'M': if (minlyap > 0.005)
1372       minlyap -= 0.005;
1373        break;
1374   case 'N': minlyap += 0.005;
1375        break;
1376   case 'p':
1377   case 'P': negative = (!negative);
1378       FlushBuffer(); redraw(exponents[frame], expind[frame], 1);
1379       break;
1380   case 'r': FlushBuffer(); redraw(exponents[frame], expind[frame], 1);
1381       break;
1382   case 'R': FlushBuffer(); Redraw(); break;
1383   case 's':
1384        spinlength=spinlength/2;
1385   case 'S': if (displayplanes > 1)
1386       Spin(canvas);
1387        spinlength=spinlength*2; break;
1388   case 'u': go_back(); break;
1389   case 'U': go_init(); break;
1390   case 'v':
1391   case 'V': print_values(); break;
1392   case 'W': if (numwheels < MAXWHEELS)
1393       numwheels++;
1394        else
1395       numwheels = 0;
1396        if (displayplanes > 1) {
1397         init_color();
1398        }
1399        break;
1400   case 'w': if (numwheels > 0)
1401       numwheels--;
1402        else
1403       numwheels = MAXWHEELS;
1404        if (displayplanes > 1) {
1405         init_color();
1406        }
1407        break;
1408   case 'x': Clear(); break;
1409   case 'X': Destroy_frame(); break;
1410   case 'z': Cycle_frames(); redraw(exponents[frame], expind[frame], 1);
1411       break;
1412   case 'Z': while (!XPending(dpy)) Cycle_frames();
1413       redraw(exponents[frame], expind[frame], 1); break;
1414   case 'q':
1415   case 'Q': exit(0); break;
1416   case '?':
1417   case 'h':
1418   case 'H': print_help(); break;
1419   default:  break;
1420   }
1421 }
1422
1423 /* Here's where we index into a color map. After the Lyapunov exponent is
1424  * calculated, it is used to determine what color to use for that point.
1425  * I suppose there are a lot of ways to do this. I used the following :
1426  * if it's non-negative then there's a reserved area at the lower range
1427  * of the color map that i index into. The ratio of some "minimum exponent
1428  * value" and the calculated value is used as a ratio of how high to index
1429  * into this reserved range. Usually these colors are dark red (see init_color).
1430  * If the exponent is negative, the same ratio (expo/minlyap) is used to index
1431  * into the remaining portion of the colormap (which is usually some light
1432  * shades of color or a rainbow wheel). The coloring scheme can actually make
1433  * a great deal of difference in the quality of the picture. Different colormaps
1434  * bring out different details of the dynamics while different indexing
1435  * algorithms also greatly effect what details are seen. Play around with this.
1436  */
1437 static int
1438 sendpoint(double expo)
1439 {
1440   static int index;
1441   static double tmpexpo;
1442
1443 #if 0
1444 /* The relationship minexp <= expo <= maxexp should always be true. This test
1445    enforces that. But maybe not enforcing it makes better pictures. */
1446   if (expo < minexp)
1447     expo = minexp;
1448   else if (expo > maxexp)
1449     expo = maxexp;
1450 #endif
1451
1452   point.x++;
1453   tmpexpo = (negative) ? expo : -1.0 * expo;
1454   if (tmpexpo > 0) {
1455     if (displayplanes >1) {
1456         index = (int)(tmpexpo*lowrange/maxexp);
1457         index = (index % lowrange) + startcolor;
1458     }
1459     else
1460         index = 0;
1461   }
1462   else {
1463     if (displayplanes >1) {
1464         index = (int)(tmpexpo*numfreecols/minexp);
1465         index = (index % numfreecols) + mincolindex;
1466     }
1467     else
1468         index = 1;
1469   }
1470     BufferPoint(dpy, canvas, index, point.x, point.y);
1471   if (save)
1472     exponents[frame][expind[frame]++] = expo;
1473   if (point.x >= width) {
1474     point.y++;
1475     point.x = 0;
1476     if (save) {
1477       b += b_inc;
1478       a = min_a;
1479     }
1480     if (point.y >= height)
1481       return FALSE;
1482     else
1483       return TRUE;
1484   }
1485   return TRUE;
1486 }
1487
1488 static void
1489 redisplay (Window w, XExposeEvent *event)
1490 {
1491   /*
1492   * Extract the exposed area from the event and copy
1493   * from the saved pixmap to the window.
1494   */
1495   XCopyArea(dpy, pixmap, canvas, Data_GC[0],
1496            event->x, event->y, event->width, event->height,
1497            event->x, event->y);
1498 }
1499
1500 static void
1501 resize(void)
1502 {
1503   Window r;
1504   int n, x, y;
1505   unsigned int bw, d, new_w, new_h;
1506
1507   XGetGeometry(dpy,canvas,&r,&x,&y,&new_w,&new_h,&bw,&d);
1508   if ((new_w == width) && (new_h == height))
1509     return;
1510   width = new_w; height = new_h;
1511   XClearWindow(dpy, canvas);
1512   if (pixmap)
1513     XFreePixmap(dpy, pixmap);
1514   pixmap = XCreatePixmap(dpy, DefaultRootWindow(dpy),
1515       width, height, DefaultDepth(dpy, screen));
1516   a_inc = a_range / (double)width;
1517   b_inc = b_range / (double)height;
1518   point.x = -1;
1519   point.y = 0;
1520   run = 1;
1521   a = rubber_data.p_min = min_a;
1522   b = rubber_data.q_min = min_b;
1523   rubber_data.p_max = max_a;
1524   rubber_data.q_max = max_b;
1525   freemem();
1526   setupmem();
1527         for (n=0;n<MAXFRAMES;n++)
1528     if ((n <= maxframe) && (n != frame))
1529         resized[n] = 1;
1530   InitBuffer();
1531   Clear();
1532   Redraw();
1533 }
1534
1535 static void
1536 redraw(double *exparray, int index, int cont)
1537 {
1538   static int i;
1539   static int x_sav, y_sav;
1540
1541   x_sav = point.x;
1542   y_sav = point.y;
1543
1544   point.x = -1;
1545   point.y = 0;
1546
1547   save=0;
1548   for (i=0;i<index;i++)
1549     sendpoint(exparray[i]);
1550   save=1;
1551
1552   if (cont) {
1553     point.x = x_sav;
1554     point.y = y_sav;
1555   }
1556   else {
1557     a = point.x * a_inc + min_a;
1558     b = point.y * b_inc + min_b;
1559   }
1560   FlushBuffer();
1561 }
1562
1563 static void
1564 Redraw(void)
1565 {
1566   FlushBuffer();
1567         point.x = -1;
1568         point.y = 0;
1569   run = 1;
1570         a = min_a;
1571         b = min_b;
1572   expind[frame] = 0;
1573   resized[frame] = 0;
1574 }
1575
1576 /* Store color pics in PPM format and monochrome in PGM */
1577 static void
1578 save_to_file(void)
1579 {
1580   FILE *outfile;
1581   unsigned char c;
1582   XImage *ximage;
1583   static int i,j;
1584   struct Colormap {
1585     unsigned char red;
1586     unsigned char green;
1587     unsigned char blue;
1588   };
1589   struct Colormap *colormap=NULL;
1590
1591   if (colormap)
1592     free(colormap);
1593   if ((colormap=
1594     (struct Colormap *)malloc(sizeof(struct Colormap)*maxcolor))
1595       == NULL) {
1596     fprintf(stderr,"Error malloc'ing colormap array\n");
1597     exit(-1);
1598   }
1599   outfile = fopen(outname,"w");
1600   if(!outfile) {
1601     perror(outname);
1602     exit(-1);
1603   }
1604
1605   ximage=XGetImage(dpy, pixmap, 0, 0, width, height, AllPlanes, XYPixmap);
1606
1607   if (displayplanes > 1) {
1608     for (i=0;i<maxcolor;i++) {
1609       colormap[i].red=(unsigned char)(Colors[i].red >> 8);
1610       colormap[i].green=(unsigned char)(Colors[i].green >> 8);
1611       colormap[i].blue =(unsigned char)(Colors[i].blue >> 8);
1612     }
1613     fprintf(outfile,"P%d %d %d\n",6,width,height);
1614   }
1615   else
1616     fprintf(outfile,"P%d %d %d\n",5,width,height);
1617   fprintf(outfile,"# settle=%d  dwell=%d start_x=%f\n",settle,dwell,
1618         start_x);
1619   fprintf(outfile,"# min_a=%f  a_rng=%f  max_a=%f\n",min_a,a_range,max_a);
1620   fprintf(outfile,"# min_b=%f  b_rng=%f  max_b=%f\n",min_b,b_range,max_b);
1621   if (Rflag)
1622     fprintf(outfile,"# pseudo-random forcing\n");
1623   else if (force) {
1624     fprintf(outfile,"# periodic forcing=");
1625     for (i=0;i<maxindex;i++) {
1626       fprintf(outfile,"%d",forcing[i]);
1627     }
1628     fprintf(outfile,"\n");
1629   }
1630   else
1631     fprintf(outfile,"# periodic forcing=01\n");
1632   if (Force) {
1633     fprintf(outfile,"# function forcing=");
1634     for (i=0;i<funcmaxindex;i++) {
1635       fprintf(outfile,"%d",Forcing[i]);
1636     }
1637     fprintf(outfile,"\n");
1638   }
1639   fprintf(outfile,"%d\n",numcolors-1);
1640
1641   for (j=0;j<height;j++)
1642       for (i=0;i<width;i++) {
1643     c = (unsigned char)XGetPixel(ximage,i,j);
1644     if (displayplanes > 1)
1645         fwrite((char *)&colormap[c],sizeof colormap[0],1,outfile);
1646     else
1647         fwrite((char *)&c,sizeof c,1,outfile);
1648       }
1649   fclose(outfile);
1650 }
1651
1652 static void
1653 recalc(void)
1654 {
1655   static int i, x, y;
1656
1657   minexp = maxexp = 0.0;
1658   x = y = 0;
1659   for (i=0;i<expind[frame];i++) {
1660     if (exponents[frame][i] < minexp)
1661       minexp = exponents[frame][i];
1662     if (exponents[frame][i] > maxexp)
1663       maxexp = exponents[frame][i];
1664   }
1665 }
1666
1667 static void
1668 Clear(void)
1669 {
1670       XClearWindow(dpy, canvas);
1671   XCopyArea(dpy, canvas, pixmap, Data_GC[0],
1672             0, 0, width, height, 0, 0);
1673   InitBuffer();
1674 }
1675
1676 static void
1677 show_defaults(void)
1678 {
1679
1680   printf("Width=%d  Height=%d  numcolors=%d  settle=%d  dwell=%d\n",
1681     width,height,numcolors,settle,dwell);
1682   printf("min_a=%f  a_range=%f  max_a=%f\n", min_a,a_range,max_a);
1683   printf("min_b=%f  b_range=%f  max_b=%f\n", min_b,b_range,max_b);
1684   printf("minlyap=%f  minexp=%f  maxexp=%f\n", minlyap,minexp,maxexp);
1685   exit(0);
1686 }
1687
1688 static void
1689 CreateXorGC(void)
1690 {
1691   XGCValues values;
1692
1693   values.foreground = foreground;
1694   values.line_style = LineSolid;
1695   values.function = GXxor;
1696   RubberGC = XCreateGC(dpy, DefaultRootWindow(dpy),
1697         GCForeground | GCBackground | GCFunction | GCLineStyle, &values);
1698 }
1699
1700 static void
1701 StartRubberBand(Window w, image_data_t *data, XEvent *event)
1702 {
1703   XPoint corners[5];
1704
1705   nostart = 0;
1706   data->rubber_band.last_x = data->rubber_band.start_x = event->xbutton.x;
1707   data->rubber_band.last_y = data->rubber_band.start_y = event->xbutton.y;
1708   SetupCorners(corners, data);
1709   XDrawLines(dpy, canvas, RubberGC,
1710       corners, sizeof(corners) / sizeof(corners[0]), CoordModeOrigin);
1711 }
1712
1713 static void
1714 SetupCorners(XPoint *corners, image_data_t *data)
1715 {
1716   corners[0].x = data->rubber_band.start_x;
1717   corners[0].y = data->rubber_band.start_y;
1718   corners[1].x = data->rubber_band.start_x;
1719   corners[1].y = data->rubber_band.last_y;
1720   corners[2].x = data->rubber_band.last_x;
1721   corners[2].y = data->rubber_band.last_y;
1722   corners[3].x = data->rubber_band.last_x;
1723   corners[3].y = data->rubber_band.start_y;
1724   corners[4] = corners[0];
1725 }
1726
1727 static void
1728 TrackRubberBand(Window w, image_data_t *data, XEvent *event)
1729 {
1730   XPoint corners[5];
1731   int xdiff, ydiff;
1732
1733   if (nostart)
1734     return;
1735   SetupCorners(corners, data);
1736   XDrawLines(dpy, canvas, RubberGC,
1737       corners, sizeof(corners) / sizeof(corners[0]), CoordModeOrigin);
1738   ydiff = event->xbutton.y - data->rubber_band.start_y;
1739   xdiff = event->xbutton.x - data->rubber_band.start_x;
1740   data->rubber_band.last_x = data->rubber_band.start_x + xdiff;
1741   data->rubber_band.last_y = data->rubber_band.start_y + ydiff;
1742   if (data->rubber_band.last_y < data->rubber_band.start_y ||
1743       data->rubber_band.last_x < data->rubber_band.start_x)
1744   {
1745     data->rubber_band.last_y = data->rubber_band.start_y;
1746     data->rubber_band.last_x = data->rubber_band.start_x;
1747   }
1748   SetupCorners(corners, data);
1749   XDrawLines(dpy, canvas, RubberGC,
1750       corners, sizeof(corners) / sizeof(corners[0]), CoordModeOrigin);
1751 }
1752
1753 static void
1754 EndRubberBand(Window w, image_data_t *data, XEvent *event)
1755 {
1756   XPoint corners[5];
1757   XPoint top, bot;
1758   double delta, diff;
1759
1760   nostart = 1;
1761   SetupCorners(corners, data);
1762   XDrawLines(dpy, canvas, RubberGC,
1763       corners, sizeof(corners) / sizeof(corners[0]), CoordModeOrigin);
1764   if (data->rubber_band.start_x >= data->rubber_band.last_x ||
1765       data->rubber_band.start_y >= data->rubber_band.last_y)
1766     return;
1767   top.x = data->rubber_band.start_x;
1768   bot.x = data->rubber_band.last_x;
1769   top.y = data->rubber_band.start_y;
1770   bot.y = data->rubber_band.last_y;
1771   diff = data->q_max - data->q_min;
1772   delta = (double)top.y / (double)height;
1773   data->q_min += diff * delta;
1774   delta = (double)(height - bot.y) / (double)height;
1775   data->q_max -= diff * delta;
1776   diff = data->p_max - data->p_min;
1777   delta = (double)top.x / (double)width;
1778   data->p_min += diff * delta;
1779   delta = (double)(width - bot.x) / (double)width;
1780   data->p_max -= diff * delta;
1781   fflush(stdout);
1782   set_new_params(w, data);
1783 }
1784
1785 static void
1786 set_new_params(Window w, image_data_t *data)
1787 {
1788   frame = (maxframe + 1) % MAXFRAMES;
1789   if (frame > maxframe)
1790     maxframe = frame;
1791   a_range = data->p_max - data->p_min;
1792   b_range = data->q_max - data->q_min;
1793         a_minimums[frame] = min_a = data->p_min;
1794         b_minimums[frame] = min_b = data->q_min;
1795         a_inc = a_range / (double)width;
1796         b_inc = b_range / (double)height;
1797         point.x = -1;
1798         point.y = 0;
1799   run = 1;
1800         a = min_a;
1801         b = min_b;
1802         a_maximums[frame] = max_a = data->p_max;
1803         b_maximums[frame] = max_b = data->q_max;
1804   expind[frame] = 0;;
1805   Clear();
1806 }
1807
1808 static void
1809 go_down(void)
1810 {
1811   frame++;
1812   if (frame > maxframe)
1813     frame = 0;
1814   jumpwin();
1815 }
1816
1817 static void
1818 go_back(void)
1819 {
1820   frame--;
1821   if (frame < 0)
1822     frame = maxframe;
1823   jumpwin();
1824 }
1825
1826 static void
1827 jumpwin(void)
1828 {
1829   rubber_data.p_min = min_a = a_minimums[frame];
1830   rubber_data.q_min = min_b = b_minimums[frame];
1831   rubber_data.p_max = max_a = a_maximums[frame];
1832   rubber_data.q_max = max_b = b_maximums[frame];
1833   a_range = max_a - min_a;
1834   b_range = max_b - min_b;
1835         a_inc = a_range / (double)width;
1836         b_inc = b_range / (double)height;
1837         point.x = -1;
1838         point.y = 0;
1839         a = min_a;
1840         b = min_b;
1841   Clear();
1842   if (resized[frame])
1843     Redraw();
1844   else
1845     redraw(exponents[frame], expind[frame], 0);
1846 }
1847
1848 static void
1849 go_init(void)
1850 {
1851   frame = 0;
1852   jumpwin();
1853 }
1854
1855 static void
1856 Destroy_frame(void)
1857 {
1858   static int i;
1859
1860   for (i=frame; i<maxframe; i++) {
1861     exponents[frame] = exponents[frame+1];
1862     expind[frame] = expind[frame+1];
1863     a_minimums[frame] = a_minimums[frame+1];
1864     b_minimums[frame] = b_minimums[frame+1];
1865     a_maximums[frame] = a_maximums[frame+1];
1866     b_maximums[frame] = b_maximums[frame+1];
1867   }
1868   maxframe--;
1869   go_back();
1870 }
1871
1872 static void
1873 InitBuffer(void)
1874 {
1875   int i;
1876
1877   for (i = 0 ; i < maxcolor; ++i)
1878     Points.npoints[i] = 0;
1879 }
1880
1881 static void
1882 BufferPoint(Display *display, Window window, int color, int x, int y)
1883 {
1884
1885 /* Guard against bogus color values. Shouldn't be necessary but paranoia
1886    is good. */
1887   if (color < 0)
1888     color = 0;
1889   else if (color >= maxcolor)
1890     color = maxcolor - 1;
1891
1892   if (Points.npoints[color] == MAXPOINTS)
1893   {
1894     XDrawPoints(display, window, Data_GC[color],
1895         Points.data[color], Points.npoints[color], CoordModeOrigin);
1896     XDrawPoints(display, pixmap, Data_GC[color],
1897         Points.data[color], Points.npoints[color], CoordModeOrigin);
1898     Points.npoints[color] = 0;
1899   }
1900   Points.data[color][Points.npoints[color]].x = x;
1901   Points.data[color][Points.npoints[color]].y = y;
1902   ++Points.npoints[color];
1903 }
1904
1905 static void
1906 FlushBuffer(void)
1907 {
1908   int color;
1909
1910   for (color = 0; color < maxcolor; ++color)
1911     if (Points.npoints[color])
1912     {
1913         XDrawPoints(dpy, canvas, Data_GC[color],
1914           Points.data[color], Points.npoints[color],
1915           CoordModeOrigin);
1916         XDrawPoints(dpy, pixmap, Data_GC[color],
1917           Points.data[color], Points.npoints[color],
1918           CoordModeOrigin);
1919         Points.npoints[color] = 0;
1920     }
1921 }
1922
1923 static void
1924 print_help(void)
1925 {
1926     printf("During run-time, interactive control can be exerted via : \n");
1927     printf("Mouse buttons allow rubber-banding of a zoom box\n");
1928     printf("< halves the 'dwell', > doubles the 'dwell'\n");
1929     printf("[ halves the 'settle', ] doubles the 'settle'\n");
1930     printf("D flushes the drawing buffer\n");
1931     printf("e or E recalculates color indices\n");
1932     printf("f or F saves exponents to a file\n");
1933     printf("h or H or ? displays this message\n");
1934     printf("i decrements, I increments the stripe interval\n");
1935     printf("KJMN increase/decrease minimum negative exponent\n");
1936     printf("m increments the map index, changing maps\n");
1937     printf("p or P reverses the colormap for negative/positive exponents\n");
1938     printf("r redraws without recalculating\n");
1939     printf("R redraws, recalculating with new dwell and settle values\n");
1940     printf("s or S spins the colorwheel\n");
1941     printf("u pops back up to the last zoom\n");
1942     printf("U pops back up to the first picture\n");
1943     printf("v or V displays the values of various settings\n");
1944     printf("w decrements, W increments the color wheel index\n");
1945     printf("x or X clears the window\n");
1946     printf("q or Q exits\n");
1947 }
1948
1949 static void
1950 print_values(void)
1951 {
1952     static int i;
1953
1954     printf("\nminlyap=%f minexp=%f maxexp=%f\n",minlyap,minexp,maxexp);
1955     printf("width=%d height=%d\n",width,height);
1956     printf("settle=%d  dwell=%d start_x=%f\n",settle,dwell, start_x);
1957     printf("min_a=%f  a_rng=%f  max_a=%f\n",min_a,a_range,max_a);
1958     printf("min_b=%f  b_rng=%f  max_b=%f\n",min_b,b_range,max_b);
1959     if (Rflag)
1960   printf("pseudo-random forcing\n");
1961     else if (force) {
1962   printf("periodic forcing=");
1963   for (i=0;i<maxindex;i++)
1964     printf("%d",forcing[i]);
1965   printf("\n");
1966     }
1967     else
1968   printf("periodic forcing=01\n");
1969     if (Force) {
1970   printf("function forcing=");
1971   for (i=0;i<funcmaxindex;i++) {
1972     printf("%d",Forcing[i]);
1973   }
1974   printf("\n");
1975     }
1976     printf("numcolors=%d\n",numcolors-1);
1977 }
1978
1979 static void
1980 freemem(void)
1981 {
1982   static int i;
1983
1984         for (i=0;i<MAXFRAMES;i++)
1985                 free(exponents[i]);
1986 }
1987
1988 static void
1989 setupmem(void)
1990 {
1991   static int i;
1992
1993         for (i=0;i<MAXFRAMES;i++) {
1994                 if((exponents[i]=
1995                     (double *)malloc(sizeof(double)*width*height))==NULL){
1996                     fprintf(stderr,"Error malloc'ing exponent array.\n");
1997                     exit(-1);
1998                 }
1999         }
2000 }
2001
2002 static void
2003 setforcing(void)
2004 {
2005   static int i;
2006   for (i=0;i<MAXINDEX;i++)
2007     forcing[i] = (random() > prob) ? 0 : 1;
2008 }