// Adapted from lyap.c by Don Yacktman
// This object embodies the lyapunov calculation code.
// Use the methods provided to set the various options,
// and use the rendering methods to render part or all
// of the fractal.  (You can render parts, because I
// eventually would like to be able to farm out image
// segments to other machines on a network, Zilla-like.)
// You can take an already-calculated fractal and apply
// a color map to it through methods provided.  Just pass
// it a color map object that you have created.

// Here's the copyright from the original UNIX program "lyap",
// from which I derived much of the engine.
// The original program and an executable and output viewer for
// the NeXT are provided in a separate directory in this distribution.
// As the notice states, you cannot redistribute this program, or
// anything derived from it without including all sources, including
// your own changes!  

/*
 * by Garrett Wollman based on a function posted to Usenet by Ed Kubaitis
 * This program is Copyright (C) 1991, Garrett A. Wollman.  This
 * program may be modified and distributed for any purpose and without
 * fee, provided that this notice remains on all such copies
 * unaltered.  Binary distributions not including this source file are
 * prohibited.  Modified versions must be marked with the name of the
 * modifier and the date of modification.
 *
 * Under no circumstances shall Garrett A. Wollman or the University
 * of Vermont and State Agricultural College be held liable for any
 * damages resulting from the use or misuse of this program, whether
 * the author is aware of such possibility or not.  This program is
 * warranted solely to occupy disk space.
 *
 * I'm sorry for all this legal garbage, but it is sadly necessary
 * these days.
 */

#include <math.h>
#include <time.h>
#include <stdlib.h>
#include <stdio.h>


/*
 * Global variables
 */

char *whoami;

extern int getopt(int,char **,const char *);
extern char *optarg;
extern int optind;
extern int opterr;
extern void perror(const char *);
extern int strlen(const char *);

#define OPTIONS "r:c:v:m:d:s:g:5pto:l:"

int rows = 512,cols = 512;
double amin,amax;
double bmin,bmax;
double aincr,curA;
double bincr;

#define nColors 256		/* don't even try to change this */
struct {
  unsigned char red;
  unsigned char green;
  unsigned char blue;
} colormap[nColors];

int Settle = 600;
int Dwell = 2000;
int *Rvec;
double minLya = -5;
int showpos = 0;
double initGuess = 0.5;

/*
 * Colormap functions
 */

/*
 * set up our initial shades of grey
 */
void init_colormap(void) {
  int i;
  for(i=0; i < nColors; i++)
    colormap[i].red = colormap[i].green = colormap[i].blue = i;
}

void read_colormap(const char *name) {
  FILE *cmap;
  char buf[255];		/* magic number */
  int a,b,c,i;
  
  cmap = fopen(name,"r");
  if(!cmap) {
    perror(name);
    exit(-1);
  }

  i = 0;
  while(!feof(cmap) && i < nColors) {
    fgets(buf,sizeof buf,cmap);
    /*
     * unparseable lines are comments
     * as is anything after the <r g b> value on a single line
     *
     * per Fractint 15.1 standard
     */
    if(sscanf(buf,"%d %d %d",&a,&b,&c) != 3)
      continue;

    colormap[i].red = a;
    colormap[i].green = b;
    colormap[i].blue = c;
    i++;
  }
  fclose(cmap);
}


/*
 * Lyapunov exponent and coloring calculation
 * original function by Ed Kubaitis, used by permission
 */

/*
 * Speedup...
 *
 * Erick B. Wong noticed that $\log_2 d(x_1)+\dots +\log_2 d(x_n)$
 * is, by a principle from high-school algebra, the same as
 * $\log_2 (d(x_1) + \dots + d(x_n))$.  We take advantage of this
 * by unrolling the dwell loop by four (hence the rounding in main())
 * and accumulating by multiplication before taking the log.  Since
 * log() is extremely expensive, this should lead to a considerable
 * speedup, while still allowing complete flexibility in selecting
 * dwell values.
 *
 * Presumably, one might unroll this loop even further, with a 
 * concomitant increase in textual complexity.
 */

int lya(double a,double b) {
  double t = 0;
  double r, lya;
  int d;
  double acc;
  double x;

  x = initGuess;		/* initialize x */

  for(d=0; d < Settle; d++) {
    r = (Rvec[d]) ? b : a;
    x = r*x*(1-x);
  }

  for(d=0; d <= Dwell; d+= 4) {
    r = (Rvec[d]) ? b : a;
    x = r * x * (1-x);
    acc = fabs(r - 2*r*x);

    r = (Rvec[d+1]) ? b : a;
    x = r * x * (1-x);
    acc *= fabs(r - 2*r*x);

    r = (Rvec[d+2]) ? b : a;
    x = r * x * (1-x);
    acc *= fabs(r - 2*r*x);

    r = (Rvec[d+3]) ? b : a;
    x = r * x * (1-x);
    t += log(acc * fabs(r - 2*r*x));

    if(fabs(x-0.5) < 1e-20) {
      d += 3;
      break;
    }
  }

  lya = (t * 1.442695041)/d;

  if(showpos)			
    lya *= -1;
  if(lya < minLya)		/* sanity check! */
    lya = minLya;

  if(lya < 0) {
    return (int)(lya / minLya * (double)nColors);
  } else {
    return 0;
  }
}


/*
 * Usage message
 */
