/*ScianDatasets.c
  Eric Pepke
  August 17, 1990
  Stuff for Datasets
*/

#include "Scian.h"
#include "ScianTypes.h"
#include "ScianLists.h"
#include "ScianIcons.h"
#include "ScianWindows.h"
#include "ScianObjWindows.h"
#include "ScianVisWindows.h"
#include "ScianButtons.h"
#include "ScianTextBoxes.h"
#include "ScianControls.h"
#include "ScianDatasets.h"
#include "ScianErrors.h"
#include "ScianColors.h"
#include "ScianDialogs.h"
#include "ScianVisObjects.h"
#include "ScianStyle.h"
#include "ScianSpaces.h"
#include "ScianIDs.h"
#include "ScianArrays.h"
#include "ScianMethods.h"
#include "ScianTimers.h"
#include "ScianDepend.h"
#include "ScianFiles.h"
#include "ScianFileSystem.h"
#include "ScianObjFunctions.h"
#include "ScianSockets.h"
#include "ScianEvents.h"
#include "ScianSciences.h"
#include "ScianGeometry.h"

#define VISMENU
#define SHOWCTLMENU

ObjPtr geometryClass;			/*Class for geometry object*/
ObjPtr datasetClass, data3DScalar, data2DScalar, data1DVector, data3DUnstructSurface;
ObjPtr iconDataset;
ObjPtr icon1DVector, icon2DVector, icon3DVector, icon4DVector;
ObjPtr icon1DScalar, icon2DScalar, icon3DScalar, icon4DScalar;
ObjPtr iconGeometry;

Bool timedDatasets = true;		/*True iff reading timed datasets*/
ObjPtr dataFormClass;
ObjPtr allDatasets = NULLOBJ;		/*All datasets*/
static ObjPtr dsLocArray = NULLOBJ;	/*Location of icon for dataset*/

Bool onePalette = false;
ObjPtr commonPalette;

DatasetBuffer curFields[MAXNCURFIELDS];

#ifdef PROTO
ObjPtr NewGeometryDataset(char *name)
#else
ObjPtr NewGeometryDataset(name)
char *name;
#endif
{
    ObjPtr retVal;
    retVal = NewObject(datasetClass, 0L);
    SetVar(retVal, DEFAULTICON, iconGeometry);
    SetVar(retVal, NAME, NewString(name));

    return retVal;
}

#ifdef PROTO
ObjPtr NewDataset(char *name, int rank, long *dims, int nComponents)
#else
ObjPtr NewDataset(name, rank, dims, nComponents)
char *name;
int rank;
long *dims;
int nComponents;
#endif
/*Makes a new structured dataset
  name		is the name of the dataset
  rank		is the rank of the dataset
  dims		is a pointer to the dimensions
  nComponents	is the number of components, 0 for scalar
*/
{
    ObjPtr retVal;
    ObjPtr dimensions;
    real *elements;
    ObjPtr data;
    int k;

    /*Make the dataset*/
    retVal = NewObject(datasetClass, 0L);

    /*Make its dimensions*/
    if (rank)
    {
	dimensions = NewRealArray(1, (long) rank);
	elements = ELEMENTS(dimensions);
	for (k = 0; k < rank; ++k)
	{
	    elements[k] = (real) dims[k];
	}
	SetVar(retVal, DIMENSIONS, dimensions);
    }

    /*Set its NCOMPONENTS*/
    if (nComponents)
    {
	SetVar(retVal, NCOMPONENTS, NewInt(nComponents));
    }

    /*Set its NAME*/
    SetVar(retVal, NAME, NewString(name));

    /*Set its DATA*/
    if (nComponents)
    {
	ObjPtr *objectElements;
	long dummy;

	/*Vector data*/
	dummy = nComponents;
	data = NewArray(AT_OBJECT, 1, &dummy);

	objectElements = (ObjPtr *) ELEMENTS(data);

	for (k = 0; k < nComponents; ++k)
	{
	    objectElements[k] = NewArray(AT_REAL, rank, dims);
	}
    }
    else
    {
	/*Scalar data*/
	data = NewArray(AT_REAL, rank, dims);
    }
    SetVar(retVal, DATA, data);

    /*Set its DEFAULTICON*/
    if (nComponents)
    {
	/*Vector*/
	switch (rank)
	{
	    case 0:	/*Really should have a value, here*/
	    case 1:     SetVar(retVal, DEFAULTICON, icon1DVector); break;
	    case 2:     SetVar(retVal, DEFAULTICON, icon2DVector); break;
	    case 3:     SetVar(retVal, DEFAULTICON, icon3DVector); break;
	    default:    SetVar(retVal, DEFAULTICON, icon4DVector); break;
	}
    }
    else
    {
	/*Scalar*/
	switch (rank)
	{
	    case 0:	/*Really should have a value, here*/
	    case 1:     SetVar(retVal, DEFAULTICON, icon1DScalar); break;
	    case 2:     SetVar(retVal, DEFAULTICON, icon2DScalar); break;
	    case 3:     SetVar(retVal, DEFAULTICON, icon3DScalar); break;
	    default:    SetVar(retVal, DEFAULTICON, icon4DScalar); break;
	}
    }

    return retVal;
}

#ifdef PROTO
ObjPtr NewSeparableDataForm(char *name, int rank, long *dims, real *scales[])
#else
ObjPtr NewSeparableDataForm(name, rank, dims, scales)
char *name;
int rank;
long *dims;
real *scales[];
#endif
/*Makes a new rectilinear data form with separable axes
  name		is the name of the dataset
  rank		is the rank of the dataset
  dims		is a pointer to the dimensions
  scales	is an array, rank rank, containing elements of size
		dims[k] giving scale along that axis
*/
{
    ObjPtr retVal;
    ObjPtr dimensions;
    ObjPtr bounds;
    real *elements;
    ObjPtr *objElements;
    ObjPtr data, axes;
    int k;
    long temp;

    /*Create the data form*/
    retVal = NewObject(dataFormClass, 0L);

    /*Make its dimensions*/
    dimensions = NewRealArray(1, (long) rank);
    elements = ELEMENTS(dimensions);
    for (k = 0; k < rank; ++k)
    {
	elements[k] = (real) dims[k];
    }
    SetVar(retVal, DIMENSIONS, dimensions);

    /*Set its NAME*/
    SetVar(retVal, NAME, NewString(name));

    /*Make its BOUNDS*/
    bounds = NewRealArray(1, 2L * (long) rank);
    if (!bounds)
    {
	return NULLOBJ;
    }
    elements = ELEMENTS(bounds);
    for (k = 0; k < rank; ++k)
    {
	int i;
	real min, max;

	min = max = scales[k][0];
	for (i = 1; i < dims[k]; ++i)
	{
	    if (scales[k][i] < min) min = scales[k][i];
	    if (scales[k][i] > max) max = scales[k][i];
	}
	elements[k * 2] = min;
	elements[k * 2 + 1] = max;
    }
    SetVar(retVal, BOUNDS, bounds);

    /*Now make a dataset which specifies the field*/
    data = NewObject(datasetClass, 0);
    if (!data)
    {
	return NULLOBJ;
    }

    /*Set its NAME*/
    SetVar(data, NAME, NewString(name));

    /*Set its NCOMPONENTS*/
    SetVar(data, NCOMPONENTS, NewInt(rank));

    /*Make an array of axes for the data*/
    temp = rank;
    axes = NewArray(AT_OBJECT, 1, &temp);
    if (!axes)
    {
	return NULLOBJ;
    }
    SetVar(data, DATA, axes);

    /*Fill the axes*/
    objElements = ELEMENTS(axes);
    for (k = 0; k < rank; ++k)
    {
	ObjPtr indices, scaleArray;
	long i;

	/*Copy over the scale*/
	scaleArray = NewRealArray(1, dims[k]);
	elements = ELEMENTS(scaleArray);
	objElements[k] = scaleArray;

	for (i = 0; i < dims[k]; ++i)
	{
	    elements[i] = scales[k][i];
	}
	
	/*Give it an INDICES array*/
	indices = NewRealArray(1, 1L);
	*(real *)ELEMENTS(indices) = k;
	SetVar(objElements[k], INDICES, indices);
    }
    SetVar(retVal, DATA, data);

    return retVal;
}

#ifdef PROTO
ObjPtr NewUnstructuredDataForm(char *name, int rank, long *dims, real *bounds, ObjPtr dataset)
#else
ObjPtr NewUnstructuredDataForm(name, rank, dims, bounds, dataset)
char *name;
int rank;
long *dims;
real *bounds;
ObjPtr dataset;
#endif
/*Makes a new unstructured data form
  name		is the name of the data form
  rank		is the rank of the dataset
		0 = just nodes
		1 = + links
		2 = + faces
		3 = + cells, etc.
  dims		is a pointer to the dimensions
		This must be one greater than rank
		0 = number of nodes, 1 = number of links, etc.
  bounds	is the bounds of the dataform, rank # components of dataform
  dataset	is the dataset that is used for the data form
*/
{
    ObjPtr retVal;
    ObjPtr dimensions;
    ObjPtr boundsArray;
    real *elements;
    ObjPtr *objElements;
    ObjPtr data, axes;
    ObjPtr var;
    int nComponents;
    int k;
    long temp;

    /*Create the data form*/
    retVal = NewObject(dataFormClass, 0L);
    SetVar(retVal, UNSTRUCTURED, ObjTrue);

    /*Make its dimensions*/
    dimensions = NewRealArray(1, (long) rank + 1);
    elements = ELEMENTS(dimensions);
    for (k = 0; k < rank + 1; ++k)
    {
	elements[k] = (real) dims[k];
    }
    SetVar(retVal, DIMENSIONS, dimensions);

    /*Set its NAME*/
    SetVar(retVal, NAME, NewString(name));

    /*Calculate its BOUNDS*/
    var = GetIntVar("NewUnstructuredDataForm", dataset, NCOMPONENTS);
    if (var)
    {
	nComponents = GetInt(var);
    }
    else
    {
	return NULLOBJ;
    }
    SetVar(retVal, NCOMPONENTS, NewInt(nComponents));
    boundsArray = NewRealArray(1, 2L * (long) nComponents);
    if (!boundsArray)
    {
	return NULLOBJ;
    }
    elements = ELEMENTS(boundsArray);
    for (k = 0; k < nComponents; ++k)
    {
	elements[k * 2] = bounds[k * 2];
	elements[k * 2 + 1] = bounds[k * 2 + 1];
    }
    SetVar(retVal, BOUNDS, boundsArray);

    SetVar(retVal, DATA, dataset);

    return retVal;
}

#ifdef PROTO
void AppendPolygonToDataset(ObjPtr dataset, long nVertices, long vertices[])
#else
void AppendPolygonToDataset(dataset, nVertices, vertices)
ObjPtr dataForm;
long nVertices;
long vertices;
#endif
/*Appends an indexed polygon to a data set*/
{
    ObjPtr geometry;

    geometry = GetVar(dataset, DATA);
    if (!geometry)
    {
	geometry = NewGeometry();
	SetVar(dataset, DATA, geometry);
    }
    AppendPolygonToGeometry(geometry, nVertices, vertices);
}

#ifdef PROTO
ObjPtr NewRegularDataForm(char *name, int rank, long *dims, real *bounds)
#else
ObjPtr NewRegularDataForm(name, rank, dims, bounds)
char *name;
int rank;
long *dims;
real *bounds;
#endif
/*Makes a new regulat data form
  name		is the name of the dataset
  rank		is the rank of the dataset
  dims		is a pointer to the dimensions
  bounds	is the bounds of the dataform, rank 1, size 2 * rank
*/
{
    ObjPtr retVal;
    ObjPtr dimensions;
    ObjPtr boundsArray;
    real *elements;
    ObjPtr *objElements;
    ObjPtr data, axes;
    int k;
    long temp;

    /*Create the data form*/
    retVal = NewObject(dataFormClass, 0L);

    /*Make its dimensions*/
    dimensions = NewRealArray(1, (long) rank);
    elements = ELEMENTS(dimensions);
    for (k = 0; k < rank; ++k)
    {
	elements[k] = (real) dims[k];
    }
    SetVar(retVal, DIMENSIONS, dimensions);

    /*Set its NAME*/
    SetVar(retVal, NAME, NewString(name));

    /*Make its BOUNDS*/
    boundsArray = NewRealArray(1, 2L * (long) rank);
    if (!boundsArray)
    {
	return NULLOBJ;
    }
    elements = ELEMENTS(boundsArray);
    for (k = 0; k < rank; ++k)
    {
	elements[k * 2] = bounds[k * 2];
	elements[k * 2 + 1] = bounds[k * 2 + 1];
    }
    SetVar(retVal, BOUNDS, boundsArray);

    return retVal;
}

#ifdef PROTO
ObjPtr NewCurvilinearDataForm(char *name, int rank, long *dims, real *bounds, ObjPtr dataset)
#else
ObjPtr NewCurvilinearDataForm(name, rank, dims, bounds, dataset)
char *name;
int rank;
long *dims;
real *bounds;
ObjPtr dataset;
#endif
/*Makes a new curvilinear data
  name		is the name of the data form
  rank		is the rank of the data form
  dims		is a pointer to the dimensions
  bounds	is the bounds of the dataform, rank # components of dataform
  dataset	is the dataset that is used for the data form
*/
{
    ObjPtr retVal;
    ObjPtr dimensions;
    ObjPtr boundsArray;
    real *elements;
    ObjPtr *objElements;
    ObjPtr data, axes;
    ObjPtr var;
    int nComponents;
    int k;
    long temp;

    /*Create the data form*/
    retVal = NewObject(dataFormClass, 0L);

    /*Make its dimensions*/
    dimensions = NewRealArray(1, (long) rank);
    elements = ELEMENTS(dimensions);
    for (k = 0; k < rank; ++k)
    {
	elements[k] = (real) dims[k];
    }
    SetVar(retVal, DIMENSIONS, dimensions);

    /*Set its NAME*/
    SetVar(retVal, NAME, NewString(name));

    /*Calculate its BOUNDS*/
    var = GetIntVar("NewCurvilinearDataForm", dataset, NCOMPONENTS);
    if (var)
    {
	nComponents = GetInt(var);
    }
    else
    {
	return NULLOBJ;
    }
    SetVar(retVal, NCOMPONENTS, NewInt(nComponents));
    boundsArray = NewRealArray(1, 2L * (long) nComponents);
    if (!boundsArray)
    {
	return NULLOBJ;
    }
    elements = ELEMENTS(boundsArray);
    for (k = 0; k < nComponents; ++k)
    {
	elements[k * 2] = bounds[k * 2];
	elements[k * 2 + 1] = bounds[k * 2 + 1];
    }
    SetVar(retVal, BOUNDS, boundsArray);

    SetVar(retVal, DATA, dataset);

    return retVal;
}

#ifdef PROTO
void SetDatasetPalette(ObjPtr dataset, ObjPtr palette)
#else
void SetDatasetPalette(dataset, palette)
ObjPtr dataset, palette;
#endif
/*Sets a dataset's palette*/
{
    FieldPaletteName(palette, dataset);
    SetVar(dataset, CPALETTE, palette);
}

#ifdef PROTO
void SetDatasetForm(ObjPtr dataset, ObjPtr dataForm)
#else
void SetDatasetForm(dataset, dataForm)
ObjPtr dataset, dataForm;
/*Sets a dataset's data form*/
#endif
{
    SetVar(dataset, DATAFORM, dataForm);
}

ObjPtr GetLongName(dataset)
ObjPtr dataset;
/*Gets the long name of dataset*/
{
    FuncTyp method;
    method = GetMethod(dataset, GETLONGNAME);
    if (method)
    {
	return (*method)(dataset);
    }
    else
    {
	return NULLOBJ;
    }
}

ObjPtr GetPlainDatasetLongName(dataset)
ObjPtr dataset;
/*Gets the long name of a plain dataset*/
{
    return GetVar(dataset, NAME);
}

