/* ContourView.m
   Contour Plot object with optional color fills.
 
     It is pretty easy to use.  If instantiated within IB, all you need to do
     to get a plot is to call the following two method to get a default
     plot.  If you understand the following method, you can use the
     ContourView object easily.

     - setCartesianGridData:(float *)f :(float)xmin :(float)xmax
				  :(float)ymin :(float)ymax
				  ofSize:(int)nx :(int)ny
				  withInterpolationTo:(int)n1x :(int)n1y;
 
    f[nx*ny] is a 1-d array containing 2-d grid data such that f[iy*nx+ix]
    contains the value at (ix, iy).
    Typically, just a few messages below will produce a contour plot with color
    fills.

    [contourView setBipolar:YES];
    [contourView setCartesianGridData: fdata :1.0 :5.0 :1.0 :5.0
				ofSize: 20 :20
				withInterpolationTo: 50 :50];
    [contourView setFillEnable:YES];
    [contourView setEnabled:YES];
    [contourView updateImageCache];
    [contourView display];
 
 * ---------------------------------------------------------------------------
 * TO DO:
 *  Optimize drawing.  Especially, do not recompute everything in drawSelf::.
 *
 * ---------------------------------------------------------------------------
 * Version History:
 * V0.99 95-01-27 Izumi Ohzawa
 *    Implemented cross-hair to provide cross-section 1-D profile.
 *    Now starts out disabled so method call [contourView setEnabled:YES];
 *    must be done before it does anything useful (see above).
 *    Also, it now caches the plot for cross-hair, so -updateImageCache
 *    must be called when the plot changes.
 * V0.98 93-10-03 Izumi Ohzawa
 *    Tailor-resistant comment insertion.
 *    Also, comments moved to top of EPS file.
 * V0.97 93-07-10 Izumi Ohzawa
 *    Odd number of contours allowed.
 *    Bug with coloring tick, axis, frame colors fixed.
 * V0.96 93-02-19 Izumi Ohzawa
 *    Dash control for negative contour and grid mesh.
 *    Control for doing or not doing contour path stroke
 * V0.95 93-01-21 Izumi Ohzawa
 *    Added tick marks.
 * V0.94 93-01-20 Izumi Ohzawa
 *    More info such as grid size, etc have been added as comments.
 *    Finally, a bug with splin2.c from Numerical Recipes has been found and
 *    fixed.  Axis and grid methods added.
 * V0.93 92-12-04 Izumi
 *    Small fixes.
 * V0.92 92-12-03 Izumi
 *    It now correctly does color fills via sorting of contour drawing order.
 * V0.91 92-05-24 Izumi
 *    Does OK color fills, printing, copying PS to paste board.
 * V0.90 92-05-18 Izumi Ohzawa, izumi@pinoko.berkeley.edu
 *    The initial version.
*/

/* Version of ContourView object.  This is entered into PS code generated after NeXT
   print package prolog.
 */
#define CV_VERSION "% ContourView.m [V0.99 95-01-24] by Izumi Ohzawa, izumi@pinoko.berkeley.edu."

#import <appkit/appkit.h>

#import "ContourView.h"
#import "contour.h"
#import "matalloc.h"
#import "splin2.h"
#import <strings.h>
#import <stdlib.h>		/* malloc */
#import <math.h>
#import <time.h>

#define N6			6	/* expand 3 points on 4 sides to close all contours */
#define N3			3
#define N2			2
#define DEFAULTNCLEVELS		16	/* Default number of contour levels */


// static float ncontdash[] = {5.0, 5.0 };	/* dash for negative contours */
static float psolid[] = {0.0 };
// static float griddash[] = {1.0, 2.0 };	/* dash for grid */


@implementation ContourView

- initFrame:(NXRect *)nf
{
    self = [super initFrame:nf];
//    ClearFlag = 0;
    enabled = NO;
    invertY = NO;
    doFill = NO;
    doContourStroke = YES;		// stroke contour path
    debug = NO;
    frameON = YES;
    axisON = NO;
    gridON = NO;
    ticksON = NO;
    frameLineWidth = 1.0;
    axisLineWidth = 0.5;
    gridLineWidth = 0.2;
    tickLineWidth = 0.8;
    tickLength = 6.0;
    axisX = axisY = 0.0;		// where to draw axis
    gridXstep = gridYstep = 1.0;
    minNumPoints = 3;
    ndata = 0;
    xd = NULL;
    yd = NULL;
    fd = NULL;
    f2 = NULL;
    fd2a = NULL;
    xt = NULL;
    yt = NULL;
    SortedCntrPtr = NULL;
    comment = NULL;
    comment012 = NULL;
    bipolar = 0;
    basevalue = 0.0;
    smallx = smally = 0.0001;
    pXmin = bounds.origin.x;
    pYmin = bounds.origin.y+0.25;
    pXmax = bounds.size.width + bounds.origin.x -0.25;
    pYmax = bounds.size.height + bounds.origin.y -0.25;
    pX_Vline = pXmin - 10.0;
    pY_Hline = pYmin - 10.0;
    lastx = lasty = -10.0;

    // Allocate cache
    imageCache = [[NXImage alloc] initSize:&bounds.size];
    [imageCache useCacheWithDepth:NX_DefaultDepth];
    cacheOK = 0;
    printCrossHair = NO;
    continuousNotify = NO;

    nclevels = DEFAULTNCLEVELS;
    cA = (CntrAttribute *)malloc((size_t)nclevels*sizeof(CntrAttribute));
    contourList = NULL;
    numContours = 0;

    // default attributes of contours
    positiveColor = NX_COLORRED;
    negativeColor = NX_COLORBLUE;
    backgroundColor = NX_COLORWHITE;
    contourLineColor = NX_COLORBLACK;
    frameColor = NX_COLORBLACK;
    axisColor = NX_COLORBLACK;
    gridColor = NX_COLORBLACK;
    tickColor = NX_COLORBLACK;
    minorContourLineWidth = 1.0;
    majorContourLineWidth = 2.0;
    griddash[0] = 1.0;		// initial grid dash (mark)
    griddash[1] = 2.0;		// (space)
    Nncdash = 2;		// Negative contour dash array size
    ncontdash = (float *)malloc((Nncdash+1)*sizeof(float));
    ncontdash[0] = 6.0;		// (mark)
    ncontdash[1] = 6.0;		// (space)

    [self setDefaultContourAttributes];
    return self;
}


- free
{
    if(imageCache) [imageCache free];
    if(f2) free_fmatrix(f2, 1, nxraw, 1, nyraw);
    if(fd2a) free_fmatrix(fd2a, 1, nxraw, 1, nyraw);
    if(xt) free((void *)xt);
    if(yt) free((void *)yt);
    if(cA) free(cA);
    if(xd) free(xd);
    if(yd) free(yd);
    if(fd) free(fd);
    [super free];
    return self;
}

- (BOOL)acceptsFirstMouse { return YES; }
- (BOOL)acceptsFirstResponder { return YES; }

