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