#ifdef PROTO
static DeleteBuckets(int whichField)
#else
static DeleteBuckets(whichField)
int whichField;
#endif
/*Deletes the buckets from within a field*/
{
    if (curFields[whichField] . bucketMinMax)
    {
	/*It has buckets, get rid of them*/
	long k;
	long nBuckets;

	/*Figure out how many buckets there are*/
	nBuckets = 1;
	for (k = 0; k < curFields[whichField] . nBucketIndices; ++k)
	{
	    if (curFields[whichField] . bucketDim[k])
	    {
		nBuckets += curFields[whichField] . bucketDim[k];
	    }
	}

	/*Empty the buckets*/
	for (k = 0; k < nBuckets; ++k)
	{
	    if (curFields[whichField] . buckets[k] . indices)
	    {
		free(curFields[whichField] . buckets[k] . indices);
	    }
	    if (curFields[whichField] . buckets[k] . samples)
	    {
		free(curFields[whichField] . buckets[k] . samples);
	    }
	}
	free(curFields[whichField] . buckets);
	curFields[whichField] . buckets = 0;

	/*Get rid of the min and max*/
	free(curFields[whichField] . bucketMinMax);
	curFields[whichField] . bucketMinMax = 0;

	/*Get rid of the dimensions*/
	free(curFields[whichField] . bucketDim);
	curFields[whichField] . bucketDim = 0;
    }
}

static void CleanCurField(k)
int k;
/*Cleans the current field number k*/
{
    DeleteBuckets(k);
    if (curFields[k] . components)
    {
	int i;
	for (i = 0; i < curFields[k] . nComponents; ++i)
	{
	    if (curFields[k] . components[i] . indices)
	    {
		free(curFields[k] . components[i] . indices);
		curFields[k] . components[i] . indices = 0;
	    }
	    if (curFields[k] . components[i] . dimensions)
	    {
		free(curFields[k] . components[i] . dimensions);
		curFields[k] . components[i] . dimensions = 0;
	    }
	    if (curFields[k] . components[i] . steps)
	    {
		free(curFields[k] . components[i] . steps);
		curFields[k] . components[i] . steps = 0;
	    }
	}
	free(curFields[k] . components);
    }
    curFields[k] . nComponents = 0;
    curFields[k] . components = (Component *) 0;
}

ObjPtr MakeDatasetCurData(dataset)
ObjPtr dataset;
/*Makes the dataset's CURDATA.  Returns ObjTrue if it had to make it*/
{
    ObjPtr data = NULLOBJ, curData = NULLOBJ;
    ObjPtr retVal;

    retVal = ObjFalse;

    data = GetVar(dataset, DATA);
    if (data)
    {
	curData = GetVar(data, CURDATA);
	if (!curData)
	{
	    /*For old-fashioned objects that don't know about CURDATA*/
	    curData = data;
	}
    }
    else
    {
	curData = data;
    }
    if (!curData)
    {
	ReportError("MakeDatasetCurData", "CURDATA of DATA is missing!");
	retVal = ObjTrue;
	return retVal;
    }
    SetVar(dataset, CURDATA, curData);
    return ObjTrue;
}

ObjPtr MakeDatasetTimesteps(dataset)
ObjPtr dataset;
/*Makes the dataset's TIMESTEPS.  Returns ObjTrue if it had to make it*/
{
    ObjPtr data = NULLOBJ, timesteps = NULLOBJ;

    data = GetVar(dataset, DATA);
    if (data)
    {
	MakeVar(data, TIMESTEPS);
	timesteps = GetVar(data, TIMESTEPS);
    }
    else
    {
	SetVar(dataset, TIMESTEPS, NULLOBJ);
	return ObjTrue;
    }
    SetVar(dataset, TIMESTEPS, timesteps);
    return ObjTrue;
}

#ifdef PROTO
real *GetComponentDataPtr(int whichField, int whichComponent)
#else
real *GetComponentDataPtr(whichField, whichComponent)
int whichField;
int whichComponent;
#endif
/*Returns a pointer to the data of component in whichField.
  This is VERY DANGEROUS and should only be called from file
  readers after having done New*Dataset and before the next time
  through the event loop!*/
{
    return curFields[whichField] . components[whichComponent] . data;
}

#ifdef PROTO
Bool SampleSpatComponent(int dataField, int formField,
	int nResultComp, real resultComp[], int nSpatIndices, 
	real spatIndices[], Bool interpolate)
#else
Bool SampleSpatComponent(dataField, formField, nResultComp, resultComp, nSpatIndices, spatIndices, interpolate)
int dataField, formField;
int nResultComp;
real resultComp[];
int nSpatIndices;
real spatIndices[];
Bool interpolate;
#endif
/*Samples a field dataField defined over formField using nSpatIndices spatIndices and puts
  the result in resultComp components*/
{
    long *topIndices;
    real *rTopIndices;
    long indexBuffer[10];
    real rIndexBuffer[10];
    int nTopIndices;
    int k;

    nTopIndices = curFields[formField] . topDim;
    if (nTopIndices != curFields[dataField] . topDim)
    {
	ReportError("SampleSpatComponent", "Topological dimension mismatch");
	return false;
    }

    if (nResultComp != curFields[dataField] . nComponents)
    {
	ReportError("SampleSpatComponent", "Wrong number of result components");
	return false;
    }

    if (nTopIndices <= 10)
    {
	topIndices = &(indexBuffer[0]);
	rTopIndices = &(rIndexBuffer[0]);
    }
    else
    {
	topIndices = malloc(nTopIndices * sizeof(long));
	rTopIndices = malloc(nTopIndices * sizeof(real));
    }
    SampleToTopIndex(formField, nTopIndices, rTopIndices, nSpatIndices, spatIndices);
    for (k = 0; k < nTopIndices; ++k)
    {
	topIndices[k] = rTopIndices[k] + 0.5;
    }

    for (k = 0; k < nResultComp; ++k)
    {
	if (interpolate)
	{
	    resultComp[k] = InterpolateFieldComponent(dataField, k, rTopIndices);
	}
	else
	{
	    resultComp[k] = SelectFieldComponent(dataField, k, topIndices);
	}
    }

    if (nTopIndices > 10)
    {
	free(topIndices);
	free(rTopIndices);
    }

    return true;
}

#ifdef PROTO
real SampleSpatScalar(int dataField, int formField,
	int nSpatIndices, 
	real spatIndices[], Bool interpolate)
#else
real SampleSpatScalar(dataField, formField, nSpatIndices, spatIndices, interpolate)
int dataField, formField;
int nSpatIndices;
real spatIndices[];
Bool interpolate;
#endif
/*Samples a scalar field dataField defined over formField using nSpatIndices
  spatIndices and returns the result
*/
{
    real retVal;
    long *topIndices;
    real *rTopIndices;
    long indexBuffer[10];
    real rIndexBuffer[10];
    int nTopIndices;
    int k;

    nTopIndices = curFields[formField] . topDim;
    if (nTopIndices != curFields[dataField] . topDim)
    {
	ReportError("SampleSpatScalar", "Topological dimension mismatch");
	return false;
    }

    if (nTopIndices <= 10)
    {
	topIndices = &(indexBuffer[0]);
	rTopIndices = &(rIndexBuffer[0]);
    }
    else
    {
	topIndices = malloc(nTopIndices * sizeof(long));
	rTopIndices = malloc(nTopIndices * sizeof(real));
    }
    SampleToTopIndex(formField, nTopIndices, rTopIndices, nSpatIndices, spatIndices);
    for (k = 0; k < nTopIndices; ++k)
    {
	topIndices[k] = rTopIndices[k] + 0.5;
    }


	if (interpolate)
	{
	    retVal = InterpolateFieldScalar(dataField, rTopIndices);
	}
	else
	{
	    retVal = SelectFieldScalar(dataField, topIndices);
	}

    if (nTopIndices > 10)
    {
	free(topIndices);
	free(rTopIndices);
    }

    return retVal;
}

#ifdef PROTO
void CalcBucketIndices(long indices[], int whichField, int nSamples, real *samples)
#else
void CalcBucketIndices(indices, whichField, nSamples, samples)
long indices[];
int whichField;
int nSamples;
real *samples;
#endif
/*Calculates the bucket indices of samples within whichField.  Indices
  must be an array of size == number of components in whichField*/
{
    long k;
    real div;

    for (k = 0; k < curFields[whichField] . nComponents; ++k)
    {
	if (k >= nSamples)
	{
	    /*Higher dimensionality, don't change offset*/
	    indices[k] = 0;
	}    
	else if (samples[k] < curFields[whichField] . bucketMinMax[k][0])
	{
	    /*Less than min, it's just zero.*/
	    indices[k] = 0;
	}
	else
	{
	    long i;

	    /*Calculate one of the divisions*/
	    div =  (curFields[whichField] . bucketMinMax[k][1] -
			curFields[whichField] . bucketMinMax[k][0]) /
			curFields[whichField] . bucketDim[k];
	    i = (samples[k] - curFields[whichField] . bucketMinMax[k][0]) / div;
	    if (i >= curFields[whichField] . bucketDim[k])
	    {
		i = curFields[whichField] . bucketDim[k] - 1;
	    }
	    indices[k] = i;
	}
    }
}

#ifdef PROTO
long CalcBucketOffset(int whichField, long indices[])
#else
long CalcBucketOffset(whichField, indices)
int whichField;
long indices[]
#endif
/*Calculates the bucket offset from indices within whichField*/
{
    long k;
    real div;
    long offset;

    offset = 0;

    for (k = 0; k < curFields[whichField] . nComponents; ++k)
    {
	if (k > 0)
	{
	    offset *= curFields[whichField] . bucketDim[k - 1];
	}
	offset += indices[k];
    }

    return offset;
}

#ifdef PROTO
static void FillBuckets(int whichField, int whichIndex, int maxNIndices, long *dimensions, long *indices)
#else
static void FillBuckets(whichField, whichIndex, maxNIndices, dimensions, indices)
int whichField;
int whichIndex;
int maxNIndices;
long *dimensions;
long *indices;
#endif
/*Fills the buckets of whichField using a
  powerset technique over indices.  Stuffs them into min and max*/
{
    if (whichIndex >= maxNIndices)
    {
	/*Bottomed out.  Spit out the sample*/
	long offset;
	long i, k;
	real *samples;
	long *bucketIndices;
	register Bucket *bucket;

	/*Calculate the offset*/
	samples = (real *) malloc(curFields[whichField] . nComponents * sizeof(real));
	bucketIndices = (long *) malloc(curFields[whichField] . nComponents * sizeof(long));

	for (k = 0; k < curFields[whichField] . nComponents; ++k)
	{
	    samples[k] = SelectFieldComponent(whichField, k, indices);
	}
	CalcBucketIndices(bucketIndices, whichField,
			curFields[whichField] . nComponents,
			samples);
	offset = CalcBucketOffset(whichField, bucketIndices);

	/*OK, got the offset.  Look at the bucket.*/
	bucket = &(curFields[whichField] . buckets[offset]);

	/*Go to the next slot on the bucket*/
	if (0 == (bucket -> nIndices % ADDBUCKETS))
	{
	    /*Must allocate some new buckets*/
	    if (bucket -> nIndices)
	    {
		bucket -> indices = (long *) realloc(bucket -> indices,
		    (ADDBUCKETS + bucket -> nIndices) *
		    sizeof(long) * maxNIndices);
		bucket -> samples = (real *) realloc(bucket -> samples,
		    (ADDBUCKETS + bucket -> nIndices) *
		    sizeof(real) * curFields[whichField] . nComponents);
	    }
	    else
	    {
		bucket -> indices = (long *) malloc(
		    ADDBUCKETS *
		    sizeof(long) * maxNIndices);
		bucket -> samples = (real *) malloc(
		    ADDBUCKETS *
		    sizeof(real) * curFields[whichField] . nComponents);
	    }
	}

	/*Insert this series of indices, and the sample*/
	i = bucket -> nIndices * maxNIndices;
	for (k = 0; k < maxNIndices; ++k)
	{
	    bucket -> indices[i] = indices[k];
	    ++i;
	}
	i = bucket -> nIndices * curFields[whichField] . nComponents;
	for (k = 0; k < curFields[whichField] . nComponents; ++k)
	{
	    bucket -> samples[i] = samples[k];
	    ++i;
	}

	/*Advance pointer*/
	++(bucket -> nIndices);

	/*No more samples needed*/
	free(bucketIndices);
	free(samples);
    }
    else
    {
	int maxIndex, k;

	maxIndex = dimensions[whichIndex];
	if (maxIndex > 0)
	{
	    /*It's real, do it*/
	    for (k = 0; k < maxIndex; ++k)
	    {
		indices[whichIndex] = k;
		FillBuckets(whichField, whichIndex + 1, maxNIndices, dimensions, indices);
	    }
	}
	else
	{
	    /*Descend nude*/
	    FillBuckets(whichField, whichIndex + 1, maxNIndices, dimensions, indices);
	}
    }
}

#ifdef PROTO
static int MaxComponentIndices(int whichField, int whichComp)
#else
static int MaxComponentIndices(whichField, whichComp)
int whichField;
int whichComp;
#endif
/*Returns the maximum number of indices in whichComp of whichField*/
{
    int k;
    int maxNIndices = -1;

    for (k = 0; k < curFields[whichField] . components[whichComp] . nIndices; ++k)
    {
	if (curFields[whichField] . components[whichComp] . indices[k] > maxNIndices)
	{
	    maxNIndices = curFields[whichField] . components[whichComp] . indices[k];
	}
    }
    return maxNIndices + 1;
}

#ifdef PROTO
static void MakeBuckets(int whichField)
#else
static void MakeBuckets(whichField)
int whichField;
#endif
/*Makes the buckets for whichField*/
{
    TwoReals *minMax;
    int whichComp;
    int k;

    if (!curFields[whichField] . bucketMinMax)
    {
	/*No buckets yet.  Make some*/
	int nTraversalDims;
	long *traversalDims;
	long *traversalIndices;

	long bucketDim;
	long nBucketDims;
	long *bucketDims;

	long k;
	long nElements;
	long nBuckets;
	Bucket *buckets;
	real sample;
	int whichComp, whichDim;

	minMax = (TwoReals *) malloc(curFields[whichField] . nComponents * sizeof(TwoReals));
	curFields[whichField] . bucketMinMax = minMax;

	nTraversalDims = CountTraversalDims(whichField);
	if (!nTraversalDims)
	{
	    return;
	}
	traversalDims = (long *) malloc(nTraversalDims * sizeof(long));
	GetTraversalDims(whichField, traversalDims);
	traversalIndices = (long *) malloc(nTraversalDims * sizeof(long));

	/*Go through all the samples, setting up the minmax*/
	for (whichComp = 0; whichComp < curFields[whichField] . nComponents; ++whichComp)
	{
	    /*Prepare minMax for expansion*/
	    minMax[whichComp][0] = PLUSINF;
	    minMax[whichComp][1] = MINUSINF;

	    /*Zero the indices*/
	    for (k = 0; k < nTraversalDims; ++k)
	    {
		traversalIndices[k] = 0;
	    }

	    do
	    {
		/*Sample the location at a point*/
		sample = SelectFieldComponent(whichField, 0, traversalIndices);

		if (sample != missingData)
		{
		    if (sample < minMax[whichComp][0]) minMax[whichComp][0] = sample;
		    if (sample > minMax[whichComp][1]) minMax[whichComp][1] = sample;
		}

		/*Advance to next point*/
		for (whichDim = 0; whichDim < nTraversalDims; ++whichDim)
		{
		    if (traversalDims[whichDim] > 0)
		    {
			if ((++traversalIndices[whichDim]) >= traversalDims[whichDim])
			{
			    traversalIndices[whichDim] = 0;
			}
			else
			{
			    break;
			}
		    }
		}
	    } while (whichDim < nTraversalDims); /*Break is based on advance*/
	}

	/*Start figuring how to traverse field.*/
	curFields[whichField] . nBucketIndices = nTraversalDims;

	/*Calculate the number of elements total*/
	nElements = 1;
	for (k = 0; k < nTraversalDims; ++k)
	{
	    if (traversalDims[k] >= 0)
	    {
		nElements *= traversalDims[k];
	    }
	}

	/*Now that we have the number of elements, calculate a good # of buckets.*/
	nBuckets = nElements * BUCKETSIZEFACTOR;

	/*If it's too many to be manageable, back it off*/
	if (nBuckets > MAXBUCKETSIZE)
	{
	    nBuckets = MAXBUCKETSIZE;
	}

	/*Find a bucket dimension*/
	bucketDim = (int) pow((double) nBuckets, 1.0 / (double) (curFields[whichField] . nComponents));
	if (bucketDim <= 0) ++bucketDim;

#ifdef DEBUG
	printf("Side dimension = %d\n", bucketDim);
	printf("Target # buckets = %d\n", nBuckets);
#endif
	
	/*Make bucket dimensions*/
	bucketDims = (long *) malloc((curFields[whichField] . nComponents) * sizeof(long));
	curFields[whichField] . bucketDim = bucketDims;
	for (k = 0; k < curFields[whichField] . nComponents; ++k)
	{
	    bucketDims[k] = bucketDim;
	}

	/*Create and zero the actual bucket array*/
	buckets = (Bucket *) malloc(((long) nBuckets) * sizeof(Bucket));
	curFields[whichField] . buckets = buckets;
	for (k = 0; k < nBuckets; ++k)
	{
	    buckets[k] . indices = 0;
	    buckets[k] . samples = 0;
	    buckets[k] . nIndices = 0;
	}

	/*Fill the buckets*/
	for (k = 0; k < nTraversalDims; ++k)
	{
	    traversalIndices[k] = 0;
	}
	FillBuckets(whichField, 0, nTraversalDims, traversalDims, traversalIndices);

#ifdef DEBUG
	/*Print out bucket statistics*/
	{
	    long nFilled, maxFilled;
	    nFilled = 0;
	    maxFilled = 0;
	    for (k = 0; k < nBuckets; ++k)
	    {
		if (buckets[k] . nIndices)
		{
		    ++nFilled;
		    if (buckets[k] . nIndices > maxFilled)
			maxFilled = buckets[k] . nIndices;
		}
	    }
	    printf("%ld out of %ld buckets filled (%lg%%)\n",
		nFilled, nBuckets, ((double) nFilled) / ((double) nBuckets) * 100.0);
	    printf("Maximum bucket size = %ld\n", maxFilled);
	}
#endif

	free(traversalIndices);
	free(traversalDims);
    }
}

