// Note: you can't yet render only parts of the fractal.
// Eventually the engine needs to be moved to a separate
// object from this view so that several engines -- each
// with it's own thread -- can be started up and all feed
// into this view object.

// Adapted from lyap.c by Don Yacktman
// This object embodies the lyapunov calculation code.
// Use the methods provided to set the various options.
// 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.
 */

/* Generated by Interface Builder */

#import "LyapunovView.h"
#import "ColorMapView.h"
#import "ColorMap.h"
#import <libc.h>
#import <stdlib.h>
#import <stdio.h>
#import <strings.h>
#import <appkit/tiff.h>
#import <appkit/color.h>
#import <appkit/graphics.h>
#import <appkit/Control.h>
#import <appkit/SavePanel.h>
#import <appkit/OpenPanel.h>
#import <appkit/Matrix.h>
#import <appkit/NXBitmapImageRep.h>
#import <appkit/NXColorWell.h>
#import <dpsclient/psops.h>
#import <appkit/Application.h>
#import <dpsclient/wraps.h>
#import <sys/file.h>
#import <streams/streams.h>
#import <sys/types.h>
#import <sys/uio.h>

#define FABS(k) (((k) < 0) ? (-(k)) : (k))	// faster than 2.1 trap

double log(x)	// to replace trap in 2.x OS that slows down '040 machines
double x;		// taken verbatim from libjv.a
{
/*
	Routine for log(x) in IEEE double precision.
	Note that we have to hack around in the exponent!
	Steps: 1 x = y*2^n, with y between sqrt(2) and 1/sqrt(2)
	       2 convert to z = (y-1)/(y+1)
	       3 ln = 2*z*sum(i=0,9, z^(2*i)/(2*i+1)) + n*ln(2)
	Routine made by J.A.M. Vermaseren 20-apr-1992.
*/
	double scr, y, z, z2, sum, sqrt();
	short n;
	if ( x < 0 ) {
		printf("log of a negative number: %d\n",x);
		return (-1);
	}
	scr = x;
	n = ( (*((short *)(&scr))) & 0x7FF0 ) - 0x3FF0;
	*((short *)(&scr)) -= n;
	n >>= 4;
	if ( scr > 1.4142135623730950488 ) {
		*((short *)(&scr)) -= 0x10;
		n++;
	}
	y = scr;
	z = (y-1.)/(y+1.);
	z2 = z*z;
	sum = 2.*z*(0.99999999999999999886+z2*(
                0.33333333333333825805+z2*(
                0.19999999999649522967+z2*(
                0.14285714380662481221+z2*(
                0.11111098494595009814+z2*(
                0.09091817553294820737+z2*(
                0.07656220982777783836+z2*(
                0.07405280893003482131))))))))
		 +n*0.6931471805599453094;
	return(sum);
}

@implementation LyapunovView

- initFrame:(const NXRect *)frm		// designated initializer for a view
{
    [super initFrame:frm];
    xPos = 0.0; yPos = 0.0;
    xScale = 4.0; yScale = 4.0;
    initial = 0.5; minExp = -5.0;
    dwell = 2000;  settle = 600;
	firstInit = YES; positive = NO;
    return self;
}

