/*ScianGaryFiles.c
  Gary Fay
  2 September 1993
*/

#define HAVE_PROTOTYPES
#include "Scian.h"
#include "ScianTypes.h"
#include "ScianArrays.h"
#include "ScianIcons.h"
#include "ScianWindows.h"
#include "ScianObjWindows.h"
#include "ScianVisWindows.h"
#include "ScianVisObjects.h"
#include "ScianControls.h"
#include "ScianColors.h"
#include "ScianDialogs.h"
#include "ScianFiles.h"
#include "ScianFileSystem.h"
#include "ScianLists.h"
#include "ScianPictures.h"
#include "ScianErrors.h"
#include "ScianTimers.h"
#include "ScianDatasets.h"
#include "ScianFilters.h"
#include "ScianTextBoxes.h"
#include "ScianTitleBoxes.h"
#include "ScianButtons.h"
#include "ScianSliders.h"
#include "ScianScripts.h"
#include "ScianIDs.h"
#include "ScianVisContours.h"
#include "ScianStyle.h"
#include "ScianSpaces.h"
#include "ScianMethods.h"
#include "ScianObjFunctions.h"

#ifdef NETCDFDEF

#undef PROTO

#include 

#define OB_MISSING 0
#define OB_CLIP    1
#define OB_USE     2
#define CDFVECTORWIDTH 250
#define CDFRADIOWIDTH 250
#define GETCOMPONENTS USER1
#define USETIMEDIM USER2

#ifdef HAVE_PROTOTYPES
static ObjPtr AddnetCDFControls(ObjPtr reader, ObjPtr panelContents)
#else
static ObjPtr AddnetCDFControls(reader, panelContents)
ObjPtr reader, panelContents;
#endif