#ifdef PROTO
void SBS2(int whichField, long powerIndices[], long centerIndices[],
	  long shell, int whichIndex, long bestIndex[],
	  real *bestGuess, real sample[])
#else
void SBS2(whichField, powerIndices, centerIndices,
	  shell, whichIndex, bestIndex,
	  bestGuess, sample)
int whichField;
long powerIndices[];
long centerIndices[];
long shell;
int whichIndex;
long bestIndex[];
real *bestGuess;
real sample[];
#endif
/*Powerset recursive helper for SearchBucketForSample.  Searches a shell 
  shell wide around centerIndices, comparing each sample with bestGuess,
  returning bestIndex if found.  whichIndex is the current index, used for
  the powerset, done through powerIndices*/
{
    long i, k;
    register int nComponents;

    nComponents = curFields[whichField] . nComponents;

    if (whichIndex >= nComponents)
    {
	real d2;	/*Distance squared*/
	real d;		/*Distance*/
	register Bucket *bucket;
	long offset;

	/*Bottomed out.  Test to see if we're in the right shell*/
	for (k = 0; k < nComponents; ++k)
	{
	    if (ABS(powerIndices[k] - centerIndices[k]) != shell)
	    {
		/*Not on this shell, so don't do anything.*/
		return;
	    }
	}

	/*Calculate the offset*/
	offset = CalcBucketOffset(whichField, powerIndices);

	/*Get the bucket*/
	bucket = &(curFields[whichField] . buckets[offset]);

	/*Search for better match*/
	for (k = 0; k < bucket -> nIndices; ++k)
	{
	    /*Calculate distance*/
	    d2 = 0.0;

	    for (i = 0; i < nComponents; ++i)
	    {
		if (bucket -> samples[k * nComponents + i] == missingData)
		{
		    break;
		}
		d = bucket -> samples[k * nComponents + i] - sample[k];
		d2 += d * d;
	    }

	    if (i >= nComponents)
	    {
		/*It completed a distance*/
		if (d2 < *bestGuess)
		{
		    /*Improved guess, copy index*/
		    *bestGuess = d2;
		    for (i = 0; i < curFields[whichField] . nBucketIndices; ++i)
		    {
			bestIndex[i] =
			    bucket -> indices[k * curFields[whichField] . nBucketIndices + i];
		    }
		}
	    }
	}
    }
    else
    {
	long min, max;

	/*Go down powerset*/
	min = centerIndices[whichIndex] - shell;
	if (min < 0)
	{
	    min = 0;
	}
	max = centerIndices[whichIndex] + shell;
	if (max >= curFields[whichField] . bucketDim[whichIndex])
	{
	    max = curFields[whichField] . bucketDim[whichIndex] - 1;
	}
	for (powerIndices[whichIndex] = min;
	     powerIndices[whichIndex] <= max;
	     ++powerIndices[whichIndex])
	{
	    SBS2(whichField, powerIndices, centerIndices,
		shell, whichIndex + 1, bestIndex,
		bestGuess, sample);
	}
    }
}

#ifdef PROTO
Bool SearchBucketForSample(int whichField, real indices[], 
	long bucketIndices[], real sample[])
#else
Bool SearchBucketForSample(whichField, indices, bucketIndices, sample)
int whichField;
real indices[];
long bucketIndices[];
real sample[];
#endif
/*Searches the buckets of whichField for sample and returns in indices.
  bucketIndices is the first guess and is assumed to be DISPOSABLE.
*/
{
    long *bestIndex;		/*Best index so far*/
    long *powerIndices;		/*Power indices*/
    long k;
    long shell, nShells;
    real bestGuess;		/*Best guess so far*/

    bestIndex = (long *) malloc(sizeof(long) * curFields[whichField] . nBucketIndices);
    powerIndices = (long *) malloc(sizeof(long) * curFields[whichField] . nComponents);

    /*Calculate the maximum number of shells to use*/
    nShells = 0;
    for (k = 0; k < curFields[whichField] . nComponents; ++k)
    {
	if (bucketIndices[k] > nShells)
	{
	    nShells = bucketIndices[k];
	}
	if (curFields[whichField] . bucketDim[k] - bucketIndices[k] > bucketIndices[k])
	{
	    nShells = curFields[whichField] . bucketDim[k] - bucketIndices[k];
	}
    }

    /*Now search for a match in progressively widening shells*/
    bestGuess = plusInf;

    for (shell = 0; shell < nShells; ++shell)
    {
	SBS2(whichField, powerIndices, bucketIndices, shell, 0, bestIndex, &bestGuess, sample);
	if (bestGuess < plusInf)
	{
	    break;
	}
    }

    /*Put the best index into indices*/
    for (k = 0; k < curFields[whichField] . nBucketIndices; ++k)
    {
	indices[k] = bestIndex[k];
    }

    free(powerIndices);
    free(bestIndex);

    return (bestGuess < plusInf) ? true : false;
}


#ifdef PROTO
Bool SampleToTopIndex(int whichField, int nIndices, real indices[], 
		int nComponents, real sample[])
#else
Bool SampleToTopIndex(whichField, nIndices, indices, nComponents, sample)
int whichField;
int nIndices;
real indices[];
int nComponents;
real sample[];
#endif
/*Finds and returns in indices the closest index in whichField to sample.
  The arrays had better be big enough*/
{
    int i, comp, k, index, myIndex;
    int indicesAssigned;
    Bool succeed;
    Bool inHere, otherFree, others;
    long *longIndices;
    long *bucketIndices;
    real retVal;

    /*Find the number of components, should be == size of sample*/
    if (nComponents < curFields[whichField] . nComponents)
    {
	ReportError("SampleToTopIndex", "Not enough components passed");
	return false;
    }

    longIndices = malloc(nIndices * sizeof(long));

    /*Make all indices free*/
    for (k = 0; k < nIndices; ++k)
    {
	indices[k] = -1.0;
	longIndices[k] = -1;
    }

    /*Number of indices which have been assigned, start at 0*/
    indicesAssigned = 0;

    do
    {
	succeed = false;
	/*Search for an index and component such that
	  1) The index is free
	  2) The index appears in the component
	  3) No other component is free*/

	for (i = 0; i < nIndices; ++i)
	{
	    /*See if it's free*/
	    if (longIndices[i] < 0)
	    {
		/*See if it appears in any components with free index*/
		for (comp = 0; comp < curFields[whichField] . nComponents; ++comp)
		{
		    inHere = others = otherFree = false;
		    for (k = 0;
			 k < curFields[whichField] . components[comp] . nIndices;
			 ++k)
		    {
			index = curFields[whichField] . components[comp] . indices[k];
			if (index == i)
			{
			    myIndex = k;
			    inHere = true;
			}
			else if (index >= 0)
			{
			    if (longIndices[index] < 0)
			    {
				otherFree = true;
			    }
			    else
			    {
				others = true;
			    }
			}
		    }

		    if (inHere && !otherFree)
		    {
			/*Found one which is free*/

			if (!others &&
			    0 == (curFields[whichField] . components[comp] . flags & (CF_MONOTONIC | CF_NOTMONOTONIC)))
			{
			    int testo;
			    real next, last;

			    testo = 0;
			    
			    /*Might be monotonic.  Let's check*/
			    longIndices[i] = 0;
			    last = SelectFieldComponent(whichField, comp, longIndices);
			    for (longIndices[i] = 1;
			         longIndices[i] < curFields[whichField] . components[comp] . dimensions[myIndex];
			         ++(longIndices[i]))
			    {
				next = SelectFieldComponent(whichField, comp, longIndices);
				if (next > last)
				{
				    testo |= 1;
				}
				else if (last > next)
				{
				    testo |= 2;
				}
				last = next;
			    }
			    if (testo < 3)
			    {
				curFields[whichField] . components[comp] . flags |= CF_MONOTONIC;
			    }
			    else
			    {
				curFields[whichField] . components[comp] . flags |= CF_NOTMONOTONIC;
			    }
			}

			/*Hey!  Found it!*/

			if (curFields[whichField] . components[comp] . flags & CF_MONOTONIC)
			{
			    long beg, end, mid;
			    real begFit, endFit, midFit, diff1, diff2;

			    beg = 0;
			    end = curFields[whichField] . components[comp] . dimensions[myIndex] - 1;
			    longIndices[i] = beg;
			    begFit = SelectFieldComponent(whichField, comp, longIndices);
			    longIndices[i] = end;
			    endFit = SelectFieldComponent(whichField, comp, longIndices);

			    if (begFit < endFit)
			    {
				/*Positive monotonic*/
				if (sample[comp] <= begFit)
				{
				    longIndices[i] = beg;
				    indices[i] = beg;
				}
				else if (sample[comp] >= endFit)
				{
				    longIndices[i] = end;
				    indices[i] = end;
				}
				else
				{
				    while (beg + 1 < end)
				    {
					mid = (beg + end) / 2;
					if (mid <= beg) mid = beg + 1;
					else if (mid >= end) mid = end - 1;
					longIndices[i] = mid;
					midFit = SelectFieldComponent(whichField, comp, longIndices);
					if (midFit > sample[comp])
					{
					    endFit = midFit;
					    end = mid;
					}
					else
					{
					    begFit = midFit;
					    beg = mid;
					}
				    }
				    diff1 = begFit - sample[comp];
				    diff1 = ABS(diff1);
				    diff2 = endFit - sample[comp];
				    diff2 = ABS(diff2);
				    longIndices[i] = (diff1 > diff2) ? end : beg;
				    indices[i] = ((beg * diff2) + (end * diff1)) / (diff1 + diff2);
				}
				succeed = true;
				++indicesAssigned;
			    }
			    else
			    {
				/*Negative monotonic*/
				if (sample[comp] >= begFit)
				{
				    longIndices[i] = beg;
				    indices[i] = beg;
				}
				else if (sample[comp] <= endFit)
				{
				    longIndices[i] = end;
				    indices[i] = end;
				}
				else
				{
				    while (beg + 1 < end)
				    {
					mid = (beg + end) / 2;
					if (mid >= beg) mid = beg + 1;
					else if (mid <= end) mid = end - 1;
					longIndices[i] = mid;
					midFit = SelectFieldComponent(whichField, comp, longIndices);
					if (midFit < sample[comp])
					{
					    endFit = midFit;
					    end = mid;
					}
					else
					{
					    begFit = midFit;
					    beg = mid;
					}
				    }
				    diff1 = begFit - sample[comp];
				    diff1 = ABS(diff1);
				    diff2 = endFit - sample[comp];
				    diff2 = ABS(diff2);
				    longIndices[i] = (diff1 > diff2) ? end : beg;
				    indices[i] = ((beg * diff2) + (end * diff1)) / (diff1 + diff2);
				}
				succeed = true;
				++indicesAssigned;
			    }
			}
			else
			{
			    /*It's not monotonic.  Do a linear search*/
			    long bestFit;
			    real bestDist;
			    real diff;

			    succeed = true;
			    ++indicesAssigned;

			    longIndices[i] = 0;
			    bestFit = 0;
			    diff = SelectFieldComponent(whichField, comp, longIndices) -
				sample[comp];
			    bestDist = ABS(diff);

			    for (longIndices[i] = 1;
				 longIndices[i] < curFields[whichField] . components[comp] . dimensions[myIndex];
				 ++(longIndices[i]))
			    {
				diff = SelectFieldComponent(whichField, comp, longIndices) -
				    sample[comp];
				diff = ABS(diff);
				if (diff < bestDist)
				{
				    bestDist = diff;
				    bestFit = longIndices[i];
				}
			    }
			    longIndices[i] = bestFit;
			    if (longIndices[i] <= 0)
			    {
				real sCenter, sPlus;
				/*At left, see if it's between 0 and 1*/

				sCenter = SelectFieldComponent(whichField, comp, longIndices);
				++longIndices[i];
				sPlus = SelectFieldComponent(whichField, comp, longIndices);
				--longIndices[i];
				if (sCenter <= sample[comp] <= sPlus)
				{
				    indices[i] =
					(sPlus - sample[comp]) * longIndices[i] +
					(sample[comp] - sCenter) * (longIndices[i] + 1);
				}
				else if (sPlus <= sample[comp] <= sCenter)
				{
				    indices[i] =
					(sCenter - sample[comp]) * longIndices[i] +
					(sample[comp] - sPlus) * (longIndices[i] + 1);
				}
				else
				{
				    indices[i] = longIndices[i];
				}
			    }
			    else if (longIndices[i] >= curFields[whichField] . components[comp] . dimensions[myIndex])
			    {
				real sCenter, sMinus;
				/*At right, see if it's between right and right - 1*/

				sCenter = SelectFieldComponent(whichField, comp, longIndices);
				--longIndices[i];
				sMinus = SelectFieldComponent(whichField, comp, longIndices);
				++longIndices[i];
				if (sCenter <= sample[comp] <= sMinus)
				{
				    indices[i] =
					(sMinus - sample[comp]) * longIndices[i] +
					(sample[comp] - sCenter) * (longIndices[i] - 1);
				}
				else if (sMinus <= sample[comp] <= sCenter)
				{
				    indices[i] =
					(sCenter - sample[comp]) * longIndices[i] +
					(sample[comp] - sMinus) * (longIndices[i] - 1);
				}
				else
				{
				    indices[i] = longIndices[i];
				}
			    }
			    else
			    {
				real sCenter, sMinus, sPlus;
				/*It's in the center*/
				sCenter = SelectFieldComponent(whichField, comp, longIndices);
				--longIndices[i];
				sMinus = SelectFieldComponent(whichField, comp, longIndices);
				++longIndices[i];
				++longIndices[i];
				sPlus = SelectFieldComponent(whichField, comp, longIndices);
				--longIndices[i];
				if (sCenter <= sample[comp] <= sPlus)
				{
				    indices[i] =
					(sPlus - sample[comp]) * longIndices[i] +
					(sample[comp] - sCenter) * (longIndices[i] + 1);
				}
				else if (sPlus <= sample[comp] <= sCenter)
				{
				    indices[i] =
					(sCenter - sample[comp]) * longIndices[i] +
					(sample[comp] - sPlus) * (longIndices[i] + 1);
				}
				else if (sCenter <= sample[comp] <= sMinus)
				{
				    indices[i] =
					(sMinus - sample[comp]) * longIndices[i] +
					(sample[comp] - sCenter) * (longIndices[i] - 1);
				}
				else if (sMinus <= sample[comp] <= sCenter)
				{
				    indices[i] =
					(sCenter - sample[comp]) * longIndices[i] +
					(sample[comp] - sMinus) * (longIndices[i] - 1);
				}
				else
				{
				    indices[i] = longIndices[i];
				}
			    }
			}
		    }
		}
	    }
	}
    } while (succeed && indicesAssigned < nIndices);

    free(longIndices);

    if (indicesAssigned >= nIndices)
    {
	/*Success*/
	return true;
    }

    /*Fell through, so have to use bucket approach*/
    MakeBuckets(whichField);

    /*Calculate initial guess*/
    bucketIndices = (long *) malloc(sizeof(long) * curFields[whichField] . nComponents);

    CalcBucketIndices(bucketIndices, whichField,
			curFields[whichField] . nComponents,
			sample);

    retVal = SearchBucketForSample(whichField, indices, bucketIndices, sample);
    free(bucketIndices);
    return retVal;
}