// ========= DRAWing code now draws into cache
- drawAll
{
int j, OK2stroke = 0;
float xxa, yya;
float xg, yg, xxg, yyg;
ContourPath *cntr;
DPSContext curContext = DPSGetCurrentContext();

    NXSetColor(backgroundColor);
    NXRectFill(&bounds);
    if(!fd) return self;	// if no data set, do nothing.

    if(contourList) [self freeContourList];		/* to start out */
    numContours = 0;

    computeContour(self, nx, ny, xd, yd, fd, nclevels, cA);

    if(debug)
	fprintf(stderr, "ContourView: # levels=%d,  # contours=%d\n", nclevels, numContours);

    if(doFill) {
	/* Sort the order of contour drawing from outermost to innermost for filling */
        SortedCntrPtr = (ContourPath **)(malloc((size_t)(sizeof(ContourPath *)*numContours)));
        sort_contourList(contourList, xdmin, xdmax, ydmin, ydmax, numContours, SortedCntrPtr);
	for(j=0; j<numContours; j++) {
	    cntr = SortedCntrPtr[j];
	    [self findInsideHighLow:cntr];		/* inside contour goin up or down */
	    [self plotContour:cntr :curContext];	/* This plots one contour curve */
	}
    }
    else {
	/* Do NOT do filling -- just go through the list */
	cntr = contourList;
	while(cntr) {
	    [self plotContour:cntr :curContext];	/* This plots one contour curve */
	    cntr = cntr->next;		/* point to the next one */
	};
    }

    PSsetdash(psolid, 0, 0.0);

    /* ==== Do GRID MESH ======================================== */
    if(gridON) {
        OK2stroke = 0;
        if(NXDrawingStatus != NX_DRAWING)
            DPSPrintf(curContext, "\n%% Grid Mesh\n");
	PSsetdash(griddash, 2, 0.0);
	NXSetColor(gridColor);
	PSsetlinewidth(gridLineWidth);
        for(xg=axisX; xg < xdmax; xg += gridXstep) {
	    xxg = (xg - xdmin)*ppxu;
	    PSmoveto(xxg, pYmin);
	    PSlineto(xxg, pYmax);
	    OK2stroke++;
	}
        for(xg=axisX; xg > xdmin; xg -= gridXstep) {
	    xxg = (xg - xdmin)*ppxu;
	    PSmoveto(xxg, pYmin);
	    PSlineto(xxg, pYmax);
	    OK2stroke++;
	}

        for(yg=axisY; yg < ydmax; yg += gridYstep) {
        if(invertY)
	    yyg = (ydmax - yg)*ppyu;
	else
	    yyg = (yg - ydmin)*ppyu;
	    PSmoveto(pXmin, yyg);
	    PSlineto(pXmax, yyg);
	    OK2stroke++;
	}
        for(yg=axisY; yg > ydmin; yg -= gridYstep) {
        if(invertY)
	    yyg = (ydmax - yg)*ppyu;
	else
	    yyg = (yg - ydmin)*ppyu;
	    PSmoveto(pXmin, yyg);
	    PSlineto(pXmax, yyg);
	    OK2stroke++;
	}
        if(OK2stroke) PSstroke();
	PSsetdash(psolid, 0, 0.0);
    }  /* end of if(gridON) {} */


    /* ==== Do TICK marks ======================================== */
    if(ticksON) {
        OK2stroke = 0;
        if(NXDrawingStatus != NX_DRAWING)
            DPSPrintf(curContext, "\n%% Tick Marks\n");
	NXSetColor(tickColor);
	PSsetlinewidth(tickLineWidth);
        for(xg=axisX; xg < xdmax; xg += gridXstep) {
	    xxg = (xg - xdmin)*ppxu;
	    PSmoveto(xxg, pYmin); PSlineto(xxg, pYmin+tickLength);
	    PSmoveto(xxg, pYmax); PSlineto(xxg, pYmax-tickLength);
	    OK2stroke++;
	}
        for(xg=axisX; xg > xdmin; xg -= gridXstep) {
	    xxg = (xg - xdmin)*ppxu;
	    PSmoveto(xxg, pYmin); PSlineto(xxg, pYmin+tickLength);
	    PSmoveto(xxg, pYmax); PSlineto(xxg, pYmax-tickLength);
	    OK2stroke++;
	}

        for(yg=axisY; yg < ydmax; yg += gridYstep) {
        if(invertY)
	    yyg = (ydmax - yg)*ppyu;
	else
	    yyg = (yg - ydmin)*ppyu;
	    PSmoveto(pXmin, yyg); PSlineto(pXmin+tickLength, yyg);
	    PSmoveto(pXmax, yyg); PSlineto(pXmax-tickLength, yyg);
	    OK2stroke++;
	}
        for(yg=axisY; yg > ydmin; yg -= gridYstep) {
        if(invertY)
	    yyg = (ydmax - yg)*ppyu;
	else
	    yyg = (yg - ydmin)*ppyu;
	    PSmoveto(pXmin, yyg); PSlineto(pXmin+tickLength, yyg);
	    PSmoveto(pXmax, yyg); PSlineto(pXmax-tickLength, yyg);
	    OK2stroke++;
	}
        if(OK2stroke) PSstroke();
    }  /* end of if(ticksON) {} */

    /* ==== Do AXIS lines ======================================== */
    if(axisON) {
        if(NXDrawingStatus != NX_DRAWING)
            DPSPrintf(curContext, "\n%% Axis\n");
	NXSetColor(axisColor);
	PSsetlinewidth(axisLineWidth);
	xxa = (axisX - xdmin)*ppxu;
        if(invertY)
	    yya = (ydmax - axisY)*ppyu;
	else
	    yya = (axisY - ydmin)*ppyu;
	if(xxa > pXmin && xxa < pXmax) {
	    PSmoveto(xxa, pYmin);
	    PSlineto(xxa, pYmax);
	    PSstroke();
	}
	if(yya > pYmin && yya < pYmax) {
	    PSmoveto(pXmin, yya);
	    PSlineto(pXmax, yya);
	    PSstroke();
	}
    } /* end of if(axisON) {} */

    /* ==== Do FRAME RECTANGLE ======================================== */
    if(frameON) {
        if(NXDrawingStatus != NX_DRAWING)
            DPSPrintf(curContext, "\n%% Frame Rectangle\n");
	// Frame rectangle
	NXSetColor(frameColor);
	PSsetlinewidth(frameLineWidth);
	PSrectstroke(pXmin,pYmin,pXmax,pYmax);	/* frame rectangle */
    }

    if(doFill && SortedCntrPtr) {
	free(SortedCntrPtr);
	SortedCntrPtr = NULL;
    }
    [self freeContourList];

    return self;
}

// This draws the plot into NXImage cache
//
- updateImageCache
{
    if([imageCache lockFocus]) {
	[self drawAll];
	[imageCache unlockFocus];
    } /* end of if([imageCache lockFocus])... */
    cacheOK = 1;
    return self;
}


// ==== MOUSE EVENT, CURSOR, MARKER HANDLING METHODS =========================================


- enablePrintCrossHair:(BOOL)yn
{
    printCrossHair = yn;
    return self;
}

- enableContinuousNotify:(BOOL)yn
{
    continuousNotify = yn;
    return self;
}

/*  This method handles a mouse down.  */
- mouseDown:(NXEvent *)event
{
NXPoint	  p;
int mknum;
    if(!enabled || !cacheOK) return self;
    [window makeFirstResponder:self];		/* give me the first responder status */
    p = event->location;
    [self convertPoint:&p fromView:nil];	/* to this view's coordinate from window's */
    if ((mknum = [self isHit:&p]) >= 0)		/* hit ! */
	[self dragMarker: mknum :event];	/* modal drag loop only if there's a hit */
    [self display];
    return self;
}

/*
*  Hit Test -- Hit test is entirely done in the client side and this is faster than
*	relying on the server for hit tests.  This is possible because all test regions
*	are rectangles and not made up of complicated paths.
*  Returns maker number (0..N) of movable object.
*  Returns a negative value if NO hit.
*/
- (int) isHit:(const NXPoint *) p
{
float tempx, tempy;
NXRect hitRect;		// for hit testing
    // On first ever mouseDown, initialize pX_Vline and pY_Hline to location
    // of mouse cursor.
    if(pY_Hline < pYmin || pX_Vline < pXmin ) {
	pX_Vline = p->x;
	pY_Hline = p->y;
	if (delegate && [delegate respondsTo:@selector(markerDidMove:::::)]) {
	    // tempx, tempy are in View's coordinate in points, convert to real units.
	    tempx = (pX_Vline - bounds.origin.x)/ppxu + xdmin;
	    tempy = (pY_Hline - bounds.origin.y)/ppyu + ydmin;
	    [delegate markerDidMove:tag :2 :tempx :tempy :1];	// for cross pt
	}
	return(-1);		// indicate no hit for this
    }
    // void NXSetRect(NXRect *aRect, NXCoord x, NXCoord y, NXCoord width, NXCoord height) 
    // BOOL NXMouseInRect(const NXPoint *aPoint, const NXRect *aRect, BOOL flipped)
    // Setup hit test rectangle for trigger thresholds, pre-trigger length lines.
    // These are done later because long lines can be easily grabbed.
    NXSetRect(&hitRect, pX_Vline-4.0, pY_Hline-4.0, 8.0, 8.0);	// tiny rect at cross pt
    if(NXMouseInRect(p, &hitRect, NO))
	return(2);
    NXSetRect(&hitRect, pXmin, pY_Hline-4.0, (pXmax-pXmin), 8.0);	// horiz line
    if(NXMouseInRect(p, &hitRect, NO))
	return(0);
    NXSetRect(&hitRect, pX_Vline-4.0, pYmin, 8.0, (pYmax-pYmin));	// vert line
    if(NXMouseInRect(p, &hitRect, NO))
	return(1);
    // if no hit and current line positions near the edge, set lines to clicked location.
    if(pY_Hline <= pYmin || pY_Hline >= pYmax || pX_Vline <= pXmin || pX_Vline >= pXmax) {
	pX_Vline = p->x;
	pY_Hline = p->y;
	if (delegate && [delegate respondsTo:@selector(markerDidMove:::::)]) {
	    // tempx, tempy are in View's coordinate in points, convert to real units.
	    tempx = (pX_Vline - bounds.origin.x)/ppxu + xdmin;
	    tempy = (pY_Hline - bounds.origin.y)/ppyu + ydmin;
	    [delegate markerDidMove:tag :2 :tempx :tempy :1];	// for cross pt
	}
    }
    return(-1);			// no hit
}

