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