source: trunk/third/xscreensaver/hacks/xlyap.c @ 20148

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