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