#ifdef PROTO
real SelectFieldComponent(int whichField, int whichComponent, long indices[])
#else
real SelectFieldComponent(whichField, whichComponent, indices)
int whichField;
int whichComponent;
long indices[];
#endif
/*Selects a field component given indices into it.  There had better be enough
  indices*/
{
    register int k;
    register long offset = 0;		/*Offset into the data*/
    register Component *component;	/*The current component*/

    /*Get the current component*/
    component = &(curFields[whichField] . components[whichComponent]);

    /*Calculate the offset*/
    for (k = 0; k < component -> nIndices; ++k)
    {
	register int whichIndex;

	whichIndex = component -> indices[k];

	if (whichIndex >= 0)
	{
	    offset += component -> steps[k] * indices[whichIndex];
	}
    }
    return component -> data[offset];
}

#ifdef PROTO
real SelectFieldScalar(int whichField, long indices[])
#else
real SelectFieldScalar(whichField, indices)
int whichField;
long indices[];
#endif
/*Selects a scalar value from whichField using indices*/
{
    if (curFields[whichField] . nComponents < 2)
    {
	return SelectFieldComponent(whichField, 0, indices);
    }
    else
    {
	real squareSum, d;
	int whichComponent;

	squareSum = 0.0;
	for (whichComponent = 0;
		whichComponent < curFields[whichField] . nComponents;
		++whichComponent)
	{
	    d = SelectFieldComponent(whichField, whichComponent, indices);

	    if (d == missingData)
	    {
		return missingData;
	    }
	    squareSum += SQUARE(d);
	}
	return rsqrt(squareSum);
    }
}

#ifdef PROTO
real InterpolateFieldComponent(int whichField, int whichComponent, real indices[])
#else
real InterpolateFieldComponent(whichField, whichComponent, indices)
int whichField;
int whichComponent;
real indices[];
#endif
/*Interpolates a component at a set of real topological indices*/
{
    register int i, k;
    register long test;
    register Component *component;	/*The current component*/
    real retVal;
    register int whichIndex;
    register real index;

    component = &(curFields[whichField] . components[whichComponent]);

    for (k = 0; k < component -> nIndices; ++k)
    {
	whichIndex = component -> indices[k];

	if (whichIndex >= 0)
	{
	    index = indices[whichIndex];

	    test = index;
	    if (index != (real) test)
	    {
		/*It's non integral*/

		if (index < -0.5 ||
		    index > ((real) component -> dimensions[whichIndex]) - 0.5)
		{
		    /*It's way outside*/
		    return missingData;
		}
		else
		{
		    /*It's between two vertices*/
		    long i1, i2;
		    real v1, v2;
		    real *tempIndices;
		    int copyIndex;

		    i1 = test;
		    i2 = i1 + 1;

		    tempIndices = (real *) malloc(MaxComponentIndices(whichField, whichComponent) * sizeof(real));

		    /*Copy all the indices*/
		    for (i = 0; i < component -> nIndices; ++i)
		    {
			copyIndex = component -> indices[i];
			if (copyIndex >= 0)
			{
			    tempIndices[copyIndex] = indices[copyIndex];
			}
		    }

		    /*Now get the two values for the two indices*/
		    tempIndices[whichIndex] = i1;
		    v1 = InterpolateFieldComponent(whichField, whichComponent, tempIndices);
		    tempIndices[whichIndex] = i2;
		    v2 = InterpolateFieldComponent(whichField, whichComponent, tempIndices);

		    free(tempIndices);

		    /*Now return an interpolation*/
		    if (v1 == missingData)
		    {
			if (v2 == missingData)
			{
			    return missingData;
			}
			else
			{
			    if (i2 - index <= index - i1)
			    {
				return v2;
			    }
			    else
			    {
				return missingData;
			    }
			}
		    }
		    else if (v2 == missingData)
		    {
			if (index - i1 <= i2 - index)
			{
			    return v1;
			}
			else
			{
			    return missingData;
			}
		    }
		    else
		    {
			return (index - i1) * v2 + (i2 - index) * v1;
		    }
		}
	    }
	}
    }

    /*If we fell through here, it must be smack dab on a vertex.  Just get it*/
    {
	long *tempIndices;

	tempIndices = (long *) malloc(MaxComponentIndices(whichField, whichComponent) * sizeof(long));
	for (k = 0; k < component -> nIndices; ++k)
	{
	    whichIndex = component -> indices[k];
	    if (whichIndex >= 0)
	    {
	        tempIndices[whichIndex] = indices[whichIndex];
	    }
	}
	retVal = SelectFieldComponent(whichField, whichComponent, tempIndices);
	free(tempIndices);
	return retVal;
    }
}

#ifdef PROTO
real InterpolateFieldScalar(int whichField, real indices[])
#else
real InterpolateFieldScalar(whichField, indices)
int whichField;
real indices[];
#endif
/*Interpolates between scalar values*/
{
    if (curFields[whichField] . nComponents < 2)
    {
	return InterpolateFieldComponent(whichField, 0, indices);
    }
    else
    {
	real squareSum;
	int k;
	real d;

	squareSum = 0.0;
	for (k = 0; k < curFields[whichField] . nComponents; ++k)
	{
	    d = InterpolateFieldComponent(whichField, k, indices);
	    if (d == missingData) return missingData;
	    squareSum += SQUARE(d);
	}

	return rsqrt(squareSum);
    }
}

long GetComponentOffset(whichField, whichComponent, indices)
int whichField;
int whichComponent;
long indices[];
/*Get the offset into a field component given indices into it.  There had better be enough
  indices*/
{
    register int k;
    register long offset = 0;		/*Offset into the data*/
    register Component *component;	/*The current component*/

    /*Get the current component*/
    if (whichComponent < 0 || 
	whichComponent > curFields[whichField] . nComponents)
    {
	return 0;
    }
    component = &(curFields[whichField] . components[whichComponent]);

    /*Calculate the offset*/
    for (k = 0; k < component -> nIndices; ++k)
    {
	register int whichIndex;
	register long index;
	
	whichIndex = component -> indices[k];

	if (whichIndex >= 0)
	{
	    index = indices[whichIndex];
	    if (index < 0 || index > component -> dimensions[k])
	    {
		return 0;
	    }
	}

	offset += component -> steps[k] * index;
    }
    return offset;
}


#ifdef PROTO
void PutFieldComponent(int whichField, int whichComponent, long indices[], real data)
#else
void PutFieldComponent(whichField, whichComponent, indices, data)
int whichField;
int whichComponent;
long indices[];
real data;
#endif
/*Puts data into a field component given indices into it.  There had better be enough
  indices*/
{
    register int k;
    register long offset = 0;		/*Offset into the data*/
    register Component *component;	/*The current component*/

    /*Get the current component*/
    if (whichComponent < 0 || 
	whichComponent > curFields[whichField] . nComponents)
    {
	return;
    }
    component = &(curFields[whichField] . components[whichComponent]);

    /*Calculate the offset*/
    for (k = 0; k < component -> nIndices; ++k)
    {
	register int whichIndex;
	register long index;
	
	whichIndex = component -> indices[k];

	if (whichIndex >= 0)
	{
	    index = indices[whichIndex];
	    if (index < 0 || index > component -> dimensions[k])
	    {
		return;
	    }

	    offset += component -> steps[k] * index;
	}
    }
    component -> data[offset] = data;
}

void StuffScalarIJSquare(dest, whichField, indices)
real dest[2][2];
int whichField;
long indices[];
/*Stuffs the square in the ij plane with its lower left corner defined by
  indices into dest.  Takes data from the specified field, which can be 
  scalar, in which case the scalar value is used, or vector, in which case
  the absolute value is used.  

  There had better be enough indices.  This does not check
  validity of the indices*/
{
    int k;
    register long offset;		/*Offset into the data*/
    register long iStep = 0, jStep = 0;	/*Step in the i and j directions*/
    register Component *component;	/*The current component*/

    if (curFields[whichField] . nComponents < 2)
    {
	/*Get the current component*/
	component = &(curFields[whichField] . components[0]);

	offset = 0;
	/*Calculate the offset*/
	for (k = 0; k < component -> nIndices; ++k)
	{
	    register int whichIndex;
	    register long index;
	    
	    whichIndex = component -> indices[k];
	    
	    if (whichIndex >= 0)
	    {
		index = indices[whichIndex];
		if (whichIndex == 0)
		{
		    /*I component*/
		    iStep = component -> steps[k];
		    offset += iStep * index;
		}
		else if (whichIndex == 1)
		{
		    /*J component*/
		    jStep = component -> steps[k];
		    offset += jStep * index;
		}
		else
		{
		    offset += component -> steps[k] * index;
		}
	    }
	}
	dest[0][0] = component -> data[offset];
	dest[1][0] = component -> data[offset + iStep];
	dest[0][1] = component -> data[offset + jStep];
	dest[1][1] = component -> data[offset + iStep + jStep];
    }
    else
    {
	int whichComponent;
	real squareSum[2][2];
	Bool missing[2][2];
	real d;

	squareSum[0][0] = 0.0;
	squareSum[0][1] = 0.0;
	squareSum[1][0] = 0.0;
	squareSum[1][1] = 0.0;

	missing[0][0] = false;
	missing[0][1] = false;
	missing[1][0] = false;
	missing[1][1] = false;

	for (whichComponent = 0;
		whichComponent < curFields[whichField] . nComponents;
		++whichComponent)
	{
	    /*Get the current component*/
	    component = &(curFields[whichField] . components[whichComponent]);

	    /*Calculate the offset*/
	    offset = 0;
	    iStep = 0;
	    jStep = 0;
	    for (k = 0; k < component -> nIndices; ++k)
	    {
		register int whichIndex;
		register long index;
	    
		whichIndex = component -> indices[k];

		if (whichIndex >= 0)
		{
		    index = indices[whichIndex];
		    if (whichIndex == 0)
		    {
			/*I component*/
			iStep = component -> steps[k];
			offset += iStep * index;
		    }
		    else if (whichIndex == 1)
		    {
			/*J component*/
			jStep = component -> steps[k];
			offset += jStep * index;
		    }
		    else
		    {
			offset += component -> steps[k] * index;
		    }
		}
	    }

	    d = component -> data[offset];
	    if (d == missingData)
	    {
		missing[0][0] = true;
	    }
	    else
	    {
		squareSum[0][0] += SQUARE(d);
	    }
	    
	    d = component -> data[offset + iStep];
	    if (d == missingData)
	    {
		missing[1][0] = true;
	    }
	    else
	    {
		squareSum[1][0] += SQUARE(d);
	    }
	    
	    d = component -> data[offset + jStep];
	    if (d == missingData)
	    {
		missing[0][1] = true;
	    }
	    else
	    {
		squareSum[0][1] += SQUARE(d);
	    }
	    
	    d = component -> data[offset + iStep + jStep];
	    if (d == missingData)
	    {
		missing[1][1] = true;
	    }
	    else
	    {
		squareSum[1][1] += SQUARE(d);
	    }
	}
	dest[0][0] = (missing[0][0]) ? (missingData) : rsqrt(squareSum[0][0]);
	dest[1][0] = (missing[1][0]) ? (missingData) : rsqrt(squareSum[1][0]);
	dest[0][1] = (missing[0][1]) ? (missingData) : rsqrt(squareSum[0][1]);
	dest[1][1] = (missing[1][1]) ? (missingData) : rsqrt(squareSum[1][1]);
   }
}

void StuffIJSquare(dest, whichField, whichComponent, indices)
real dest[2][2];
int whichField;
int whichComponent;
long indices[];
/*Stuffs the square in the ij plane with its lower left corner defined by
  indices into dest.  Takes data from the specified field, which had better
  have component.  There had better be enough indices.  This does not check
  validity of the indices*/
{
    int k;
    register long offset = 0;		/*Offset into the data*/
    register long iStep = 0, jStep = 0;	/*Step in the i and j directions*/
    register Component *component;	/*The current component*/
    register real *data;

    /*Get the current component*/
    component = &(curFields[whichField] . components[whichComponent]);

    /*Calculate the offset*/
    for (k = 0; k < component -> nIndices; ++k)
    {
	register int whichIndex;
	register long index;
	
	whichIndex = component -> indices[k];

	if (whichIndex >= 0)
	{
	    index = indices[whichIndex];
	    if (whichIndex == 0)
	    {
		/*I component*/
		iStep = component -> steps[k];
		offset += iStep * index;
	    }
	    else if (whichIndex == 1)
	    {
		/*J component*/
		jStep = component -> steps[k];
		offset += jStep * index;
	    }
	    else
	    {
		offset += component -> steps[k] * index;
	    }
	}
    }
    dest[0][0] = component -> data[offset];
    dest[1][0] = component -> data[offset + iStep];
    dest[0][1] = component -> data[offset + jStep];
    dest[1][1] = component -> data[offset + iStep + jStep];
}

Bool SetCurField(whichField, field)
int whichField;
ObjPtr field;
/*Sets current field whichField to field.  Returns true iff it did*/
{
    FuncTyp method;

    /*Clean the field before setting it*/
    CleanCurField(whichField);

    curFields[whichField] . topDim = GetTopDim(field);
    curFields[whichField] . fieldObj = field;

    method = GetMethodSurely("SetCurField", field, REGISTERFIELD);
    if (method)
    {
	ObjPtr result;
	result = (*method)(field, whichField);
	if (IsTrue(result))
	{
	    /*See if it's appropriate to set the idiot flag*/
	    long whichComponent;
	    register Component *component;
	    int k;

	    curFields[whichField] . idiotFlag = true;

	    for (whichComponent = 0;
		 whichComponent < curFields[whichField] . nComponents;
		 ++whichComponent)
	    {
		component = &(curFields[whichField] . components[whichComponent]);
		if (curFields[whichField] . topDim != 1)
		{
		    curFields[whichField] . idiotFlag = false;
		    break;
		}

		if (curFields[whichField] . topDim !=
		    component -> nIndices)
		{
		    curFields[whichField] . idiotFlag = false;
		    break;
		}

		for (k = component -> nIndices - 1; k >= 0; --k)
		{
		    if (component -> indices[k] != k)
		    {
			curFields[whichField] . idiotFlag = false;
			break;
		    }
		    if (component -> steps[k] !=
			((k == component -> nIndices - 1) ?
			    1 :
			    component -> dimensions[k + 1] * component -> steps[k + 1]))
		    {
			curFields[whichField] . idiotFlag = false;
			break;
		    }
		}
	    }

	    return true;
	}
	else
	{
	    return false;
	}
    }
    else
    {
	return false;
    }
}

Bool SetCurForm(whichField, field)
int whichField;
ObjPtr field;
/*Sets current field whichField to the data form of field*/
{
    FuncTyp method;

    CleanCurField(whichField);

    curFields[whichField] . topDim = GetTopDim(field);
    curFields[whichField] . fieldObj = GetVar(field, DATAFORM);

    method = GetMethodSurely("SetCurField", field, REGISTERFORM);
    if (method)
    {
	ObjPtr result;
	result = (*method)(field, whichField);
	if (IsTrue(result))
	{
	    /*See if it's appropriate to set the idiot flag*/
	    long whichComponent;
	    register Component *component;
	    int k;

	    curFields[whichField] . idiotFlag = true;

	    for (whichComponent = 0;
		 whichComponent < curFields[whichField] . nComponents;
		 ++whichComponent)
	    {
		component = &(curFields[whichField] . components[whichComponent]);
		if (curFields[whichField] . topDim !=
		    component -> nIndices)
		{
		    curFields[whichField] . idiotFlag = false;
		    break;
		}

		if (curFields[whichField] . topDim != 1)
		{
		    curFields[whichField] . idiotFlag = false;
		    break;
		}

		for (k = component -> nIndices - 1; k >= 0; --k)
		{
		    if (component -> indices[k] != k)
		    {
			curFields[whichField] . idiotFlag = false;
			break;
		    }
		    if (component -> steps[k] !=
			((k == component -> nIndices - 1) ?
			    1 :
			    component -> dimensions[k + 1] * component -> steps[k + 1]))
		    {
			curFields[whichField] . idiotFlag = false;
			break;
		    }
		}
	    }

	    return true;
	}
	else
	{
	    return false;
	}
     }
    return false;
}

