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