// Draging of Lines -- called only when there's a MouseDown hit on something
// mk = marker or line index that got MouseDown.
- dragMarker:(int)mk :(NXEvent *)event
{
int      old_mask, drag_happened=0;
NXPoint	 p;
float	 tempx=0.0, tempy=0.0;
    lastx = lasty = 0.0;
    old_mask = [window addToEventMask:NX_MOUSEUPMASK|NX_MOUSEDRAGGEDMASK];
    event = [NXApp getNextEvent:NX_MOUSEUPMASK|NX_MOUSEDRAGGEDMASK];

    while (event->type != NX_MOUSEUP) {
	if (event->type == NX_MOUSEDRAGGED) {
	    p = event->location;		
	    [self convertPoint:&p fromView:nil];
	    if(!drag_happened) {
		drag_happened = 1;		// indicate that at least 1 drag event came
	    }
	}
	switch(mk) {
	    case 0:	/* horiz line */
		tempy = p.y;
	        if(tempy > pYmax) tempy = bounds.origin.y+bounds.size.height;
		if(tempy < pYmin) tempy = bounds.origin.y;
		pY_Hline = tempy;
		tempx = 0.0;
		break;
	    case 1:	/* vert line */
		tempx = p.x;
		if(tempx > pXmax) tempx = bounds.origin.x+bounds.size.width;
		if(tempx < pXmin) tempx = bounds.origin.x;
		pX_Vline = tempx;
		tempy = 0.0;
		break;
	    case 2:
		tempx = p.x;
		if(tempx > pXmax) tempx = bounds.origin.x+bounds.size.width;
		if(tempx < pXmin) tempx = bounds.origin.x;
		pX_Vline = tempx;
		tempy = p.y;
	        if(tempy > pYmax) tempy = bounds.origin.y+bounds.size.height;
		if(tempy < pYmin) tempy = bounds.origin.y;
		pY_Hline = tempy;
		break;

	    default:		// shouldn't happen
		break;
	}  /* end of switch(mk) */

	// Update display only if values have changed since last time
	if(tempx != lastx || tempy != lasty) {
	    lastx = tempx; lasty = tempy;
	    [self display];

	    // Report change in the marker position to delegate.
	    // Do it for drag.
	    if (drag_happened && delegate && continuousNotify
		&& [delegate respondsTo:@selector(markerDidMove:::::)]) {
		// tempx, tempy are in View's coordinate in points, convert to real units.
		tempx = (tempx - bounds.origin.x)/ppxu + xdmin;
		tempy = (tempy - bounds.origin.y)/ppyu + ydmin;
		[delegate markerDidMove:tag :mk :tempx :tempy :0];	// for drag
	    }
	    // NXPing();
	}


	event = [NXApp getNextEvent:NX_MOUSEUPMASK|NX_MOUSEDRAGGEDMASK
				waitFor:1000  threshold:NX_BASETHRESHOLD];
    }  /* end of while( <not mouse-up>) loop */

    // ### MouseUp detected, do necessary things and exit drag loop  ###
    [window setEventMask:old_mask];	// restore original event mask

    // Report change in the marker position to delegate.
    // Do it for mouse up event ater drag.
    if (drag_happened && delegate && [delegate respondsTo:@selector(markerDidMove:::::)]) {
	// tempx, tempy should already be correct and converted
	tempx = (pX_Vline - bounds.origin.x)/ppxu + xdmin;
	tempy = (pY_Hline - bounds.origin.y)/ppyu + ydmin;
    	[delegate markerDidMove:tag :mk :tempx :tempy :1];	// mouse up
    }
    return self;
}

// fakes a delegage message without mouse action using current marker positions.
- fakeMakerMovedFor:(int)mk
{
float tempx, tempy;
    if(!enabled || !cacheOK) return self;
    if (delegate && [delegate respondsTo:@selector(markerDidMove:::::)]) {
	tempx = (pX_Vline - bounds.origin.x)/ppxu + xdmin;
	tempy = (pY_Hline - bounds.origin.y)/ppyu + ydmin;
    	[delegate markerDidMove:tag :mk :tempx :tempy :1];	// mouse up
    }
    return self;
}

// Programmatically sets markers to specified position
// (float)pos in real XY unit, not in points.
- setMarkerPositionFor:(int)mk to:(float)pos
{
    if(mk)
        pX_Vline = (pos - xdmin)*ppxu;		// actually shoud add bounds.origin.x
    else
        pY_Hline = (pos - ydmin)*ppyu;
    [self fakeMakerMovedFor:mk];
    [self display];
    return self;
}


// This will change scale of plotting, but will not touch internally
// stored grid data: xd[], yd[] arrays.  May be used for zooming?
//
- setScaleXY
{
   ppxu = bounds.size.width/(xdmax - xdmin);	// scaling factor points per unit of X
   ppyu = bounds.size.height/(ydmax - ydmin);	// scaling factor points per unit of Y

   smallx = fabs(xdmax - xdmin) * 0.0001;
   smally = fabs(ydmax - ydmin) * 0.0001;
   return self;
}

- setEnabled:(BOOL)yn
{
    enabled = yn;
    return self;
}


// Set a comment string to be inserted into PostScript stream when printed or
// copied onto pasteboard as EPS.  The comment string must already be formatted
// as valid PostScript comments: Each new line must start with '%'.
//
// The comment does not appear at the top of EPS file.  It will be inserted
// immmediately AFTER the NeXT printPackage.ps, but before all the drawing
// commands for the contour appear.
// This is used to identify in detail the data used to create the contour
// plot, so we can look inside the EPS file later and relate to the raw data.
//
- setComment:(char *)cstr :(int)len
{
char *ptr, *ptr2, ch;
int  lfcount;
    if(comment) {
	free((void *)comment);
	comment = NULL;
    }
    comment = (char *)malloc((size_t)(len*sizeof(char)+2));
    strncpy(comment, cstr, len);
    comment[len] = '\0';	/* string terminator */

    // We need a duplicate comment with '\n' translated to '\012' char sequence
    if(comment012) {
	free((void *)comment012);
	comment012 = NULL;
    }
    lfcount = 0;
    ptr = comment;
    while( (ch = *ptr++) != '\0')
	if(ch == '\n') lfcount++;
    comment012 = (char *)malloc((size_t)(len + lfcount*4 + 2));
    ptr = comment;
    ptr2 = comment012;
    while( (ch = *ptr++) != '\0') {
	if(ch == '\n') {
	    *ptr2++ = '\\';
	    *ptr2++ = '0'; *ptr2++ = '1'; *ptr2++ = '2'; 
	}
	else
	    *ptr2++ = ch;
    }
    *ptr2 = '\0';
    return self;
}

// Given array x[] with npts elements, returns min and max values in the array.
//
- findMinMax:(float *)x :(int)npts :(float *)vmin :(float *)vmax
{
int i;
float fmin= MAXFLOAT;
float fmax= MINFLOAT;
    for(i=0; i<npts; i++)
    {
	if(x[i] < fmin) fmin = x[i];
	if(x[i] > fmax) fmax = x[i];
    }
    *vmin = fmin;
    *vmax = fmax;
    return self;
}