- appDidInit:sender				// final initialization
{        	// set up default parameters on control panel
    activeMap = [activeMapView colorMap];
	
	// set up default position
    [[posMatrix findCellWithTag:0] setFloatValue:xPos];
    [[posMatrix findCellWithTag:1] setFloatValue:yPos];
    
    // set up default scale
    [[scaleMatrix findCellWithTag:0] setFloatValue:xScale];
    [[scaleMatrix findCellWithTag:1] setFloatValue:yScale];
    
    // set up upper r.h. corner
    [[sizeMatrix findCellWithTag:0] setFloatValue:(xPos + xScale)];
    [[sizeMatrix findCellWithTag:1] setFloatValue:(yPos + yScale)];
    [[windowSizeMatrix findCellWithTag:0] setFloatValue:(XSIZE)];
    [[windowSizeMatrix findCellWithTag:1] setFloatValue:(YSIZE)];
    
    // set up default initial value
    [[boundsMatrix findCellWithTag:1] setFloatValue:initial];
    [[boundsMatrix findCellWithTag:0] setFloatValue:minExp];
    
    // set up default depth
    [[depthMatrix findCellWithTag:0] setIntValue:dwell];
    [[depthMatrix findCellWithTag:1] setIntValue:settle];
    
	// set up polarity
	[polarityMatrix selectCellWithTag:positive];

    // set up default pattern
    if (firstInit) [[pattern findCellWithTag:0] setStringValue:"aaab"];
    firstInit = NO;
	[saveMenuButton setEnabled:NO];
	[reApplyButton setEnabled:NO];
    return self;
}



// 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 = initial;		/* initialize x */

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

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

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

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

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

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

	lya = (t * 1.442695041)/d;

	if (positive) lya *= -1;
	if (lya < minExp) lya = minExp;

	if (lya < 0) {
		return (int)(lya / minExp * (double)NUMCOLORS);
	} else {
		return 0;
	}
}

- go:sender				// begin calculation of image
{
    register int pix, y, x;
	double aa, bb;
    int r, g, b;
	unsigned char c;
	NXRect pixel, line;
    
    pixel.size.width = 1.0;
    pixel.size.height = 1.0;
    line.size.height = 1.0;
    line.size.width = XSIZE;
    line.origin.x = 0.0;
    
    [self newParam:self];	// be sure we're up to date
    
	// main calc loop: calc line & draw it.
	[self lockFocus];
	for (y=0; y<YSIZE; y++) {
		bb = yPos + y*yScale/YSIZE; //  get physical y-coord
		for (x=0; x<XSIZE; x++) {
			aa = xPos + x*xScale/XSIZE; //  get physical x-coord
			space[x][y] = [self lya:aa :bb]; // Lyapunov exponent approx.
			// put into pixel map w/color	    
			pix = (x+(YSIZE-y-1)*XSIZE)*RGB;
			c = (unsigned char)space[x][y];
			[activeMap colorFor:c :&r :&g :&b];
			pixels[pix] = r;  // pix+RED = pix because RED=0
			pixels[pix+GREEN] = g;
			pixels[pix+BLUE] = b;
		}
		// plot a line of it
		line.origin.y = y;
		NXImageBitmap(&line, XSIZE, 1, 8, 3,
			NX_MESHED, NX_COLORMASK,
			&pixels[(int)((YSIZE-y-1)*XSIZE*RGB+0.5)],
		NULL, NULL, NULL, NULL);
		[window flushWindow];
		NXPing();
	}
	[self unlockFocus];
	[saveMenuButton setEnabled:YES];
	[reApplyButton setEnabled:YES];
	return self;
}

- reApply:sender		// re do coloring on plot
{
    register int x, y, pix;
	register unsigned char c;
	int r, g, b;
    NXRect pixel, line;
    
    [self newParam:self];	// be sure we're up to date
    // set up plotting bounds
    pixel.size.width = 1; pixel.size.height = 1;
    line.size.height = 1; line.size.width = XSIZE;
    line.origin.x = 0;
    // re-apply colors
    [self lockFocus];
    for (y=0; y<YSIZE; y++) {
		for (x=0; x<XSIZE; x++) {
			pix = (x+(YSIZE-y-1)*XSIZE)*RGB;
			c = (unsigned char)space[x][y];
			[activeMap colorFor:c :&r :&g :&b];
			pixels[pix] = r;  // pix+RED = pix because RED=0
			pixels[pix+GREEN] = g;
			pixels[pix+BLUE] = b;
		}
		// plot a line of it
		line.origin.y = y;
		NXImageBitmap(&line, XSIZE, 1, 8, 3,
			NX_MESHED, NX_COLORMASK,
			&pixels[(int)((YSIZE-y-1)*XSIZE*RGB+0.5)],
			NULL, NULL, NULL, NULL);
		[window flushWindow];
		NXPing();
	} 
	[self unlockFocus];   
	[saveMenuButton setEnabled:YES];
	[reApplyButton setEnabled:YES];
	return self;
}

