home *** CD-ROM | disk | FTP | other *** search
/ Education Sampler 1992 [NeXTSTEP] / Education_1992_Sampler.iso / Programming / Source / Lyapunov / lyap.m < prev    next >
Encoding:
Text File  |  1992-07-31  |  9.2 KB  |  404 lines
  1. // Adapted from lyap.c by Don Yacktman
  2. // This object embodies the lyapunov calculation code.
  3. // Use the methods provided to set the various options,
  4. // and use the rendering methods to render part or all
  5. // of the fractal.  (You can render parts, because I
  6. // eventually would like to be able to farm out image
  7. // segments to other machines on a network, Zilla-like.)
  8. // You can take an already-calculated fractal and apply
  9. // a color map to it through methods provided.  Just pass
  10. // it a color map object that you have created.
  11.  
  12. // Here's the copyright from the original UNIX program "lyap",
  13. // from which I derived much of the engine.
  14. // The original program and an executable and output viewer for
  15. // the NeXT are provided in a separate directory in this distribution.
  16. // As the notice states, you cannot redistribute this program, or
  17. // anything derived from it without including all sources, including
  18. // your own changes!  
  19.  
  20. /*
  21.  * by Garrett Wollman based on a function posted to Usenet by Ed Kubaitis
  22.  * This program is Copyright (C) 1991, Garrett A. Wollman.  This
  23.  * program may be modified and distributed for any purpose and without
  24.  * fee, provided that this notice remains on all such copies
  25.  * unaltered.  Binary distributions not including this source file are
  26.  * prohibited.  Modified versions must be marked with the name of the
  27.  * modifier and the date of modification.
  28.  *
  29.  * Under no circumstances shall Garrett A. Wollman or the University
  30.  * of Vermont and State Agricultural College be held liable for any
  31.  * damages resulting from the use or misuse of this program, whether
  32.  * the author is aware of such possibility or not.  This program is
  33.  * warranted solely to occupy disk space.
  34.  *
  35.  * I'm sorry for all this legal garbage, but it is sadly necessary
  36.  * these days.
  37.  */
  38.  
  39. #include <math.h>
  40. #include <time.h>
  41. #include <stdlib.h>
  42. #include <stdio.h>
  43.  
  44.  
  45. /*
  46.  * Global variables
  47.  */
  48.  
  49. char *whoami;
  50.  
  51. extern int getopt(int,char **,const char *);
  52. extern char *optarg;
  53. extern int optind;
  54. extern int opterr;
  55. extern void perror(const char *);
  56. extern int strlen(const char *);
  57.  
  58. #define OPTIONS "r:c:v:m:d:s:g:5pto:l:"
  59.  
  60. int rows = 512,cols = 512;
  61. double amin,amax;
  62. double bmin,bmax;
  63. double aincr,curA;
  64. double bincr;
  65.  
  66. #define nColors 256        /* don't even try to change this */
  67. struct {
  68.   unsigned char red;
  69.   unsigned char green;
  70.   unsigned char blue;
  71. } colormap[nColors];
  72.  
  73. int Settle = 600;
  74. int Dwell = 2000;
  75. int *Rvec;
  76. double minLya = -5;
  77. int showpos = 0;
  78. double initGuess = 0.5;
  79.  
  80. /*
  81.  * Colormap functions
  82.  */
  83.  
  84. /*
  85.  * set up our initial shades of grey
  86.  */
  87. void init_colormap(void) {
  88.   int i;
  89.   for(i=0; i < nColors; i++)
  90.     colormap[i].red = colormap[i].green = colormap[i].blue = i;
  91. }
  92.  
  93. void read_colormap(const char *name) {
  94.   FILE *cmap;
  95.   char buf[255];        /* magic number */
  96.   int a,b,c,i;
  97.   
  98.   cmap = fopen(name,"r");
  99.   if(!cmap) {
  100.     perror(name);
  101.     exit(-1);
  102.   }
  103.  
  104.   i = 0;
  105.   while(!feof(cmap) && i < nColors) {
  106.     fgets(buf,sizeof buf,cmap);
  107.     /*
  108.      * unparseable lines are comments
  109.      * as is anything after the <r g b> value on a single line
  110.      *
  111.      * per Fractint 15.1 standard
  112.      */
  113.     if(sscanf(buf,"%d %d %d",&a,&b,&c) != 3)
  114.       continue;
  115.  
  116.     colormap[i].red = a;
  117.     colormap[i].green = b;
  118.     colormap[i].blue = c;
  119.     i++;
  120.   }
  121.   fclose(cmap);
  122. }
  123.  
  124.  
  125. /*
  126.  * Lyapunov exponent and coloring calculation
  127.  * original function by Ed Kubaitis, used by permission
  128.  */
  129.  
  130. /*
  131.  * Speedup...
  132.  *
  133.  * Erick B. Wong noticed that $\log_2 d(x_1)+\dots +\log_2 d(x_n)$
  134.  * is, by a principle from high-school algebra, the same as
  135.  * $\log_2 (d(x_1) + \dots + d(x_n))$.  We take advantage of this
  136.  * by unrolling the dwell loop by four (hence the rounding in main())
  137.  * and accumulating by multiplication before taking the log.  Since
  138.  * log() is extremely expensive, this should lead to a considerable
  139.  * speedup, while still allowing complete flexibility in selecting
  140.  * dwell values.
  141.  *
  142.  * Presumably, one might unroll this loop even further, with a 
  143.  * concomitant increase in textual complexity.
  144.  */
  145.  
  146. int lya(double a,double b) {
  147.   double t = 0;
  148.   double r, lya;
  149.   int d;
  150.   double acc;
  151.   double x;
  152.  
  153.   x = initGuess;        /* initialize x */
  154.  
  155.   for(d=0; d < Settle; d++) {
  156.     r = (Rvec[d]) ? b : a;
  157.     x = r*x*(1-x);
  158.   }
  159.  
  160.   for(d=0; d <= Dwell; d+= 4) {
  161.     r = (Rvec[d]) ? b : a;
  162.     x = r * x * (1-x);
  163.     acc = fabs(r - 2*r*x);
  164.  
  165.     r = (Rvec[d+1]) ? b : a;
  166.     x = r * x * (1-x);
  167.     acc *= fabs(r - 2*r*x);
  168.  
  169.     r = (Rvec[d+2]) ? b : a;
  170.     x = r * x * (1-x);
  171.     acc *= fabs(r - 2*r*x);
  172.  
  173.     r = (Rvec[d+3]) ? b : a;
  174.     x = r * x * (1-x);
  175.     t += log(acc * fabs(r - 2*r*x));
  176.  
  177.     if(fabs(x-0.5) < 1e-20) {
  178.       d += 3;
  179.       break;
  180.     }
  181.   }
  182.  
  183.   lya = (t * 1.442695041)/d;
  184.  
  185.   if(showpos)            
  186.     lya *= -1;
  187.   if(lya < minLya)        /* sanity check! */
  188.     lya = minLya;
  189.  
  190.   if(lya < 0) {
  191.     return (int)(lya / minLya * (double)nColors);
  192.   } else {
  193.     return 0;
  194.   }
  195. }
  196.  
  197.  
  198. /*
  199.  * Usage message
  200.  */
  201. #ifdef __GNUC__
  202. volatile
  203. #endif
  204. void usage(void) {
  205.   fprintf(stderr,
  206.       "%s: usage:\n\n"
  207.       "%s [-r rows] [-c cols] [-v program] [-p] [-t]\n\t"
  208.       "[-o file] [-l colormap]\n\t"
  209.       "[-m minLya] [-d Dwell] [-s Settle] [-g initGuess]\n\t"
  210.       "[-5] amin amax bmin bmax bits\n",whoami,whoami);
  211.   exit(-1);
  212. }
  213.  
  214. int main(int argc,char **argv) {
  215.   int c,i,j;
  216.   char *outname = NULL;
  217.   FILE *outfile;
  218.  
  219.   char *bits;
  220.  
  221.   double curB;
  222.   int onRow,onCol;
  223.   int *oneRow;
  224.  
  225.   int outtext = 0;
  226.  
  227.   init_colormap();
  228.   opterr = 0;
  229.  
  230.   whoami = *argv;
  231.   while((c = getopt(argc,argv,OPTIONS)) != EOF) {
  232.     switch(c) {
  233.     case 'r':
  234.       rows = atoi(optarg);
  235.       if(!rows)
  236.         rows = 512;
  237.       break;
  238.  
  239.     case 'c':
  240.       cols = atoi(optarg);
  241.       if(!cols)
  242.         cols = 512;
  243.       break;
  244.  
  245.     case 'o':
  246.       outname = optarg;
  247.       break;
  248.  
  249.     case 'l':
  250.       read_colormap(optarg);
  251.       break;
  252.  
  253.     case 'm':
  254.       if(!sscanf(optarg,"%lf",&minLya))
  255.     usage();
  256.       break;
  257.  
  258.     case 'd':
  259.       Dwell = atoi(optarg);
  260.       if(!Dwell)
  261.     usage();
  262.       /*
  263.        * By forcing the Dwell to be a multiple of four, we can use the
  264.        * speedup noted in lya().  We *could* just do this silently,
  265.        * by if it might make a difference to someone...
  266.        */
  267.       if (Dwell % 4) {
  268.         fprintf(stderr,"Dwell value of %d being rounded off to %d.\n",
  269.           Dwell, Dwell += 4 - (Dwell % 4));
  270.       }
  271.       break;
  272.  
  273.     case 's':
  274.       Settle = atoi(optarg);
  275.       if(!Settle)        /* for settle close to 0, use 1 */
  276.         usage();
  277.       break;
  278.  
  279.     case 'g':
  280.       if(!sscanf(optarg,"%lf",&initGuess))
  281.         usage();
  282.       break;
  283.  
  284.     case 'p':
  285.       showpos = 1;
  286.       break;
  287.  
  288.     case 't':
  289.       outtext = 1;
  290.       break;
  291.  
  292.     default:            /* includes '?' (and '#' when !MULTIPROC) */
  293.       usage();
  294.     }
  295.   }
  296.  
  297.   if(!argv[optind])
  298.     usage();
  299.  
  300.   /*
  301.    * Note to language nit-pickers... The code below is legal because
  302.    * the || operator is a sequence point.
  303.    */
  304.   if(!sscanf(argv[optind++],"%lf",&amin) || !argv[optind] ||
  305.      !sscanf(argv[optind++],"%lf",&amax) || !argv[optind] ||
  306.      !sscanf(argv[optind++],"%lf",&bmin) || !argv[optind] ||
  307.      !sscanf(argv[optind++],"%lf",&bmax) || !argv[optind])
  308.     usage();
  309.  
  310.   bits = argv[optind];
  311.  
  312.   if(Dwell < Settle)
  313.     Dwell = Settle;        /* a bit of sanity */
  314.  
  315.   j = strlen(bits);
  316.  
  317.   Rvec = (int *)malloc((Dwell+1) * sizeof *Rvec);
  318.   if(!Rvec) {
  319.     fputs("Cannot allocate space for Markus vector!\n",stderr);
  320.     exit(-1);
  321.   }
  322.  
  323.   for(i=0; i <= Dwell; i++) {
  324.     if(bits[i % j] == 'a')
  325.       Rvec[i] = 0;
  326.     else if(bits[i % j] == 'b')
  327.       Rvec[i] = 1;
  328.     else {
  329.       fprintf(stderr,"Invalid character \\%o ('%c') in bit string\n",
  330.           (int)(unsigned char)bits[i % j], bits[i % j]);
  331.       exit(-1);
  332.     }
  333.   }
  334.  
  335.   oneRow = (int *)malloc(cols * sizeof *oneRow);
  336.   if(!oneRow) {
  337.     fputs("Cannot allocate space for single-row buffer!\n",stderr);
  338.     exit(-1);
  339.   }
  340.  
  341.   aincr = (amax - amin) / rows;
  342.   bincr = (bmax - bmin) / cols;
  343.  
  344.   if(outname) {
  345.     outfile = fopen(outname,"w");
  346.     if(!outfile) {
  347.       perror(outname);
  348.       exit(-1);
  349.     }
  350.   } else {
  351.     outfile = stdout;
  352.   }
  353.  
  354.   /*
  355.    * Print output header... same as DKBTrace
  356.    */
  357.   if (outtext) {
  358.       fprintf(outfile, "P3 %d %d", cols, rows);
  359.   } else {
  360.   putc('\0',outfile);  putc('\0',outfile);
  361.   putc('\0',outfile);  putc('\0',outfile);
  362.   putc('\0',outfile);  putc('\0',outfile);
  363.   putc('\0',outfile);  putc('\0',outfile);
  364.   putc('\0',outfile);  putc('\0',outfile);
  365.   putc('\0',outfile);  putc('\0',outfile);
  366.   fprintf(stderr, "Sizes: %d %d %d %d\n",cols%256,cols/256,rows%256,rows/256);
  367.   putc((unsigned char)(cols%256),outfile);
  368.   putc((unsigned char)(cols/256),outfile);
  369.   putc((unsigned char)(rows%256),outfile);
  370.   putc((unsigned char)(rows/256),outfile);
  371.   putc('\0',outfile);  putc('\0',outfile);
  372.   }
  373.  
  374.  
  375.   for(onRow = 0, curA = amax; onRow < rows; curA -= aincr, onRow++) {
  376.     for(onCol = cols - 1, curB = bmax; onCol; curB -= bincr, onCol--) {
  377.       oneRow[onCol] = lya(curA,curB);
  378.     }
  379.     
  380.     for(i=0; i < cols; i++) {
  381.       unsigned char c;
  382.       c = (unsigned char)oneRow[i];
  383.       if(outtext) {
  384.           fprintf(outfile,"%d %d %d ",colormap[c].red,colormap[c].green,
  385.           colormap[c].blue);
  386.       } else {
  387.       putc(colormap[c].red,outfile);
  388.       putc(colormap[c].green,outfile);
  389.       putc(colormap[c].blue,outfile);
  390.       }
  391.     }
  392.     if(outtext)
  393.       fputc('\n',outfile);
  394.     if(!(onRow % 16)) {
  395.       fprintf(stderr,"Done row %d/%d \n",onRow,rows);
  396.     }
  397.   }
  398.  
  399.   if(outname)
  400.     fclose(outfile);
  401.  
  402.   return 0;
  403. }
  404.