ObjPtr RegisterDatasetForm(dataset, whichField)
int whichField;
ObjPtr dataset;
/*Registers the form of a dataset in whichField*/
{
    ObjPtr dataForm;

    dataForm = GetVar(dataset, DATAFORM);
    if (dataForm)
    {
	return SetCurField(whichField, dataForm) ? ObjTrue : ObjFalse;
    }
    ReportError("RegisterDatasetForm", "No DATAFORM");
    return ObjFalse;
}

Bool RegisterComponent(whichField, whichComponent, field)
int whichField, whichComponent;
ObjPtr field;
/*Registers field as whichComponent of whichField*/
{
    FuncTyp method;

    if (whichComponent < 0 || whichComponent >= curFields[whichField] . nComponents)
    {
	ReportError("RegisterComponent", "Component number out of range");
	return false;
    }

    method = GetMethodSurely("RegisterComponent", field, REGISTERCOMP);
    if (method)
    {
	ObjPtr result;
	result = (*method)(field, whichField, whichComponent);
	return IsTrue(result) ? true : false;
    }
    return false;
}

static ObjPtr globalDataset;		/*Global dataser for ForAllVisWindows*/
static void ChangeWindowForDataset(window)
WinInfoPtr window;
/*Changes a given window upon the knowledge that the dataset has changed*/
{
    ObjPtr space;
    space = FindSpace(window);
    if (space)
    {
	ObjPtr clock;
	ImInvalid((ObjPtr) window);
	clock = GetVar(space, CLOCK);
	if (clock)
	{
	    WinInfoPtr dialogExists;
	    dialogExists = DialogExists((WinInfoPtr) clock, NewString("Clock"));
	    if (dialogExists)
	    {
		ObjPtr timeControl;
		timeControl = GetVar((ObjPtr) dialogExists, TIMECONTROL);
		if (timeControl)
		{
		    RecalcScroll(timeControl);
		}
		ImInvalid((ObjPtr) dialogExists);
	    }
	}
    }
}

void DatasetChanged(dataset)
ObjPtr dataset;
/*Alerts the system that a dataset has changed*/
{
    globalDataset = dataset;
    ForAllVisWindows(ChangeWindowForDataset);
}

static ObjPtr ChangeDatasetInterpolate(object)
ObjPtr object;
/*Changes, according to the value of object, whether a dataset is interpolated*/
{
    ObjPtr repObj, value, dataObj;

    repObj = GetObjectVar("ChangeDatasetInterpolate", object, REPOBJ);
    if (!repObj)
    {
	return ObjFalse;
    }

    value = GetValue(object);

    SetVar(repObj, INTERPOLATEP, IsTrue(value) ? ObjTrue : ObjFalse);

    dataObj = GetVar(repObj, DATA);
    if (dataObj)
    {
	SetVar(dataObj, INTERPOLATEP, IsTrue(value) ? ObjTrue : ObjFalse);
    }

    DatasetChanged(repObj);
    return ObjTrue;
}

static ObjPtr HideDatasetWindow(window)
ObjPtr window;
/*Hide a dataset window, just return OK*/
{
    return ObjTrue;
}

static void NotPaletteAlert()
{
    WinInfoPtr errWindow;
    errWindow = AlertUser(UIERRORALERT, (WinInfoPtr) 0, "Only palettes can be dropped here.", 0, 0, "");
    SetVar((ObjPtr) errWindow, HELPSTRING,
	NewString("You have tried to drop an object which is not a palette into the palette \
box of a dataset.  Try again using a real palette."));
}

static ObjPtr DropInPaletteCorral(corral, object, x, y)
ObjPtr corral, object;
int x, y;
/*Drops an icon in a palette corral*/
{
    ObjPtr dataset;
    ObjPtr palette;
    ObjPtr icon;
    ObjPtr name;
    ObjPtr defaultIcon;
    ObjPtr contents;
    ObjPtr button;

    /*Find the visualization object*/
    dataset = GetObjectVar("DropInPaletteCorral", corral, REPOBJ);
    if (!dataset)
    {
	return ObjFalse;
    }

    /*Get the field object*/
    palette = GetVar(object, REPOBJ);
    if (!palette)
    {
	return ObjFalse;
    }
    if (!IsPalette(palette))
    {
	DoTask(NotPaletteAlert);
	return ObjFalse;
    }

    /*Create an icon for it*/
    name = GetStringVar("DropInPaletteCorral", palette, NAME);
    if (!name)
    {
	return ObjFalse;
    }

    defaultIcon = GetVar(palette, DEFAULTICON);
    if (defaultIcon)
    {
	ObjPtr locArray;
	real loc[2];
	icon = NewObject(defaultIcon, 0);
	SetVar(icon, NAME, name);
	loc[0] = x;
	loc[1] = y;
	locArray = NewRealArray(1, 2L);
	CArray2Array(locArray, loc);
	SetVar(icon, ICONLOC, locArray);
    }
    else
    {
	icon = NewIcon(x, y, ICONQUESTION, GetString(name));
    }
    
    /*Make the icon point to the palette*/
    SetVar(icon, REPOBJ, palette);

    /*Make it the only icon in the corral*/
    contents = NewList();
    PrefixList(contents, icon);
    SetVar(corral, CONTENTS, contents);
    SetVar(contents, PARENT, corral);
    SetVar(icon, CORRAL, corral);

    SetVar(icon, PARENT, corral);
    RecalcScroll(corral);

    /*Make this object the colored object*/
    SetVar(dataset, CPALETTE, palette);

    ImInvalid(corral);
    DatasetChanged(dataset);

    return ObjTrue;
}

static ObjPtr AddDatasetControls(object, contents)
ObjPtr object, contents;
/*Adds controls for a dataset*/
{
    ObjPtr var;
    ObjPtr corral;
    long info;
    ObjPtr textBox, checkBox, icon, palette, name, button;
    int left, right, bottom, top;
    char *s;

    /*Add the info description*/
    left = MINORBORDER;
    right = CWINWIDTH - 2 * CORRALBORDER - CWINCORRALWIDTH - MINORBORDER;
    top = CWINHEIGHT - MINORBORDER;
    bottom = MINORBORDER;

    info = GetDatasetInfo(object);

    s = tempStr;
    sprintf(s, "Data: ");
    while (*s) ++s;
    if (info & DS_TIMEDEPENDENT)
    {
	sprintf(s, "Time-dependent ");
    }
    else
    {
	sprintf(s, "Static ");
    }

    while (*s) ++s;
    if (info &  (DS_HASGEOMETRY | DS_HASNEWGEOMETRY))
    {
	sprintf(s, "geometric data\n");
    }
    else
    {
	SetCurField(FIELD1, object);
	if (info & DS_VECTOR)
	{
	    long nComponents;
	    nComponents = GetNComponents(FIELD1);
	    sprintf(s, "%d-vector %s\n", nComponents,
		info & DS_HASFORM ? "field" : "data");
	}
	else
	{
	    sprintf(s, "scalar %s\n", info & DS_HASFORM ? "field" : "data");
	}
    }
    while (*s) ++s;
    sprintf(s, "\nGrid: ");
    while (*s) ++s;
    if (info & DS_HASFORM)
    {
	char *t;
	long topDim;
	ObjPtr formObj;

	sprintf(s, "%d-dimensional %s grid\n", 
	    topDim = GetTopDim(object),
	    info & DS_UNSTRUCTURED ? "unstructured" : "structured");

	if (0 == (info & DS_UNSTRUCTURED))
	{
	    ObjPtr dimsArray;
	    int k;
	    real *meat;
	    t = s;
	    while (*t) ++t;

	    MakeVar(object, DATAFORM);
	    formObj = GetVar(object, DATAFORM);

	    dimsArray = GetFixedArrayVar("AddDatasetControls", formObj, DIMENSIONS, 1, topDim);
	    if (dimsArray)
	    {
		meat = ArrayMeat(dimsArray);
		for (k = 0; k < topDim; ++k)
		{
		    sprintf(t, "%g", meat[k]);
		    while (*t) ++t;
		    if (k < topDim - 1)
		    {
			sprintf(t, " by ");
			while (*t) ++t;
		    }
		}
		sprintf(t, "\n");
	    }
	}
    }
    else
    {
	sprintf(s, "none");
    }
    while (*s) ++s;
    bottom = top - DCINFOHEIGHT;
    textBox = NewTextBox(left, right, 
		     bottom, top,
		     0, "Info Text", tempStr);
    PrefixList(contents, textBox);
    SetVar(textBox, PARENT, contents);
    SetVar(textBox, TEXTFONT, NewString(DCFONT));
    SetVar(textBox, TEXTSIZE, NewInt(DCFONTSIZE));
    top = bottom - MAJORBORDER;

    /*Add the interpolate check box*/
    if (info & DS_TIMEDEPENDENT)
    {
	bottom = top - CHECKBOXHEIGHT;
	/*Create the Display Light check box*/
	checkBox = NewCheckBox(left, right, bottom, top,
			       "Interpolate in time",
			       GetPredicate(object, INTERPOLATEP));
	SetVar(checkBox, HELPSTRING,
	    NewString("If this box is checked, the field will automatically be \
interpolated in time.  If it is not checked, the time step nearest the current \
time will always be used.\n"));
	PrefixList(contents, checkBox);
	SetVar(checkBox, PARENT, contents);
	SetVar(checkBox, REPOBJ, object);
	SetMethod(checkBox, CHANGEDVALUE, ChangeDatasetInterpolate);
    }

    /*Add in the colors corral*/
    bottom = MAJORBORDER + BUTTONHEIGHT + MINORBORDER + TEXTBOXHEIGHT + TEXTBOXSEP;
    top = bottom + ONECORRALHEIGHT;
    right = left + ONECORRALWIDTH;
    corral = NewIconCorral(NULLOBJ,
		       left, right, bottom, top,
		       0);
    SetVar(corral, SINGLECORRAL, ObjTrue);
    SetVar(corral, TOPDOWN, ObjTrue);
    SetVar(corral, NAME, NewString("Color Palette"));
    PrefixList(contents, corral);
    SetVar(corral, HELPSTRING,
    NewString("This corral shows the color palette of the data in the dataset.  This \
may be another dataset.  Right now, you cannot drag any icons \
into the corral."));
    SetVar(corral, PARENT, contents);
    SetVar(corral, REPOBJ, object);
    SetMethod(corral, DROPINCONTENTS, DropInPaletteCorral);

    /*Create the source text box*/
    textBox = NewTextBox(left, right, 
			    bottom - TEXTBOXSEP - TEXTBOXHEIGHT, bottom - TEXTBOXSEP,
			    0, "Colors Text", "Color Palette");
    PrefixList(contents, textBox);
    SetVar(textBox, PARENT, contents);
    SetTextAlign(textBox, CENTERALIGN);

    /*Make a palette icon*/
    MakeVar(object, CPALETTE);
    palette = GetVar(object, CPALETTE);
    if (palette)
    {
	MakeVar(palette, NAME);
	name = GetVar(palette, NAME);
	icon = GetVar(palette, DEFAULTICON);
	if (icon)
	{
	    icon = NewObject(icon, 0);
	    SetVar(icon, NAME, name);
	}
	else
	{
	    icon = NewIcon(0, 0, ICONQUESTION, GetString(name));
	}
	SetVar(icon, ICONLOC, NULLOBJ);
	SetVar(icon, REPOBJ, palette);
	SetVar(icon, CORRAL, corral);
	DropIconInCorral(corral, icon);
    }

    /*Make a show info button*/
    button = NewFunctionButton((WinInfoPtr) GetVar(contents, OWNERWINDOW),
	    left, right, 
	    MINORBORDER, MINORBORDER + BUTTONHEIGHT, OF_SHOW_CONTROLS); 
    if (button)
    {
	SetVar(button, PARENT, contents);
	SetVar(button, STICKINESS, NewInt(STICKYBOTTOM + FLOATINGRIGHT + STICKYLEFT));
	PrefixList(contents, button);
    }

    return ObjTrue;
}

#ifdef PROTO
static int AddControlButton(ObjPtr object, WinInfoPtr controlWindow, ObjPtr controlField, ObjPtr panel, int top)
#else
static int AddControlButton(object, controlWindow, controlField, panel, top)
ObjPtr object;
WinInfoPtr controlWindow;
ObjPtr controlField;
ObjPtr panel;
int top;
#endif
/*Adds a number of control buttons from object into contents in window starting at
  top.  Returns a new top.*/
{
    ObjPtr curObj;
    ObjPtr icon;
    ObjPtr contents;
    Bool firstTime = true;

    contents = GetVar(controlField, CONTENTS);

    /*Fill the control field up with buttons*/
    MakeVar(object, CONTROLICON);

    curObj = object;
    icon = NULLOBJ;
    while (curObj)
    {
	if (!icon)
	{
	    icon = Get1Var(curObj, CONTROLICON);
	}
	if (icon)
	{
	    ObjPtr button;
	    ObjPtr panelContents;
	    FuncTyp method;
	    int whichIcon;
	    char *name;
	    ObjPtr var;

	    /*Give the button a chance to add controls*/
	    method = Get1Method(curObj, ADDCONTROLS);
	    if (method)
	    {
		var = GetIntVar("ShowVisControls", icon, WHICHICON);
		if (var)
		{
		    whichIcon = GetInt(var);
		}
		else
		{
		    whichIcon = ICONQUESTION;
		}

		var = GetStringVar("ShowVisControls", icon, NAME);
		if (var)
		{
		    name = GetString(var);
		}
		else
		{
		    name = "Unknown";
		}

		button = NewIconLabeledButton(0, CWINCCWIDTH, top - CWINICONBUTHEIGHT, top,
			whichIcon, UIYELLOW, name, BS_PITTED);
		SetMethod(button, ICONEXTRADRAW, GetMethod(icon, ICONEXTRADRAW));

		SetMethod(button, CHANGEDVALUE, ChangeControlPanelButton);

		/*Make a new panel contents just for this button*/
		panelContents = NewList();
		PrefixList(panelContents, controlField);
		SetVar(panelContents, PARENT, panel);
		SetVar(button, PANELCONTENTS, panelContents);
		SetVar(button, PANEL, panel);
		SetVar(panelContents, OWNERWINDOW, (ObjPtr) controlWindow);

		icon = NULLOBJ;
		(*method)(object, panelContents);
		PostfixList(contents, button);
		SetVar(button, PARENT, controlField);
		top -= CWINICONBUTHEIGHT + MINORBORDER;

		if (firstTime)
		{
		    ObjPtr mainDataset, notes;


		    /*Do the notes*/
		    notes = GetVar(object, NOTES);
		    if (notes)
		    {
			ObjPtr infoField;
			ThingListPtr runner;
		
			runner = LISTOF(notes);
			while (runner)
			{
			    int l, r, b, t;
			    char objName[200];
			    ObjPtr fieldContents;

			    name = GetString(runner -> thing);
			    runner = runner -> next;
		
			    if (runner)
			    {
				/*Make the info*/
				ObjPtr textBox;

				button = NewIconLabeledButton(0, CWINCCWIDTH, top - CWINICONBUTHEIGHT, top,
					    ICONNOTE, UIYELLOW, name, BS_PITTED);
		    
				SetMethod(button, CHANGEDVALUE, ChangeControlPanelButton);
		    
				/*Make a new panel contents just for this button*/
				panelContents = NewList();
				PrefixList(panelContents, controlField);
				SetVar(panelContents, PARENT, panel);
				SetVar(button, PANELCONTENTS, panelContents);
				SetVar(button, PARENT, panelContents);
				SetVar(button, PANEL, panel);
				SetVar(panelContents, OWNERWINDOW, (ObjPtr) controlWindow);
	
				strcpy(objName, name);
				strcat(objName, " info field");
				infoField = NewControlField(MINORBORDER, 
					    CWINWIDTH - CORRALBORDER - CWINCORRALWIDTH - MINORBORDER,
					    MINORBORDER, CWINHEIGHT - MINORBORDER,
					    objName, OBJECTSFROMTOP | BARRIGHT);	
				PrefixList(panelContents, infoField);
				SetVar(infoField, PARENT, panel);
				SetVar(infoField, BACKGROUND, NewInt(UIGRAY75));
				SetVar(infoField, BORDERTYPE, NewInt(1));
				fieldContents = GetVar(infoField, CONTENTS);
    
				PostfixList(contents, button);
				SetVar(button, PARENT, controlField);
				top -= CWINICONBUTHEIGHT + MINORBORDER;
    
				/*Now make the text box inside the window*/
				Get2DIntBounds(infoField, &l, &r, &b, &t);
				strcpy(objName, name);
				strcat(objName, " info text");
				textBox = NewTextBox(MINORBORDER, r - l - 2 * MINORBORDER,
					     -MINORBORDER - 100, /*Will change later*/
					     -MINORBORDER,
					     0, objName, GetString(runner -> thing));
				Set2DIntBounds(textBox,
					    MINORBORDER, r - l - 2 * MINORBORDER,
					    -MINORBORDER - TextHeight(textBox),
					    -MINORBORDER);
				PrefixList(fieldContents, textBox);
				SetVar(textBox, PARENT, infoField);
				RecalcScroll(infoField);

				runner = runner -> next;
			    }
			}
		    }

		    firstTime = false;
		    mainDataset = GetVar(object, MAINDATASET);
		    if (mainDataset)
		    {
			top = AddControlButton(mainDataset, controlWindow, controlField, panel, top);
		    }
		}
	    }
	}

	curObj = ClassOf(curObj);
    }

    return top;
}