#ifdef __GNUC__
volatile
#endif
void usage(void) {
  fprintf(stderr,
	  "%s: usage:\n\n"
	  "%s [-r rows] [-c cols] [-v program] [-p] [-t]\n\t"
	  "[-o file] [-l colormap]\n\t"
	  "[-m minLya] [-d Dwell] [-s Settle] [-g initGuess]\n\t"
	  "[-5] amin amax bmin bmax bits\n",whoami,whoami);
  exit(-1);
}

int main(int argc,char **argv) {
  int c,i,j;
  char *outname = NULL;
  FILE *outfile;

  char *bits;

  double curB;
  int onRow,onCol;
  int *oneRow;

  int outtext = 0;

  init_colormap();
  opterr = 0;

  whoami = *argv;
  while((c = getopt(argc,argv,OPTIONS)) != EOF) {
    switch(c) {
    case 'r':
      rows = atoi(optarg);
      if(!rows)
	    rows = 512;
      break;

    case 'c':
      cols = atoi(optarg);
      if(!cols)
	    cols = 512;
      break;

    case 'o':
      outname = optarg;
      break;

    case 'l':
      read_colormap(optarg);
      break;

    case 'm':
      if(!sscanf(optarg,"%lf",&minLya))
	usage();
      break;

    case 'd':
      Dwell = atoi(optarg);
      if(!Dwell)
	usage();
      /*
       * By forcing the Dwell to be a multiple of four, we can use the
       * speedup noted in lya().  We *could* just do this silently,
       * by if it might make a difference to someone...
       */
      if (Dwell % 4) {
	    fprintf(stderr,"Dwell value of %d being rounded off to %d.\n",
		  Dwell, Dwell += 4 - (Dwell % 4));
      }
      break;

    case 's':
      Settle = atoi(optarg);
      if(!Settle)		/* for settle close to 0, use 1 */
	    usage();
      break;

    case 'g':
      if(!sscanf(optarg,"%lf",&initGuess))
	    usage();
      break;

    case 'p':
      showpos = 1;
      break;

    case 't':
      outtext = 1;
      break;

    default:			/* includes '?' (and '#' when !MULTIPROC) */
      usage();
    }
  }

  if(!argv[optind])
    usage();

  /*
   * Note to language nit-pickers... The code below is legal because
   * the || operator is a sequence point.
   */
  if(!sscanf(argv[optind++],"%lf",&amin) || !argv[optind] ||
     !sscanf(argv[optind++],"%lf",&amax) || !argv[optind] ||
     !sscanf(argv[optind++],"%lf",&bmin) || !argv[optind] ||
     !sscanf(argv[optind++],"%lf",&bmax) || !argv[optind])
    usage();

  bits = argv[optind];

  if(Dwell < Settle)
    Dwell = Settle;		/* a bit of sanity */

  j = strlen(bits);

  Rvec = (int *)malloc((Dwell+1) * sizeof *Rvec);
  if(!Rvec) {
    fputs("Cannot allocate space for Markus vector!\n",stderr);
    exit(-1);
  }

  for(i=0; i <= Dwell; i++) {
    if(bits[i % j] == 'a')
      Rvec[i] = 0;
    else if(bits[i % j] == 'b')
      Rvec[i] = 1;
    else {
      fprintf(stderr,"Invalid character \\%o ('%c') in bit string\n",
	      (int)(unsigned char)bits[i % j], bits[i % j]);
      exit(-1);
    }
  }

  oneRow = (int *)malloc(cols * sizeof *oneRow);
  if(!oneRow) {
    fputs("Cannot allocate space for single-row buffer!\n",stderr);
    exit(-1);
  }

  aincr = (amax - amin) / rows;
  bincr = (bmax - bmin) / cols;

  if(outname) {
    outfile = fopen(outname,"w");
    if(!outfile) {
      perror(outname);
      exit(-1);
    }
  } else {
    outfile = stdout;
  }

  /*
   * Print output header... same as DKBTrace
   */
  if (outtext) {
  	fprintf(outfile, "P3 %d %d", cols, rows);
  } else {
  putc('\0',outfile);  putc('\0',outfile);
  putc('\0',outfile);  putc('\0',outfile);
  putc('\0',outfile);  putc('\0',outfile);
  putc('\0',outfile);  putc('\0',outfile);
  putc('\0',outfile);  putc('\0',outfile);
  putc('\0',outfile);  putc('\0',outfile);
  fprintf(stderr, "Sizes: %d %d %d %d\n",cols%256,cols/256,rows%256,rows/256);
  putc((unsigned char)(cols%256),outfile);
  putc((unsigned char)(cols/256),outfile);
  putc((unsigned char)(rows%256),outfile);
  putc((unsigned char)(rows/256),outfile);
  putc('\0',outfile);  putc('\0',outfile);
  }


  for(onRow = 0, curA = amax; onRow < rows; curA -= aincr, onRow++) {
    for(onCol = cols - 1, curB = bmax; onCol; curB -= bincr, onCol--) {
      oneRow[onCol] = lya(curA,curB);
    }
    
    for(i=0; i < cols; i++) {
      unsigned char c;
      c = (unsigned char)oneRow[i];
      if(outtext) {
		  fprintf(outfile,"%d %d %d ",colormap[c].red,colormap[c].green,
		  colormap[c].blue);
      } else {
	  putc(colormap[c].red,outfile);
	  putc(colormap[c].green,outfile);
	  putc(colormap[c].blue,outfile);
      }
    }
    if(outtext)
	  fputc('\n',outfile);
    if(!(onRow % 16)) {
      fprintf(stderr,"Done row %d/%d \n",onRow,rows);
    }
  }

  if(outname)
    fclose(outfile);

  return 0;
}