- setInvertY:(BOOL)invy
{
    invertY = invy;
    return self;
}

- setDebugEnable:(BOOL)de
{
    debug = de;
    return self;
}

- setBipolar:(BOOL)yn
{
    if(yn) bipolar = 1;
    else   bipolar = 0;
    [self setDefaultContourAttributes];
    return self;
}

- setFillEnable:(BOOL)fe
{
    doFill = fe;
    return self;
}

- setContourStrokeEnable:(BOOL)fe
{
    doContourStroke = fe;
    return self;
}

// Set the number of contour levels.
// It will be forced to an EVEN number if not already.
// It will automatically set default contour level values.
// If customized contour levels are needed,
// Call - setContourLevelArray:(float *)caa :(int)nl instead.
//
- setNumberOfContourLevels:(int)nl
{    
    if(nl == nclevels) return self;	/* no change */
    if(cA) free((void *)cA);
//    nclevels = nl/2*2;		// force it to be even
    nclevels = nl;		// force it to be even
    cA = (CntrAttribute *)malloc((size_t)nclevels*sizeof(CntrAttribute));
    [self setDefaultContourAttributes];
    return self;
}

// Set the range of contour level values used.
// This can be used to override default autoscaling.
// For real flexibility use -setContourLevelArray:: method.
//
- setMinMaxOfContourLevels:(float)min :(float)max
{
    fdmin = min;
    fdmax = max;
    [self setDefaultContourAttributes];
    return self;
}

- setContourLineColor:(NXColor)clc
{
    contourLineColor = clc;
    [self setDefaultContourAttributes];
    return self;
}

- setFrameColor:(NXColor)fc
{
    frameColor = fc;
    return self;
}

- setTickColor:(NXColor)fc
{
    tickColor = fc;
    return self;
}

- setGridColor:(NXColor)fc
{
    gridColor = fc;
    return self;
}

- setAxisColor:(NXColor)fc
{
    axisColor = fc;
    return self;
}

- setBackgroundColor:(NXColor)bc
{
    backgroundColor = bc;
    [self setDefaultContourAttributes];
    return self;
}

- setFillColors:(NXColor)pe :(NXColor)ne
{
    positiveColor = pe;
    negativeColor = ne;
    [self setDefaultContourAttributes];
    return self;
}


// This method allows customization of contour attributes to use
// overriding the default ones.
//
- setContourAttributeArray:(CntrAttribute *)caa :(int)nl;
{
// FIXME #### copy caa to cA.
    return self;
}

// Setup contour levels and other attributes to use based on the max and min of
// values in the data.  This will be called automatically when grid
// data are set.
//
- setDefaultContourAttributes
{
int i, ncl2;
float absfmax, step;
float rbase, gbase, bbase;	// RGB of background
float rp, gp, bp;		// RGB of positive extreme
float rn, gn, bn;		// RGB of negative extreme
float rpstep, gpstep, bpstep;
float rnstep, gnstep, bnstep;
float r, g, b;

    NXConvertColorToRGB(positiveColor, &rp, &gp, &bp);
    NXConvertColorToRGB(negativeColor, &rn, &gn, &bn);
    NXConvertColorToRGB(backgroundColor, &rbase, &gbase, &bbase);

    if(fabs(fdmax) > fabs(fdmin))
	absfmax = fabs(fdmax);
    else
	absfmax = fabs(fdmin);
//    if(fdmin<0.0 && fdmax>0.0)
    if(bipolar)
    {
	/* bipolar data */
        ncl2 = nclevels/2;
        nclevels = ncl2 * 2;	// Quietly change nclevels -- this may not be great.
//	bipolar = 1;		/* instance var flag */
	step = absfmax /(float)ncl2;
        rpstep = (rp - rbase)/(float)ncl2;
	gpstep = (gp - gbase)/(float)ncl2;
	bpstep = (bp - bbase)/(float)ncl2;
        rnstep = (rn - rbase)/(float)ncl2;
	gnstep = (gn - gbase)/(float)ncl2;
	bnstep = (bn - bbase)/(float)ncl2;
	for(i=0; i<ncl2; i++)
	{
	    /* positive */
	    cA[ncl2+i].dash = 0;
	    cA[ncl2+i].level = step*i + step/2.0;
	    r = rpstep*(i+1) + rbase;
	    g = gpstep*(i+1) + gbase;
	    b = bpstep*(i+1) + bbase;
	    cA[ncl2+i].fillcolor_hi = NXConvertRGBToColor(r, g, b);
	    if(i) cA[ncl2+i].fillcolor_lo = cA[ncl2+i-1].fillcolor_hi;
	    else  cA[ncl2+i].fillcolor_lo = backgroundColor;
	    cA[ncl2+i].linecolor = contourLineColor;
	    cA[ncl2+i].linewidth = minorContourLineWidth;
	    /* negative */
	    cA[ncl2-1-i].dash = 1;
	    cA[ncl2-1-i].level = -(step*i + step/2.0);
	    r = rnstep*(i+1) + rbase;
	    g = gnstep*(i+1) + gbase;
	    b = bnstep*(i+1) + bbase;
	    cA[ncl2-1-i].fillcolor_lo = NXConvertRGBToColor(r, g, b);
	    if(i) cA[ncl2-1-i].fillcolor_hi = cA[ncl2-i].fillcolor_lo;
	    else  cA[ncl2-1-i].fillcolor_hi = backgroundColor;
	    cA[ncl2-1-i].linecolor = contourLineColor;
	    cA[ncl2-1-i].linewidth = minorContourLineWidth;
	}
    }
    else
    {
	/* monopolar data */
//	bipolar = 0;		/* store flag in instance variable */
        if( fdmin < 0.0 )
	    fdmin = 0.0;		// This is a bug. for all negative data.
	step = (fdmax - fdmin)/(float)(nclevels+1);
        rpstep = (rp - rbase)/(float)nclevels;
	gpstep = (gp - gbase)/(float)nclevels;
	bpstep = (bp - bbase)/(float)nclevels;
	for(i=0; i<nclevels; i++) {
	    cA[i].dash = 0;
	    cA[i].level = fdmin + step*(i+1);
	    r = rpstep*(i+1) + rbase;
	    g = gpstep*(i+1) + gbase;
	    b = bpstep*(i+1) + bbase;
	    cA[i].fillcolor_hi = NXConvertRGBToColor(r, g, b);
	    if(i) cA[i].fillcolor_lo = cA[i-1].fillcolor_hi;
	    else  cA[i].fillcolor_lo = backgroundColor;
	    cA[i].linecolor = contourLineColor;
	    cA[i].linewidth = minorContourLineWidth;
	}
    }
    return self;
}

/* The following two methods must be used AFTER the method
 * - setCartesianGridData:(float *)f :(float)xmin :(float)xmax
 *				  :(float)ymin :(float)ymax
 *				  ofSize:(int)nnx :(int)nny
 *				  withInterpolationTo:(int)n1x :(int)n1y
 * has been called with
 * active interpolation.
*/

- horizCrossSectionAtY:(float)fy :(int)np forXs:(float *)fx :(float *)fcurve
{
int i;
    for(i=0; i<np; i++)
	splin2(xt-1, yt-1, f2, fd2a, nxraw, nyraw, fx[i], fy, &fcurve[i]);
    return self;
}

- vertCrossSectionAtX:(float)fx :(int)np forYs:(float *)fy :(float *)fcurve
{
int i;
    for(i=0; i<np; i++)
	splin2(xt-1, yt-1, f2, fd2a, nxraw, nyraw, fx, fy[i], &fcurve[i]);
    return self;
}

- xMinMax:(float *)xmin :(float *)xmax;
{
    *xmin = xdmin;
    *xmax = xdmax;
    return self;
}

- yMinMax:(float *)ymin :(float *)ymax
{
    *ymin = ydmin;
    *ymax = ydmax;
    return self;
}

- zMinMax:(float *)zmin :(float *)zmax
{
    *zmin = fdmin;
    *zmax = fdmax;
    return self;
}