static ObjPtr ShowDatasetControls(object, windowName)
ObjPtr object;
char *windowName;
/*Makes a new control window to control dataset object*/
{
    WinInfoPtr controlWindow;
    ObjPtr var;
    ObjPtr panel;
    ObjPtr controlField;
    ObjPtr contents;
    int left, right, bottom, top, width;
    WinInfoPtr dialogExists;
    ObjPtr icon;

    dialogExists = DialogExists((WinInfoPtr) object, NewString("Controls"));
    controlWindow = GetDialog((WinInfoPtr) object, NewString("Controls"), windowName, 
	CWINWIDTH, CWINHEIGHT, CWINWIDTH, CWINHEIGHT, WINFIXEDSIZE + WINDBUF);
    
    if (!dialogExists)
    {
	SetVar((ObjPtr) controlWindow, REPOBJ, object);

	/*Add in help string*/
	SetVar((ObjPtr) controlWindow, HELPSTRING,
	    NewString("This window shows controls for a dataset object.  \
At the right is an icon corral showing a series of icons.  Each icon represents a \
set of attributes of the visualization object.  On the left are the controls for \
the selected set of attributes.  \
Use Help In Context and click on the various controls to find out what they do.  \
Click on a different icon to choose a different set of attributes."));

	/*Add in a panel*/
	panel = NewPanel(greyPanelClass, 0, CWINWIDTH, 0, CWINHEIGHT);
	if (!panel)
	{
	    return ObjFalse;
	}
	contents = GetVar((ObjPtr) controlWindow, CONTENTS);
	PrefixList(contents, panel);
	SetVar(panel, PARENT, (ObjPtr) controlWindow);

	contents = GetVar(panel, CONTENTS);

	/*Add in a control field*/
	controlField = NewControlField(CWINWIDTH - CORRALBORDER - CWINCORRALWIDTH,
			       CWINWIDTH - CORRALBORDER,
			       CORRALBORDER,
			       CWINHEIGHT - CORRALBORDER,
				"Dataset Object Attributes", OBJECTSFROMTOP | BARRIGHT);
	SetVar(controlField, HELPSTRING,
		NewString("This icon button group shows sets of attributes of the dataset \
object.  The left side of the panel shows controls for the \
attribute given by the selected icon button.  To show another set of \
attributes, press another button."));
	SetVar(controlField, PARENT, panel);
	PrefixList(contents, controlField);

	contents = GetVar(controlField, CONTENTS);

	top = -MAJORBORDER;
	AddControlButton(object, controlWindow, controlField, panel, top);

	/*Adjust the scroll bars*/
	RecalcScroll(controlField);

	if (LISTOF(contents))
	{
	    SetValue(LISTOF(contents) -> thing, NewInt(1));
	    SetVar(controlField, BUTTON, LISTOF(contents) -> thing);
	}
    }

    return (ObjPtr) controlWindow;
}

ObjPtr FindDatasetByName(name)
ObjPtr name;
/*Finds a dataset by its name, returns it or nothing*/
{
    ThingListPtr runner;
    runner = LISTOF(allDatasets);
    while (runner)
    {
	if (Eql(GetVar(runner -> thing, NAME), name))
	{
	    return runner -> thing;
	}
	runner = runner -> next;
    }
    return NULLOBJ;
}

#ifdef PROTO
Bool InsertDatasetTimeStep(ObjPtr masterDataset, ObjPtr timeSlice, real time)
#else
Bool InsertDatasetTimeStep(masterDataset, timeSlice, time)
ObjPtr masterDataset, timeSlice;
real time;
#endif
/*Inserts a dataset time step into a master dataset*/
{
    /*Add on to an existing dataset*/
    InsertTimeSlice(GetVar(masterDataset, DATA), NewReal(time), GetVar(timeSlice, DATA));

    /*Register the dataset as changed*/
    DatasetChanged(masterDataset);
    return true;
}

#ifdef PROTO
Bool RegisterTimeDataset(ObjPtr dataset, real time)
#else
Bool RegisterTimeDataset(dataset, time)
ObjPtr dataset;
real time;
#endif
/*Registers a timeslice dataset, creating a new one or appending on to an old one
  based on time*/
{
    ObjPtr var;
    ObjPtr priorDataset;

    /*Find out if there already exists a time object like this*/
    var = GetVar(dataset, NAME);
    priorDataset = FindDatasetByName(var);
    if (priorDataset)
    {
	InsertDatasetTimeStep(priorDataset, dataset, time);
	/*Make the prior dataset the dataset in use for the icon test*/
	dataset = priorDataset;
    }
    else
    {
	ObjPtr timeSteps, timeData;
	/*Make a new dataset*/

	timeSteps = NewList();
	timeData = NewList();
	PrefixList(timeSteps, NewReal(time));
	PrefixList(timeData, GetVar(dataset, DATA));

	SetVar(dataset, DATA, NewTimedObject(timeSteps, timeData));
	RegisterNewDataset(dataset);
    }
    return true;
}

#ifdef PROTO
void SetDatasetTimeFormat(ObjPtr dataset, int format)
#else
void SetDatasetTimeFormat(dataset, format)
ObjPtr dataset;
int format;
#endif
/*Sets a dataset's time format*/
{
    SetVar(dataset, TIMEFORMAT, NewInt(format));
}

#ifdef PROTO
Bool RegisterNewDataset(ObjPtr dataset)
#else
Bool RegisterNewDataset(dataset)
ObjPtr dataset;
#endif
/*Registers a brand new dataset, doesn't look for name or anything*/
{
    ObjPtr corral, contents;
    ThingListPtr runner;
    WinInfoPtr datasetsWindow;
    int k;
    char *name;
    ObjPtr icon;
    ObjPtr defaultIcon;
    ObjPtr var;

    datasetsWindow = DatasetsWindow();
    PopWindow(datasetsWindow);

    IdleAllWindows();

    LongOperation();

    /*Tell dataset about its filename and directory*/
    if (curFileName)
    {
	SetVar(dataset, FILENAME, NewString(curFileName));
    }
    getcwd(tempStr, TEMPSTRSIZE);
    SetVar(dataset, DIRECTORY, NewString(tempStr));


    corral = GetObjectVar("RegisterNewDataset", (ObjPtr) datasetsWindow, CORRAL);
    var = GetStringVar("RegisterNewDataset", dataset, NAME);
    if (!var) return false;
    name = GetString(var);

    k = 1;
    while (FindDatasetByName(var))
    {
	++k;
	sprintf(tempStr, "%s (%d)", name, k);
	var = NewString(tempStr); 
    }
    if (k > 1)
    {
	name = GetString(var);
	SetVar(dataset, NAME, var);
    }

    /*Search for an icon*/
    contents = GetListVar("RegisterNewDataset", corral, CONTENTS);
    if (!contents) return false;
    runner = LISTOF(contents);
    while (runner)
    {
	ObjPtr repObj;
	repObj = GetVar(runner -> thing, REPOBJ);
	if (repObj == dataset)
	{
	    break;
	}
	runner = runner -> next;
    }

    if (!runner)
    {
	/*Must create a new icon*/
	defaultIcon = GetVar(dataset, DEFAULTICON);
	if (defaultIcon)
	{
	    icon = NewObject(defaultIcon, 0);
	    SetVar(icon, NAME, NewString(name));
	}
	else
	{
	    icon = NewIcon(0, 0, ICONQUESTION, name);
	    SetVar(icon, HELPSTRING,
		NewString("This icon represents a dataset.  You can select it and use the controls \
in this window to modify or visualize it, or you can drag it into the \
corral of an existing visualization window to visualize it there."));

	}
	SetVar(icon, REPOBJ, dataset);
	SetVar(icon, CORRAL, corral);
	SetVar(icon, ICONLOC, dsLocArray);
	dsLocArray = NULLOBJ;
	DropIconInCorral(corral, icon);
    }

    PrefixList(allDatasets, dataset);
    return true;
}

Bool RegisterDataset(dataset)
ObjPtr dataset;
/*Registers dataset as a dataset*/
{
    ObjPtr var;
    char *name;
    char *sPtr;

    var = GetStringVar("RegisteDataset", dataset, NAME);
    if (!var) return false;
    name = GetString(var);

    if (timedDatasets)
    {
	/*See if there is a time name*/
	sPtr = name;
	while (*sPtr)
	{
	    if (*sPtr == '@')
	    {
		/*Yes, it's a time slice!*/
		strncpy(tempStr, name, TEMPSTRSIZE);
		if (sPtr - name < TEMPSTRSIZE)
		{
		    tempStr[sPtr - name] = 0;
		}
		else
		{
		    tempStr[TEMPSTRSIZE] = 0;
		}
    
		/*Edit the dataset's name*/
		var = NewString(tempStr);
		SetVar(dataset, NAME, var);
    
		/*Copy the time*/
		strncpy(tempStr, sPtr + 1, TEMPSTRSIZE);
		tempStr[TEMPSTRSIZE] = 0;
    
		/*Get the edited name back into name*/
		name = GetString(var);
		break;
	    }
	    ++sPtr;
	}
    
	if (*sPtr == '@')
	{
	    /*It must be a time slice*/
	    real time;
	    int format;
    
	    /*Go through and parse the time*/
	    if (ParseTime(&time, &format, tempStr) <= 0)
	    {
		ReportError("RegisterDataset","Bad time format");
	    }
	    SetDatasetTimeFormat(dataset, format);
	    return RegisterTimeDataset(dataset, time);
	}
	else
	{
	    return RegisterNewDataset(dataset);
	}
    }
    else
    {
	return RegisterNewDataset(dataset);
    }
}

void ReportFileFormatError()
/*Task to report a bad read file error*/
{
    WinInfoPtr errWindow;
    errWindow = AlertUser(UIERRORALERT, DatasetsWindow(),
	"SciAn can only read data files with defined formats.  Set the file formats with Set Format and try again.",
	0, 0, "Yeah, right.");
    SetVar((ObjPtr) errWindow, HELPSTRING,
	NewString("You have tried to open a file cannot be opened because it \
is not a SciAn data file or its file format has not been set.  If you know what \
the proper format is, you can select the icon representing the file and press the \
Set Format button.  You will be presented with a series of choices.  Select the \
proper format, and try this operation again."));
}

void ReportDirError()
/*Task to report a directory tried to drag into datasets*/
{
    WinInfoPtr errWindow;
    errWindow = AlertUser(UIERRORALERT, DatasetsWindow(),
	"Directories cannot be opened as datasets.",
	0, 0, "Yeah, right.");
    SetVar((ObjPtr) errWindow, HELPSTRING,
	NewString("You have tried to open an entire directory as a dataset.  \
This does not work.  If you want to open all the files contained in the \
directory, open the directory, select all the files, and open them."));
}

static ObjPtr DropInDatasetsCorral(corral, icon, x, y)
ObjPtr corral, icon;
int x, y;
/*Drops an icon in a vis corral*/
{
    ObjPtr repObj;
    repObj = GetVar(icon, REPOBJ);
    if (repObj && IsFile(repObj))
    {
	ObjPtr fileName;		/*Name of the file*/
	ObjPtr fileType;		/*Type of the file*/
	ObjPtr fileReader;		/*Reader of the file*/
	FuncTyp method;
	real loc[2];
	loc[0] = x;
	loc[1] = y;

	fileName = GetVar(repObj, NAME);
	fileType = GetVar(repObj, FILETYPE);
	fileReader = GetVar(repObj, READER);

	if (GetInt(fileType) == DIRFILE)
	{
	    DoUniqueTask(ReportDirError);
	    return ObjFalse;
	}
	if (!fileName || !fileReader)
	{
	    DoUniqueTask(ReportFileFormatError);
	    return ObjFalse;
	}
	dsLocArray = NewRealArray(1, 2L);
	CArray2Array(dsLocArray, loc);

	method = GetMethod(repObj, OPEN);
	if (method)
	{
	    (*method)(repObj);
	}
	dsLocArray = NULLOBJ;
	return ObjTrue;
    }
    else
    {
	return ObjFalse;
    }
}

WinInfoPtr NewDatasetsWindow(void)
/*Create a new Datasets window*/
{
    WinInfoPtr objWin;
    ObjPtr panel, contents, corral, button;
    int l, r, b, t;
    int bw;

    /*Create the window*/
    objWin = NewObjWindow(NULLOBJ, "Datasets", WINDBUF, DSWINWIDTH, DSWINHEIGHT, SCRWIDTH, SCRHEIGHT);
    ContentsExpectWindowSize(objWin, DSWINWIDTH, DSWINHEIGHT);
    l = 0; r = DSWINWIDTH; b = 0; t = DSWINHEIGHT;

    /*Set a null but successful HIDE routine*/
    SetMethod((ObjPtr) objWin, HIDE, HideDatasetWindow);

    /*Add a help string*/
    SetVar((ObjPtr) objWin, HELPSTRING, 
	NewString("This window shows all the datasets that have been opened in \
Scian."));

    /*Put in a panel*/
    panel = NewPanel(greyPanelClass, 0, DSWINWIDTH, 0, DSWINHEIGHT);
    SetVar(panel, STICKINESS, NewInt(STICKYLEFT + STICKYRIGHT +
				     STICKYBOTTOM + STICKYTOP));

    contents = GetVar((ObjPtr) objWin, CONTENTS);
    PrefixList(contents, panel);
    SetVar(panel, PARENT, (ObjPtr) objWin);

    /*Put in buttons and an icon corral*/
    contents = GetListVar("NewDatasetsWindow", panel, CONTENTS);
    if (!contents)
    {
	return 0;
    }
    /*Make an icon corral*/
    corral = NewIconCorral(NULLOBJ, l + MINORBORDER, r - MINORBORDER, b + 3 * MINORBORDER + 2 * BUTTONHEIGHT, t - MINORBORDER, BARRIGHT + BARBOTTOM);
    SetVar(corral, STICKINESS, NewInt(STICKYLEFT + STICKYRIGHT +
				     STICKYBOTTOM + STICKYTOP));
    SetVar(corral, TOPDOWN, ObjTrue);
    SetVar((ObjPtr) objWin, CORRAL, corral);
    SetVar(corral, NAME, NewString("Datasets Corral"));
    SetVar(corral, HELPSTRING,
	NewString("This corral contains icons for all the datasets that have \
been read into SciAn.  You can visualize, show the controls of, or modify the datasets by selecting \
some of them and pressing the buttons at the bottom of the window.  You can delete \
datasets by choosing Delete from the Object menu.  You can also open new datasets by dragging icons from \
the file window into this corral."));
    SetMethod(corral, DROPINCONTENTS, DropInDatasetsCorral);
    PrefixList(contents, corral);
    SetVar(corral, PARENT, panel);

    l += MINORBORDER;
    r -= MINORBORDER;
    b += 2 * MINORBORDER + BUTTONHEIGHT;
    t = b + BUTTONHEIGHT;
    bw = (r - l - MINORBORDER) / 2;

    /*Make a visualize button*/
    button = NewFunctionButton(objWin,
		l, l + bw,
		b, b + BUTTONHEIGHT, OF_VISUALIZE); 
    if (button)
    {
	SetVar(button, PARENT, panel);
	SetVar(button, STICKINESS, NewInt(STICKYBOTTOM + STICKYLEFT + FLOATINGRIGHT));
	PrefixList(contents, button);
    }

    /*Make a visualize as... button*/
    button = NewFunctionButton(objWin,
		r - bw, r,
		b, b + BUTTONHEIGHT, OF_VISUALIZE_AS); 
    if (button)
    {
	SetVar(button, PARENT, panel);
	SetVar(button, STICKINESS, NewInt(STICKYBOTTOM + FLOATINGLEFT + STICKYRIGHT));
	PrefixList(contents, button);
    }

    t = b - MINORBORDER;
    b = t - BUTTONHEIGHT;
    /*Make a show info button*/
    button = NewFunctionButton(objWin,
		l, l + bw, 
		b, b + BUTTONHEIGHT, OF_SHOW_CONTROLS); 
    if (button)
    {
	SetVar(button, PARENT, panel);
	SetVar(button, STICKINESS, NewInt(STICKYBOTTOM + STICKYLEFT + FLOATINGRIGHT));
	PrefixList(contents, button);
    }

    /*Make a modify button*/
    button = NewFunctionButton(objWin, 
		r - bw, r,
		b, b + BUTTONHEIGHT, OF_MODIFY);
    if (button)
    {
	SetVar(button, PARENT, panel);
	SetVar(button, STICKINESS, NewInt(STICKYBOTTOM + FLOATINGLEFT + STICKYRIGHT));
	PrefixList(contents, button);
    }

    return objWin;
}

