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