// ==========================================================================
// Set surface data for regular Cartesian grid with optional interpolation.
// This is probably the most important method of this object.
// Use this if you want to plot data on uniform XY grid.
// Input data f[] of size (nx,ny) will be interpolated into (n1x, n1y) for
//     smoother contour plot if desired.  If no interpolation is needed,
//     make (n1x, n1y) same as (nx, ny).  This will turn off interpolation.
// 
//     (int)  nx, ny            : index sizes of input data
//     (int) n1x, n1y		: size of interpolated grid
//     for (i = 0; i < nx*ny; i++) (float)f[i]   : f value at (x,y)
//     f[i] is a 1-D array.  f(x,y) is in f[x + nx*y]
//
//     (xmin,ymin) and (xmax,ymax) define the (X,Y) domain of the plot.
//
// We extend the final (interpolated) data 3 points on each border to
// let the contour algorithm close all contours.  1-point wide region copy
// the data of the original border, and 2 points margins outside that are
// set to the base value.  This way, open contour lines that would have
// terminated at a border will get closed outside the original domain.
// This extra region for contour closure outside the original domain is
// outside the clip path, thus will not show up in the final plot.
// ---------------------------------------------------------------------------
//
- setCartesianGridData:(float *)f :(float)xmin :(float)xmax
				  :(float)ymin :(float)ymax
				  ofSize:(int)nnx :(int)nny
				  withInterpolationTo:(int)n1x :(int)n1y
{
int i, ix, iy;
int oldnx, oldny;
int redoAlloc = 0;
float xstep, ystep, fxtemp, fytemp;
float *fd1;

    xdmin = xmin;
    xdmax = xmax;
    ydmin = ymin;
    ydmax = ymax;
    if(nnx != nxraw || nny != nyraw)
	redoAlloc = 1;			// redo allocation of arrays
    oldnx = nxraw;			// save for use in free_fmatrix()
    oldny = nyraw;
    nxraw = nnx;
    nyraw = nny;
    // do reallocation only if # of data points has changed.
    nx = n1x+N6;
    ny = n1y+N6;
    nx1 = n1x;
    ny1 = n1y;
    if(ndata != ((n1x+N6)*(n1y+N6))) {
        ndata = (n1x+N6)*(n1y+N6);		/* new grid size */
    	if(xd) free(xd);
	if(yd) free(yd);
	if(fd) free(fd);
    	xd = (float *)malloc((size_t)ndata*sizeof(float));
    	yd = (float *)malloc((size_t)ndata*sizeof(float));
    	fd = (float *)malloc((size_t)ndata*sizeof(float));
    }

// ===== Interpolation Code =============================
    /* ######## Interpolate (nx, ny) to new size grid (n1x, n1y) ######## */
    // Allocate temporary input arrays for interpolation
    if(redoAlloc) {
	if(f2) free_fmatrix(f2, 1, oldnx, 1, oldny);
	if(fd2a) free_fmatrix(fd2a, 1, oldnx, 1, oldny);
	if(xt) free((void *)xt);
	if(yt) free((void *)yt);
	f2 = NULL;
	fd2a = NULL;
	xt = NULL;
	yt = NULL;
    }
    if(!xt) xt = (float *)malloc((size_t)nxraw*sizeof(float));
    if(!yt) yt = (float *)malloc((size_t)nyraw*sizeof(float));
    if(!f2) f2 = fmatrix(1, nxraw, 1, nyraw);
    if(!fd2a) fd2a = fmatrix(1, nxraw, 1, nyraw);

    // This will be freed within this method.
    fd1 =  (float *)malloc((size_t)nx1*ny1*sizeof(float)); /* for interpolation */

    xstep = (xdmax-xdmin)/(float)(nxraw-1);
    ystep = (ydmax-ydmin)/(float)(nyraw-1);
    for(ix=0; ix<nxraw; ix++)
	xt[ix] = xstep*ix + xdmin;
    for(iy=0; iy<nyraw; iy++)
	yt[iy] = ystep*iy + ydmin;

    /* copy original data to 2-D array for interpolation */
    for(ix=0; ix <nxraw; ix++)
	for(iy=0; iy <nyraw; iy++)
	    f2[ix+1][iy+1] = f[iy*nxraw+ix];

    /* precompute second-derivative arrays for splin2()
    * 2-nd derivative is returned in fd2a[][].
    * (xt-1) and (yt-1) are for passing C's zero-offset array to
    * splie2() which expects 1-offset arrays. */
    splie2(xt-1, yt-1, f2, nxraw, nyraw, fd2a);

    // Interpolate f(nnx,nny) to finer array fd1(n1x,n1y)
    xstep = (xdmax-xdmin)/(float)(nx1-1);
    ystep = (ydmax-ydmin)/(float)(ny1-1);
    for(ix=0; ix <nx1; ix++) {
	for(iy=0; iy <ny1; iy++) {
	    fxtemp = xstep*ix+xdmin;
	    fytemp = ystep*iy+ydmin;
	    splin2(xt-1, yt-1, f2, fd2a, nxraw, nyraw, fxtemp, fytemp, &fd1[iy*nx1+ix]);
    	}
    }

    /* Copy interpolated array to expanded array while setting border values. */
    for(i=0; i<ndata; i++) {
	ix = i % nx;
	iy = i / nx;
	if(ix >= N3 && ix < (nx-N3) && iy >= N3 && iy < (ny-N3))
	    fd[i] = fd1[(iy-N3)*n1x+(ix-N3)];		/* copy interpolated data */
	else if(iy < N3 && ix >= N3 && ix <(nx-N3))    	/* bottom */
	    fd[i] = fd1[ix-N3];
	else if(iy >= (ny-N3) && ix >= N3 && ix <(nx-N3))  	/* top */
	    fd[i] = fd1[(n1y-1)*n1x+(ix-N3)];
	else if(ix < N3 && iy >= N3 && iy <(ny-N3))		/* left */
	    fd[i] = fd1[(iy-N3)*n1x];
	else if(ix >=(nx-N3) && iy >= N3 && iy <(ny-N3))	/* right */
	    fd[i] = fd1[(iy-N3)*n1x + n1x-1];
	else if(ix < N3 && iy <N3)		/* lower left */
	    fd[i] = fd1[0];
	else if(ix < N3 && iy >= (ny-N3))	/* uppper left */
	    fd[i] = fd1[(n1y-1)*n1x];
	else if(ix >= (nx-N3) && iy < N3)	/* lower right */
	    fd[i] = fd1[n1x-1];
	else if(ix >= (nx-N3) && iy >= (ny-N3))
	    fd[i] = fd1[(n1y-1)*n1x + n1x-1];

	/* wipe 2 pixel border to base level */
	if(ix < N2 || ix >= (nx-N2) || iy <N2 || iy >= (ny-N2))
	    fd[i] = basevalue;		/* totally out */
    }  /* end for(i=0; i< ndata...) */


    // Free temporary input array.
    // if(f2) free_fmatrix(f2, 1, nxraw, 1, nyraw);
    // if(fd2a) free_fmatrix(fd2a, 1, nxraw, 1, nyraw);
    // if(xt) free((void *)xt);
    // if(yt) free((void *)yt);

    if(fd1) free((void *)fd1);

// ===== Interpolation Code =============================

    // Generate proper xd[], yd[] array for countour_().
    xstep = (xdmax-xdmin)/(float)(nx1-1);
    ystep = (ydmax-ydmin)/(float)(ny1-1);
    for(iy=0; iy<ny; iy++) {
	for(ix=0; ix<nx; ix++) {
	    xd[ix+iy*nx] = xstep*(ix-N3) + xdmin;
	    yd[ix+iy*nx] = ystep*(iy-N3) + ydmin;
	}
    }
    // Do appropriate scaling
    [self setScaleXY];

    // Get min and max fd[]
//    [self findMinMax:f :nxraw*nyraw :&fdmin :&fdmax];	// get it from original
    [self findMinMax:fd :ndata :&fdmin :&fdmax];	// get it from interpolated

    [self setDefaultContourAttributes];

    return self;
}



