home *** CD-ROM | disk | FTP | other *** search
/ Skunkware 5 / Skunkware 5.iso / src / X11 / lyap / lyap.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-05-09  |  52.6 KB  |  1,671 lines
  1. /*
  2.  *    @(#) lyap.c 12.1 95/05/09 SCOINC
  3.  */
  4. /*************************************************************************
  5.  *                                                                       *
  6.  *  Copyright (c) 1992, 1993 Ronald Joe Record                           *
  7.  *                                                                       *
  8.  *  All rights reserved. No part of this program or publication may be   *
  9.  *  reproduced, transmitted, transcribed, stored in a retrieval system,  *
  10.  *  or translated into any language or computer language, in any form or *
  11.  *  by any means, electronic, mechanical, magnetic, optical, chemical,   *
  12.  *  biological, or otherwise, without the prior written permission of:   *
  13.  *                                                                       *
  14.  *      Ronald Joe Record (408) 458-3718                                 *
  15.  *      212 Owen St., Santa Cruz, California 95062 USA                   *
  16.  *                                                                       *
  17.  *************************************************************************/
  18.  
  19. /* Lyap - calculate and display Lyapunov exponents */
  20.  
  21. /* Written by Ronald Record (rr@sco.com) 03 Sep 1991 */
  22.  
  23. /* The idea here is to calculate the Lyapunov exponent for a periodically
  24.  * forced logistic map (later i added several other nonlinear maps of the unit
  25.  * interval). In order to turn the 1-dimensional parameter space of the
  26.  * logistic map into a 2-dimensional parameter space, select two parameter
  27.  * values ('a' and 'b') then alternate the iterations of the logistic map using
  28.  * first 'a' then 'b' as the parameter. This program accepts an argument to 
  29.  * specify a forcing function, so instead of just alternating 'a' and 'b', you
  30.  * can use 'a' as the parameter for say 6 iterations, then 'b' for 6 iterations
  31.  * and so on. An interesting forcing function to look at is abbabaab (the
  32.  * Morse-Thue sequence, an aperiodic self-similar, self-generating sequence).
  33.  * Anyway, step through all the values of 'a' and 'b' in the ranges you want,
  34.  * calculating the Lyapunov exponent for each pair of values. The exponent
  35.  * is calculated by iterating out a ways (specified by the variable "settle")
  36.  * then on subsequent iterations calculating an average of the logarithm of
  37.  * the absolute value of the derivative at that point. Points in parameter
  38.  * space with a negative Lyapunov exponent are colored one way (using the
  39.  * value of the exponent to index into a color map) while points with a
  40.  * non-negative exponent are colored differently. 
  41.  * 
  42.  * The algorithm was taken from the September 1991 Scientific American article
  43.  * by A. K. Dewdney who gives credit to Mario Markus of the Max Planck Institute
  44.  * for its creation. Additional information and ideas were gleaned from the
  45.  * discussion on alt.fractals involving Stephen Hall, Ed Kubaitis, Dave Platt
  46.  * and Baback Moghaddam. Assistance with colormaps and spinning color wheels
  47.  * and X was gleaned from Hiram Clawson. Rubber banding code was adapted from
  48.  * an existing Mandelbrot program written by Stacey Campbell.
  49.  */
  50.  
  51. #include "lyap.h"
  52.  
  53. static char *lyap_c_id = "@(#) lyap.c 12.1 95/05/09 SCOINC";
  54. static char *version = LYAP_VERSION;
  55.  
  56. extern void BufferPoint(), InitBuffer(), FlushBuffer();
  57.  
  58. void
  59. Cleanup() {
  60.     XCloseDisplay(dpy);
  61. }
  62.  
  63. main(ac, av)
  64.     int ac;
  65.     char **av;
  66. {
  67.     int i, j; 
  68.     int x_center=0, y_center=0;
  69.     XSizeHints hint;
  70.     extern void init_canvas(), init_data(), init_color(), parseargs();
  71.     extern void Clear(), restor_picture(); 
  72.     extern void center_horizontal(), center_vertical();
  73.  
  74.     parseargs(ac, av);
  75.     dpy = XOpenDisplay("");
  76.     screen = DefaultScreen(dpy);
  77.     background = BlackPixel(dpy, screen);
  78.     if (width < 1)
  79.         width = XDisplayWidth(dpy, screen);
  80.     if (height < 1)
  81.         height = XDisplayHeight(dpy, screen);
  82.     setupmem();
  83.     init_data();
  84.     if (displayplanes > 1)
  85.         foreground = startcolor;
  86.     else
  87.         foreground = WhitePixel(dpy,XDefaultScreen(dpy));
  88.     if (xposition < 0) {
  89.         x_center = 1;
  90.         xposition = Max(0,(XDisplayWidth(dpy, screen) - width)/2);
  91.     }
  92.     if (yposition < 0) {
  93.         y_center = 1;
  94.         yposition = Max(0,(XDisplayHeight(dpy, screen) - height)/2);
  95.     }
  96.     hint.x = xposition;
  97.     hint.y = yposition;
  98.     hint.width = width;
  99.     hint.height = height;
  100.     hint.flags = PPosition | PSize;
  101.     /*
  102.      * Create the window to display the Lyapunov exponents 
  103.      */
  104.     canvas = XCreateSimpleWindow(dpy, DefaultRootWindow(dpy),
  105.         hint.x, hint.y, hint.width, hint.height,
  106.         5, foreground, background);
  107.     XSetStandardProperties(dpy, canvas, "Lyap by Ron Record", 
  108.         "Lyap", None, av, ac, &hint);
  109.     XChangeProperty( dpy, canvas, XA_WM_CLASS, XA_STRING, 8, PropModeReplace, 
  110.                     "lyap", strlen("lyap"));
  111.     hint.x = 0;
  112.     hint.y = XDisplayHeight(dpy, screen) / 2;
  113.     hint.width = XDisplayWidth(dpy, screen) / 3;
  114.     hint.height = hint.y / 2;
  115.     info = XCreateSimpleWindow(dpy, DefaultRootWindow(dpy),
  116.                 hint.x, hint.y, hint.width, hint.height, 5, 0, 1);
  117.     XSetWindowBackground(dpy, info, BlackPixel(dpy, screen));
  118.     /* Title */
  119.     XSetStandardProperties(dpy,info,"Info","Info",None,av,ac,&hint);
  120.     hint.x = XDisplayWidth(dpy, screen) / 2;
  121.     hint.y = 0;
  122.     hint.width = hint.x;
  123.     hint.height = 4 * XDisplayHeight(dpy, screen) / 6;
  124.     help = XCreateSimpleWindow(dpy, DefaultRootWindow(dpy),
  125.                 hint.x, hint.y, hint.width, hint.height, 5, 0, 1);
  126.     XSetWindowBackground(dpy, help, BlackPixel(dpy, screen));
  127.     /* Title */
  128.     XSetStandardProperties(dpy,help,"Help","Help",None,av,ac,&hint);
  129.     init_canvas();
  130.     XMapRaised(dpy, canvas);
  131.     if (x_center)
  132.         center_horizontal(dpy, canvas, width);
  133.     if (y_center)
  134.         center_vertical(dpy, canvas, height);
  135.     /* Try to write into a new color map */
  136.     cmap = XCreateColormap(dpy, canvas, DefaultVisual(dpy, screen), AllocAll);
  137.     if (displayplanes > 1)
  138.         init_color(dpy, canvas, cmap, Colors, startcolor, mincolindex, 
  139.                     numcolors, numwheels, "lyap", "Lyap", 0);
  140.     else
  141.         XQueryColors(dpy, DefaultColormap(dpy, DefaultScreen(dpy)), 
  142.                 Colors, numcolors);
  143.     /* install new color map */
  144.     XSetWindowColormap(dpy, canvas, cmap);
  145.     XSetWindowColormap(dpy, help, cmap);
  146.     XSetWindowColormap(dpy, info, cmap);
  147.     pixmap = XCreatePixmap(dpy, DefaultRootWindow(dpy), width, height, 
  148.                            DefaultDepth(dpy, screen));
  149.     rubber_data.band_cursor=XCreateFontCursor(dpy, XC_hand2);
  150.     CreateXorGC();
  151.     Clear();
  152.     if (restfile)
  153.         restor_picture();
  154.     if (demo)
  155.         for(;;) {
  156.             if (!run) {
  157.                 DemoSpin(dpy, cmap, Colors, startcolor, maxcolor, delay, 2);
  158.                 for (i=0; i<=MAXWHEELS; i++) {
  159.                     init_color(dpy, canvas, cmap, Colors, startcolor, 
  160.                             mincolindex, numcolors, i, "lyap", "Lyap", 0);
  161.                     sleep(1);
  162.                 }
  163.                 if (!run) {
  164.                     Cleanup();
  165.                     exit(0);
  166.                 } 
  167.             }
  168.         }
  169.     else {
  170.         XSelectInput(dpy,help,KeyPressMask|ExposureMask);
  171.         XSelectInput(dpy,info,KeyPressMask|ExposureMask);
  172.         XSelectInput(dpy,canvas,KeyPressMask|ButtonPressMask|ButtonMotionMask|
  173.                 ButtonReleaseMask|ExposureMask|StructureNotifyMask);
  174.         for(;;)
  175.             main_event();
  176.     }
  177. }
  178.  
  179. main_event()
  180. {
  181.     int n;
  182.     XEvent event;
  183.  
  184.     if (complyap() == TRUE)
  185.         run=0;
  186.     n = XEventsQueued(dpy, QueuedAfterFlush);
  187.     while (n--) {
  188.         XNextEvent(dpy, &event);
  189.         switch(event.type) {
  190.             case KeyPress: Getkey(&event); break;
  191.             case Expose: redisplay(canvas, &event); break;
  192.             case ConfigureNotify: resize(); break;
  193.             case ButtonPress: StartRubberBand(canvas, &rubber_data, &event);
  194.                  break;
  195.             case MotionNotify: TrackRubberBand(canvas, &rubber_data, &event);
  196.                  break;
  197.             case ButtonRelease: EndRubberBand(canvas, &rubber_data, &event);
  198.                  break;
  199.         }
  200.     }
  201. }
  202.  
  203. /* complyap() is the guts of the program. This is where the Lyapunov exponent
  204.  * is calculated. For each iteration (past some large number of iterations)
  205.  * calculate the logarithm of the absolute value of the derivative at that
  206.  * point. Then average them over some large number of iterations. Some small
  207.  * speed up is achieved by utilizing the fact that log(a*b) = log(a) + log(b).
  208.  */
  209. complyap()
  210. {
  211.     register i, bindex, findex;
  212.     double total, prod, x, r;
  213.     extern void save_to_file(), store_to_file();
  214.  
  215.     if (!run)
  216.         return TRUE;
  217.     if (a >= max_a)
  218.         if (sendpoint(lyapunov) == TRUE)
  219.             return FALSE;
  220.         else {
  221.             run=0;
  222.             FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points, 0, maxcolor);
  223.             if (savefile)
  224.                 save_to_file();
  225.             else if (storefile)
  226.                 store_to_file();
  227.             return TRUE;
  228.         }
  229.     if (b >= max_b) {
  230.         run=0;
  231.         FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points, 0, maxcolor);
  232.         if (savefile)
  233.             save_to_file();
  234.         else if (storefile)
  235.             store_to_file();
  236.         return TRUE;
  237.     }
  238.     prod = 1.0;
  239.     total = 0.0;
  240.     bindex = 0;
  241.     x = start_x;
  242.     r = (forcing[bindex]) ? b : a;
  243. #ifdef MAPS
  244.     findex = 0;
  245.     map = Maps[Forcing[findex]];
  246. #endif
  247.     for (i=0;i<settle;i++) {       /* Here's where we let the thing */
  248.         x = (*map)(x, r);       /* "settle down". There is usually */
  249.         if (++bindex >= maxindex) { /* some initial "noise" in the */
  250.             bindex = 0;       /* iterations. How can we optimize */
  251.             if (Rflag)        /* the value of settle ??? */
  252.                 setforcing();
  253.         }
  254.         r = (forcing[bindex]) ? b : a;
  255. #ifdef MAPS
  256.         if (++findex >= funcmaxindex)
  257.             findex = 0;
  258.         map = Maps[Forcing[findex]];
  259. #endif
  260.     }
  261. #ifdef MAPS
  262.     deriv = Derivs[Forcing[findex]];
  263. #endif
  264.     if (useprod) {            /* using log(a*b) */
  265.         for (i=0;i<dwell;i++) {
  266.             x = (*map)(x, r);
  267.             prod *= ABS((*deriv)(x, r));
  268.             /* we need to prevent overflow and underflow */
  269.             if ((prod > 1.0e12) || (prod < 1.0e-12)) {
  270.                 if (prod == 0.0)
  271.                     total += LN_MINDOUBLE;
  272.                 else if (prod < 0.0)
  273.                     total += LN_MAXDOUBLE;
  274.                 else
  275.                     total += log(prod);
  276.                 prod = 1.0;
  277.             }
  278.             if (++bindex >= maxindex) {
  279.                 bindex = 0;
  280.                 if (Rflag)
  281.                     setforcing();
  282.             }
  283.             r = (forcing[bindex]) ? b : a;
  284. #ifdef MAPS
  285.             if (++findex >= funcmaxindex)
  286.                 findex = 0;
  287.             map = Maps[Forcing[findex]];
  288.             deriv = Derivs[Forcing[findex]];
  289. #endif
  290.         }
  291.         if (prod == 0.0)
  292.             total += LN_MINDOUBLE;
  293.         else if (prod < 0.0)
  294.             total += LN_MAXDOUBLE;
  295.         else
  296.             total += log(prod);
  297.         lyapunov = (total * M_LOG2E) / (double)dwell;
  298.     }
  299.     else {                /* use log(a) + log(b) */
  300.         for (i=0;i<dwell;i++) {
  301.             x = (*map)(x, r);
  302.             total += log(ABS((*deriv)(x, r)));
  303.             if (++bindex >= maxindex) {
  304.                 bindex = 0;
  305.                 if (Rflag)
  306.                     setforcing();
  307.             }
  308.             r = (forcing[bindex]) ? b : a;
  309. #ifdef MAPS
  310.             if (++findex >= funcmaxindex)
  311.                 findex = 0;
  312.             map = Maps[Forcing[findex]];
  313.             deriv = Derivs[Forcing[findex]];
  314. #endif
  315.         }
  316.         lyapunov = (total * M_LOG2E) / (double)dwell;
  317.     }
  318.     if (sendpoint(lyapunov) == TRUE)
  319.         return FALSE;
  320.     else {
  321.         run=0;
  322.         FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points, 0, maxcolor);
  323.         if (savefile)
  324.             save_to_file();
  325.         else if (storefile)
  326.             store_to_file();
  327.         return TRUE;
  328.     }
  329. }
  330.  
  331. double
  332. logistic(x, r)            /* the familiar logistic map */
  333. double x, r;
  334. {
  335.     return(r * x * (1.0 - x));
  336. }
  337.  
  338. double
  339. dlogistic(x, r)            /* the derivative of logistic map */
  340. double x, r;
  341. {
  342.     return(r - (2.0 * r * x));
  343. }
  344.  
  345. double
  346. circle(x, r)            /* sin() hump or sorta like the circle map */
  347. double x, r;
  348. {
  349.     extern double sin();
  350.  
  351.     return(r * sin(M_PI * x));
  352. }
  353.  
  354. double
  355. dcircle(x, r)            /* derivative of the "sin() hump" */
  356. double x, r;
  357. {
  358.     extern double cos();
  359.  
  360.     return(r * M_PI * cos(M_PI * x));
  361. }
  362.  
  363. double
  364. leftlog(x, r)            /* left skewed logistic */
  365. double x, r;
  366. {
  367.     double d;
  368.  
  369.     d = 1.0 - x;
  370.     return(r * x * d * d);
  371. }
  372.  
  373. double
  374. dleftlog(x, r)            /* derivative of the left skewed logistic */
  375. double x, r;
  376. {
  377.     return(r * (1.0 - (4.0 * x) + (3.0 * x * x)));
  378. }
  379.  
  380. double
  381. rightlog(x, r)            /* right skewed logistic */
  382. double x, r;
  383. {
  384.     return(r * x * x * (1.0 - x));
  385. }
  386.  
  387. double
  388. drightlog(x, r)            /* derivative of the right skewed logistic */
  389. double x, r;
  390. {
  391.     return(r * ((2.0 * x) - (3.0 * x * x)));
  392. }
  393.  
  394. double
  395. doublelog(x, r)            /* double logistic */
  396. double x, r;
  397. {
  398.     double d;
  399.  
  400.     d = 1.0 - x;
  401.     return(r * x * x * d * d);
  402. }
  403.  
  404. double
  405. ddoublelog(x, r)        /* derivative of the double logistic */
  406. double x, r;
  407. {
  408.     double d;
  409.  
  410.     d = x * x;
  411.     return(r * ((2.0 * x) - (6.0 * d) + (4.0 * x * d)));
  412. }
  413.  
  414. void
  415. init_data()
  416. {
  417.     static int i;
  418.  
  419.     numcolors = XDisplayCells(dpy, XDefaultScreen(dpy));
  420.     displayplanes = DisplayPlanes(dpy, XDefaultScreen(dpy));
  421.     if (numcolors > maxcolor)
  422.         numcolors = maxcolor;
  423.     if (numcolors <= 16) {
  424.         maxcolor = 16; startcolor = 0; 
  425.         color_offset = 0; mincolindex = 1;
  426.     }
  427.     numfreecols = numcolors - mincolindex;
  428.     lowrange = mincolindex - startcolor;
  429.     a_inc = a_range / (double)width;
  430.     b_inc = b_range / (double)height;
  431.     if (!restfile) {
  432.         point.x = point.y = 0;
  433.         a = min_a; b = min_b;
  434.     }
  435.     rubber_data.p_min = min_a;
  436.     rubber_data.q_min = min_b;
  437.     rubber_data.p_max = max_a;
  438.     rubber_data.q_max = max_b;
  439.     if (show)
  440.         show_defaults();
  441.     InitBuffer(&Points, maxcolor);
  442.     srand48(time(0));
  443. }
  444.  
  445. void
  446. init_canvas()
  447. {
  448.     static int i;
  449.  
  450.     /*
  451.      * create default, writable, graphics contexts for the canvas.
  452.      */
  453.     for (i=0; i<maxcolor; i++) {
  454.         Data_GC[i] = XCreateGC(dpy, DefaultRootWindow(dpy),
  455.                 (unsigned long) NULL, (XGCValues *) NULL);
  456.         /* set the background to black */
  457.         XSetBackground(dpy,Data_GC[i],BlackPixel(dpy,XDefaultScreen(dpy)));
  458.         /* set the foreground of the ith context to i */
  459.         XSetForeground(dpy, Data_GC[i], i);
  460.     }
  461.     XSetForeground(dpy,Data_GC[0],BlackPixel(dpy,XDefaultScreen(dpy)));
  462.     XSetForeground(dpy,Data_GC[1],WhitePixel(dpy,XDefaultScreen(dpy)));
  463. }
  464.  
  465. void
  466. parseargs(ac, av)
  467. int ac;
  468. char **av;
  469. {
  470.     static int c;
  471.     static int i;
  472.     int bindex=0, findex;
  473.     char *ch;
  474.     extern int optind;
  475.     extern char *optarg;
  476.     extern double atof();
  477.     extern void check_params(), usage(), restor_params();
  478.  
  479.     map = Maps[0];
  480.     deriv = Derivs[0];
  481.     maxexp=minlyap; minexp= -1.0 * minlyap;
  482.  
  483.     while ((c=getopt(ac,av,
  484.         "Ldpuvc:i:m:C:W:H:I:M:N:O:R:S:a:b:D:F:f:o:w:h:s:x:y:"))!=EOF){
  485.         switch (c) {
  486.         case 'C':    mincolindex=atoi(optarg); break;
  487.         case 'D':    if (strcmp(optarg, "elay"))
  488.                         dwell=atoi(optarg); 
  489.                      else
  490.                         delay = atoi(av[optind++]);
  491.                      break;
  492. #ifdef MAPS
  493.         case 'F':    funcmaxindex = strlen(optarg);
  494.                 if (funcmaxindex > FUNCMAXINDEX)
  495.                     usage();
  496.                 ch = optarg;
  497.                 Force++;
  498.                 for (findex=0;findex<funcmaxindex;findex++) {
  499.                     Forcing[findex] = (int)(*ch++ - '0');
  500.                     if (Forcing[findex] >= NUMMAPS)
  501.                         usage();
  502.                 }
  503.                 break;
  504. #endif
  505.         case 'H':    height=atoi(optarg); break;
  506.         case 'I':    start_x=atof(optarg); break;
  507.         case 'L':    useprod=0; break;
  508.         case 'M':    minlyap=ABS(atof(optarg)); 
  509.                 maxexp=minlyap; minexp= -1.0 * minlyap; break;
  510.         case 'N':    maxcolor=ABS(atoi(optarg)); 
  511.                 if ((maxcolor - startcolor) <= 0)
  512.                     startcolor = 0;
  513.                 if ((maxcolor - mincolindex) <= 0) {
  514.                     mincolindex = 1;
  515.                     color_offset = 0;
  516.                 }
  517.                 break;
  518.         case 'O':    color_offset=atoi(optarg); break;
  519.         case 'R':    prob=atof(optarg); Rflag++; setforcing(); break;
  520.         case 'S':    settle=atoi(optarg); break;
  521.         case 'W':    width=atoi(optarg); break;
  522.         case 'a':    min_a=atof(optarg); aflag++; break;
  523.         case 'b':    min_b=atof(optarg); bflag++; break;
  524.         case 'c':    numwheels=atoi(optarg); break;
  525.         case 'd':    demo++; break;
  526.         case 'f':    maxindex = strlen(optarg);
  527.                 if (maxindex > MAXINDEX)
  528.                     usage();
  529.                 ch = optarg;
  530.                 force++;
  531.                 while (bindex < maxindex) {
  532.                     if (*ch == 'a')
  533.                         forcing[bindex++] = 0;
  534.                     else if (*ch == 'b')
  535.                         forcing[bindex++] = 1;
  536.                     else
  537.                         usage();
  538.                     ch++;
  539.                 }
  540.                 break;
  541.         case 'h':    b_range=atof(optarg); hflag++; break;
  542.         case 'i':    restfile++; inname=optarg; break;
  543.         case 'm':    mapindex=atoi(optarg); 
  544.                 if ((mapindex >= NUMMAPS) || (mapindex < 0))
  545.                     usage();
  546.                 map = Maps[mapindex];
  547.                 deriv = Derivs[mapindex];
  548.                 if (!aflag)
  549.                     min_a = amins[mapindex];
  550.                 if (!wflag)
  551.                     a_range = aranges[mapindex];
  552.                 if (!bflag)
  553.                     min_b = bmins[mapindex];
  554.                 if (!hflag)
  555.                     b_range = branges[mapindex];
  556.                 if (!Force)
  557.                     for (i=0;i<FUNCMAXINDEX;i++)
  558.                         Forcing[i] = mapindex;
  559.                 break;
  560.         case 'o':    savefile++; outname=optarg; break;
  561.         case 'p':    negative--; break;
  562.         case 's':    storefile++; savname=optarg; break;
  563.         case 'u':    usage(); break;
  564.         case 'v':    show=1; break;
  565.         case 'w':    a_range=atof(optarg); wflag++; break;
  566.         case 'x':    xposition=atoi(optarg); break;
  567.         case 'y':    yposition=atoi(optarg); break;
  568.         case '?':    usage(); break;
  569.         }
  570.     }
  571.     if (restfile)
  572.         restor_params();
  573.     else {
  574.         max_a = min_a + a_range;
  575.         max_b = min_b + b_range;
  576.     }
  577.     a_minimums[0] = min_a; b_minimums[0] = min_b;
  578.     a_maximums[0] = max_a; b_maximums[0] = max_b;
  579.     if (Force)
  580.         if (maxindex == funcmaxindex)
  581.             for (findex=0;findex<funcmaxindex;findex++)
  582.                 check_params(Forcing[findex],forcing[findex]);
  583.         else
  584.             fprintf(stderr, "Warning! Unable to check parameters\n");
  585.     else
  586.         check_params(mapindex,2);
  587. }
  588.  
  589. void
  590. print_help() 
  591. {
  592.     static char str[80];
  593.     static int x_str, y_str, spacing;
  594.     static int ascent, descent, dir;
  595.     static XCharStruct overall;
  596.     static GC gc;
  597.  
  598.     gc = Data_GC[1];
  599.     XClearWindow(dpy, help);
  600.     x_str = 10;
  601.     y_str = 20;
  602.     sprintf(str,"During run-time, interactive control can be exerted via : ");
  603.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  604.     XQueryTextExtents(dpy,(XID)XGContextFromGC(gc),"Hey!",
  605.             4,&dir,&ascent,&descent,&overall);
  606.     spacing = ascent + descent + 5;
  607.     y_str += 2 * spacing;
  608.     sprintf(str,"Mouse buttons allow rubber-banding of a zoom box");
  609.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  610.     y_str += spacing;
  611.     sprintf(str,"< halves the 'dwell', > doubles the 'dwell'");
  612.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  613.     y_str += spacing;
  614.     sprintf(str,"[ halves the 'settle', ] doubles the 'settle'");
  615.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  616.     y_str += spacing;
  617.     sprintf(str,"D flushes the drawing buffer");
  618.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  619.     y_str += spacing;
  620.     sprintf(str,"d descends a ladder of windows created via the u or U keys");
  621.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  622.     y_str += spacing;
  623.     sprintf(str,"e or E recalculates color indices");
  624.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  625.     y_str += spacing;
  626.     sprintf(str,"f saves program state to a file");
  627.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  628.     y_str += spacing;
  629.     sprintf(str,"F saves window to a PPM file");
  630.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  631.     y_str += spacing;
  632.     sprintf(str,"h or H or ? displays this message");
  633.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  634.     sprintf(str,"i decrements, I increments the stripe interval");
  635.     y_str += spacing;
  636.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  637.     sprintf(str,"KJMN increase/decrease minimum negative exponent");
  638.     y_str += spacing;
  639.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  640.     sprintf(str,"m increments the map index, changing maps");
  641.     y_str += spacing;
  642.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  643.     sprintf(str,"p or P reverses the colormap for negative/positive exponents");
  644.     y_str += spacing;
  645.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  646.     sprintf(str,"r redraws without recalculating");
  647.     y_str += spacing;
  648.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  649.     sprintf(str,"R redraws, recalculating with new dwell and settle values");
  650.     y_str += spacing;
  651.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  652.     sprintf(str,"s or S spins the colorwheel");
  653.     y_str += spacing;
  654.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  655.     sprintf(str,"u pops back up to the last zoom");
  656.     y_str += spacing;
  657.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  658.     sprintf(str,"U pops back up to the first picture");
  659.     y_str += spacing;
  660.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  661.     sprintf(str,"v or V displays the values of various settings");
  662.     y_str += spacing;
  663.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  664.     sprintf(str,"w decrements, W increments the color wheel index");
  665.     y_str += spacing;
  666.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  667.     sprintf(str,"x or X clears the window");
  668.     y_str += spacing;
  669.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  670.     sprintf(str,"q or Q exits");
  671.     y_str += 2*spacing;
  672.     XDrawImageString(dpy,help,gc,x_str,y_str,str,strlen(str));
  673. }
  674.  
  675. void
  676. check_params(mapnum, parnum)
  677. int mapnum;
  678. int parnum;
  679. {
  680.  
  681.     if (parnum != 1) {
  682.         if ((max_a > pmaxs[mapnum]) || (min_a < pmins[mapnum])) {
  683.             fprintf(stderr, "Warning! Parameter 'a' out of range.\n");
  684.             fprintf(stderr, "You have requested a range of (%lf,%lf).\n",
  685.                     min_a,max_a);
  686.             fprintf(stderr, "Valid range is (%lf,%lf).\n",
  687.                     pmins[mapnum],pmaxs[mapnum]);
  688.         }
  689.     }
  690.     if (parnum != 0) {
  691.         if ((max_b > pmaxs[mapnum]) || (min_b < pmins[mapnum])) {
  692.             fprintf(stderr, "Warning! Parameter 'b' out of range.\n");
  693.             fprintf(stderr, "You have requested a range of (%lf,%lf).\n",
  694.                     min_b,max_b);
  695.             fprintf(stderr, "Valid range is (%lf,%lf).\n",
  696.                     pmins[mapnum],pmaxs[mapnum]);
  697.         }
  698.     }
  699. }
  700.  
  701. void
  702. usage()
  703. {
  704.     fprintf(stderr,"lyap [-Ls][-W#][-H#][-a#][-b#][-w#][-h#][-x xstart]\n");
  705.     fprintf(stderr,"\t[-M#][-S#][-D#][-f string][-r#][-O#][-C#][-c#][-m#]\n");
  706. #ifdef MAPS
  707.     fprintf(stderr,"\t[-F string]\n");
  708. #endif
  709.     fprintf(stderr,"\tWhere: -C# specifies the minimum color index\n");
  710.     fprintf(stderr,"\t       -r# specifies the maxzimum rgb value\n");
  711.     fprintf(stderr,"\t       -u displays this message\n");
  712.     fprintf(stderr,"\t       -a# specifies the minimum horizontal parameter\n");
  713.     fprintf(stderr,"\t       -b# specifies the minimum vertical parameter\n");
  714.     fprintf(stderr,"\t       -w# specifies the horizontal parameter range\n");
  715.     fprintf(stderr,"\t       -h# specifies the vertical parameter range\n");
  716.     fprintf(stderr,"\t       -D# specifies the dwell\n");
  717.     fprintf(stderr,"\t       -S# specifies the settle\n");
  718.     fprintf(stderr,"\t       -H# specifies the initial window height\n");
  719.     fprintf(stderr,"\t       -W# specifies the initial window width\n");
  720.     fprintf(stderr,"\t       -O# specifies the color offset\n");
  721.     fprintf(stderr,"\t       -c# specifies the desired color wheel\n");
  722.     fprintf(stderr,"\t       -m# specifies the desired map (0-4)\n");
  723.     fprintf(stderr,"\t       -f aabbb specifies a forcing function of 00111\n");
  724. #ifdef MAPS
  725.     fprintf(stderr,"\t       -F 00111 specifies the function forcing function\n");
  726. #endif
  727.     fprintf(stderr,"\t       -L indicates use log(x)+log(y) rather than log(xy)\n");
  728.     fprintf(stderr,"\tDuring display :\n");
  729.     fprintf(stderr,"\t     Use the mouse to zoom in on an area\n");
  730.     fprintf(stderr,"\t     e or E recalculates color indices\n");
  731.     fprintf(stderr,"\t     f or F saves exponents to a file\n");
  732.     fprintf(stderr,"\t     KJmn increase/decrease minimum negative exponent\n");
  733.     fprintf(stderr,"\t     r or R redraws\n");
  734.     fprintf(stderr,"\t     s or S spins the colorwheel\n");
  735.     fprintf(stderr,"\t     w or W changes the color wheel\n");
  736.     fprintf(stderr,"\t     x or X clears the window\n");
  737.     fprintf(stderr,"\t     q or Q exits\n");
  738.     exit(1);
  739. }
  740.  
  741. Cycle_frames()
  742. {
  743.     static int i;
  744.     extern void redraw();
  745.  
  746.     for (i=0;i<=maxframe;i++)
  747.         redraw(exponents[i], expind[i], 1); 
  748. }
  749.  
  750. void
  751. print_info() 
  752. {
  753.     Window r, child;
  754.     int x, y;
  755.     unsigned int bw, d, new_w, new_h;
  756.     static int i;
  757.     static char str[80];
  758.     static int x_str, y_str, spacing;
  759.     static int ascent, descent, dir;
  760.     static XCharStruct overall;
  761.     static GC gc;
  762.     static XWindowAttributes attr;
  763.     static char *Mapname[NUMMAPS] = { "logistic", "circle", "leftlog",
  764.             "rightlog", "doublelog" };
  765.  
  766.     XGetWindowAttributes(dpy, info, &attr);
  767.     if (attr.map_state == IsUnmapped)
  768.         return;
  769.     if (XGetGeometry(dpy,canvas,&r,&x,&y,&new_w,&new_h,&bw,&d) != 0)
  770.         XTranslateCoordinates(dpy, canvas, r, x, y, &x, &y, &child);
  771.  
  772.     gc = Data_GC[1];
  773.     XClearWindow(dpy, info);
  774.     x_str = 10;
  775.     y_str = 20;
  776.     sprintf(str,"Current values are as follows : ");
  777.     XDrawImageString(dpy,info,gc,x_str,y_str,str,strlen(str));
  778.     XQueryTextExtents(dpy,(XID)XGContextFromGC(gc),"Hey!",
  779.             4,&dir,&ascent,&descent,&overall);
  780.     spacing = ascent + descent + 5;
  781.     y_str += 2 * spacing;
  782.     sprintf(str,"Using the %s map with geometry = %dx%d+%d+%d", 
  783.             Mapname[mapindex], new_w, new_h, x, y);
  784.     XDrawImageString(dpy,info,gc,x_str,y_str,str,strlen(str));
  785.     sprintf(str,"width = %d height = %d run = %d",width, height, run);
  786.     y_str += spacing;
  787.     XDrawImageString(dpy,info,gc,x_str,y_str,str,strlen(str));
  788.     sprintf(str,"minlyap = %lf minexp = %lf maxexp = %lf",
  789.             minlyap,minexp,maxexp);
  790.     y_str += spacing;
  791.     XDrawImageString(dpy,info,gc,x_str,y_str,str,strlen(str));
  792.     sprintf(str,"settle = %d  dwell = %d start_x = %lf",settle,dwell, start_x);
  793.     y_str += spacing;
  794.     XDrawImageString(dpy,info,gc,x_str,y_str,str,strlen(str));
  795.     sprintf(str,"min_a = %lf  a_rng = %lf  max_a = %lf",min_a,a_range,max_a);
  796.     y_str += spacing;
  797.     XDrawImageString(dpy,info,gc,x_str,y_str,str,strlen(str));
  798.     sprintf(str,"min_b = %lf  b_rng = %lf  max_b = %lf",min_b,b_range,max_b);
  799.     y_str += spacing;
  800.     XDrawImageString(dpy,info,gc,x_str,y_str,str,strlen(str));
  801.     y_str += spacing;
  802.     if (Rflag) {
  803.         sprintf(str,"pseudo-random forcing");
  804.         XDrawImageString(dpy,info,gc,x_str,y_str,str,strlen(str));
  805.     }
  806.     else if (force) {
  807.         sprintf(str,"periodic forcing = ");
  808.         XDrawImageString(dpy,info,gc,x_str,y_str,str,strlen(str));
  809.         for (i=0;i<maxindex;i++) {
  810.             x_str += XTextWidth(XQueryFont(dpy, XGContextFromGC(gc)), 
  811.                             str, strlen(str));
  812.             sprintf(str,"%d",forcing[i]);
  813.             XDrawImageString(dpy,info,gc,x_str,y_str,str,strlen(str));
  814.         }
  815.         x_str = 10;
  816.     }
  817.     else {
  818.         sprintf(str,"periodic forcing = 01");
  819.         XDrawImageString(dpy,info,gc,x_str,y_str,str,strlen(str));
  820.     }
  821.     if (Force) {
  822.         sprintf(str,"function forcing = ");
  823.         y_str += spacing;
  824.         XDrawImageString(dpy,info,gc,x_str,y_str,str,strlen(str));
  825.         for (i=0;i<funcmaxindex;i++) {
  826.             x_str += XTextWidth(XQueryFont(dpy, XGContextFromGC(gc)), 
  827.                             str, strlen(str));
  828.             sprintf(str,"%d",Forcing[i]);
  829.             XDrawImageString(dpy,info,gc,x_str,y_str,str,strlen(str));
  830.         }
  831.         x_str = 10;
  832.     }
  833.     sprintf(str,"numcolors = %d delay = %d stripe interval = %d",
  834.             numcolors, delay, stripe_interval);
  835.     y_str += spacing;
  836.     XDrawImageString(dpy,info,gc,x_str,y_str,str,strlen(str));
  837. }
  838.  
  839. Getkey(event)
  840. XKeyEvent *event;
  841. {
  842.     unsigned char key;
  843.     static int i, running, spindir=0, spinning=0;
  844.     static XWindowAttributes attr;
  845.     extern void init_color(), recalc(), print_help();
  846.     extern void go_init(), go_back(), go_down(), Clear(), write_cmap();
  847.     extern void save_to_file(), Redraw(), redraw(), store_to_file();
  848.  
  849.     if (XLookupString(event, (char *)&key, sizeof(key), (KeySym *)0,
  850.             (XComposeStatus *) 0) > 0)
  851.         switch (key) {
  852.             case '\015': /* write out current colormap to $HOME/.<prog>map */
  853.                     write_cmap(dpy,cmap,Colors,maxcolor,"lyap","Lyap");
  854.                     break;
  855.             case '(': delay -= 25; if (delay < 0) delay = 0; break;
  856.             case ')': delay += 25; break;
  857.             case '<': dwell /= 2; if (dwell < 1) dwell = 1; break;
  858.             case '>': dwell *= 2; break;
  859.             case '[': settle /= 2; if (settle < 1) settle = 0; break;
  860.             case ']': settle *= 2; if (settle < 1) settle = 1; break;
  861.             case 'd': go_down(); break;
  862.             case 'D': FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points, 
  863.                                 0, maxcolor); 
  864.                     break;
  865.             case 'e':
  866.             case 'E': FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points, 
  867.                                 0, maxcolor);
  868.                   dorecalc = (!dorecalc);
  869.                   if (dorecalc)
  870.                       recalc(); 
  871.                   else {
  872.                       maxexp = minlyap; minexp = -1.0 * minlyap;
  873.                   }
  874.                   redraw(exponents[frame], expind[frame], 1);
  875.                   break;
  876.             case 'f': FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points, 
  877.                         0, maxcolor);
  878.                       store_to_file(); break;
  879.             case 'F': FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points, 
  880.                         0, maxcolor);
  881.                       save_to_file(); break;
  882.             case 'i': 
  883.                   if (stripe_interval > 0) {
  884.                       stripe_interval--;
  885.                       if (displayplanes > 1) {
  886.                           running = run; run = 0;
  887.                           init_color(dpy,canvas,cmap,Colors,startcolor,
  888.                                mincolindex,numcolors,numwheels,"lyap","Lyap",0);
  889.                           run = running;
  890.                       }
  891.                   }
  892.                   break;
  893.             case 'I': stripe_interval++;
  894.                   if (displayplanes > 1) {
  895.                       running = run; run = 0;
  896.                       init_color(dpy,canvas,cmap,Colors,startcolor,
  897.                             mincolindex,numcolors,numwheels,"lyap","Lyap",0);
  898.                       run = running;
  899.                   }
  900.                   break;
  901.             case 'K': if (minlyap > 0.05)
  902.                           minlyap -= 0.05;
  903.                    break;
  904.             case 'J': minlyap += 0.05; 
  905.                    break;
  906.             case 'm': mapindex++;
  907.                   if (mapindex >= NUMMAPS)
  908.                       mapindex=0;
  909.                   map = Maps[mapindex];
  910.                   deriv = Derivs[mapindex];
  911.                   if (!aflag)
  912.                       min_a = amins[mapindex];
  913.                   if (!wflag)
  914.                       a_range = aranges[mapindex];
  915.                   if (!bflag)
  916.                       min_b = bmins[mapindex];
  917.                   if (!hflag)
  918.                       b_range = branges[mapindex];
  919.                   if (!Force)
  920.                       for (i=0;i<FUNCMAXINDEX;i++)
  921.                            Forcing[i] = mapindex;
  922.                   max_a = min_a + a_range;
  923.                   max_b = min_b + b_range;
  924.                   a_minimums[0] = min_a; b_minimums[0] = min_b;
  925.                   a_maximums[0] = max_a; b_maximums[0] = max_b;
  926.                   a_inc = a_range / (double)width;
  927.                   b_inc = b_range / (double)height;
  928.                   point.x = 0;
  929.                   point.y = 0;
  930.                   a = rubber_data.p_min = min_a;
  931.                   b = rubber_data.q_min = min_b;
  932.                   rubber_data.p_max = max_a;
  933.                   rubber_data.q_max = max_b;
  934.                   Clear();
  935.                   run = 1;
  936.                   break;
  937.             case 'M': if (minlyap > 0.005)
  938.                       minlyap -= 0.005;
  939.                   break;
  940.             case 'N': minlyap += 0.005;
  941.                   break;
  942.             case 'p':
  943.             case 'P': negative = (!negative); 
  944.                   FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points, 
  945.                                 0, maxcolor); 
  946.                   redraw(exponents[frame], expind[frame], 1); 
  947.                   break;
  948.             case 'r': FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points, 
  949.                             0, maxcolor); 
  950.                   redraw(exponents[frame],expind[frame],1); 
  951.                   break;
  952.             case 'R': FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points, 
  953.                             0, maxcolor); 
  954.                   Redraw(); 
  955.                   break;
  956.             case 'S': spinning=0; break;
  957.             case 's': spinning=1; spindir=(!spindir); 
  958.                   Spin(dpy,cmap,Colors,startcolor,numcolors,delay,spindir); 
  959.                   break;
  960.             case 'u': go_back(); break;
  961.             case 'U': go_init(); break;
  962.             case 'v':
  963.             case 'V': XGetWindowAttributes(dpy, info, &attr);
  964.                     if (attr.map_state != IsUnmapped)
  965.                         XUnmapWindow(dpy, info);
  966.                     else {
  967.                         XMapRaised(dpy, info);
  968.                         print_info(); 
  969.                     }
  970.                     break;
  971.             case '\027': /* (ctrl-W) read color palette from $HOME/.lyapmap */
  972.                   numwheels = 0;
  973.                   if (displayplanes > 1) {
  974.                       running = run; run = 0;
  975.                       init_color(dpy,canvas,cmap,Colors,startcolor,
  976.                             mincolindex,numcolors,numwheels,"lyap","Lyap",0);
  977.                       run = running;
  978.                   }
  979.                   break;
  980.             case 'W': if (numwheels < MAXWHEELS)
  981.                           numwheels++;
  982.                       else
  983.                           numwheels = 0;
  984.                       if (displayplanes > 1) {
  985.                           running = run; run = 0;
  986.                             init_color(dpy,canvas,cmap,Colors,startcolor,
  987.                                mincolindex,numcolors,numwheels,"lyap","Lyap",0);
  988.                           run = running;
  989.                       }
  990.                       break;
  991.             case 'w': if (numwheels > 0)
  992.                           numwheels--;
  993.                       else
  994.                           numwheels = MAXWHEELS;
  995.                       if (displayplanes > 1) {
  996.                           running = run; run = 0;
  997.                             init_color(dpy,canvas,cmap,Colors,startcolor,
  998.                                mincolindex,numcolors,numwheels,"lyap","Lyap",0);
  999.                           run = running;
  1000.                       }
  1001.                       break;
  1002.             case 'x': Clear(); break;
  1003.             case 'X': Destroy_frame(); break;
  1004.             case 'z': running = run; run = 0;
  1005.                   Cycle_frames(); 
  1006.                   run = running; redraw(exponents[frame], expind[frame], 1);
  1007.                   break;
  1008.             case 'Z': running = run; run = 0;
  1009.                   while (!XPending(dpy)) Cycle_frames(); 
  1010.                   run = running; redraw(exponents[frame], expind[frame], 1); 
  1011.                   break;
  1012.             case 'q':
  1013.             case 'Q': Cleanup(); exit(0); break;
  1014.             case '?':
  1015.             case 'h':
  1016.             case 'H': XGetWindowAttributes(dpy, help, &attr);
  1017.                     if (attr.map_state != IsUnmapped)
  1018.                         XUnmapWindow(dpy, help);
  1019.                     else {
  1020.                         XMapRaised(dpy, help);
  1021.                         print_help(); 
  1022.                     }
  1023.                     break;
  1024.             default:  break;
  1025.         }
  1026.         if (spinning)
  1027.             Spin(dpy,cmap,Colors,startcolor,numcolors,delay,spindir); 
  1028.         print_info();
  1029. }
  1030.  
  1031. /* Here's where we index into a color map. After the Lyapunov exponent is
  1032.  * calculated, it is used to determine what color to use for that point.
  1033.  * I suppose there are a lot of ways to do this. I used the following :
  1034.  * if it's non-negative then there's a reserved area at the lower range
  1035.  * of the color map that i index into. The ratio of some "minimum exponent
  1036.  * value" and the calculated value is used as a ratio of how high to index
  1037.  * into this reserved range. Usually these colors are dark red (see init_color).
  1038.  * If the exponent is negative, the same ratio (expo/minlyap) is used to index
  1039.  * into the remaining portion of the colormap (which is usually some light
  1040.  * shades of color or a rainbow wheel). The coloring scheme can actually make
  1041.  * a great deal of difference in the quality of the picture. Different colormaps
  1042.  * bring out different details of the dynamics while different indexing
  1043.  * algorithms also greatly effect what details are seen. Play around with this.
  1044.  */
  1045. sendpoint(expo)
  1046. double expo;
  1047. {
  1048.     static int index;
  1049.     static double tmpexpo;
  1050.     extern double exp();
  1051.  
  1052.     tmpexpo = (negative) ? expo : -1.0 * expo;
  1053.     if (tmpexpo > 0) {
  1054.         if (displayplanes >1) {
  1055.             index = (int)(tmpexpo*lowrange/maxexp);
  1056.             index = (index % lowrange) + startcolor;
  1057.         }
  1058.         else
  1059.             index = 0;
  1060.     }
  1061.     else {
  1062.         if (displayplanes >1) {
  1063.             index = (int)(tmpexpo*numfreecols/minexp);
  1064.             index = (index % numfreecols) + mincolindex;
  1065.         }
  1066.         else
  1067.             index = 1;
  1068.     }
  1069.     BufferPoint(dpy, canvas, pixmap, Data_GC, &Points, index, point.x, point.y);
  1070.     point.x++;
  1071.     a += a_inc;
  1072.     if (save)
  1073.         exponents[frame][expind[frame]++] = expo;
  1074.     if (point.x >= width) {
  1075.         point.y++;
  1076.         point.x = 0;
  1077.         if (save) {
  1078.             b += b_inc;
  1079.             a = min_a;
  1080.         }
  1081.         if (point.y >= height)
  1082.             return FALSE;
  1083.         else
  1084.             return TRUE;
  1085.     }
  1086.     return TRUE;
  1087. }
  1088.  
  1089. void 
  1090. redisplay (w, event)
  1091. Window          w;
  1092. XExposeEvent    *event;
  1093. {
  1094.     if (event->window == help)
  1095.         print_help();
  1096.     else if (event->window == info)
  1097.         print_info();
  1098.     else {
  1099.         /*
  1100.          * Extract the exposed area from the event and copy
  1101.          * from the saved pixmap to the window.
  1102.          */
  1103.         XCopyArea(dpy, pixmap, canvas, Data_GC[0], 
  1104.                event->x,event->y,event->width,event->height,event->x, event->y);
  1105.     }
  1106. }
  1107.  
  1108. void
  1109. resize()
  1110. {
  1111.     Window r;
  1112.     int n, x, y;
  1113.     unsigned int bw, d, new_w, new_h;
  1114.     static XWindowAttributes attr;
  1115.     extern void Clear(), Redraw();
  1116.  
  1117.     if (restfile && demo)
  1118.         return;
  1119.     XGetGeometry(dpy,canvas,&r,&x,&y,&new_w,&new_h,&bw,&d);
  1120.     if ((new_w == width) && (new_h == height)) {
  1121.         print_info();
  1122.         return;
  1123.     }
  1124.     freemem();
  1125.     width = new_w; height = new_h;
  1126.     XClearWindow(dpy, canvas);
  1127.     if (pixmap)
  1128.         XFreePixmap(dpy, pixmap);
  1129.     pixmap = XCreatePixmap(dpy, DefaultRootWindow(dpy), 
  1130.             width, height, DefaultDepth(dpy, screen));
  1131.     a_inc = a_range / (double)width;
  1132.     b_inc = b_range / (double)height;
  1133.     point.x = 0;
  1134.     point.y = 0;
  1135.     run = 1;
  1136.     a = rubber_data.p_min = min_a;
  1137.     b = rubber_data.q_min = min_b;
  1138.     rubber_data.p_max = max_a;
  1139.     rubber_data.q_max = max_b;
  1140.     setupmem();
  1141.     for (n=0;n<MAXFRAMES;n++)
  1142.         if ((n <= maxframe) && (n != frame))
  1143.             resized[n] = 1;
  1144.     InitBuffer(&Points, maxcolor);
  1145.     Clear();
  1146.     print_info();
  1147.     Redraw();
  1148. }
  1149.  
  1150. void
  1151. redraw(exparray, index, cont)
  1152. double *exparray;
  1153. int index, cont;
  1154. {
  1155.     static int i;
  1156.     static int x_sav, y_sav;
  1157.  
  1158.     x_sav = point.x;
  1159.     y_sav = point.y;
  1160.  
  1161.     point.x = 0;
  1162.     point.y = 0;
  1163.  
  1164.     save=0;
  1165.     for (i=0;i<index;i++)
  1166.         sendpoint(exparray[i]);
  1167.     save=1;
  1168.     
  1169.     if (cont) {
  1170.         point.x = x_sav;
  1171.         point.y = y_sav;
  1172.     }
  1173.     else {
  1174.         a = point.x * a_inc + min_a;
  1175.         b = point.y * b_inc + min_b;
  1176.     }
  1177.     FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points, 0, maxcolor);
  1178. }
  1179.  
  1180. void
  1181. Redraw() 
  1182. {
  1183.     FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points, 0, maxcolor);
  1184.     point.x = 0;
  1185.     point.y = 0;
  1186.     run = 1;
  1187.     a = min_a;
  1188.     b = min_b;
  1189.     expind[frame] = 0;
  1190.     resized[frame] = 0;
  1191. }
  1192.  
  1193. /* Store colormap indices in file so we can read them in later */
  1194. void
  1195. store_to_file() 
  1196. {
  1197.     FILE *outfile;
  1198.     unsigned char c;
  1199.     XImage *ximage;
  1200.     static int i,j;
  1201.  
  1202.     outfile = fopen(savname,"w");
  1203.     if(!outfile) {
  1204.         perror(savname);
  1205.         Cleanup();
  1206.         exit(-1);
  1207.     }
  1208.     fprintf(outfile,"# width=%d height=%d \n",width,height);
  1209.     fprintf(outfile,"# settle=%d dwell=%d start_x=%lf \n",settle,dwell,start_x);
  1210.     fprintf(outfile,"# min_a=%lf a_rng=%lf max_a=%lf \n",min_a,a_range,max_a);
  1211.     fprintf(outfile,"# min_b=%lf b_rng=%lf max_b=%lf \n",min_b,b_range,max_b);
  1212.     fprintf(outfile,"# run=%d point.x=%d point.y=%d \n",run,point.x,point.y);
  1213.     fprintf(outfile,"# negative=%d mincolindex=%d startcolor=%d \n",
  1214.             negative,mincolindex,startcolor);
  1215.     fprintf(outfile,"# minexp=%lf maxexp=%lf a=%lf b=%lf \n",
  1216.             minexp,maxexp,a,b);
  1217.     fprintf(outfile,"# Force=%d force=%d useprod=%d \n",Force,force,useprod);
  1218.     fprintf(outfile,"# mapindex=%d maxindex=%d \n", mapindex, maxindex);
  1219.     fprintf(outfile,"# numwheels=%d Rflag=%d color_offset=%d \n",
  1220.             numwheels,Rflag,color_offset);
  1221.     fprintf(outfile,"# maxcolor=%d start_x=%lf minlyap=%lf prob=%lf\n",
  1222.             maxcolor,start_x,minlyap,prob);
  1223.     fprintf(outfile,"# maxindex=%d periodic forcing=", maxindex);
  1224.         for (i=0;i<maxindex;i++) {
  1225.             fprintf(outfile," %d",forcing[i]);
  1226.         }
  1227.     fprintf(outfile," \n");
  1228.     fprintf(outfile,"# funcmaxindex=%d function forcing=", funcmaxindex);
  1229.     for (i=0;i<funcmaxindex;i++)
  1230.          fprintf(outfile," %d",Forcing[i]);
  1231.     fprintf(outfile," \n");
  1232.     fprintf(outfile,"# numcolors=%d\n",numcolors);
  1233.  
  1234.     ximage=XGetImage(dpy, pixmap, 0, 0, 
  1235.             (unsigned int)width, (unsigned int)height, AllPlanes, ZPixmap);
  1236.  
  1237.     for (j=0;j<=point.y;j++)
  1238.         for (i=0;i<width;i++) {
  1239.             c = (unsigned char)XGetPixel(ximage,i,j);
  1240.             fwrite((char *)&c,sizeof c,1,outfile);
  1241.         }
  1242.     fclose(outfile);
  1243. }
  1244.  
  1245. /* Store color pics in PPM format and monochrome in PGM */
  1246. void
  1247. save_to_file() 
  1248. {
  1249.     FILE *outfile;
  1250.     unsigned char c;
  1251.     XImage *ximage;
  1252.     static int i,j;
  1253.     struct Colormap {
  1254.         unsigned char red;
  1255.         unsigned char green;
  1256.         unsigned char blue;
  1257.     };
  1258.     struct Colormap *colormap=NULL;
  1259.  
  1260.     if (colormap)
  1261.         free(colormap);
  1262.     if ((colormap=
  1263.         (struct Colormap *)malloc(sizeof(struct Colormap)*maxcolor))
  1264.             == NULL) {
  1265.         fprintf(stderr,"Error malloc'ing colormap array\n");
  1266.         Cleanup();
  1267.         exit(-1);
  1268.     }
  1269.     outfile = fopen(outname,"w");
  1270.     if(!outfile) {
  1271.         perror(outname);
  1272.         Cleanup();
  1273.         exit(-1);
  1274.     }
  1275.  
  1276.     ximage=XGetImage(dpy, pixmap, 0, 0, width, height, AllPlanes, ZPixmap);
  1277.  
  1278.     if (displayplanes > 1) {
  1279.         for (i=0;i<maxcolor;i++) {
  1280.             colormap[i].red=(unsigned char)(Colors[i].red >> 8);
  1281.             colormap[i].green=(unsigned char)(Colors[i].green >> 8);
  1282.             colormap[i].blue =(unsigned char)(Colors[i].blue >> 8);
  1283.         }
  1284.         fprintf(outfile,"P%d %d %d\n",6,width,height);
  1285.     }
  1286.     else
  1287.         fprintf(outfile,"P%d %d %d\n",5,width,height);
  1288.     fprintf(outfile,"# settle=%d  dwell=%d start_x=%lf\n",settle,dwell,start_x);
  1289.     fprintf(outfile,"# min_a=%lf  a_rng=%lf  max_a=%lf\n",min_a,a_range,max_a);
  1290.     fprintf(outfile,"# min_b=%lf  b_rng=%lf  max_b=%lf\n",min_b,b_range,max_b);
  1291.     if (Rflag)
  1292.         fprintf(outfile,"# pseudo-random forcing\n");
  1293.     else if (force) {
  1294.         fprintf(outfile,"# periodic forcing=");
  1295.         for (i=0;i<maxindex;i++) {
  1296.             fprintf(outfile,"%d",forcing[i]);
  1297.         }
  1298.         fprintf(outfile,"\n");
  1299.     }
  1300.     else
  1301.          fprintf(outfile,"# periodic forcing=01\n");
  1302.     if (Force) {
  1303.          fprintf(outfile,"# function forcing=");
  1304.          for (i=0;i<funcmaxindex;i++)
  1305.              fprintf(outfile,"%d",Forcing[i]);
  1306.          fprintf(outfile,"\n");
  1307.     }
  1308.     fprintf(outfile,"%d\n",numcolors-1);
  1309.  
  1310.     for (j=0;j<height;j++)
  1311.         for (i=0;i<width;i++) {
  1312.             c = (unsigned char)XGetPixel(ximage,i,j);
  1313.             if (displayplanes > 1)
  1314.                 fwrite((char *)&colormap[c],sizeof colormap[0],1,outfile);
  1315.             else
  1316.                 fwrite((char *)&c,sizeof c,1,outfile);
  1317.         }
  1318.     fclose(outfile);
  1319. }
  1320.  
  1321. /* Read saved file with parameters */
  1322. void
  1323. restor_params() 
  1324. {
  1325.     unsigned char s[16];
  1326.     static int i;
  1327.  
  1328.     infile = fopen(inname,"r");
  1329.     if(!infile) {
  1330.         perror(inname);
  1331.         Cleanup();
  1332.         exit(-1);
  1333.     }
  1334.  
  1335.     fscanf(infile,"%1s%6s%d%7s%d",s,s,&width,s,&height);
  1336.     fscanf(infile,"%1s%7s%d%6s%d%8s%lf",s,s,&settle,s,&dwell,s,&start_x);
  1337.     fscanf(infile,"%1s%6s%lf%6s%lf%6s%lf", s,s,&min_a,s,&a_range,s,&max_a);
  1338.     fscanf(infile,"%1s%6s%lf%6s%lf%6s%lf", s,s,&min_b,s,&b_range,s,&max_b);
  1339.     fscanf(infile,"%1s%4s%d%8s%d%8s%d",s,s,&run,s,&point.x,s,&point.y);
  1340.     fscanf(infile,"%1s%9s%d%12s%d%11s%d",
  1341.             s,s,&negative,s,&mincolindex,s,&startcolor);
  1342.     fscanf(infile,"%1s%7s%lf%7s%lf%2s%lf%2s%lf",
  1343.             s,s,&minexp,s,&maxexp,s,&a,s,&b);
  1344.     fscanf(infile,"%1s%6s%d%6s%d%8s%d",s,s,&Force,s,&force,s,&useprod);
  1345.     fscanf(infile,"%1s%9s%d%9s%d",
  1346.             s,s,&mapindex,s,&maxindex);
  1347.     fscanf(infile,"%1s%10s%d%6s%d%13s%d",
  1348.             s,s,&numwheels,s,&Rflag,s,&color_offset);
  1349.     fscanf(infile,"%1s%9s%d%8s%lf%8s%lf%5s%lf",
  1350.             s,s,&maxcolor,s,&start_x,s,&minlyap,s,&prob);
  1351.     fscanf(infile,"%1s%9s%d%s%s", s,s,&maxindex,s,s);
  1352.     for(i=0;i<maxindex;i++) {
  1353.         fscanf(infile, "%d", &forcing[i]);
  1354.     }
  1355.     fscanf(infile,"%1s%13s%d%s%s", s,s,&funcmaxindex,s,s);
  1356.     for(i=0;i<funcmaxindex;i++)
  1357.         fscanf(infile, "%d", &Forcing[i]);
  1358.     fscanf(infile, "%1s%10s%d", s,s,&numcolors);
  1359.     fread(s, sizeof s[0], 1, infile);
  1360.     map = Maps[mapindex];
  1361.     deriv = Derivs[mapindex];
  1362. }
  1363.  
  1364. void
  1365. restor_picture() {
  1366.     static int i, j, k;
  1367.     unsigned char c;
  1368.  
  1369.     for (j=0;j<point.y;j++)
  1370.         for (i=0;i<width;i++) {
  1371.             k = fread(&c, sizeof c, 1, infile);
  1372.             if (k) {
  1373.                 BufferPoint(dpy,canvas,pixmap,Data_GC,&Points,(int)c,i,j);
  1374.                 if (k < mincolindex)
  1375.                     exponents[frame][expind[frame]++] = 
  1376.                         -(double)k/(double)mincolindex;
  1377.                 else
  1378.                     exponents[frame][expind[frame]++] = 
  1379.                         (double)k/(double)numfreecols;
  1380.             }
  1381.             else
  1382.                 break;
  1383.         }
  1384.     for (j=0;j<point.x;j++) {
  1385.         k = fread(&c, sizeof c, 1, infile);
  1386.         if (k) {
  1387.             BufferPoint(dpy,canvas,pixmap,Data_GC,&Points,(int)c,j,point.y);
  1388.             if (k < mincolindex)
  1389.                 exponents[frame][expind[frame]++] = 
  1390.                     -(double)k/(double)mincolindex;
  1391.             else
  1392.                 exponents[frame][expind[frame]++] = 
  1393.                     (double)k/(double)numfreecols;
  1394.         }
  1395.         else
  1396.             break;
  1397.     }
  1398.     fclose(infile);
  1399.     FlushBuffer(dpy, canvas, pixmap, Data_GC, &Points, 0, maxcolor);
  1400.     XCopyArea(dpy, pixmap, canvas, Data_GC[0], 0, 0, width, height, 0, 0);
  1401. }
  1402.  
  1403. void
  1404. recalc() 
  1405. {
  1406.     static int i, index, x, y;
  1407.     
  1408.     minexp = maxexp = 0.0;
  1409.     x = y = 0;
  1410.     for (i=0;i<expind[frame];i++) {
  1411.         if (exponents[frame][i] < minexp)
  1412.             minexp = exponents[frame][i];
  1413.         if (exponents[frame][i] > maxexp)
  1414.             maxexp = exponents[frame][i];
  1415.     }
  1416. }
  1417.  
  1418. void
  1419. Clear() 
  1420. {
  1421.     XFillRectangle(dpy, pixmap, Data_GC[0], 0, 0, width, height);
  1422.     XCopyArea(dpy, pixmap, canvas, Data_GC[0], 
  1423.                 0, 0, width, height, 0, 0);
  1424.     InitBuffer(&Points, maxcolor);
  1425. }
  1426.  
  1427. void
  1428. show_defaults() 
  1429. {
  1430.  
  1431.     printf("Width=%d  Height=%d  numcolors=%d  settle=%d  dwell=%d\n",
  1432.         width,height,numcolors,settle,dwell);
  1433.     printf("min_a=%lf  a_range=%lf  max_a=%lf\n", min_a,a_range,max_a);
  1434.     printf("min_b=%lf  b_range=%lf  max_b=%lf\n", min_b,b_range,max_b);
  1435.     printf("minlyap=%lf  minexp=%lf  maxexp=%lf\n", minlyap,minexp,maxexp);
  1436.     exit(0);
  1437. }
  1438.  
  1439. void
  1440. CreateXorGC()
  1441. {
  1442.     XGCValues values;
  1443.  
  1444.     values.foreground = foreground;
  1445.     values.line_style = LineSolid;
  1446.     values.function = GXxor;
  1447.     RubberGC = XCreateGC(dpy, DefaultRootWindow(dpy),
  1448.           GCForeground | GCBackground | GCFunction | GCLineStyle, &values);
  1449. }
  1450.  
  1451. void 
  1452. StartRubberBand(w, data, event)
  1453. Window w;
  1454. image_data_t *data;
  1455. XEvent *event;
  1456. {
  1457.     XPoint corners[5];
  1458.     extern void SetupCorners();
  1459.  
  1460.     nostart = 0;
  1461.     data->rubber_band.last_x = data->rubber_band.start_x = event->xbutton.x;
  1462.     data->rubber_band.last_y = data->rubber_band.start_y = event->xbutton.y;
  1463.     SetupCorners(corners, data);
  1464.     XDrawLines(dpy, canvas, RubberGC,
  1465.         corners, sizeof(corners) / sizeof(corners[0]), CoordModeOrigin);
  1466. }
  1467.  
  1468. void
  1469. SetupCorners(corners, data)
  1470. XPoint *corners;
  1471. image_data_t *data;
  1472. {
  1473.     corners[0].x = data->rubber_band.start_x;
  1474.     corners[0].y = data->rubber_band.start_y;
  1475.     corners[1].x = data->rubber_band.start_x;
  1476.     corners[1].y = data->rubber_band.last_y;
  1477.     corners[2].x = data->rubber_band.last_x;
  1478.     corners[2].y = data->rubber_band.last_y;
  1479.     corners[3].x = data->rubber_band.last_x;
  1480.     corners[3].y = data->rubber_band.start_y;
  1481.     corners[4] = corners[0];
  1482. }
  1483.  
  1484. void 
  1485. TrackRubberBand(w, data, event)
  1486. Window w;
  1487. image_data_t *data;
  1488. XEvent *event;
  1489. {
  1490.     XPoint corners[5];
  1491.     extern void SetupCorners();
  1492.  
  1493.     if (nostart)
  1494.         return;
  1495.     SetupCorners(corners, data);
  1496.     XDrawLines(dpy, canvas, RubberGC, corners, 
  1497.                sizeof(corners) / sizeof(corners[0]), CoordModeOrigin);
  1498.     data->rubber_band.last_x = Min(event->xbutton.x, width);
  1499.     data->rubber_band.last_y = Min(event->xbutton.y, height);
  1500.     if (data->rubber_band.last_y < data->rubber_band.start_y ||
  1501.         data->rubber_band.last_x < data->rubber_band.start_x) {
  1502.             data->rubber_band.last_y = data->rubber_band.start_y;
  1503.             data->rubber_band.last_x = data->rubber_band.start_x;
  1504.     }
  1505.     SetupCorners(corners, data);
  1506.     XDrawLines(dpy, canvas, RubberGC, corners, 
  1507.                sizeof(corners) / sizeof(corners[0]), CoordModeOrigin);
  1508. }
  1509.  
  1510. void 
  1511. EndRubberBand(w, data, event)
  1512. Window w;
  1513. image_data_t *data;
  1514. XEvent *event;
  1515. {
  1516.     XPoint corners[5];
  1517.     XPoint top, bot;
  1518.     double delta, diff;
  1519.     extern void set_new_params(), SetupCorners();
  1520.  
  1521.     nostart = 1;
  1522.     SetupCorners(corners, data);
  1523.     XDrawLines(dpy, canvas, RubberGC,
  1524.         corners, sizeof(corners) / sizeof(corners[0]), CoordModeOrigin);
  1525.     if (data->rubber_band.start_x >= data->rubber_band.last_x ||
  1526.         data->rubber_band.start_y >= data->rubber_band.last_y)
  1527.         return;
  1528.     top.x = data->rubber_band.start_x;
  1529.     bot.x = data->rubber_band.last_x;
  1530.     top.y = data->rubber_band.start_y;
  1531.     bot.y = data->rubber_band.last_y;
  1532.     diff = data->q_max - data->q_min;
  1533.     delta = (double)top.y / (double)height;
  1534.     data->q_min += diff * delta;
  1535.     delta = (double)(height - bot.y) / (double)height;
  1536.     data->q_max -= diff * delta;
  1537.     diff = data->p_max - data->p_min;
  1538.     delta = (double)top.x / (double)width;
  1539.     data->p_min += diff * delta;
  1540.     delta = (double)(width - bot.x) / (double)width;
  1541.     data->p_max -= diff * delta;
  1542.     fflush(stdout);
  1543.     set_new_params(w, data);
  1544. }
  1545.  
  1546. void
  1547. set_new_params(w, data)
  1548. Window w;
  1549. image_data_t *data;
  1550. {
  1551.     static XWindowAttributes attr;
  1552.     extern void Clear();
  1553.  
  1554.     frame = (maxframe + 1) % MAXFRAMES;
  1555.     if (frame > maxframe)
  1556.         maxframe = frame;
  1557.     a_range = data->p_max - data->p_min;
  1558.     b_range = data->q_max - data->q_min;
  1559.     a_minimums[frame] = min_a = data->p_min;
  1560.     b_minimums[frame] = min_b = data->q_min;
  1561.     a_inc = a_range / (double)width;
  1562.     b_inc = b_range / (double)height;
  1563.     point.x = 0;
  1564.     point.y = 0;
  1565.     run = 1;
  1566.     a = min_a;
  1567.     b = min_b;
  1568.     a_maximums[frame] = max_a = data->p_max;
  1569.     b_maximums[frame] = max_b = data->q_max;
  1570.     expind[frame] = 0;
  1571.     Clear();
  1572.     print_info();
  1573. }
  1574.  
  1575. void
  1576. go_down() 
  1577. {
  1578.     static int i;
  1579.     
  1580.     frame++;
  1581.     if (frame > maxframe)
  1582.         frame = 0;
  1583.     jumpwin();
  1584. }
  1585.  
  1586. void
  1587. go_back() 
  1588. {
  1589.     static int i;
  1590.     
  1591.     frame--;
  1592.     if (frame < 0)
  1593.         frame = maxframe;
  1594.     jumpwin();
  1595. }
  1596.  
  1597. jumpwin()
  1598. {
  1599.     rubber_data.p_min = min_a = a_minimums[frame];
  1600.     rubber_data.q_min = min_b = b_minimums[frame];
  1601.     rubber_data.p_max = max_a = a_maximums[frame];
  1602.     rubber_data.q_max = max_b = b_maximums[frame];
  1603.     a_range = max_a - min_a;
  1604.     b_range = max_b - min_b;
  1605.     a_inc = a_range / (double)width;
  1606.     b_inc = b_range / (double)height;
  1607.     point.x = 0;
  1608.     point.y = 0;
  1609.     a = min_a;
  1610.     b = min_b;
  1611.     Clear();
  1612.     if (resized[frame])
  1613.         Redraw();
  1614.     else
  1615.         redraw(exponents[frame], expind[frame], 0);
  1616. }
  1617.  
  1618. void
  1619. go_init() 
  1620. {
  1621.     static int i;
  1622.     
  1623.     frame = 0;
  1624.     jumpwin();
  1625. }
  1626.  
  1627. Destroy_frame()
  1628. {
  1629.     static int i;
  1630.  
  1631.     for (i=frame; i<maxframe; i++) {
  1632.         exponents[frame] = exponents[frame+1];
  1633.         expind[frame] = expind[frame+1];
  1634.         a_minimums[frame] = a_minimums[frame+1];
  1635.         b_minimums[frame] = b_minimums[frame+1];
  1636.         a_maximums[frame] = a_maximums[frame+1];
  1637.         b_maximums[frame] = b_maximums[frame+1];
  1638.     }
  1639.     maxframe--;
  1640.     go_back();
  1641. }
  1642.  
  1643. freemem()
  1644. {
  1645.     static int i;
  1646.     for (i=0;i<MAXFRAMES;i++)
  1647.         free(exponents[i]);
  1648. }
  1649.  
  1650. setupmem()
  1651. {
  1652.     static int i;
  1653.  
  1654.     for (i=0;i<MAXFRAMES;i++) {
  1655.         if((exponents[i]=
  1656.             (double *)malloc(sizeof(double)*width*height))==NULL){
  1657.                 fprintf(stderr,"Error malloc'ing exponent array.\n");
  1658.                 exit(-1);
  1659.         }
  1660.     }
  1661. }
  1662.  
  1663. setforcing()
  1664. {
  1665.     static int i;
  1666.     extern double drand48();
  1667.  
  1668.     for (i=0;i<MAXINDEX;i++)
  1669.         forcing[i] = (drand48() > prob) ? 0 : 1;
  1670. }
  1671.