{  
    int left, right, bottom, top;
    ObjPtr titleBox, radio, checkBox, button;
    ObjPtr var;

    left = MAJORBORDER;
    bottom = MAJORBORDER;

    right = left + CDFVECTORWIDTH;
    top = bottom + CHECKBOXHEIGHT;

    checkBox = NewCheckBox(left, right, bottom, top, "Read unlimited time dimension",
        GetPredicate(reader, USETIMEDIM));
    PrefixList(panelContents, checkBox);
    SetVar(checkBox,PARENT,panelContents);
    SetVar(checkBox,HELPSTRING, NewString("This controls whether the NetCDF \
file reader will use the unlimited dimension as time.  NetCDF allows each dataset to have several limited dimensions and one limited \
dimension.  To use the unlimited dimension at time, \
define a variable with the same name as the unlimited dimension and store \
the time values.  Individual timesteps do not have to be in order or regularly spaced.  \
You can store all the time steps of one dataset into a single file or store individual \
time steps in separate files.\n\n\
If this box is not checked, then the unlimited dimension will just be interpreted as \
another spatial dimension."));
    AssocDirectControlWithVar(checkBox, reader, USETIMEDIM);

    bottom = top + MINORBORDER;
    top = bottom + CHECKBOXHEIGHT;

    checkBox = NewCheckBox(left, right, bottom, top, "Read 2- and 3- Vector Data",
        GetPredicate(reader, READVECTOR));
    PrefixList(panelContents, checkBox);
    SetVar(checkBox,PARENT,panelContents);
    SetVar(checkBox,HELPSTRING, NewString("This controls whether the netCDF file \
reader will read 2 and 3 dimensioned vector data.  If it is checked, whenever \
the first dimension or the last dimension of the dataset is 2 or 3, then the \
data will be interpreted as vector data.  The dataset is stored in terms of \
its X-component, Y-Component, and Z-Component.  If the box is not checked, \
data will be read as scalar data.\n\n\
For example, say you have a 3-D velocity field (40 by 50 by 60), where the \
velocity has three cartesian components.  Store it in netCDF as a variable \
with dimensions 60 by 50 by 40 by 3 in FORTRAN and 3 by 40 by 50 by 60 in C.  \
(An alternative approach would by to store the variable with dimensions 3 by \
60 by 50 by 40 in FORTRAN and 40 by 50 by 60 by 3 in C.)"));
    AssocDirectControlWithVar(checkBox, reader, READVECTOR);

    bottom = top + MINORBORDER;
    top = bottom + CHECKBOXHEIGHT;

    checkBox = NewCheckBox(left, right, bottom, top, "Include Vector Component Datasets",
        GetPredicate(reader, GETCOMPONENTS));
    PrefixList(panelContents, checkBox);
    SetVar(checkBox, PARENT, panelContents);
    SetVar(checkBox, HELPSTRING, NewString("This controls whether the netCDF \
file reader seperates the vector datasets into separate components.  If the \
read vector data box and this box are both checked, then each vector dataset \
will also have two or three other datasets for each of its components.  These \
component datasets can be used to visualize only one component.  If this control is not selected, then \
you may only view the vector field."));
    AssocDirectControlWithVar(checkBox, reader, GETCOMPONENTS);

    bottom = top + MINORBORDER;
    top = bottom + CHECKBOXHEIGHT;

    checkBox = NewCheckBox(left,right,bottom,top,"Wrap Polar Coordinates",
        GetPredicate(reader,WRAPPOLAR));
    if (!GetVar(reader,WRAPPOLAR)) SetVar(reader,WRAPPOLAR,ObjTrue);
    PrefixList(panelContents, checkBox);
    SetVar(checkBox,PARENT,panelContents);
    SetVar(checkBox,HELPSTRING,NewString("This controls whether the netCDF reader will automatically wrap scales in polar or spherical coordinates when reading the dataset.  When it is off, the scales along the axes will be used as is.  We recommend that you always leave thes check box on.  \n\n\
This option is needed because most SciAn visualization objects expect the data to be defined over a topologically continuous grid and convert spherical coordinates to Mercator projection by default.  If SciAn sees a sequence along a scale such as 178, 179, -180, -179, there is a discontinuity between 179 and -180.  Even though this sequence is continuous in spherical coordinates, it is discontinuous in the Mercator projection.  When the check is on, a sequence like this will automatically be converted to 178, 179, 180, 181.  When SciAn is enhanced to do spherical coordinates natively, both ways will work.\n\n\
In order for this check box to have any effect, you must create attributes for each polar dimension named \"Coordinates\" and \"polarUnit\". Coordinates must be set to \"polar\" or \"spherical\" and \"polarUnit\" must be set to \"degrees\", \"radians\", or \"gradians\" for the button to work.  The attribute names are case sensitive, but the settings aren't."));
    AssocDirectControlWithVar(checkBox,reader,WRAPPOLAR);

    bottom = top + MINORBORDER;
    top = bottom + CHECKBOXHEIGHT;

    checkBox = NewCheckBox(left, right, bottom, top, "Compress Data",
        GetPredicate(reader, COMPRESSDATA));
    if (!GetVar(reader, COMPRESSDATA)) SetVar(reader, COMPRESSDATA, ObjFalse);
    PrefixList(panelContents, checkBox);
    SetVar(checkBox, PARENT, panelContents);
    SetVar(checkBox, HELPSTRING, NewString("This controls whether the netCDF \
reader will store the data in compressed form.  In compressed form, only one \
byte is used to store each numerical value.  Compressed form takes only 1/4 \
as much space as normal, with an associated loss in precision.  If the \
minimum and maximum are set within the file, they are used to determine the \
compression mapping.  Otherwise, the minimum and maximum are determined from \
looking at the data."));
    AssocDirectControlWithVar(checkBox, reader, COMPRESSDATA);

    bottom = top + MINORBORDER;
    right = left + 2 * MINORBORDER + CDFRADIOWIDTH;
    top = bottom + 2 * MINORBORDER + 2 * CHECKBOXSPACING + 3 * CHECKBOXHEIGHT + TITLEBOXTOP;
    titleBox = NewTitleBox(left, right, bottom, top, "Handle Out-Of-Bounds Data");
    PrefixList(panelContents, titleBox);
    SetVar(titleBox, PARENT, panelContents);
    radio = NewRadioButtonGroup("Out of Bounds Radio");
    SetVar(radio, PARENT, panelContents);
    SetVar(radio, REPOBJ, reader);
    PrefixList(panelContents, radio);

    left += MINORBORDER;
    right -= MINORBORDER;
    top -= TITLEBOXTOP + MINORBORDER;
    bottom = top - CHECKBOXHEIGHT;
    button = NewRadioButton(left, right, bottom, top, "Treat as missing");
    AddRadioButton(radio, button);
    SetVar(button, HELPSTRING, NewString("When this button is down, values in the dataset which are out of bounds are treated as missing data. \
This provides another way of specifying missing data beyond the NetCDF \"missing_value\" variable."));

    top = bottom - CHECKBOXSPACING;
    bottom = top - CHECKBOXHEIGHT;
    button = NewRadioButton(left, right, bottom, top, "Clip to bounds");
    AddRadioButton(radio, button);
    SetVar(button, HELPSTRING, NewString("When this button is down, values in the dataset which are out of bounds are clipped to the bounds."));

    top = bottom - CHECKBOXSPACING;
    bottom = top - CHECKBOXHEIGHT;
    button = NewRadioButton(left, right, bottom, top, "Use as is");
    AddRadioButton(radio, button);
    SetVar(button, HELPSTRING, NewString("When this buttom is down, values in the dataset which are out of bounds are used as ordinary data values."));

    var = GetIntVar("AddnetCDFControls", reader, OUTOFBOUNDS);
    if (!var) SetVar(reader, OUTOFBOUNDS, NewInt(0));
    AssocDirectControlWithVar(radio, reader, OUTOFBOUNDS);
    SetVar(radio, HELPSTRING, NewString("This radio button group controls how \
values in the dataset which are out of bounds are treated.  A data value is \
considered out of bounds if it is outside the range given by the minimum and \
maximum, as specified in a call to NCATG or ncattget to the attribute \
valid_range.  If no minimum and maximum have been set, this radio button \
group has no effect."));

    return ObjTrue;
}

#ifdef HAVE_PROTOTYPES
ObjPtr ReadnetCDFFile(ObjPtr reader, char *fileName)
#else
ObjPtr ReadnetCDFFile(reader, fileName)
ObjPtr reader;
char *fileName;
#endif