// Set XY grid points and values at each grid points.
// Grid does not have to be Cartesian or regular.
// Use this routine if you need flexibility.  No interpolation of
// grid is provided.  If you want interpolation,
// you have to do that yourself external to this object.
//     (int)  nj, nk                       : index sizes of data
//     for (i = 0; i < nj*nk; i++) (float)x[i]   : x coordinate
//     for (i = 0; i < nj*nk; i++) (float)y[i]   : y coordinate
//     for (i = 0; i < nj*nk; i++) (float)f[i]   : f value at (x[i], y[i])}
//     f[i] is a 1-D array.  f(j,k) is in f[j + nj*k].  Same for x[], y[].
//
// ***  FIXME ***
// NOTE: This methond does not support contour closure by expanding the domain.
//	 THIS METHOD PROBABLY DOES NOT WORK ANYMORE.
//
- setGridAndValueData:(float *)x :(float *)y :(float *)f ofSize:(int)nj :(int)nk
{
int i;
    // do reallocation only if # of data points has changed.
    nx = nj;
    ny = nk;
    if(ndata != (nj*nk))
    {
        ndata = nj*nk;
    	if(xd) free(xd);
	if(yd) free(yd);
	if(fd) free(fd);
    	xd = (float *)malloc((size_t)ndata*sizeof(float));
    	yd = (float *)malloc((size_t)ndata*sizeof(float));
    	fd = (float *)malloc((size_t)ndata*sizeof(float));
    }
    for(i=0; i<ndata; i++)	// copy data into instance variables
    {
	xd[i] = x[i];
	yd[i] = y[i];
	fd[i] = f[i];
    }

    // Do appropriate scaling
    // Get min and max of xd[], yd[], and fd[]
    [self findMinMax:xd :ndata :&xdmin :&xdmax];
    [self findMinMax:yd :ndata :&ydmin :&ydmax];
    [self findMinMax:fd :ndata :&fdmin :&fdmax];
    [self setScaleXY];
    [self setDefaultContourAttributes];

    return self;
}



// Controls width and whether a rectangular frame is drawn.
//
- setFrameEnable:(BOOL)fflag lineWidth:(float)fw
{

    frameON= fflag;	// Draw bounds frame
    frameLineWidth=fw;
    return self;
}

- setAxisEnable:(BOOL)fflag
{
    axisON = fflag;
    return self;
}

- setGridEnable:(BOOL)fflag
{
    gridON = fflag;
    return self;
}

- setTicksEnable:(BOOL)fflag
{
    ticksON = fflag;
    return self;
}

- setAxisLineWidth:(float)alw position:(float)x :(float)y
{
    axisLineWidth = alw;
    axisX = x;
    axisY = y;
    return self;
}

- setGridLineWidth:(float)glw spacing:(float)xs :(float)ys
{
    gridLineWidth = glw;
    gridXstep = xs;
    gridYstep = ys;
    return self;
}


- setNumXgrid:(int)numXgrid
{
    gridXstep = (xdmax - xdmin)/(float)(numXgrid-1);
    return self;
}


- setTickLineWidth:(float)tlw length:(float)tll
{
    tickLineWidth = tlw;
    tickLength = tll;
    return self;
}

- setContourLineWidthMinor:(float)clw andMajor:(float)clwm
{
    minorContourLineWidth = clw;
    majorContourLineWidth = clwm;
    [self setDefaultContourAttributes];
    return self;
}

- setGridDash:(float)lmark: (float)lspace
{
    griddash[0] = lmark;
    griddash[1] = lspace;
    return self;
}

- setContourDash:(float *)darray :(int)nd
{
int i;
    if(nd != Nncdash) {
	if(ncontdash) free(ncontdash);
        Nncdash = nd;
	ncontdash = (float *)malloc((Nncdash+1)*sizeof(float));
    }
    for(i=0; i<Nncdash; i++)
	ncontdash[i] = darray[i];
    return self;
}




// Contours with number of points less than nmp are not plotted.
// This is used to eliminate tiny "speckles" in contour plot.
//
- setMinNumberOfPointsPerContour:(int)mnp
{
    minNumPoints = mnp;
    return self;
}

- clear:sender
{
/*    ClearFlag = 1;
    [self display];
    ClearFlag = 0;
*/
    return self;
}

- copy:sender
{
    [self copyPScode:sender];
    return self;
}

// Copies the current view into the Pasteboard as PostScript.
//
- copyPScode:sender
{
NXStream *psStream;
id        pb;
char     *data;
int       dataLen, maxDataLen;

  /* Open a stream on memory where we will collect the PostScript */
    psStream = NXOpenMemory(NULL, 0, NX_WRITEONLY);
    if (!psStream)
      return self;

  /* Tell the Pasteboard we're going to copy PostScript */
    pb = [Pasteboard new];
    [pb declareTypes:&NXPostScriptPboardType num:1 owner:self];

  /* writes the PostScript for the whole plot as EPS into the stream */
    [self copyPSCodeInside:&bounds to:psStream];

  /* get the buffered up PostScript out of the stream */
    NXGetMemoryBuffer(psStream, &data, &dataLen, &maxDataLen);

  /* put the buffer in the Pasteboard, free the stream (and the buffer) */
    [pb writeType:NXPostScriptPboardType data:data length:dataLen];
    NXCloseMemory(psStream, NX_FREEBUFFER);
    return self;
}


// Save EPS for the view to a fileName.
//
- (int)saveEPStoFile:(const char *)fileName
{
BOOL saveOK;
int fdes;     // File descriptor
int retval = 0;
NXStream *stream = NULL;

    if(saveOK = 
	  (((fdes = open(fileName, O_WRONLY|O_CREAT|O_TRUNC, 0644)) != -1) &&
            (stream = NXOpenFile(fdes, NX_WRITEONLY)))) {
        [self copyPSCodeInside:&bounds to:stream];
    }
    if(stream) NXClose(stream);
    if(fdes != -1) close(fdes);

    if(!saveOK) retval = -1;
    return retval;
}



// Call back method messaged from computeContour() function.
// The method accumulates contour path information into a linked list
// each time it is called from the function.
// This method is called back from conputeContour() function for
// each contour path it finds.  ContourView object passes the id of
// itself as an argument to computeContour() so it can call back.
// icont = level index for contour.
// np    = number of points in contour path
// x, y  = arrays containing x and y values of each point on contour path.
//
- accumContour:(int)icont :(int)np :(float *)x :(float *)y
{
ContourPath *newcntr;	/* scratch pad pointer to new contourList items */
int c_closed = 1;

    if(np < minNumPoints) return self;		/* too few points for contour */

    if (fabs(x[0] - x[np-1]) < smallx && fabs(y[0] - y[np-1]) < smally) {
	c_closed = 1;
	x[np-1] = x[0];  y[np-1] = y[0];	/* guarantee closure */
    }
    else
	c_closed = 0;			/* contour not closed */

    /* ---- Save new contour in a linked list ----------------------------------- */
    /* Allocate new item for contourList */
    newcntr = (ContourPath *)malloc((size_t)sizeof(ContourPath));
    newcntr->num_pts = np;		/* # of (x,y) points in contour */
    newcntr->closed = c_closed;		/* contour is closed flag */
    newcntr->levelindex = icont;	/* levelindex-th contour (indx to contour level) */
    newcntr->level = cA[icont].level;	/* value of contour */
    /* Allocate arrays for X and Y coordinates. */
    newcntr->x = (float *)malloc((size_t)(np*sizeof(float)));
    newcntr->y = (float *)malloc((size_t)(np*sizeof(float)));
    /* copy (x,y) coordinates of contour path */
    memcpy((void *)newcntr->x, x, np * sizeof(float));
    memcpy((void *)newcntr->y, y, np * sizeof(float));
    /* insert new item into the list */
    newcntr->next = contourList;	/* point to previous one as next */
    contourList = newcntr;		/* for next time */
    numContours++;			/* increment number of contours counter */
    /* ---- New contour saved in a linked list ----------------------------------- */
    return self;
}