WinInfoPtr DatasetsWindow()
/*Returns or creates a datasets window*/
{
    WinInfoPtr retVal;

    retVal = GetWinFromTitle("Datasets");
    if (!retVal)
    {
	retVal = NewDatasetsWindow();
    }
    return retVal;
}

static ObjPtr GetDataMinMax(object)
ObjPtr object;
/*Gets and returns the min and max of the data in object for a dataset*/
{
    ObjPtr retVal;
    retVal = GetVar(object, MINMAX);
    if (retVal && IsRealArray(retVal) && RANK(retVal) == 1 && DIMS(retVal)[0] == 2)
    {
	return retVal;
    }
    else
    {
	if (GetDatasetInfo(object) &  (DS_HASGEOMETRY | DS_HASNEWGEOMETRY))
	{
	    real minmax[2];
	    minmax[0] = 0.0;
	    minmax[1] = 1.0;
	    retVal = NewRealArray(1, 2L);
	    CArray2Array(retVal, minmax);
	    SetVar(object, MINMAX, retVal);
	}
	else
	{
	    /*Must be a field*/
	    int nTraversalDims;		/*Number of dimensions to traverse*/
	    int whichDim;		/*Dimension counter*/
	    long *traversalDims = 0;	/*Dimensions of the dataset to traverse*/
	    long *index = 0;		/*Counting index*/
	    real sample;
	    real minmax[2];

	    SetCurField(FIELD5, object);

	    /*Get the information on traversing the dataset*/
	    nTraversalDims = CountTraversalDims(FIELD5);
	    if (nTraversalDims)
	    {
		traversalDims = (long *) malloc(sizeof(long) * nTraversalDims);
		index = (long *) malloc(sizeof(long) * nTraversalDims);
	    }
	    else
	    {
		index = (long *) malloc(sizeof(long));
	    }
	    GetTraversalDims(FIELD5, traversalDims);

	    /*Zero the index*/
	    for (whichDim = 0; whichDim < nTraversalDims; ++whichDim)
	    {
		index[whichDim] = 0;
	    }

	    sample = SelectFieldScalar(FIELD5, index);
	    minmax[0] = minmax[1] = sample;

	    /*Traverse all the points*/
	    do
	    {
		/*Sample the location at a point*/
		sample = SelectFieldScalar(FIELD5, index);
		minmax[0] = MIN(minmax[0], sample);
		minmax[1] = MAX(minmax[1], sample);

		/*Advance to next point*/
		for (whichDim = 0; whichDim < nTraversalDims; ++whichDim)
		{
		    if (traversalDims[whichDim] > 0)
		    {
			if ((++index[whichDim]) >= traversalDims[whichDim])
			{
			    index[whichDim] = 0;
			}
		        else
			{
			    break;
			}
		    }
		}
	    } while (whichDim < nTraversalDims); /*Break is based on advance*/

	    /*Free up temporary storage*/
	    SAFEFREE(traversalDims);
	    SAFEFREE(index);

	    retVal = NewRealArray(1, 2L);
	    CArray2Array(retVal, minmax);
	    SetVar(object, MINMAX, retVal);

	    return retVal;
	}
    }
}

static ObjPtr MakeDatasetPalette(object)
ObjPtr object;
/*Gets a palette for a dataset object*/
{
    ObjPtr palette;
    palette = GetVar(object, CPALETTE);
    if (!palette)
    {
	ObjPtr minMax, directory, name;
	real min, max, *elements;
	ObjPtr var;

	MakeVar(object, MINMAX);
	minMax = GetVar(object, MINMAX);
	if (minMax)
	{
	    elements = ELEMENTS(minMax);
	    min = elements[0];
	    max = elements[1];
	}
	else
	{
	    min = 0.0;
	    max = 1.0;
	}

	if (onePalette)
	{
	    if (commonPalette)
	    {
		min = MIN(min, ((PPtr) commonPalette) -> min);
		max = MAX(max, ((PPtr) commonPalette) -> max);
		palette = commonPalette;
	    }
	}
	if (!palette)
	{
	    var = GetVar(object, UNITSNAME);
	    if (var)
	    {
		palette = UnitsNameToPalette(GetString(var), min, max);
	    }
	    else
	    {
		palette = UnitsNameToPalette((char *) 0, min, max);
	    }
	    FieldPaletteName(palette, object);
	}

	/*If there is a saved palette, use it.*/
	directory = GetVar(object, DIRECTORY);
	if (directory)
	{
	    ReadObjectControls(palette, directory);
	}

	SetVar(object, CPALETTE, palette);
	if (onePalette)
	{
	    commonPalette = palette;
	}
    }
    return ObjTrue;
}

static ObjPtr GetPlainDatasetInfo(dataset)
ObjPtr dataset;
/*Returns a sum of flags defined in ScianDatasets.h giving info on the dataset*/
{
    long retVal;
    ObjPtr data, form, var;

    retVal = 0;
    if (form = GetVar(dataset, DATAFORM))
    {
	retVal |= DS_HASFORM;
	if (GetPredicate(form, UNSTRUCTURED))
	{
	    retVal |= DS_UNSTRUCTURED;
	}
    }
    if (data = GetVar(dataset, DATA))
    {
	if (IsTimedObj(data))
	{
	    ObjPtr timeData;

	    retVal |= DS_TIMEDEPENDENT;
	    timeData = GetVar(data, TIMEDATA);
	    if (timeData)
	    {
		long dimension = 0;
		data = GetObjectElement(timeData, &dimension);
		if (IsPicture(data))
		{
		    retVal |= DS_HASGEOMETRY;
		}
		else if (IsGeometry(data))
		{
		    retVal |= DS_HASNEWGEOMETRY;
		}
		else
		{
		    retVal |= DS_HASFIELD;
		}
	    }
	    else
	    {
		return 0;
	    }
	}
	else
	{
	    if (IsPicture(data))
	    {
		retVal |= DS_HASGEOMETRY;
	    }
	    else if (IsGeometry(data))
	    {
		retVal |= DS_HASNEWGEOMETRY;
	    }
	    else
	    {
		retVal |= DS_HASFIELD;
	    }
	}
    }
    if (var = GetVar(dataset, NCOMPONENTS))
    {
	retVal |= DS_VECTOR;
    }
    return NewInt(retVal);
}

int GetTopDim(dataset)
ObjPtr dataset;
/*Returns the topological or computational dimensionality of the dataset
  Returns 0 if it has no topological dimensions*/
{
    FuncTyp method;
    method = GetMethodSurely("GetTopDim", dataset, GETTOPDIM);
    if (method)
    {
	ObjPtr result;
	result = (*method)(dataset);
	if (result)
	{
	    return GetInt(result);
	}
	else
	{
	    ReportError("GetTopDim", "GETTOPDIM did not return dimension");
	    return 0;
	}
    }
    else
    {
	return 0;
    }
}

int GetSpatialDim(dataset)
ObjPtr dataset;
/*Returns the spatial dimensionality of the dataset
  Returns 0 if it has no topological dimensions*/
{
    FuncTyp method;
    method = GetMethodSurely("GetSpatialDim", dataset, GETSPATIALDIM);
    if (method)
    {
	ObjPtr result;
	result = (*method)(dataset);
	if (result)
	{
	    return GetInt(result);
	}
	else
	{
	    return 1;
	}
    }
    else
    {
	return 1;
    }
}

long GetDatasetInfo(dataset)
ObjPtr dataset;
/*Returns the spatial dimensionality of the dataset
  Returns 0 if it has no topological dimensions*/
{
    FuncTyp method;
    method = GetMethodSurely("GetDatasetInfo", dataset, GETDATASETINFO);
    if (method)
    {
	ObjPtr result;
	result = (*method)(dataset);
	if (result)
	{
	    return GetInt(result);
	}
	else
	{
	    ReportError("GetDatasetInfo", "No info returned from GETDATASETINFO");
	    return 0;
	}
    }
    else
    {
	return 0;
    }
}

ObjPtr GetDataFormTopDim(dataForm)
ObjPtr dataForm;
/*Returns the topological dimension of dataForm*/
{
    ObjPtr dims;

    dims = GetArrayVar("GetDataFormTopDim", dataForm, DIMENSIONS);
    if (dims)
    {
	if (GetPredicate(dataForm, UNSTRUCTURED))
	{
	    return NewInt(DIMS(dims)[0] - 1);
	}
	else
	{
	    return NewInt(DIMS(dims)[0]);
	}
    }
    else
    {
	return NewInt(0);
    }
}

static ObjPtr GetDatasetTopDim(dataset)
ObjPtr dataset;
/*Returns the topological or computational dimensionality of the dataset
  Returns 0 if it has no topological dimensions*/
{
    ObjPtr dataForm;
    FuncTyp method;

    dataForm = GetVar(dataset, DATAFORM);

    if (dataForm)
    {
	method = GetMethod(dataForm, GETTOPDIM);
	if (method)
	{
	    return (*method)(dataForm);
	}
    }
    else
    {
	ObjPtr data;
	/*No data form, must get the top dim of its DATA*/
	data = GetVar(dataset, DATA);
	if (data)
	{
	    method = GetMethod(data, GETTOPDIM);
	    if (method)
	    {
		return (*method)(data);
	    }
	}
    }
	
    return NewInt(0);
}

static ObjPtr GetDatasetSpatialDim(dataset)
ObjPtr dataset;
/*Returns the spatial dimensionality of the dataset*/
{
    ObjPtr dataForm;
    ObjPtr nComponents, bounds;

    dataForm = GetVar(dataset, DATAFORM);
    if (!dataForm)
    {
	return NewInt(0);
    }

    nComponents = GetVar(dataForm, NCOMPONENTS);
    if (nComponents)
    {
	return nComponents;
    }

    bounds = GetVar(dataForm, BOUNDS);
    if (bounds)
    {
	return NewInt(DIMS(bounds)[0] / 2);
    }
    return ObjFalse;
}

ObjPtr GetDatasetFormBounds(dataset)
ObjPtr dataset;
/*Returns the bounds of the dataset's data form*/
{
    ObjPtr dataForm;
    ObjPtr dims;
    
    dataForm = GetVar(dataset, DATAFORM);
    if (!dataForm)
    {
	return 0;
    }

    dims = GetArrayVar("GetDatasetFormBounds", dataForm, BOUNDS);
    return dims;
}

ObjPtr GetDatasetFormDims(dataset)
ObjPtr dataset;
/*Returns the dimensions of the dataset's data form*/
{
    FuncTyp method;
    method = GetMethodSurely("GetDatasetFormDims", dataset, GETFORMDIMS);

    if (method)
    {
	return (*method)(dataset);
    }
    else
    {
	return NULLOBJ;
    }
}

ObjPtr GetPlainDatasetFormDims(dataset)
ObjPtr dataset;
/*Returns the dimensions of the dataset's data form*/
{
    ObjPtr dataForm;
    ObjPtr dims;
    
    dataForm = GetVar(dataset, DATAFORM);
    if (!dataForm)
    {
	return 0;
    }

    dims = GetArrayVar("GetPlainDatasetFormDims", dataForm, DIMENSIONS);
    return dims;
}

ObjPtr GetDatasetKEdges(dataset, k)
ObjPtr dataset;
int k;
/*Returns a list of k-edges of the dataset.  Only works for unstructured 
  data forms.
*/
{
    ObjPtr dataForm;
    dataForm = GetVar(dataset, DATAFORM);
    if (!dataForm)
    {
	return NULLOBJ;
    }
    if (k == 1)
    {
	return GetVar(dataForm, EDGES);
    }
    else if (k == 2)
    {
	return GetVar(dataForm, CELLS);
    }
    else
    {
	ReportError("GetDatasetKEdges", "Edges beyond 2 are not defined");
	return NULLOBJ;
    }
}

static ObjPtr RegisterDatasetField(dataset, whichField)
ObjPtr dataset;
int whichField;
/*Registers datastet's data as whichField.*/
{
    ObjPtr data;
    FuncTyp method;

    MakeVar(dataset, CURDATA);
    data = GetVar(dataset, CURDATA);

    if (!data)
    {
	ReportError("RegisterDatasetField", "No CURDATA");
	return ObjFalse;
    }

    /*It's an array or vector.  Stuff it recursively*/
    method = GetMethodSurely("RegisterDatasetField", data, REGISTERFIELD);
    if (method)
    {
	return (*method)(data, whichField);
    }

    return ObjFalse;
}

int GetNComponents(whichField)
int whichField;
/*Returns the number of components of whichField*/
{
    return curFields[whichField] . nComponents;
}

int CountComponentDims(whichField, whichComponent)
int whichField, whichComponent;
/*Returns a count of whichField's component dimensions*/
{
    return curFields[whichField] . components[whichComponent] . nIndices;
}

long *GetComponentDims(whichField, whichComponent)
int whichField, whichComponent;
/*Returns a count of whichField's component dimensions*/
{
    return curFields[whichField] . components[whichComponent] . dimensions;
}

#ifdef PROTO
Bool CompareReals(real *r1, real *r2, long nReals)
#else
Bool CompareReals(r1, r2, lnReals)
real *r1;
real *r2;
long nReals;
#endif
/*Compares a bunch of reals, with any luck fast.  Uses indices, not pointers,
  to take advantage of MIPS architecture*/
{
    register long k, max;
    register real *rc1, *rc2;

    rc1 = r1;
    rc2 = r2;
    max = nReals;

    for (k = 0; k < max; ++k)
    {
	if (rc1[k] != rc2[k]) return false;
    }
    return true;
}

#ifdef PROTO
Bool IdenticalFields(int field1, int field2)
#else
Bool IdenticalFields(field1, field2)
int field1, field2;
#endif
/*Returns true if the two fields are identical.  Tries to be quick about it.
  May return false falsely.
*/
{
    register int k;

    /*Compare gross quantities*/
    if (curFields[field1] . fieldObj == curFields[field2] . fieldObj)
    {
	return true;
    }

    if (curFields[field1] . topDim !=
	curFields[field2] . topDim) return false;
    if (curFields[field1] . nComponents !=
	curFields[field2] . nComponents) return false;

    for (k = 0; k < curFields[field1] . nComponents; ++k)
    {
	register Component *comp1, *comp2;
	register int i;

	comp1 = &(curFields[field1] . components[k]);
	comp2 = &(curFields[field2] . components[k]);

	if (comp1 -> dataSize != comp2 -> dataSize) return false;
	if (!CompareReals(comp1 -> data, comp2 -> data, comp1 -> dataSize))
	{
	    return false;
	}
    }

    /*Have to assume they're identical*/
    return true;
}

