home *** CD-ROM | disk | FTP | other *** search
/ Hackers Handbook - Millenium Edition / Hackers Handbook.iso / files / c_scripts / chaos.c < prev    next >
Encoding:
C/C++ Source or Header  |  1997-09-25  |  6.3 KB  |  312 lines
  1. #if 0
  2.     cc -o chaos chaos.c -lm -lX11
  3.     exit
  4. #endif
  5.  
  6. /*
  7.  * Chaos program.  Copyright 1991 Ken Shirriff
  8.  * Send comments, etc. to shirriff@sprite.Berkeley.EDU
  9.  *
  10.  * This program draws the "chaotic attractor".
  11.  * Usage: chaos [-r min max] [-i ilow ihigh] [-s step]\n");
  12.  * -r gives the minimum and maximum values of the logistic equation parameter.
  13.  * -i gives the iteration range to be plotted and used for the Lyapunov
  14.  * exponent computation.  (The iterations below ilow are not used, so they
  15.  * can be used to allow the function to settle.
  16.  * -s gives the step size (i.e. how big the steps are in logistic equation
  17.  * parameter space.  1=1 pixel.
  18.  *
  19.  * After drawing, left click prints the X position of the mouse in terms
  20.  * of the logistic equation parameter.
  21.  * Right click ends the program.
  22.  *
  23.  * Compile program with "sh chaos.c"
  24.  *
  25.  */
  26.  
  27.  
  28. #include <X11/Xlib.h>
  29. #include <stdio.h>
  30. #include <math.h>
  31.  
  32. struct Point {
  33.     double x,y;
  34. };
  35. typedef struct Point pos_t;
  36.  
  37. /*
  38.  * SIZE is the size of the X window to use.
  39.  */
  40. #define SIZE 600
  41.  
  42. /*
  43.  * The logistic equation.
  44.  */
  45. double logis(x,s)
  46. double x,s;
  47. {
  48.     double xp;
  49.     xp = s*x*(1-x);
  50.     return xp;
  51. }
  52.  
  53. main(argc,argv)
  54. int argc;
  55. char **argv;
  56. {
  57.     double s;    /* Our logistics equation driving parameter. */
  58.     int i;
  59.     double x;    /* Our iterated x value. */
  60.     pos_t p;
  61.     float lya;    /* Lyapunov exponent. */
  62.     int N0=10, NT=100; /* Low and high iterations to plot. */
  63.     float S0=0, S1=4.0; /* Low and high parameter values. */
  64.     float div=1;    /* Plotting resolution in pixels. */
  65. #define abs(x) ((x)>0?(x):-(x))
  66.  
  67.     if (argc < 1) {
  68. usage:
  69.     fprintf(stderr,"Usage: [-r min max] [-i ilow ihigh] [-s step]\n");
  70.     exit(1);
  71.     }
  72.     while (1) {
  73.     if (argc==1) {
  74.         break;
  75.     } else if (strncmp("-h",argv[1],2)==0) {
  76.         goto usage;
  77.     } else if (strcmp("-r",argv[1])==0 && argc>3) {
  78.         sscanf(argv[2],"%f", &S0);
  79.         sscanf(argv[3],"%f", &S1);
  80.         if (S0>=S1) {
  81.         fprintf(stderr,"Invalid range\n");
  82.         exit(-1);
  83.         }
  84.         argc -= 3;
  85.         argv += 3;
  86.     } else if (strcmp("-i",argv[1])==0 && argc>3) {
  87.         sscanf(argv[2],"%d", &N0);
  88.         sscanf(argv[3],"%d", &NT);
  89.         argc -= 3;
  90.         argv += 3;
  91.     } else if (strcmp("-s",argv[1])==0 && argc>2) {
  92.         sscanf(argv[2],"%f", &div);
  93.         argc -= 2;
  94.         argv += 2;
  95.     } else if (strcmp("-i",argv[1])==0 && argc>3) {
  96.         sscanf(argv[2],"%d", &N0);
  97.         sscanf(argv[3],"%d", &NT);
  98.         argc -= 3;
  99.         argv += 3;
  100.     } else {
  101.         goto usage;
  102.     }
  103.     }
  104.  
  105.     xinit();
  106.     range(S0,-2.,S1,1.5);
  107.  
  108.     /*
  109.      * Loop across the different 'C' values (in the variable s)
  110.      */
  111.     for (s=S0;s<S1;s+=(S1-S0)/SIZE*div) {
  112.     x = .5;
  113.     lya = 0;
  114.     /*
  115.      * Do NT iterations of the logistics equation.
  116.      */
  117.     for (i=0;i<NT;i++) {
  118.         x = logis(x,s);
  119.         if (i>N0) {
  120.         /* 
  121.          * Plot the point.
  122.          */
  123.         p.x = s;
  124.         p.y = x;
  125.         point(&p);
  126.         /*
  127.          * Sum up the Lyapunov exponent.
  128.          */
  129.         lya += log(abs(s-2*x*s))/(NT-N0);
  130.         }
  131.     }
  132.     /*
  133.      * Plot the Lyapunov exponent and 0 axis.
  134.      */
  135.     p.x = s;
  136.     p.y = lya/2-1; /* Scale lya to lower part of window. */
  137.     point(&p);
  138.     p.x = s;
  139.     p.y = 0./2-1;
  140.     point(&p);
  141.     }
  142.  
  143.     /*
  144.      * Flush the drawing buffer.
  145.      */
  146.     drawflush();
  147.     printf("Done\n");
  148.     click();
  149.     endit();
  150. }
  151.  
  152. /*
  153.  * -------------------------------------------------------------
  154.  * X drivers
  155.  */
  156.  
  157.  
  158. Display *dp;
  159. Window w;
  160. GC gc;
  161.  
  162. float Xxmin, Xxmax, Xymin, Xymax;
  163.  
  164. #define SCALEX(v) ((int)(((v)-(Xxmin))*SIZE/(Xxmax-(Xxmin))+.5))
  165. #define ISCALEX(n) (Xxmin+(n)*(Xxmax-Xxmin)/SIZE)
  166. #define SCALEY(v) (SIZE-((int)(((v)-(Xymin))*SIZE/(Xymax-(Xymin))+.5)))
  167. #define ISCALEY(n) (Xymin+(SIZE-n)*(Xymax-Xymin)/SIZE)
  168.  
  169. /*
  170.  * Initialize the X window.
  171.  */
  172.  
  173. xinit()
  174. {
  175.     Window root;
  176.     Screen *sc;
  177.     int dscreen;
  178.     XGCValues gcvals;
  179.     XSetWindowAttributes watt;
  180.     int depth;
  181.     XEvent xevent;
  182.     int done=0;
  183.  
  184.     if ((dp=XOpenDisplay(NULL))==NULL) {
  185.     fprintf(stderr,"Could not open Display!0\n");
  186.     exit(1);
  187.     }
  188.     XSynchronize(dp, True);
  189.     dscreen = XDefaultScreen(dp);
  190.     sc = ScreenOfDisplay(dp,dscreen);
  191.     depth = DefaultDepth(dp,dscreen);
  192.     root = DefaultRootWindow(dp);
  193.     watt.background_pixel = WhitePixelOfScreen(sc);
  194.     watt.bit_gravity = StaticGravity;
  195.     w = XCreateWindow(dp,root,0,0,SIZE,SIZE,0,
  196.        depth,InputOutput,CopyFromParent,CWBackPixel|CWBitGravity, &watt);
  197.     XStoreName(dp,w,"Test");
  198.     gcvals.background = BlackPixelOfScreen(sc);
  199.     gcvals.foreground = BlackPixelOfScreen(sc);
  200.     gc = XCreateGC(dp,w, GCForeground|GCBackground, &gcvals);
  201.     XSetFunction(dp,gc,GXcopy);
  202.     XSelectInput(dp,w,ExposureMask|EnterWindowMask|PointerMotionMask|
  203.     ButtonPressMask|StructureNotifyMask|SubstructureNotifyMask);
  204.     XMapWindow(dp,w);
  205.     XSync(dp,False);
  206.     XClearArea(dp,w,0,0,SIZE, SIZE, False);
  207.     while (!done) {
  208.     XNextEvent(dp, &xevent);
  209.     switch (((XAnyEvent *)&xevent)->type) {
  210.         case ButtonPress:
  211.         if ( ((XButtonEvent*)&xevent)->button == Button3) {
  212.             done = 1;
  213.         }
  214.         break;
  215.         case Expose:
  216.         done = 1;
  217.         break;
  218.         default:
  219.         break;
  220.     }
  221.     }
  222. }
  223.  
  224. /*
  225.  * Set the plotting range.
  226.  */
  227. range(pxmin, pymin, pxmax, pymax)
  228. float pxmin, pxmax, pymin, pymax;
  229. {
  230.     Xxmin = pxmin;
  231.     Xxmax = pxmax;
  232.     Xymin = pymin;
  233.     Xymax = pymax;
  234. }
  235.  
  236. /*
  237.  * Plot the point.  Points are buffered and plotted in large bunches
  238.  * to make it run faster.
  239.  */
  240. #define BUF 1000
  241. XPoint buf[BUF];
  242. int bufp=0;
  243. point(p)
  244. pos_t *p;
  245. {
  246.     if (bufp>=BUF) {
  247.     XDrawPoints(dp,w,gc,buf,bufp,CoordModeOrigin);
  248.     bufp=0;
  249.     }
  250.     buf[bufp].x = SCALEX(p->x);
  251.     buf[bufp].y = SCALEY(p->y);
  252.     bufp++;
  253. }
  254.  
  255. /*
  256.  * Flush the point drawing buffer.
  257.  */
  258. drawflush()
  259. {
  260.     if (bufp>0) {
  261.     XDrawPoints(dp,w,gc,buf,bufp,CoordModeOrigin);
  262.     bufp=0;
  263.     }
  264. }
  265.  
  266. /*
  267.  * Wait for mouse click.
  268.  */
  269. click()
  270. {
  271.     XEvent xevent;
  272.     int done=0;
  273.     while (!done) {
  274.     XNextEvent(dp, &xevent);
  275.     switch (((XAnyEvent *)&xevent)->type) {
  276.         case ButtonPress:
  277.         if ( ((XButtonEvent*)&xevent)->button == Button3) {
  278.             done = 1;
  279.         } else if ( ((XButtonEvent*)&xevent)->button == Button1) {
  280.             /*
  281.              * If we get a click, print the mouse position.
  282.              */
  283.             int x,y;
  284.             float fy;
  285.             x =  ((XButtonEvent *)&xevent)->x;
  286.             y =  ((XButtonEvent *)&xevent)->y;
  287.             fy = ISCALEY(y);
  288.             /*
  289.              * Rescale the Lyapunov window coordinates.
  290.              */
  291.             if (fy<-.2) {
  292.             printf("fy=%f\n", fy);
  293.             fy = (fy+1)*2;
  294.             }
  295.             printf("(%f,%f)\n", ISCALEX(x), fy);
  296.         }
  297.         break;
  298.         default:
  299.         break;
  300.     }
  301.     }
  302. }
  303.  
  304. /*
  305.  * Close the window.
  306.  */
  307. endit()
  308. {
  309.     XDestroyWindow(dp, w);
  310.     XCloseDisplay(dp);
  311. }
  312.