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