// This generates comments block that survives editing by Tailor.app
// for keeping track of the data source and parameters for the plot.
//
- beginTailorGroupWithComments
{
int j;
float rr, gg, bb;
long timenow;
char timestr[32];
DPSContext curContext = DPSGetCurrentContext();

    // Do it only if printing or copying to pasteboard.
    if(NXDrawingStatus != NX_DRAWING) {
	DPSPrintf(curContext, "\n%% Comments below should survive Tailor editing.\n\n");
        DPSPrintf(curContext, "(%s\\012", CV_VERSION);
	if(comment012)
	    DPSPrintf(curContext, "%s\\012", comment012);
        DPSPrintf(curContext, "%% X (H) range: %g - %g\\012", xdmin, xdmax);
        DPSPrintf(curContext, "%% Y (V) range: %g - %g\\012", ydmin, ydmax);
	if(gridON)
	   DPSPrintf(curContext, "%% Grid interval (xintvl, yintvl) = (%g, %g)\\012",
				gridXstep, gridYstep);
	if(axisON)
	   DPSPrintf(curContext, "%% Axes at (x, y) = (%g, %g)\\012", axisX, axisY);
	DPSPrintf(curContext, "%% Raw grid (%d, %d) interpolated to (%d, %d)\\012",
				nxraw, nyraw, nx1, ny1);
	DPSPrintf(curContext, "%% Number of contour levels: %d\\012", nclevels);
	DPSPrintf(curContext, "%% Number of contour paths: %d\\012", numContours);
	DPSPrintf(curContext, "%% Z (value) range: %g -- %g\\012", fdmin, fdmax);
	DPSPrintf(curContext, "%% Removed paths with fewer than %d points.\\012",
				minNumPoints);
	DPSPrintf(curContext, "%% Cross-hair at: (%.3g, %.3g)\\012",
			(pX_Vline - bounds.origin.x)/ppxu + xdmin,
			(pY_Hline - bounds.origin.y)/ppyu + ydmin);
	DPSPrintf(curContext, "%% ---- Graphic Attributes ----\\012");
	DPSPrintf(curContext, "%% Contour LineWidth: %g (minor), %g (major)\\012",
					minorContourLineWidth, majorContourLineWidth);
	DPSPrintf(curContext, "%% Dash (Neg Contour [mark space ...]): [");
	    for(j=0; j<Nncdash; j++)
		DPSPrintf(curContext, "%g ", ncontdash[j]);
        DPSPrintf(curContext, "]\\012");
	if(gridON) {
	    DPSPrintf(curContext, "%% Grid Dash [mark space] = [%g %g]\\012",
				griddash[0], griddash[1]);
	}		
	NXConvertColorToRGB(positiveColor, &rr, &gg, &bb);
	DPSPrintf(curContext, "%% Positive Color: %g, %g, %g\\012", rr, gg, bb);
	NXConvertColorToRGB(negativeColor, &rr, &gg, &bb);
	DPSPrintf(curContext, "%% Negative Color: %g, %g, %g\\012", rr, gg, bb);
	NXConvertColorToRGB(backgroundColor, &rr, &gg, &bb);
	DPSPrintf(curContext, "%% BackGND Color: %g, %g, %g\\012", rr, gg, bb);
	NXConvertColorToRGB(contourLineColor, &rr, &gg, &bb);
	DPSPrintf(curContext, "%% ContourLine Color: %g, %g, %g\\012", rr, gg, bb);
	DPSPrintf(curContext, ") TailorGroupBegin\n\n");

	if(comment) {
	  DPSPrintf(curContext, "%% For Acrobat Distiller to generate a Note with comments.\n\n");

	  DPSPrintf(curContext, "[ /Rect [0 0 300 150]\n");
	  DPSPrintf(curContext, "/Contents (");
	  DPSPrintf(curContext, "%s)\n", comment);
	    time(&timenow); strcpy(timestr, ctime(&timenow)); timestr[24] = '\0';
	  DPSPrintf(curContext, "/Title (ContourView plot: %s)\n", timestr);
	  DPSPrintf(curContext, "/Open false\n");   /* initialize for iconized appearance */
	  DPSPrintf(curContext, "/ANN pdfmark\n\n");
	}
    }
    return self;
}

- endTailorGroup
{
DPSContext curContext = DPSGetCurrentContext();
    if(NXDrawingStatus != NX_DRAWING) {
        DPSPrintf(curContext, "\n%% End of ContourView drawing.\n");
	DPSPrintf(curContext, "TailorGroupEnd\n\n");
    }
    return self;
}

// This puts the data trace info right after the %% comments.
// %%EndComments etc moved to after this, which is not good, but OK.
//
- endHeaderComments
{
int j;
float rr, gg, bb;
DPSContext curContext = DPSGetCurrentContext();
    if(NXDrawingStatus != NX_DRAWING) {
        DPSPrintf(curContext, "\n\n%s\n", CV_VERSION);
	if(comment)
	    DPSPrintf(curContext, "%s\n", comment);
        DPSPrintf(curContext, "%% X (horiz) range: %g - %g\n", xdmin, xdmax);
        DPSPrintf(curContext, "%% Y (vert) range: %g - %g\n", ydmin, ydmax);
	if(gridON)
	   DPSPrintf(curContext, "%% Grid interval (xintvl, yintvl) = (%g, %g)\n",
				gridXstep, gridYstep);
	if(axisON)
	   DPSPrintf(curContext, "%% Axes at (x, y) = (%g, %g)\n", axisX, axisY);
	DPSPrintf(curContext, "%% Raw grid (%d, %d) interpolated to (%d, %d)\n",
				nxraw, nyraw, nx1, ny1);
	DPSPrintf(curContext, "%% Number of contour levels: %d\n", nclevels);
	DPSPrintf(curContext, "%% Number of contour paths: %d\n", numContours);
	DPSPrintf(curContext, "%% Z (value) range: %g -- %g\n", fdmin, fdmax);
	DPSPrintf(curContext, "%% Removed paths with fewer than %d points.\n",
				minNumPoints);
	DPSPrintf(curContext, "%% Cross-hair at: (%.3g, %.3g)\n",
			(pX_Vline - bounds.origin.x)/ppxu + xdmin,
			(pY_Hline - bounds.origin.y)/ppyu + ydmin);
	DPSPrintf(curContext, "%% ---- Graphic Attributes ----\n");
	DPSPrintf(curContext, "%% Contour LineWidth: %g (minor), %g (major)\n",
					minorContourLineWidth, majorContourLineWidth);
	DPSPrintf(curContext, "%% Dash (Neg Contour [mark space ...]): [");
	    for(j=0; j<Nncdash; j++)
		DPSPrintf(curContext, "%g ", ncontdash[j]);
        DPSPrintf(curContext, "]\n");
	if(gridON) {
	    DPSPrintf(curContext, "%% Grid Dash [mark space] = [%g %g]\n",
				griddash[0], griddash[1]);
	}		
	NXConvertColorToRGB(positiveColor, &rr, &gg, &bb);
	DPSPrintf(curContext, "%% Positive Color (RGB): %g, %g, %g\n", rr, gg, bb);
	NXConvertColorToRGB(negativeColor, &rr, &gg, &bb);
	DPSPrintf(curContext, "%% Negative Color (RGB): %g, %g, %g\n", rr, gg, bb);
	NXConvertColorToRGB(backgroundColor, &rr, &gg, &bb);
	DPSPrintf(curContext, "%% BackGND Color (RGB): %g, %g, %g\n", rr, gg, bb);
	NXConvertColorToRGB(contourLineColor, &rr, &gg, &bb);
	DPSPrintf(curContext, "%% ContourLine Color (RGB): %g, %g, %g\n", rr, gg, bb);
	DPSPrintf(curContext, "%%\n\n");
    }
    return [super endHeaderComments];
}

// This is inserted for embedding comments that can persist across
// editing by Tailor.app (FirstClass).  Tailor strips regular comments
// from its output.
//
- endPrologue
{
DPSContext curContext = DPSGetCurrentContext();
    if(NXDrawingStatus != NX_DRAWING) {
	DPSPrintf(curContext, "\n\n%%%%BeginResource: procset (Tailor group constructs)\n");
	DPSPrintf(curContext, "/TailorGroupBegin {pop} def /TailorGroupEnd {} def\n");
	DPSPrintf(curContext, "/pdfmark where\n");
	DPSPrintf(curContext, "{pop} {userdict /pdfmark /cleartomark load put} ifelse\n");
	DPSPrintf(curContext, "%%%%EndResource\n\n");
    }
    return [super endPrologue];
}