{
    ObjPtr var;
    int ncid;                        /*netCDF file ID                          */
    int nDatasets;                   /*Number of possible datasets             */
    int recdim;                      /*ID of unlimited time dim if there is one*/
    int recvar;                      /*ID of unlimited time variable           */
    int r;                           /*Random counter                          */
    int whichSet;                    /*Index into current dataset              */
    long recSize;                    /*Size of time dimension                  */
    static long start[1] = {
	0    };       /*Beginning of scales                     */
    static long count[1];            /*Length of scales                        */
    real *recScale;                  /*Holds time scale                        */
    real missingValue = missingData; /*User set value used to represent missing
    				    data in netCDF file                     */

    char recName[MAX_NC_NAME];       /*Name of time dimension                  */
    char errMes[255];                /*Used to return error messages           */
    Bool missingSet = false;         /*True if a missing value has been set for
    				    the netCDF file                         */

    /*The following are read from file reader 
				    control panel:                          */
    Bool useTimeDim = false;         /*True iff the unlimited dimension was    
    				    used for time                           */
    Bool readVector = false;         /*True iff vector data is to be read      */
    Bool getComponents = true;       /*True iff vector data is also to be split
                                        into its corresponding components       */
    Bool wrapPolar = true;           /*True iff wrapping polar scales          */
    Bool compressData = true;        /*True iff data is to be compressed       */
    int outOfBoundsAction;           /*Action to take on out of bounds data    */
    int timeFormat = 0;		  /*Format for time			    */
    int stringLength;		  /*Length of string                        */
    nc_type testType;		  /*Type to test			    */
    char *timeString;		  /*Time string                             */
    char *s;			  /*Pointer into time string		    */

    ncopts = 0; /*Turn netCDF fatal error checking off*/
    /*Open netCDF file with read access only*/
    ncid = ncopen(fileName,NC_NOWRITE);
    if (ncid == -1)
    {
	sprintf(errMes,"Cannot read netCDF file %s.  The file is missing or is \
not a netCDF file.\n",fileName);
	FileFormatError("ReadnetCDFFile",errMes);
	return NULLOBJ;
    }
    /*Get number of possible datasets and dimensions*/
    if (-1 == ncinquire(ncid, (int *) 0, &nDatasets, (int *) 0, &recdim))
    {
	FileFormatError("ReadnetCDFFile", "Can't get netCDF file information");
	return NULLOBJ;
    }

    /*Get fileReader settings*/
    useTimeDim = GetPredicate(reader, USETIMEDIM);
    readVector = GetPredicate(reader, READVECTOR);
    getComponents = GetPredicate(reader,GETCOMPONENTS);
    wrapPolar = GetPredicate(reader,WRAPPOLAR);
    compressData = GetPredicate(reader, COMPRESSDATA);
    if ((var = GetIntVar("ReadnetCDFFile",reader,OUTOFBOUNDS)) == NULLOBJ)
	outOfBoundsAction = OB_MISSING;
    else
	outOfBoundsAction = GetInt(var);

    if (useTimeDim && (recdim != -1))
    {
	if (-1 == ncdiminq(ncid,recdim,recName,&recSize))
	{
	    FileFormatError("ReadnetCDFFile","Can't get information for the unlimited time dimension.");
	    return NULLOBJ;
	}
	recvar = ncvarid(ncid,recName);
	if (!(recScale = (real*) Alloc(sizeof(real) * recSize)))
	{
	    OMErr;
	    return NULLOBJ;
	}
	count[0] = recSize;
	if (recvar == -1 || (-1 == ncvarget(ncid,recvar,start,count,(void *) recScale)))
	{
	    sprintf(errMes,"Can't read in time scale.  Guessing time scale from 0 to %d",recSize);
	    FileFormatError("ReadnetCDFFile",errMes);
	    for (r = 0; r < recSize; ++r)
		recScale[r] = (real) r;
	}

	/*Try to get the format for time*/
	if (-1 != ncattinq(ncid, recvar, "time_format", &testType, &stringLength) &&
	    stringLength > 0)
	{
	    timeString = Alloc(stringLength + 1);
	    timeString[stringLength] = 0;
	    if (-1 == ncattget(ncid, recvar, "time_format", timeString))
	    {
		FileFormatError("ReadnetCDFFile", "Cannot get time format");
		timeFormat = TF_SECONDS + TF_SUBSECONDS;
	    }
	    else
	    {
		s = timeString;

		while (*s)
		{
		    switch(tolower(*s))
		    {
		    case 'h':		/*Hours*/
			timeFormat |= TF_HOURS;
			break;
		    case 'm':		/*Minutes*/
			timeFormat |= TF_MINUTES;
			break;
		    case 's':		/*Seconds*/
			timeFormat |= TF_SECONDS;
			break;
		    case '.':		/*Subseconds*/
			timeFormat |= TF_SUBSECONDS;
			break;
		    }
		    ++s;
		}

		if (!timeFormat) timeFormat = TF_SECONDS + TF_SUBSECONDS;
	    }
	    Free(timeString);
	}
	else
	{
	    timeFormat = TF_SECONDS + TF_SUBSECONDS;
	}
    }
    else
	recSize = 1;

    /*Register each set one at a time*/
    for (whichSet = 0; whichSet < nDatasets; ++whichSet)
    {
	int varid;                           /*netCDF ID for current dataset     */
	int initRank;                        /*Rank of the dataset               */
	int Rank;                            /*Adjusted for time dimension       */
	int adjustedRank;                    /*Adjusted for vectors and time     */
	int whichDim;                        /*Which dimension is being worked on*/
	int whichTime;                       /*Which timeset is being worked on  */
	int whichComponent;                  /*Which component is being worked on*/
	int vectorDim;                       /*Dimension of the vector components*/
	int i, k;                            /*Random Counters                   */
	int nComponents;                     /*# of components                   */
	int initDimPermutations[MAX_NC_DIMS];/*Permutations of dimensions, i.e., 
					   for each dimension, which spatial 
						   dimension it corresponds to       */
	int *DimPermutations;                /*Adjusted for time dimension       */
	int *adjustedDimPermutations;        /*Adjusted for vectors and time     */
	int *cdfToTopological = 0;           /*netCDF dimension corresponding to 
						   each topological dimension        */
	int DimIDs[MAX_NC_DIMS];             /*netCDF ID for each dimension      */
	long initDimSize[MAX_NC_DIMS];       /*Array holding the size of each dim*/
	long *DimSize;                       /*Adjusted for time dimension       */
	long *adjustedDimSize;               /*Adjusted for vectors and time     */
	long *permutedDimSize = 0;           /*Permuted to spatial dimensions    */
	long *dataStart;                     /*Start Corner to read netCDF data  */
	long *imap;                          /*How to map netCDF array to memory */
	long arraySize;                      /*Size of netCDF array              */
	long dataSize;                       /*Size of timeset in netCDF array   */
	long componentSize;                  /*Size of component in netCDF array */
	long dataIndex;                      /*Index to next netCDF value to 
						   write to scian next               */
	long *index;                         /*Indices of data value             */
	long *permutedIndex;                 /*Permuted to spatial dimensions    */
	float min = PLUSINF;                 /*Current minimum and maximium      */
	float max = MINUSINF;
	real VminMax[2];                     /*Array to hold vector min & max    */
	real minMax[2];                      /*Array to hold min and max         */
	real cTable[256];                    /*Compression Table                 */
	float *tempData;                     /*Holds data read from netCDF file  */
	float *data;                         /*Holds data for one timeset        */
	float *dataChunk;                    /*Holds data for one component      */
	real **initScales;                   /*Scale for axes                    */
	real **Scales;                       /*Adjusted for time dimension       */
	real **adjustedScales;               /*Adjusted for vectors and time     */
	real **permutedScales = 0;           /*Permuted to spatial dimensions    */
	unsigned char cd;                    /*Used for data compression         */
	char initDimName[MAX_NC_DIMS][MAX_NC_NAME];/*Array to dim names          */
	char dimsLabel[256];                 /*Label of a dimension, used to 
						   determine which axis              */
	char datasetName[256];               /*Name of current dataset           */
	char dataCoord[256];                 /*Set to type of coord system       */
	char dimsUnit[256];                  /*Units of polar dimension          */
	char *blanklessName;                 /*Dataset name w/o leading blanks   */
	char *tempName;                      /*temporary character pointer       */
	char *componentName;                 /*Name of component for scian       */
	Bool polarData = false;              /*True iff data is polar            */
	Bool minMaxSet = false;              /*True iff min and max were read in */
	Bool isLeftHanded = false;           /*True iff axes could be left-handed*/
	Bool dimUsed[MAX_NC_DIMS];           /*Flag to see if dim used           */
	ObjPtr xName = NULLOBJ;              /*Holds the names for each of the   */
	ObjPtr yName = NULLOBJ;              /*three spatial dimensions          */
	ObjPtr zName = NULLOBJ;
	ObjPtr curField;                     /*Field being assembled             */
	ObjPtr componentField = NULLOBJ;     /*Component field being assembled   */
	ObjPtr dataForm;                     /*Scian dataset dataform            */
	ObjPtr componentDataForm;            /*Component dataform                */
	ObjPtr minMaxArray;                  /*Holds min and max for dataset     */
	ObjPtr VminMaxArray;                 /*Holds min & max for vector dataset*/
	/*Get datset name, rank and the netCDF ID for each dimension*/
	if (-1 == ncvarinq(ncid, whichSet, datasetName, (nc_type *) 0, &initRank, DimIDs, (int *) 0))
	{
	    FileFormatError("ReadnetCDFFile", "Can't get the dataset information.");
	    return NULLOBJ;
	}

	if (initRank < 1)
	{
#if 0
	    /*It's a scalar value.  See if it's a reference to a vector dataset.*/
	    int len;
	    nc_type attType;
	    char *xName = 0, *yName = 0, *zName = 0;
	    char *finalName;
	    ObjPtr xDataset = NULLOBJ, yDataset = NULLOBJ, zDataset = NULLOBJ;
	    ObjPtr compDatasets, vecDataset;

	    /*Look for the x component variable*/
	    if (-1 != ncattinq(ncid, whichSet, "Xcomp", &attType, &len))
	    {
		if (attType != NC_CHAR)
		{
		    char errMes[255];
		    sprintf(errMes, "Xcomp attribute of variable %s must be characters",
			datasetName);
		    FileFormatError("ReadnetCDFFile",errMes);
		}
		else if (len <= 0)
		{
		    char errMes[255];
		    sprintf(errMes, "Xcomp attribute of variable %s must contain some characters",
			datasetName);
		    FileFormatError("ReadnetCDFFile",errMes);
		}
		else
		{
		    xName = Alloc(len + 1);
		    ncattget(ncid, whichSet, "Xcomp", xName);
		    xName[len] = 0;
		    
		    xDataset = FindDatasetByName(xName);
		    if (!xDataset)
		    {
			char errMes[255];
			sprintf(errMes, "Vector component dataset %s not found.",
			    xName);
			FileFormatError("ReadnetCDFFile",errMes);
		    }
		}
	    }

	    /*Look for the y component variable*/
	    if (-1 != ncattinq(ncid, whichSet, "Ycomp", &attType, &len))
	    {
		if (attType != NC_CHAR)
		{
		    char errMes[255];
		    sprintf(errMes, "Ycomp attribute of variable %s must be characters",
			datasetName);
		    FileFormatError("ReadnetCDFFile",errMes);
		}
		else if (len <= 0)
		{
		    char errMes[255];
		    sprintf(errMes, "Ycomp attribute of variable %s must contain some characters",
			datasetName);
		    FileFormatError("ReadnetCDFFile",errMes);
		}
		else
		{
		    yName = Alloc(len + 1);
		    ncattget(ncid, whichSet, "Ycomp", yName);
		    yName[len] = 0;

		    yDataset = FindDatasetByName(yName);
		    if (!yDataset)
		    {
			char errMes[255];
			sprintf(errMes, "Vector component dataset %s not found.",
			    yName);
			FileFormatError("ReadnetCDFFile",errMes);
		    }
		}
	    }

	    /*Look for the z component variable*/
	    if (-1 != ncattinq(ncid, whichSet, "Zcomp", &attType, &len))
	    {
		if (attType != NC_CHAR)
		{
		    char errMes[255];
		    sprintf(errMes, "Zcomp attribute of variable %s must be characters",
			datasetName);
		    FileFormatError("ReadnetCDFFile",errMes);
		}
		else if (len <= 0)
		{
		    char errMes[255];
		    sprintf(errMes, "Zcomp attribute of variable %s must contain some characters",
			datasetName);
		    FileFormatError("ReadnetCDFFile",errMes);
		}
		else
		{
		    zName = Alloc(len + 1);
		    ncattget(ncid, whichSet, "Zcomp", zName);
		    zName[len] = 0;

		    zDataset = FindDatasetByName(zName);
		    if (!zDataset)
		    {
			char errMes[255];
			sprintf(errMes, "Vector component dataset %s not found.",
			    zName);
			FileFormatError("ReadnetCDFFile",errMes);
		    }
		}
	    }
	    SAFEFREE(xName);
	    SAFEFREE(yName);
	    SAFEFREE(zName);

	    compDatasets = NewList();
	    if (xDataset) PostfixList(compDatasets, xDataset);
	    if (yDataset) PostfixList(compDatasets, yDataset);
	    if (zDataset) PostfixList(compDatasets, zDataset);

	    vecDataset = NewFilter(vectorJoinerClass, compDatasets, 
		-1, -1, -1, -1);

	    if (vecDataset)
	    {
		SetVar(vecDataset, NAME, NewString(datasetName));
		if (NULLOBJ == FindDatasetByName(datasetName))
		{
		    RegisterDataset(vecDataset);
		}
	    }
#endif
	}
	else
	{
	    /*Check if this is not a dimensional variable*/
	    if (-1 == ncdimid(ncid,datasetName))/*It is not a cooridinate variable*/
	    {
		curField = 0;/*Start with null field*/
		/*Assume it is x,y,z*/
		xName = NewString("X");
		yName = NewString("Y");
		zName = NewString("Z");
    
		/*Set up a place for the scales*/
		if (!(initScales = Alloc(sizeof(real *) * initRank)))
		{
		    OMErr();
		    return NULLOBJ;
		}
		for (whichDim = 0; whichDim < initRank; ++whichDim)
		{
		    initScales[whichDim] = 0;
		    /*Get dimension names and sizes*/
		    if (-1 == ncdiminq(ncid, DimIDs[whichDim], initDimName[whichDim], &initDimSize[whichDim]))
		    {
			sprintf(errMes,"Cannot get dimension %d information for dataset %s.",whichDim + 1,
			    datasetName);
			FileFormatError("ReadnetCDFFile", errMes);
			return NULLOBJ;
		    }
		}
    
		/*See if time dimension id used and adjust*/
		if (useTimeDim && (recdim != -1))
		{
		    Rank            = initRank - 1;
		    DimSize         = &(initDimSize[1]);
		    Scales          = &(initScales[1]);
		    DimPermutations = &(initDimPermutations[1]);
		}
		else
		{
		    Rank            = initRank;
		    DimSize         = initDimSize;
		    Scales          = initScales;
		    DimPermutations = initDimPermutations;
		}
    
		/*See if its a vector or scalar dataset*/
		if (readVector && (DimSize[0] == 2 || DimSize[0] == 3))
		{/*Its vector with the vector dimension at 0*/
		    vectorDim               = 0;
		    nComponents             = DimSize[0];
		    adjustedRank            = Rank - 1;
		    adjustedDimSize         = &(DimSize[1]);
		    adjustedScales          = &(Scales[1]);
		    adjustedDimPermutations = &(DimPermutations[1]);
		}
		else if (readVector && (initDimSize[initRank-1] == 2 || DimSize[Rank-1] == 3))
		{/*Its vector withthe vector dimension at initrank-1*/
		    vectorDim               = initRank - 1;
		    nComponents             = DimSize[Rank - 1];
		    adjustedRank            = Rank - 1;
		    adjustedDimSize         = DimSize;
		    adjustedScales          = Scales;
		    adjustedDimPermutations = DimPermutations;
		}
		else
		/*Its scalar*/
		{
		    vectorDim               = -1;
		    nComponents             = 1;
		    adjustedRank            = Rank;
		    adjustedDimSize         = DimSize;
		    adjustedScales          = Scales;
		    adjustedDimPermutations = DimPermutations;
		}
    
		for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
		{
		    real *tempScale;
    
		    count[0] = adjustedDimSize[whichDim];
    
		    if (vectorDim == initRank -1)
			varid = ncvarid(ncid, initDimName[whichDim+initRank-Rank]);
		    else
			varid = ncvarid(ncid, initDimName[whichDim+initRank-adjustedRank]);
    
		    if (wrapPolar && (varid != -1))
		    {
			if ((-1 == ncattget(ncid,varid,"Coordinates",dataCoord)) || ((0 != strcmp2(dataCoord,
			    "polar")) && (0 != strcmp2(dataCoord, "spherical"))))
			    /*It's not polar*/
			    polarData = false;
			else
			    /*It's polar*/
			    polarData = true;
		    }
    
		    /*Get scales from the file*/
		    adjustedScales[whichDim] = (real *) Alloc(adjustedDimSize[whichDim] * sizeof(real));
		    if (!adjustedScales[whichDim]) {
			OMErr(); 
			return NULLOBJ;
		    }
		    tempScale = (real *) Alloc(adjustedDimSize[whichDim] * sizeof(real));
		    if (!tempScale) {
			OMErr(); 
			return NULLOBJ;
		    }
    
		    if ((varid != -1) && (-1 != ncvarget(ncid, varid, start, count, (void *) tempScale)))
		    {
			for (i = 0; i < adjustedDimSize[whichDim]; ++i)
			    adjustedScales[whichDim][i] = tempScale[i];
			if (wrapPolar && polarData)
			{
			    double circ = 0.0;
    
			    if (-1 != ncattget(ncid,varid,"polarUnit",dimsUnit))
			    {
				if ((strcmp(dimsUnit, "degrees") == 0)) circ=360.0;
				else if ((strcmp(dimsUnit, "radians") == 0)) circ=2.0*M_PI;
				else if ((strcmp(dimsUnit, "gradians") == 0)) circ=400.0;
				if (circ > 0.0)
				{
				    if (adjustedScales[whichDim][1] < adjustedScales[whichDim][0])
				    {
					for (i = 2; i < adjustedDimSize[whichDim]; ++i)
					    while (adjustedScales[whichDim][i] > adjustedScales[whichDim][i-1])
						adjustedScales[whichDim][i] -= circ;
				    }
				    else
				    {
					for (i = 2; i < adjustedDimSize[whichDim]; ++i)
					    while (adjustedScales[whichDim][i] < adjustedScales[whichDim][i-1])
						adjustedScales[whichDim][i] += circ;
				    }
				}
			    }
			}
		    }
		    else
		    {/*Guess the scale*/
			sprintf(tempStr, "Guessing scale %d from %d to %d", whichDim, 0, initDimSize[whichDim]
			    - 1);
			FileFormatError("ReadnetCDFFile", tempStr);
			for (i = 0; i < adjustedDimSize[whichDim]; ++i)
			    adjustedScales[whichDim][i] = (real) i;
		    }
		    Free(tempScale);
    
		    /*Get the strings for this dimension*/
		    adjustedDimPermutations[whichDim] = -1;
    
		    /*Do stuff with the axis names*/
		    {
			if (vectorDim == initRank - 1)
			    strcpy(dimsLabel, initDimName[whichDim+initRank-Rank]);
			else
			    strcpy(dimsLabel, initDimName[whichDim+initRank-adjustedRank]);
			/*See if the dimension name corresponds to an axis*/
			adjustedDimPermutations[whichDim] = FindAxisDimension(dimsLabel);
			if (adjustedDimPermutations[whichDim] < 0)
			{
			    sprintf(tempStr,"Warning: Unrecognized axis '%s'", dimsLabel);
			    FileFormatError("ReadnetCDFFile", tempStr);
			}
			else
			{
			    switch(adjustedDimPermutations[whichDim])
			    {
			    case 0:
				xName = NewString(dimsLabel);
				break;
			    case 1:
				yName = NewString(dimsLabel);
				break;
			    case 2:
				zName = NewString(dimsLabel);
				break;
			    default:
				break;
			    }
			}
		    }
		}
    
		/*Fill in unpermuted axes*/
		for (k = 0; k < MAX_NC_DIMS; ++k)
		    dimUsed[k] = false;
    
		for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
		    if (adjustedDimPermutations[whichDim] >= 0)
			if (dimUsed[adjustedDimPermutations[whichDim]])
			{
			    sprintf(tempStr,"Dimension %d is used more than once", whichDim);
			    FileFormatError("ReadnetCDFFile", tempStr);
			    adjustedDimPermutations[whichDim] = -1;
			}
			else
			    dimUsed[adjustedDimPermutations[whichDim]] = true;
    
		k = 0;
		for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
		    if (adjustedDimPermutations[whichDim] < 0)
		    {
			while (dimUsed[k]) ++k;
			adjustedDimPermutations[whichDim] = k;
			dimUsed[k] = true;
			++k;
		    }
    
		/*Make netCDF to topological dimension mapping*/
		cdfToTopological = (int *) Alloc(sizeof(int) * adjustedRank);
		for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
		{
		    cdfToTopological[whichDim] = 0;
		    for (k = 0; k < adjustedDimPermutations[whichDim]; ++k)
			if (dimUsed[k])
			    ++cdfToTopological[whichDim];
		}
    
		/*Make permuted dim size*/
		permutedDimSize = (long *) Alloc(sizeof(long) * adjustedRank);
		for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
		    permutedDimSize[cdfToTopological[whichDim]] = adjustedDimSize[whichDim];
    
		/*Make permuted scales*/
		permutedScales = (real **) Alloc(sizeof(real *) * adjustedRank);
		for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
		    permutedScales[cdfToTopological[whichDim]] = adjustedScales[whichDim];
    
		/*See if it is left-handed or not*/
		isLeftHanded = false;
		for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
		    if (adjustedScales[whichDim][1] < adjustedScales[whichDim][0])
			isLeftHanded = isLeftHanded ? false : true;
    
		/*Calculate the size of the array and set up the starting indices*/
		if (!(dataStart = (long *) Alloc(sizeof(long) * initRank)))
		{
		    OMErr;
		    return NULLOBJ;
		}
		if (!(imap = (long *) Alloc(sizeof(long) * initRank)))
		{
		    OMErr;
		    return NULLOBJ;
		}
    
		arraySize = 1;
		for (i = 0; i < initRank; ++i)
		{
		    arraySize *= initDimSize[i];
		    dataStart[i] = 0;
		}
    
		if (vectorDim == initRank - 1)
		{
		    imap [initRank - 2] = sizeof(real);
		    i = initRank - 3;
		    while (i >= 0)
		    {
			if (i == 0 && (useTimeDim && recdim != -1))
			{
			    imap[vectorDim] = initDimSize[1] * imap[1];
			    imap[0]         = initDimSize[vectorDim] * imap[vectorDim];
			}
			else
			{
			    imap[i] = initDimSize[i+1] * imap[i+1];
			    if (i == 0)
				imap[vectorDim] = initDimSize[0] * imap[0];
			}
			--i;
		    }
		}
    
		/*Read in the data and get size of components if it is vector dataset*/
		tempData = (real *) Alloc(arraySize * sizeof(real));
    
		dataSize = arraySize / recSize;
    
		componentSize = dataSize / nComponents;
    
		if (vectorDim == initRank - 1)
		{
		    if (-1 == ncvargetg(ncid,whichSet,dataStart,initDimSize,0,imap,(void *) tempData))
		    {
			FileFormatError("ReadnetCDFFile", "Can't read field data.");
			return NULLOBJ;
		    }
		}
		else 
		{
		    if (-1 == ncvarget(ncid, whichSet, dataStart, initDimSize, (void *) tempData))
		    {
			FileFormatError("ReadnetCDFFile", "Can't read field data.");
			return NULLOBJ;
		    }
		}
    
		/*See if their is missing data*/
		if (-1 != ncattget(ncid, whichSet, "missing_value", &missingValue))
		    missingSet = true;
		else
		    missingSet = false;
    
		/*See if there is a min and max set, otherwise calculate it*/
		if (-1 != ncattinq(ncid, whichSet, "valid_range", (nc_type *) 0, (int *) 0))
		{
		    ncattget(ncid, whichSet, "valid_range", minMax);
		    min = minMax[0];
		    max = minMax[1];
		    if (min > max)
		    {
			minMax[1] = min;
			minMax[0] = max;
			min = minMax[0];
			max = minMax[1];
		    }
		    minMaxSet = true;
		}
		else
		{
		    min = PLUSINF;
		    minMax[0] = min;
		    max = MINUSINF;
		    minMax[1] = max;
		    minMaxSet = false;
		}
    
		/*Set missing Data*/
		for (dataIndex = 0; dataIndex < arraySize; ++dataIndex)
		    if (missingSet && (tempData[dataIndex] == missingValue))
		    {
			tempData[dataIndex] = missingData;
		    }
    
		/*Deal with out of bounds data*/
		if (minMaxSet && (outOfBoundsAction != OB_USE))
		{
		    for (dataIndex = 0; dataIndex < arraySize; ++ dataIndex)
			if (tempData[dataIndex] != missingData)
			{
			    if (tempData[dataIndex] < min)
				tempData[dataIndex] = (outOfBoundsAction == OB_CLIP) ? min : missingData;
			    else if (tempData[dataIndex] > max)
				tempData[dataIndex] = (outOfBoundsAction == OB_CLIP) ? max : missingData;
			}
		}
    
		if (!minMaxSet)
		{
		    for (dataIndex = 0; dataIndex < arraySize; ++ dataIndex)
			if (tempData[dataIndex] != missingData)
			    if (minMaxSet)
			    {
				min = MIN(min, tempData[dataIndex]);
				max = MAX(max, tempData[dataIndex]);
			    }
			    else
			    {
				min = max = tempData[dataIndex];
				minMaxSet = true;
			    }
		    minMax[0] = min;
		    minMax[1] = max;
		}
    
		/*Set min and max*/
		if (minMaxSet)
		    if (vectorDim < 0)
		    {
			minMaxArray = NewRealArray(1, 2L);
			CArray2Array(minMaxArray, minMax);
		    }
		    else
		    {
			VminMax[0] = 0.0;
			VminMax[1] = MAX(max, -min);
			VminMaxArray = NewRealArray(1, 2L);
			CArray2Array(VminMaxArray, VminMax);
			if (getComponents && vectorDim >= 0)
			{
			    minMaxArray = NewRealArray(1, 2L);
			    CArray2Array(minMaxArray, minMax);
			}
		    }
    
		/*Make compression table*/
		if (compressData)
		{
		    cTable[0] = missingData;
		    for (k = 1; k < 256; ++k)
			cTable[k] = min + (max - min) * (k - 1) / 245.0;
		}
    
		/*Loop through and register each time dataset and component*/
		for (whichTime = 0; whichTime < recSize; ++whichTime)
		{
		    data = tempData + whichTime * dataSize;
    
		    /*Set up the new dataform, dataset, and current field for scian*/
		    blanklessName = datasetName;
		    while (*blanklessName && *blanklessName == ' ') ++blanklessName;
		    dataForm = NewSeparableDataForm(blanklessName, adjustedRank, permutedDimSize, permutedScales);
		    if (compressData)
			curField = NewCompressedDataset(blanklessName, adjustedRank, permutedDimSize, vectorDim
			    >= 0 ? nComponents : 0, cTable);
		    else
			curField = NewDataset(blanklessName, adjustedRank, permutedDimSize, vectorDim >=
			    0 ? nComponents : 0);
		    if (isLeftHanded)
			SetVar(curField, ISLEFTHANDED, ObjTrue);
		    if (xName)
			SetVar(curField, XNAME, xName);
		    if (yName)
			SetVar(curField, YNAME, yName);
		    if (zName)
			SetVar(curField, ZNAME, zName);
		    if (vectorDim < 0)
			SetVar(curField,MINMAX,minMaxArray);
		    else
			SetVar(curField,MINMAX,VminMaxArray);
		    SetDatasetForm(curField, dataForm);
		    SetCurField(FIELD1, curField);
    
		    for (whichComponent = 0; whichComponent < nComponents; ++whichComponent)
		    {
			dataChunk = data + whichComponent * componentSize;
    
			/*If it is vector datset, set up dataform, dataset, and field for each component as well*/
			if (getComponents && vectorDim >= 0)
			{
			    tempName = blanklessName;
			    i = 0;
			    while (*tempName != '@' && *tempName != '\0')
			    {
				tempStr[i] = *tempName;
				++i;
				++tempName;
			    }
			    tempStr[i] = '.';
			    ++i;
			    if (whichComponent == 0)
				tempStr[i] = 'X';
			    else if (whichComponent == 1)
				tempStr[i] = 'Y';
			    else
				tempStr[i] = 'Z';
			    ++i;
			    while (*tempName != '\0')
			    {
				tempStr[i] = *tempName;
				++i;
				++tempName;
			    }
			    tempStr[i] = '\0';
			    componentName = tempStr;
			    componentDataForm = NewSeparableDataForm(componentName, adjustedRank, permutedDimSize,
				permutedScales);
			    if (compressData)
				componentField = NewCompressedDataset(componentName, adjustedRank,
				    permutedDimSize, 0, cTable);
			    else
				componentField = NewDataset(componentName, adjustedRank, permutedDimSize,
				    0);
			    if (isLeftHanded)
				SetVar(componentField, ISLEFTHANDED, ObjTrue);
			    if (xName)
				SetVar(componentField, XNAME, xName);
			    if (yName)
				SetVar(componentField, YNAME, yName);
			    if (zName)
				SetVar(componentField, ZNAME, zName);
			    SetVar(componentField,MINMAX,minMaxArray);
			    SetDatasetForm(componentField, componentDataForm);
			    SetCurField(FIELD2, componentField);
			}
    
			/*Set up scian indices*/
			index = (long *) Alloc(sizeof(long) * adjustedRank);
			permutedIndex = (long *) Alloc(sizeof(long) * adjustedRank);
    
			/*Copy tempData to the Scian dataset*/
			for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
			    index[whichDim] = 0;
			for (dataIndex = 0; dataIndex < componentSize; ++ dataIndex)
			{
			    for (whichDim = 0; whichDim < adjustedRank; ++whichDim)
				permutedIndex[cdfToTopological[whichDim]] = index[whichDim];
			    if (compressData)
			    {
				if (dataChunk[dataIndex] == missingData)
				    cd = 0;
				else
				{
				    cd = 1.5 + (dataChunk[dataIndex]-min) / (max-min)*254.0;
				    if (cd < 1) cd = 1;
				    if (cd > 255) cd = 255;
				}
				PutCompressedFieldComponent(FIELD1,whichComponent,permutedIndex,cd);
				if (getComponents &&(vectorDim >= 0))
				    PutCompressedFieldComponent(FIELD2,0,permutedIndex,cd);
			    }
			    else
			    {
				PutFieldComponent(FIELD1,whichComponent,permutedIndex,dataChunk[dataIndex]);
				if (getComponents && (vectorDim >= 0))
				    PutFieldComponent(FIELD2,0,permutedIndex,dataChunk[dataIndex]);
			    }
			    for (whichDim = adjustedRank - 1; whichDim >= 0; --whichDim)
			    {
				if ((++index[whichDim]) >= adjustedDimSize[whichDim])
				    index[whichDim] = 0;
				else
				    break;
			    }
			}
			SAFEFREE(index);
			SAFEFREE(permutedIndex);
    
			/*register component field*/
    
			SetDatasetTimeFormat(componentField, timeFormat);
    
			if (getComponents && vectorDim >= 0)
			    if (useTimeDim && recdim != -1)
				RegisterTimeDataset(componentField,recScale[whichTime]);
			    else
				RegisterDataset(componentField);
		    }
		    /*register dataset field*/
    
		    SetDatasetTimeFormat(curField, timeFormat);
    
		    if (useTimeDim && recdim != -1)
			RegisterTimeDataset(curField,recScale[whichTime]);
		    else
			RegisterDataset(curField);
		}
    
		/*Free temp data*/
		SAFEFREE(dataStart);
		SAFEFREE(imap);
		SAFEFREE(tempData);
		SAFEFREE(permutedDimSize);
		SAFEFREE(cdfToTopological);
		for (i = 0; i < initRank; ++i)
		    if (initScales[i])  Free(initScales[i]);
		SAFEFREE(initScales);
		SAFEFREE(permutedScales);
	    }
	}
    }
    return NULLOBJ;
}