- newParam:sender			// set up new parameters for calc
{
    int count, len;
    
    // get x-y coordinates
    xPos = [[posMatrix findCellWithTag:0] floatValue];
    yPos = [[posMatrix findCellWithTag:1] floatValue];
    
    // get x-y scaling
    xScale = [[scaleMatrix findCellWithTag:0] floatValue];
    yScale = [[scaleMatrix findCellWithTag:1] floatValue];
    
    // get initial value for each point
    minExp  = [[boundsMatrix findCellWithTag:0] floatValue];
    initial = [[boundsMatrix findCellWithTag:1] floatValue];
    
    // get depth
    settle = [[depthMatrix findCellWithTag:1] intValue];
    dwell = [[depthMatrix findCellWithTag:0] intValue];
	if (dwell % UNRAVEL) dwell += UNRAVEL - (dwell % UNRAVEL); // round off
	if (dwell < settle) {	// sanity check
		if (settle % UNRAVEL) settle += UNRAVEL - (settle % UNRAVEL);
		dwell = settle;
	}
    
    // get pattern
    patternString = [[pattern findCellWithTag:0] stringValue];
	len = strlen(patternString);
	if (patternNum) free(patternNum);
	patternNum = (int *)malloc((dwell + 1) * sizeof (int));
	if(!patternNum) {
		fprintf(stderr, "Cannot allocate space for Markus vector!\n");
	}
    for(count=0; count <= dwell; count++) {
        patternNum[count] = ((patternString[count % len] == 'a') ||
			(patternString[count % len] == 'x')) ? 0 : 1;
    }
    patternLength = count;
	
	positive = [[polarityMatrix selectedCell] tag];
	
	[self appDidInit:self];
    return self;
}

- drawSelf:(NXRect *)rects :(int)rectCount	// redraws the screen.
{
    PSsetgray(NX_DKGRAY);
    NXRectFill(&bounds);
    return self;
}