// Draws objects (lines/markers) that are manipulated by the user.
- drawMarkers
{
DPSContext curContext = DPSGetCurrentContext();

    if(pY_Hline > pYmin && pY_Hline < pYmax && pX_Vline > pXmin && pX_Vline < pXmax) {
        if(NXDrawingStatus != NX_DRAWING) {
            DPSPrintf(curContext, "\n%% Cross-hair\n");
        }
	PSmoveto(pXmin, pY_Hline);
	PSrlineto(bounds.size.width, 0.0);
	PSmoveto(pX_Vline, pYmin);
	PSrlineto(0.0, bounds.size.height);
	PSsetgray(0.66667);			// stroke with lightgray, 0 linewidth
	PSsetlinewidth(0.0);
	PSstroke();
    }
    return self;
}


// Mandatory method for any View that does something useful.
// This will override the default -drawSelf::, which does no real drawing.
// Should never be called directly.  It will be called by [self display];
- drawSelf:(const NXRect *)r:(int)count
{    
    if(!enabled || !cacheOK) return self;
    if(NXDrawingStatus == NX_DRAWING) {
	[imageCache composite: NX_SOVER toPoint: &bounds.origin];
	[self drawMarkers];		// vert/horiz lines
    }
    else {
	[self beginTailorGroupWithComments];
	[self drawAll];			// this does actual PS generation
	if(printCrossHair)
	    [self drawMarkers];		// vert/horiz lines
	[self endTailorGroup];
    }
    return self;
}


// This plots a single contour path with optional color fill.
//
- plotContour:(ContourPath *)cntr :(DPSContext)cContext
{
int i, np;
float xp, yp;
static char *ocmsg[] = {"open", "closed" };
static char *hlmsg[] = {"low", "high" };
    np = cntr->num_pts;		/* # points in this contour */
    if(NXDrawingStatus != NX_DRAWING)
        DPSPrintf(cContext,
	        "\n%% ### Contour: level index=%d, value=%g, #pts=%d, [%s], inside=%s\n",
		cntr->levelindex, cntr->level, np,
		ocmsg[(cntr->closed)?1:0], hlmsg[(cntr->hi_inside)?1:0]);

    for(i=0; i<np; i++) {
	xp = (cntr->x[i] - xdmin)*ppxu;
        if(invertY)
	    yp = (ydmax - cntr->y[i])*ppyu;
	else	
	    yp = (cntr->y[i] - ydmin)*ppyu;
	if(i==0)
	    PSmoveto(xp, yp);
	else
	    PSlineto(xp, yp);
    }
    if(cntr->closed)
	    PSclosepath();

    if(doFill) {
	if(cntr->hi_inside)
	    NXSetColor(cA[cntr->levelindex].fillcolor_hi);
	else
	    NXSetColor(cA[cntr->levelindex].fillcolor_lo);
	if(doContourStroke) {
	    PSgsave();
	    PSfill();
	    PSgrestore();	// save-restore contour path for later stroke
	}
	else
	    PSfill();		// Fill only, no stroke.
    }

    if(doContourStroke || !doFill) {
	if(cA[cntr->levelindex].dash)
	    PSsetdash(ncontdash, 2, 0.0);
	    /* fprintf(fpp, "[4 4] 0 setdash\n"); */ /* dashed curve */
	else
	    PSsetdash(psolid, 0, 0.0);
	    /* fprintf(fpp, "[] 0 setdash\n");	*/	/* solid curve */
	NXSetColor(cA[cntr->levelindex].linecolor);
	PSsetlinewidth(cA[cntr->levelindex].linewidth);
	PSstroke();
    }
    return self;
}

// Finds if inside of contour path is going up or going down
// (higher or lower than the contour level).
// This is critical for filling contours with correct color.
// Contour lines are drawn at a certain height.  Colors that should
// fill the inside depends on whether the inside is higher than the level
// or lower.  This method tries to figure this out for each contour.
//
- findInsideHighLow:(ContourPath *)cntr
{
int i, i3, j, ii=0, ix=0, iy=0, npp, inside=0, gpfound=0, jx=0, jy=0;
float xstep, ystep, xg=0.0, yg=0.0, fg= 0.0;
static int ipx[] = { 0, -1,  1,  0, -1,  1,  0, -1,  1 };
static int ipy[] = { 0,  0,  0, -1, -1, -1,  1,  1,  1 };

    xstep = (xdmax-xdmin)/(float)(nx1-1);
    ystep = (ydmax-ydmin)/(float)(ny1-1);
    npp = cntr->num_pts;

    // Typical setting for hi_inside flag
    if(cntr->level > 0.0) cntr->hi_inside = 1;
    else		  cntr->hi_inside = 0;

    // Now try to determine it exactly.  Start with i=3 in an attempt to get
    // a good point that is not on the edge.
    for(i3=3; i3<(npp+3); i3++) {
	i = i3 % npp;		/* i = [3, .. (npp-1), 0, 1, 2] */
	if( ![self pointInDomain:cntr->x[i] :cntr->y[i]] )
		continue;		// this contour point is outside domain
	ix  = (cntr->x[i] - xdmin)/xstep + 0.5;	// indices of closest grid point
	iy  = (cntr->y[i] - ydmin)/ystep + 0.5;
	// find indices (jx, jy) which is inside contour and inside domain
	// by checking 3x3 points centered on (ix, iy)
	for(j=0; j<9; j++) {
	    jx = ix + ipx[j];	
	    jy = iy + ipy[j];
	    xg = xstep*(float)jx + xdmin;
	    yg = ystep*(float)jy + ydmin;
	    if([self pointInDomain:xg :yg] && non_zero_winding(xg, yg, cntr->x, cntr->y, npp)==2)
	    {
		inside = 1;
		break;		// out of for(j..)
	    }
	}
	if(!inside) continue;	// try another point on contour

	ii = i;
	gpfound = 1;
	break;				// good point found.
    }  /* end of for(i=0...) */

    // If inside point not found, use closed point that is outside
    if(!inside) {
        for(i3=3; i3<(npp+3); i3++) {
	    i = i3 % npp;		/* i = [3, .. (npp-1), 0, 1, 2] */
	    if( ![self pointInDomain:cntr->x[i] :cntr->y[i]] )
		    continue;		// this contour point is outside domain
	    jx  = (cntr->x[i] - xdmin)/xstep + 0.5;	// indices of closest grid point
	    jy  = (cntr->y[i] - ydmin)/ystep + 0.5;
	    xg = xstep*(float)jx + xdmin;
	    yg = ystep*(float)jy + ydmin;
	    if( ![self pointInDomain:xg :yg])	// at least it must be inside domain
		continue;

	    ii = i;
	    gpfound = 1;
	    if(non_zero_winding(xg, yg, cntr->x, cntr->y, npp) == 2) inside = 1;
	    break;
	}
    }

    fg = fd[(jy+N3)*nx+jx+N3];		// grid value
    if(inside) {
	if(fg > cntr->level) cntr->hi_inside = 1;	// inside is high
	else		     cntr->hi_inside = 0;
    } else {
	if(fg > cntr->level) cntr->hi_inside = 0;	// outside is high
	else		     cntr->hi_inside = 1;
    }

    if(debug && !gpfound)
	fprintf(stderr,
		"ContourView: Can't determine hi_inside: idx=%d, val=%g, #pt=%d, hi_ins=%d\n",
			cntr->levelindex, cntr->level, npp, cntr->hi_inside);
    return self;
}

- (BOOL)pointInDomain:(float)xx :(float)yy
{
    if(xx >= xdmin && xx <= xdmax && yy >= ydmin && yy <= ydmax)
	return YES;
    else
	return NO;
}


- freeContourList;
{
ContourPath *cntr;
	while((cntr = contourList))
	{
	    if(cntr->x) free(cntr->x);
	    if(cntr->y) free(cntr->y);
	    contourList = cntr->next;
	    free(cntr);
	};
        contourList = NULL;
	return self;
}

- setDelegate:anObject
{
    delegate = anObject;
    return self;
}

- delegate
{
    return delegate;
}

- setTag:(int)t
{
    tag = t;
    return self;
}
- (int)tag
{
    return tag;
}

@end