#ifdef HAVE_PROTOTYPES
void  InitGaryFiles(void)
#else
void  InitGaryFiles()
#endif
{
    ObjPtr  fileReader;

    fileReader = NewFileReader("NetCDF");
    SetVar(fileReader, EXTENSION, NewString("nc"));
    SetVar(fileReader, OUTOFBOUNDS, NewInt(1));
    AddSnapVar(fileReader, OUTOFBOUNDS);
    SetVar(fileReader, HELPSTRING, NewString("This file reader reads \
scientific data using the netCDF data format produced by Unidata.  Libraries \
for writing and reading netCDF files and documentation are available via \
anonymous ftp from unidata.ucar.edu from pub/netcdf/cetcdf.tar.Z."));
    SetVar(fileReader, USETIMEDIM, ObjTrue);
    AddSnapVar(fileReader, USETIMEDIM);
    SetVar(fileReader, READVECTOR, ObjTrue);
    AddSnapVar(fileReader, READVECTOR);
    SetVar(fileReader, GETCOMPONENTS, ObjFalse);
    AddSnapVar(fileReader, GETCOMPONENTS);
    SetVar(fileReader,WRAPPOLAR, ObjTrue);
    AddSnapVar(fileReader, WRAPPOLAR);
    SetVar(fileReader, COMPRESSDATA, ObjFalse);
    AddSnapVar(fileReader, COMPRESSDATA);
    SetMethod(fileReader,READALL,ReadnetCDFFile);
    SetMethod(fileReader,ADDCONTROLS,AddnetCDFControls);
    ApplySavedSettings(fileReader);
}

#ifdef HAVE_PROTOTYPES
void   KillGaryFiles(void)
#else
void   KillGaryFiles()
#endif
{
}
#endif