- mouseDown:(NXEvent *)e	// stolen from MandelView.m
/*
 * We implement the mouseDown method so the user can sweep out a section of
 * the view and select that as the current x,y,scale coordinates.
 */
{
  int looping = YES;
  int oldMask;
  NXRect bbox;
  NXPoint startPoint, currPoint;
  double newX, newY, newDX, newDY;
  NXCoord aspectRatio;
  float mvX = xPos;
  float mvY = yPos;
  float mvDX = xScale;
  float mvDY = yScale;
  
  aspectRatio = bounds.size.height / bounds.size.width;

  oldMask =  [window addToEventMask:NX_MOUSEDRAGGEDMASK];
  startPoint = e->location;
  [self convertPoint:&startPoint fromView:nil];
  NXSetRect(&bbox,startPoint.x,startPoint.y,0.0,0.0);
  [self lockFocus];
  while (looping) {
    e=[NXApp getNextEvent:NX_MOUSEUPMASK | NX_MOUSEDRAGGEDMASK];
    currPoint = e->location;
    [self convertPoint:&currPoint fromView:nil];
    bbox.size.width = 2*(currPoint.x - startPoint.x);
    bbox.size.height = 2*(currPoint.y - startPoint.y);
    /* Normalize bbox to always have positive width and height */
    if (bbox.size.width < 0) bbox.size.width = -bbox.size.width;
    if (bbox.size.height < 0) bbox.size.height = -bbox.size.height;
    /*
     * constrain the box to have the aspect ratio of the view.  Choose
     * whichever dimension is closer to the desired result.
     */
    if ((bbox.size.height/bbox.size.width) > aspectRatio) {
      bbox.size.height = bbox.size.width * aspectRatio;
    } else {
      bbox.size.width = bbox.size.height / aspectRatio;
    }
    /* The startPoint is always at the center of the bbox */
    bbox.origin.x = startPoint.x - .5 * bbox.size.width;
    bbox.origin.y = startPoint.y - .5 * bbox.size.height;
    PSnewinstance();
    if (looping = (e->type == NX_MOUSEDRAGGED)) {
      PSsetinstance(YES);
      PSsetgray(NX_WHITE);
      NXFrameRect(&bbox);
      PSsetinstance(NO);
    }
  }
  [self unlockFocus];
  [window setEventMask:oldMask];
  if ((bbox.size.width > 0) && (bbox.size.height > 0)) {
    /*
     * At this point, bbox is in window coordinates.  Set new controller
     * parameters based on this view's coordinates and the bounding box.
     */
    newDX = (bbox.size.width*mvDX)/bounds.size.width;
    newDY = (bbox.size.height*mvDY)/bounds.size.height;
    newX = ((startPoint.x + bounds.origin.x) / bounds.size.width);
    newX = newX * mvDX + mvX - newDX/2;
    newY = ((startPoint.y + bounds.origin.y) / bounds.size.height);
    newY = newY * mvDY + mvY - newDY/2;

    /*
     * Note that we only update the text fields -- we
     * don't update the internal instance variables.  If we were to update
     * the internal ivs, then the user would only get one chance with the
     * mouse down method because subsequent mouse-downs would work in terms
     * of the new coordinate system.
     */
    [[posMatrix findCellWithTag:0] setFloatValue:newX];
    [[posMatrix findCellWithTag:1] setFloatValue:newY];
    [[scaleMatrix findCellWithTag:0] setFloatValue:newDX];
    [[scaleMatrix findCellWithTag:1] setFloatValue:newDY];
  }
  return self;
}

- saveAsTIFF:sender;			// save image in a .tiff file
{
    if ([[SavePanel new] runModalForDirectory:[[OpenPanel new] directory]
    		file:"UNTITLED.tiff"]) {
        [self saveTIFF:[[SavePanel new] filename]];
    }
    return self;
}

- (BOOL)saveTIFF:(const char *)name
{
	id imageRep;
	int fd;
	NXStream *stream;

	imageRep = [[NXBitmapImageRep alloc] initData:pixels
		pixelsWide:(int)(XSIZE) pixelsHigh:(int)(YSIZE)
		bitsPerSample:8 samplesPerPixel:3
		hasAlpha:NO isPlanar:NO colorSpace:NX_RGBColorSpace
		bytesPerRow:(int)(XSIZE*3) bitsPerPixel:0];
	if (imageRep == nil) {
		fprintf(stderr, "Couldn't create ImageRep.");
		return NO;
	}
     
    fd = open(name, O_CREAT | O_WRONLY | O_TRUNC, 0666);
    if (!fd) { printf("Error opening save file %s.\n", name); }
    stream = NXOpenFile(fd, NX_WRITEONLY);
    if (stream == NULL) {
		fprintf(stderr, "Error opening save stream.\n");
		return NO;
	}
    [imageRep writeTIFF:stream usingCompression:NX_TIFF_COMPRESSION_LZW];
    NXFlush(stream);
    NXClose(stream);
    close(fd);
    return YES;
}

- windowWillResize:sender toSize:(NXSize *)frameSize	// watch our window
{
    if (frameSize->width > (MAXXSIZE + 2)) frameSize->width = MAXXSIZE + 2;
    if (frameSize->height > (MAXYSIZE + 32)) frameSize->height = MAXYSIZE + 32;
    if (frameSize->width < 66) frameSize->width = 66;
    if (frameSize->height < 96) frameSize->height = 96;
	return self;
}

- windowDidResize:sender								// watch our window
{
	[saveMenuButton setEnabled:NO];
	[reApplyButton setEnabled:NO];
	return self;
}



@end