#ifdef PROTO
int CountTraversalDims(int whichField)
#else
int CountTraversalDims(whichField)
int whichField;
#endif
/*Counts the traversal dimensions in whichField*/
{
    int k, curIndex;
    int maxDim = 0, temp;

    for (k = 0; k < curFields[whichField] . nComponents; ++k)
    {
	temp = MaxComponentIndices(whichField, k);
	if (temp > maxDim) maxDim = temp;
    }
    return maxDim;
}

#ifdef PROTO
void GetTraversalDims(int whichField, long dims[])
#else
void GetTraversalDims(whichField, dims)
int whichField;
long dims[];
#endif
/*Gets the traversal dims in whichField into dims, which had better be big enough*/
{
    int whichComponent, whichIndex, index, indIndex;
    int nIndices;

    nIndices = CountTraversalDims(whichField);

    for (indIndex = 0; indIndex < nIndices; ++indIndex)
    {
	dims[indIndex] = -1;
    }

    for (whichComponent = 0; whichComponent < curFields[whichField] . nComponents; ++whichComponent)
    {
	for (whichIndex = 0; whichIndex < curFields[whichField] . components[whichComponent] . nIndices; ++whichIndex)
	{
	    indIndex = curFields[whichField] . components[whichComponent] . indices[whichIndex];
	    if (indIndex >= 0)
	    {
		dims[indIndex] = curFields[whichField] . components[whichComponent] . dimensions[whichIndex];
	    }
	}
    }
}

static ObjPtr RegisterDataFormField(form, whichField)
ObjPtr form;
int whichField;
/*Registers data form form in field whichField*/
{
    ObjPtr dataset;			/*Dataset that may define this dataform*/

    MakeVar(form, DATA);
    dataset = GetVar(form, DATA);
    if (dataset)
    {
	Bool test;
	/*It's based on a dataset*/

	MakeVar(dataset, CURDATA);
	test = SetCurField(whichField, dataset);
	curFields[whichField] . topDim = GetTopDim(form);
	return test ? ObjTrue : ObjFalse;
    }

    if (GetPredicate(form, UNSTRUCTURED))
    {
	/*Unstructured grid.  This is an error, as all unstructured grids
	  must have a dataset in their form*/
	ReportError("RegisterDataFormField", "Unstructured grid without data");
    }
    else
    {
	ObjPtr dimsArray, boundsArray;
	real *dimsPtr, *boundsPtr;
	Component *components;
	int whichComponent;
	long topDim;
	long dataSize;
	Bool omErr = false;		/*True if om error occurred*/

	dimsArray = GetArrayVar("RegisterDataFormField", form, DIMENSIONS);
	boundsArray = GetArrayVar("RegisterDataFormField", form, BOUNDS);
	if (!dimsArray || !boundsArray)
	{
	    return ObjFalse;
	}

	topDim = DIMS(dimsArray)[0];

	/*Topological dimension is now topDim  Create that many components*/
	components = (Component *) malloc(sizeof(Component) * topDim);

	dimsPtr = ELEMENTS(dimsArray);
	boundsPtr = ELEMENTS(boundsArray);

	for (whichComponent = 0; whichComponent < topDim; ++whichComponent)
	{
	    int *indices = 0;
	    long *dimensions = 0;
	    long *steps = 0;
	    ObjPtr dataArray;

	    /*Set up one index, the current component*/
	    components[whichComponent] . nIndices = 1;
	    indices = malloc(sizeof(int));
	    if (indices)
	    {
		*indices = whichComponent;
	    }
	    else
	    {
		omErr = true;
		OMErr();
	    }

	    /*One set of dimensions*/
	    dimensions = malloc(sizeof(int));
	    if (dimensions)
	    {
		*dimensions = *dimsPtr;
	    }
	    else
	    {
		omErr = true;
		OMErr();
	    }

	    /*One step, unity*/
	    steps = malloc(sizeof(int));
	    if (steps)
	    {
		*steps = 1;
	    }
	    else
	    {
		omErr = true;
		OMErr();
	    }

	    /*Fill in a data array*/
	    dataArray = NewRealArray(1, (long) *dimsPtr);
	    if (dataArray)
	    {
		real min, max, *dp;
		long k;
		dp = ELEMENTS(dataArray);
		min = *boundsPtr;
		max = *(boundsPtr + 1);

		if (*dimsPtr < 1.5)
		{
		    *dp = (min + max) * 0.5;
		}
		else for (k = 0; k < (long) *dimsPtr; ++k)
		{
		    *dp = min + (max - min) * (k / (*dimsPtr - 1));
		    ++dp;
		}
	    }
	    else
	    {
		omErr = true;
		OMErr();
	    }
	    components[whichComponent] . flags = 0;
	    components[whichComponent] . data = ELEMENTS(dataArray);
	    components[whichComponent] . dataSize = (long) *dimsPtr;
	    components[whichComponent] . indices = indices;
	    components[whichComponent] . dimensions = dimensions;
	    components[whichComponent] . steps = steps;

	    ++dimsPtr;
	    boundsPtr += 2;
	}
	curFields[whichField] . components = components;
	curFields[whichField] . nComponents = topDim;
    }
    return ObjTrue;
}

static ObjPtr DeleteDatasetIcon(icon)
ObjPtr icon;
/*Delete repObj from dataset icon within a corral.  Always OK.*/
{
    DeleteFromList(allDatasets, GetVar(icon, REPOBJ));
    return ObjTrue;
}

static ObjPtr CleanupDataset(object)
ObjPtr object;
/*Cleans up a dataset*/
{
    return ObjTrue;
}

static ObjPtr MakeDatasetIconHelp(object, class)
ObjPtr object;
ObjPtr class;
/*Makes a dataset icon help string*/
{
    long info;
    long topDim;
    ObjPtr repObj;
    char *s;

    repObj = GetVar(object, REPOBJ);
    if (!repObj)
    {
	SetVar(class, HELPSTRING,
	   NewString("This icon is supposed to represent a dataset.  Unfortunately, \
it does not appear to be connected to anything right now.  See the section on Reporting Bugs \
to find out how to report this bug."));
	return ObjTrue;
    }
    info = GetDatasetInfo(repObj);
    if (info & DS_HASFORM)
    {
	topDim = GetTopDim(repObj);
    }
    else
    {
	topDim = 0;
    }
    s = tempStr;
    sprintf(s,
	"This icon represents a dataset, which contains");
    while (*s) ++s;
    if (info & (DS_HASGEOMETRY | DS_HASNEWGEOMETRY))
    {
	sprintf(s, "%s geometry.  ",
	    info & DS_TIMEDEPENDENT ? " time-dependent" : "");
    }
    else
    {
	if (topDim)
	{
	    sprintf(s, " a %d-D%s %s field.  ",
		topDim,
		info & DS_TIMEDEPENDENT ? " time-dependent" : "",
		info & DS_VECTOR ? "vector" : "scalar");
 	}
	else
	{
	    sprintf(s, " a%s %s field.  ",
		info & DS_TIMEDEPENDENT ? " time-dependent" : "",
		info & DS_VECTOR ? "vector" : "scalar");
	}
    }
    while (*s) ++s;
    strcpy(s, "You can select it and use the controls \
in this window to modify or visualize it, or you can drag it into the \
corral of an existing visualization window to visualize it there.");

    SetVar(class, HELPSTRING, NewString(tempStr));
    return ObjTrue;
}

ObjPtr DatasetPaletteSet(dataset)
ObjPtr dataset;
/*Sets PALETTESET on dataset*/
{
    SetVar(dataset, PALETTESET, ObjTrue);
    return ObjTrue;
}

static ObjPtr VisualizeDataset(object)
ObjPtr object;
/*Visualizes a single dataset in the current window*/
{
    if (object)
    {
	IdleAllWindows();
	AddObjToSpace(object, FindSpace(selWinInfo), GetVar((ObjPtr) selWinInfo, CORRAL), NULLOBJ, NULLOBJ);
	return ObjTrue;
    }
    return ObjFalse;
}

ObjPtr VisualizeDatasetAs(object)
ObjPtr object;
/*Visualizes an object as...*/
{
    AddToVisualizeAsList(object);
    return ObjTrue;
}

ObjPtr MakeDatasetControlIcon(dataset)
ObjPtr dataset;
/*Makes a dataset's control icon*/
{
    SetVar(dataset, CONTROLICON, GetVar(dataset, DEFAULTICON));
    return ObjTrue;
}

static ObjPtr MakeCurDataForm(dataset)
ObjPtr dataset;
/*Makes a dataset's or a data form's CURDATAFORM dummy variable*/
{
    SetVar(dataset, CURDATAFORM, ObjTrue);
    return ObjTrue;
}

void InitDatasets()
/*Initializes all the kinds of data sets*/
{
    int k;
    for (k = 0; k < MAXNCURFIELDS; ++k)
    {
	curFields[k] . nComponents = 0;
	curFields[k] . components = (Component *) 0;
    }

    iconDataset = NewIcon(0, 0, ICONQUESTION, "Dataset");
    AddToReferenceList(iconDataset);
    SetMethod(iconDataset, MAKE1HELPSTRING, MakeDatasetIconHelp);
    SetMethod(iconDataset, DELETEICON, DeleteDatasetIcon);

    icon4DScalar = NewObject(iconDataset, 0);
    SetVar(icon4DScalar, WHICHICON, NewInt(ICON4DSCALAR));
    AddToReferenceList(icon4DScalar);

    icon3DScalar = NewObject(iconDataset, 0);
    SetVar(icon3DScalar, WHICHICON, NewInt(ICON3DSCALAR));
    AddToReferenceList(icon3DScalar);

    icon2DScalar = NewObject(iconDataset, 0);
    SetVar(icon2DScalar, WHICHICON, NewInt(ICON2DSCALAR));
    AddToReferenceList(icon2DScalar);

    icon1DScalar = NewObject(iconDataset, 0);
    SetVar(icon1DScalar, WHICHICON, NewInt(ICON1DSCALAR));
    AddToReferenceList(icon1DScalar);


    icon4DVector = NewObject(iconDataset, 0);
    SetVar(icon4DVector, WHICHICON, NewInt(ICON4DVECTOR));
    AddToReferenceList(icon4DVector);

    icon3DVector = NewObject(iconDataset, 0);
    SetVar(icon3DVector, WHICHICON, NewInt(ICON3DVECTOR));
    AddToReferenceList(icon3DVector);

    icon2DVector = NewObject(iconDataset, 0);
    SetVar(icon2DVector, WHICHICON, NewInt(ICON2DVECTOR));
    AddToReferenceList(icon2DVector);

    icon1DVector = NewObject(iconDataset, 0);
    SetVar(icon1DVector, WHICHICON, NewInt(ICONVECTORS));
    AddToReferenceList(icon1DVector);

    datasetClass = NewObject(advertiseableClass, 0);
    AddToReferenceList(datasetClass);
    SetVar(datasetClass, CLASSID, NewInt(CLASS_DATASET));
    SetMethod(datasetClass, GETLONGNAME, GetPlainDatasetLongName);
    SetVar(datasetClass, DEFAULTICON, iconDataset);
    SetVar(datasetClass, NAME, NewString("Dataset"));
    SetVar(datasetClass, DOUBLECLICK, NewString(OF_SHOW_CONTROLS));
#ifdef OLDCODE
    SetMethod(datasetClass, LOCALCOPY, MakeLocalCopy);
#endif
    SetMethod(datasetClass, CLEANUP, CleanupDataset);
    SetMethod(datasetClass, GETDATASETINFO, GetPlainDatasetInfo);
    SetMethod(datasetClass, MINMAX, GetDataMinMax);
    SetMethod(datasetClass, NEWCTLWINDOW, ShowDatasetControls);
    SetMethod(datasetClass, SHOWCONTROLS, NewControlWindow);
    SetMethod(datasetClass, REGISTERFIELD, RegisterDatasetField);
    SetMethod(datasetClass, REGISTERFORM, RegisterDatasetForm);
    SetMethod(datasetClass, GETTOPDIM, GetDatasetTopDim);
    SetMethod(datasetClass, GETSPATIALDIM, GetDatasetSpatialDim);
    SetMethod(datasetClass, GETFORMDIMS, GetPlainDatasetFormDims);
    SetMethod(datasetClass, DELETE, DeleteObject);
    SetMethod(datasetClass, ADDCONTROLS, AddDatasetControls);
    SetMethod(datasetClass, CONTROLICON, MakeDatasetControlIcon);
    DeclareIndirectDependency(datasetClass, PALETTESET, CPALETTE, CHANGED);
    DeclareIndirectDependency(datasetClass, PALETTESET, CPALETTE, JUSTCOLORCHANGE);
    SetMethod(datasetClass, CPALETTE, MakeDatasetPalette);
    SetMethod(datasetClass, PALETTESET, DatasetPaletteSet);
    DeclareDependency(datasetClass, CURDATA, DATA);
    DeclareIndirectDependency(datasetClass, CURDATA, DATA, CURDATA);
    SetMethod(datasetClass, CURDATA, MakeDatasetCurData);
    DeclareDependency(datasetClass, TIMESTEPS, DATA);
    DeclareIndirectDependency(datasetClass, TIMESTEPS, DATA, TIMESTEPS);
    SetMethod(datasetClass, TIMESTEPS, MakeDatasetTimesteps);
    DeclareDependency(datasetClass, CURDATA, TIME);
    SetMethod(datasetClass, VISUALIZE, VisualizeDataset);
    SetMethod(datasetClass, VISUALIZEAS, VisualizeDatasetAs);
    DeclareIndirectDependency(datasetClass, CURDATAFORM, DATAFORM, CURDATAFORM);
    SetMethod(datasetClass, CURDATAFORM, MakeCurDataForm);

    dataFormClass = NewObject(NULLOBJ, 0);
    AddToReferenceList(dataFormClass);
    SetVar(dataFormClass, NCOMPONENTS, NewInt(3));
    SetMethod(dataFormClass, REGISTERFIELD, RegisterDataFormField);
    SetMethod(dataFormClass, GETTOPDIM, GetDataFormTopDim);
    DeclareIndirectDependency(dataFormClass, CURDATAFORM, DATA, CURDATA);
    SetMethod(dataFormClass, CURDATAFORM, MakeCurDataForm);
    /*SetMethod(dataFormClass, NORMALS, MakeDataFormNormals);*/

    data3DScalar = NewObject(datasetClass, 0);
    SetVar(data3DScalar, DEFAULTICON, icon3DScalar);
    AddToReferenceList(data3DScalar);

    data2DScalar = NewObject(datasetClass, 0);
    SetVar(data2DScalar, DEFAULTICON, icon2DScalar);
    AddToReferenceList(data2DScalar);

    data1DVector = NewObject(datasetClass, 0);
    SetVar(data1DVector, DEFAULTICON, icon1DVector);
    AddToReferenceList(data1DVector);

    data3DUnstructSurface = NewObject(datasetClass, 0);
    SetVar(data3DUnstructSurface, DEFAULTICON, icon2DScalar);
    AddToReferenceList(data3DUnstructSurface);

    /*Class for a geometry file*/
    iconGeometry = NewIcon(0, 0, ICONGEOMETRY, "Geometry");
    SetMethod(iconGeometry, DELETEICON, DeleteDatasetIcon);
    geometryClass = NewObject(datasetClass, 0);
    AddToReferenceList(geometryClass);
    SetVar(geometryClass, NAME, NewString("Geometry"));
    SetVar(geometryClass, DEFAULTICON, iconGeometry);
    SetMethod(geometryClass, DELETE, DeleteObject);

    allDatasets = NewList();
#if 0
    SetMethod(allDatasets, MARK, (FuncTyp) 0);
#endif
    AddToReferenceList(allDatasets);
}

void KillDatasets()
/*Kills the data sets*/
{
    int k;
    for (k = 0; k < MAXNCURFIELDS; ++k)
    {
	CleanCurField(k);
    }
    DeleteThing(allDatasets);
    DeleteThing(geometryClass);
    DeleteThing(data3DUnstructSurface);
    DeleteThing(data1DVector);
    DeleteThing(data3DScalar);
    DeleteThing(data2DScalar);
    DeleteThing(dataFormClass);
    DeleteThing(datasetClass);

    DeleteThing(icon1DVector);
    DeleteThing(icon2DVector);
    DeleteThing(icon3DVector);
    DeleteThing(icon4DVector);

    DeleteThing(icon1DScalar);
    DeleteThing(icon2DScalar);
    DeleteThing(icon3DScalar);
    DeleteThing(icon4DScalar);
    DeleteThing(iconDataset);
}