/*******************************************************************************/
/********************************STF READER ************************************/
/*******************************************************************************/
/*
  10/26/1992
  ScianHwuFiles.c
  ReadSTFFile routine

  Tzong-Yow Hwu

  The file format is:
  
      RANK <rank>
      DIMENSIONS <xdim> <ydim> ....
      SCALAR|VECTOR <nComponents>
      [ BOUNDS     <xmin> <xmax> <ymin> <ymax> .... ] {required when no grid}
      [ INTERLACED|NONINTERLACED ]  { default is INTERLACED }  
      [ ORDER COLUMN|ROW ]          { defalut is COLUMN major }
      [ NAME <datasetName> ]
      [ TIME <time> ]
      DATA
      <values> ......
      .
      .
      .
      VECTOR <nComponents>
      [INTERLACED|NONINTERLACED] { takes previous spec as default }
      [ ORDER COLUMN|ROW ]       { takes previous spec as default }
      GRID
      <values> ......
      .
      .
      .
      END    { the default is reset to INTERLACED and COLUMN major when 
	       anther datasets begin }
      .
      .
      .
      .

  Grid is optional for any dataset, if a grid is not present then regular
  dataform is assumed.  The order of grid and data for a dataset is 
  arbitrary in the data file.  However, rank and dimensions are required 
  before either value is processed and the rank must precede dimensions. 
  If no grid is specified, then bounds is required.  
  If no order is specified for a set of values, then the default of 
  column major is assumed.
  If no dataset name is specified for the dataset, then it inherits from the
  previous dataset in the file.  If this happens to the first dataset, then
  it takes the name of the file.
*/
#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"
#include "ScianHwuFiles.h"

#define STF_MAXHEADER 9
#define STF_MAXNAMELEN 256

/* field type */
#define STF_SCALAR 0
#define STF_VECTOR 1

/* vector field data permutation type */
#define STF_NONINTERLACED 0
#define STF_INTERLACED 1

#define STF_COLUMNMAJOR 0
#define STF_ROWMAJOR 1

/* for datasetType */
#define STF_ERROR -1
#define STF_END  0
#define STF_DATA 1
#define STF_GRID 2

/* for specFlag */
#define STF_ISTYPE 1
#define STF_ISDIMS 2
#define STF_ISBOUNDS 4
#define STF_ISNAME 8
#define STF_ISORDER 16
#define STF_ISRANK 32 
#define STF_ISPERMUTATION 64 
#define STF_ISTIME 128

typedef struct  {
     int fieldType;   /* SCALAR or VECTOR */
     long *dimensions; /* xdim, ydim, zdim .... */
     real *bounds;  /* bounds array */
     int nBounds; /* number of bounds provided by the user if any */
     char datasetName[STF_MAXNAMELEN];
     int order; /* COLUMNMAJOR(default) or ROWMAJOR */
     int rank;
     int permutation; /* INTERLACED or NONINTERLACED */
     int nComponents;  /* number of components for vector data set */
     real time;
     } STF_spec;

#ifdef HAVE_PROTOTYPES
static int parseHeader(FILE *inFile, int *specFlag, STF_spec *datasetSpec)
#else
static int parseHeader(inFile, specFlag, datasetSpec)
FILE *inFile;
int *specFlag;
STF_spec *datasetSpec;
#endif
{
  int datasetType = STF_END; /* DATA, GRID, ERROR, or END */
  int headerLine = 0;
  char dirStr[360]; 
  int counter;
  int c;

  for (SkipBlanks(inFile); (c = getc(inFile)) >= 0 && headerLine < STF_MAXHEADER;
       SkipBlanks(inFile))
  {
    ungetc(c, inFile);

    if (1 == fscanf(inFile, "%s", dirStr))
      {
        if (dirStr[0] == '#')
          {
            /* a comment */
            ReadLn(inFile);
            continue;
          }

        if (0 == strcmp2(dirStr, "RANK"))
          {
            if ( *specFlag & STF_ISRANK )
              {
                FileFormatError("ReadSTFFile", "Error reading header: double RANK specifications.");
                datasetType = STF_ERROR;
                break;
              }
            SkipBlanks(inFile);
            if ( 1 != fscanf( inFile, "%d", &(datasetSpec->rank) ) )
              {
                FileFormatError("ReadSTFFile", "Error reading header: RANK value format error.");
                datasetType = STF_ERROR;
                break;
              }
            if ( datasetSpec->rank < 1  )
              /* rank must be greater than or equal to 1 */
              {
                FileFormatError("ReadSTFFile", "Error reading header: invalid RANK value.");
                datasetType = STF_ERROR;
                break;
              }
             ReadLn(inFile);
             headerLine++;
             *specFlag |= STF_ISRANK;
             continue; 
           }

        if (0 == strcmp2(dirStr, "TIME"))
          {
            if ( *specFlag & STF_ISTIME )
              {
                FileFormatError("ReadSTFFile", "Error reading header: double TIME specifications.");
                datasetType = STF_ERROR;
                break;
              }
            SkipBlanks(inFile);
            if ( 1 != fscanf( inFile, "%g", (float *) &(datasetSpec->time) ) )
              {
                FileFormatError("ReadSTFFile", "Error reading header: invalid TIME value.");
                datasetType = STF_ERROR;
                break;
              }

             ReadLn(inFile);
             headerLine++;
             *specFlag |= STF_ISTIME;
             continue; 
           }

         if (0 == strcmp2(dirStr, "SCALAR"))
           {
             if ( *specFlag & STF_ISTYPE )
               {
                 FileFormatError("ReadSTFFile", "Error reading header: double specifications of SCALAR or VECTOR header line.");
                 datasetType = STF_ERROR;
                 break;
               }
             headerLine++;
             ReadLn(inFile);
             datasetSpec->fieldType = STF_SCALAR;
             datasetSpec->nComponents = 0;
             *specFlag |= STF_ISTYPE;     
             continue;
           }

         if (0 == strcmp2(dirStr, "VECTOR"))
           {
             if ( *specFlag & STF_ISTYPE )
               {
                 FileFormatError("ReadSTFFile", "Error reading header: double specifications of SCALAR or VECTOR header line.");
                 datasetType = STF_ERROR;
                 break;
               }

             SkipBlanks(inFile);
             if ( 1 != fscanf(inFile, "%d", &(datasetSpec->nComponents)) )
               {
                 FileFormatError("ReadSTFFile", "Error reading header: VECTOR component number format error.");
                 datasetType = STF_ERROR;
                 break;
               }
             if ( datasetSpec->nComponents < 1 )
               {
                 FileFormatError("ReadSTFFile", "Error reading header: invalid VECTOR component value.");
                 datasetType = STF_ERROR;
                 break;
               }

             headerLine++;

             ReadLn(inFile);
             datasetSpec->fieldType = STF_VECTOR;
             *specFlag |= STF_ISTYPE;     
             continue;
           }

         if (0 == strcmp2(dirStr, "INTERLACED"))
           {
             if ( *specFlag & STF_ISPERMUTATION )
               {
                 FileFormatError("ReadSTFFile", "Error reading header: double INTERLACED/NONINTERLACED specifications.");
                 datasetType = STF_ERROR;
                 break;
               }

             headerLine++;
             ReadLn(inFile);
             datasetSpec->permutation = STF_INTERLACED;
             *specFlag |= STF_ISPERMUTATION;
             continue;
           }

         if (0 == strcmp2(dirStr, "NONINTERLACED"))
           {
             if ( *specFlag & STF_ISPERMUTATION )
               {
                 FileFormatError("ReadSTFFile", "Error reading header: double INTERLACED/NONINTERLACED specifications.");
                 datasetType = STF_ERROR;
                 break;
               }

              headerLine++;
              ReadLn(inFile);
              datasetSpec->permutation = STF_NONINTERLACED;
              *specFlag |= STF_ISPERMUTATION;
              continue;
           }

         if (0 == strcmp2(dirStr, "NAME"))
           {
             if ( *specFlag & STF_ISNAME )
               {
                 FileFormatError("ReadSTFFile", "Error reading header: double NAME specifications.");
                 datasetType = STF_ERROR;
                 break;
               }

             headerLine++;
             SkipBlanks(inFile);

             if ( 1 != fscanf(inFile, "%s", datasetSpec->datasetName))
               {
                 FileFormatError("ReadSTFFile", "Error reading header: NAME specificaton format error.");
                 datasetType = STF_ERROR;
                 break;
               }
             ReadLn(inFile);
             *specFlag |= STF_ISNAME;
             continue;
           }

         if (0 == strcmp2(dirStr, "ORDER"))
           {
             if ( *specFlag & STF_ISORDER )
               {
                 FileFormatError("ReadSTFFile", "Error reading header: double ORDER specifications.");
                 datasetType = STF_ERROR;
                 break;
               }

             headerLine++;
             SkipBlanks(inFile);

             if ( 1 != fscanf(inFile, "%s", dirStr))
               {
                 FileFormatError("ReadSTFFile", "Error reading header: ORDER specification format error.");
                 datasetType = STF_ERROR;
                 break;
               }
             if ( 0 == strcmp2(dirStr, "COLUMN") )
               {
                 datasetSpec->order = STF_COLUMNMAJOR;
               }
             else if ( 0 == strcmp2(dirStr, "ROW") )
               {
                 datasetSpec->order = STF_ROWMAJOR;
               }
             else
               {
                 FileFormatError("ReadSTFFile", "Error reading header: invalid ORDER specification.  It should be COLUMN nor ROW.");
                 datasetType = STF_ERROR;
                 break;
               }

             ReadLn(inFile);
             *specFlag |= STF_ISORDER;
             continue;
           }

         if (0 == strcmp2(dirStr, "DIMENSIONS"))
           {
             if ( !(*specFlag & STF_ISRANK) )
               {
                 FileFormatError("ReadSTFFile", "Error reading header: RANK missing befor DIMENSIONS");
                 datasetType = STF_ERROR;
                 break;
               }
                 
             if ( *specFlag & STF_ISDIMS )
               {
                 FileFormatError("ReadSTFFile", "Error reading header: double DIMENSIONS specifications.");
                 datasetType = STF_ERROR;
                 break;
               }

             headerLine++;

            /* allocate storage for dimensions*/
            if ( !(datasetSpec->dimensions = (long *) Alloc(datasetSpec->rank * sizeof(long))) )
              {
                FileFormatError("ReadSTFFile", "Unable to allocate memory");
                datasetType = STF_ERROR;
                break;
              }

             *specFlag |= STF_ISDIMS;

             /* read all the dimensions */
             for(counter = 0; counter < datasetSpec->rank; counter++)
                {
                  SkipBlanks(inFile);
                  if ( 1 != fscanf(inFile, "%ld", datasetSpec->dimensions + counter) )
                    {
                      FileFormatError("ReadSTFFile", "Error reading header: DIMENSIONS value format error. ");
                      datasetType = STF_ERROR;
                      break;
                    }
                }

             if (datasetType == STF_ERROR)
               break;

             ReadLn(inFile);
             continue;
           }

         if (0 == strcmp2(dirStr, "BOUNDS"))
           {
             if ( *specFlag & STF_ISBOUNDS )
               {
                 FileFormatError("ReadSTFFile", "Error reading header: double BOUNDS specifications.");
                 datasetType = STF_ERROR;
                 break;
               }

             headerLine++;

             /* allocate storage for bounds*/
             if ( !(datasetSpec->bounds = (real *) Alloc(sizeof(real))) )
               {
                 FileFormatError("ReadSTFFile", "Unable to allocate memory.");
                 datasetType = STF_ERROR;
                 break;
               }

             *specFlag |= STF_ISBOUNDS;

             /* read all the bounds */
             for( counter = 0, SkipBlanks(inFile);
                  (c = getc(inFile))!= '\n' && c >= 0;
                  counter++, SkipBlanks(inFile)
                )
                {
                  ungetc(c, inFile);
                  if ( !( datasetSpec->bounds = (real *) Realloc(datasetSpec->bounds, (counter+1)*sizeof(real))) )
                    {
                      FileFormatError("ReadSTFFile", "Unable to allocate memory.");
                      datasetType = STF_ERROR;
                      break;
                    }
                  if ( 1 != fscanf(inFile, "%f", (float *) (datasetSpec->bounds + counter)) )
                    {
                      FileFormatError("ReadSTFFile", "Error reading header: BOUNDS value format error.");
                      datasetType = STF_ERROR;
                      break;
                    }
                 }
             
             if (datasetType == STF_ERROR)
               /* this break gets out of the outer for loop */
               break;

             /* the number of bounds must be an even number */
             if ( counter == 0 || counter % 2 != 0 )
               {
                 FileFormatError("ReadSTFFile", "Error reading header: BOUNDS value format error.");
                 datasetType = STF_ERROR;
                 break;
               }

             datasetSpec->nBounds = counter;
             continue;
           }

         if (0 == strcmp2(dirStr, "DATA"))
           {
             headerLine++;
             ReadLn(inFile);
             datasetType = STF_DATA;
             break;
           }

         if (0 == strcmp2(dirStr, "GRID"))
           {
             headerLine++;
             ReadLn(inFile);
             datasetType = STF_GRID;
             break;
           }

         if (0 == strcmp2(dirStr, "END"))
           {
             ReadLn(inFile);
             break;
           }

         FileFormatError("ReadSTFFile, unknown header line specification", dirStr);
         datasetType = STF_ERROR;
         break;
       }
     else
       {
         /* a blank line*/
         ReadLn(inFile);
         continue;
       }
  }

  if ( datasetType == STF_END )
    {
      if ( headerLine != 0 )
        {
          FileFormatError("ReadSTFFile",
          "Error reading header: Incomplete header specification for dataset.");
          datasetType = STF_ERROR;
        }
      else
        {
          return datasetType; 
        }
    }

  if ( 
       (datasetType == STF_ERROR) || 

       /* rank and dimensions are always required */
       ( !(*specFlag & STF_ISRANK) || !(*specFlag & STF_ISDIMS) ) ||

       /* for grid data set, fieldtype is required */
       /* and its fieldtype must be vector */
       ( (datasetType == STF_GRID) && (!(*specFlag & STF_ISTYPE) || (datasetSpec->fieldType != STF_VECTOR) ) ) ||
       
       /* for field data set, fieldtype is required */
       /* and for vector dataset, permutation is required */
       ( (datasetType == STF_DATA) && (!(*specFlag & STF_ISTYPE)) ) 
     )
    {
      switch ( datasetType )
        {
          case STF_ERROR:
            /* error message was broadcasted */
            break;
          case STF_DATA:
            if ( !(*specFlag & STF_ISRANK) )
              {
                FileFormatError("ReadSTFFile",
                "Error reading header: RANK missing before DATA.");
              }
            else if ( !(*specFlag & STF_ISDIMS) )
              {
                FileFormatError("ReadSTFFile",
                "Error reading header: DIMENSIONS missing before DATA.");
              }
            else if ( !(*specFlag & STF_ISTYPE) )
              {
                FileFormatError("ReadSTFFile",
                "Error reading header: SCALAR/VECTOR missing before DATA.");
              }
            break;
          case STF_GRID:
            if ( !(*specFlag & STF_ISRANK) )
              {
                FileFormatError("ReadSTFFile",
                "Error reading header: RANK missing before GRID.");
              }
            else if ( !(*specFlag & STF_ISDIMS) )
              {
                FileFormatError("ReadSTFFile",
                "Error reading header: DIMENSIONS missing before GRID.");
              }
            else if ( !(*specFlag & STF_ISTYPE) )
              {
                FileFormatError("ReadSTFFile",
                "Error reading header: VECTOR missing before GRID.");
              }
            else if ( datasetSpec->fieldType != STF_VECTOR )
              {
                FileFormatError("ReadSTFFile",
                "Error reading header: GRID must be vector type");
              }
            break;
          default:
              FileFormatError("ReadSTFFile",
              "Error reading header: header format error");
            break;
        }
      datasetType = STF_ERROR;
    }
  return datasetType;
} /* parseHeader */

#ifdef HAVE_PROTOTYPES
static int getStr(FILE *inFile, char *token)
#else
static int getStr(inFile, token) 
FILE *inFile;
char *token;
#endif
{
  int i;
  int c;

  /* skip comment lines and blanks */
  do 
    {
      while((c = fgetc(inFile)) != EOF && (c == ' ' || c == '\t'))
         ;
      if (c == '#')
        while((c = fgetc(inFile)) != EOF && c != '\n')
           ;
    } while(c == '\n');

  /* get a token */
  for (i = 0; c != EOF && c != ' ' && c != '\t' && c != '\n'; i++)
    {
      token[i] = c;
      c = fgetc(inFile);
    }

  if (c != EOF)
    {
      ungetc(c, inFile);
    }
      
  token[i] = '\0';
  return i;
} /* getStr */ 


#ifdef HAVE_PROTOTYPES
static real *readValue(FILE *inFile, STF_spec *datasetSpec, ObjPtr dataset)
#else
static real *readValue(inFile, datasetSpec, dataset)
FILE *inFile;
STF_spec *datasetSpec;
ObjPtr dataset;
#endif
/* return pointer to bounds of the dataset just read */
/* return null if an error occurred */
{

  real *bounds;
  long *index;
  int counter;
  int nComponents;

  
  /* determine number of components first */
  if ( datasetSpec->fieldType == STF_SCALAR )
    {
      nComponents = 1;
    }
  else
    {
      nComponents = datasetSpec->nComponents;
    }


  /* alloc bounds array and init it to contain 0 values */
  if ( !(bounds = (real *) Alloc( 2*nComponents*sizeof(real) )) )
    {
      FileFormatError("ReadSTFFile", "Unable to allocate memory");
      return NULL;
    }
  for ( counter = 0; counter < 2*nComponents; counter++)
     {
       bounds[counter] = 0.0;
     }


  /* alloc index array and init it to contain 0 values */
  if ( !(index = (long *) Alloc( datasetSpec->rank*sizeof(long) )) )
    {
      FileFormatError("ReadSTFFile", "Unable to allocate memory");
      return 0;
    }
  for ( counter = 0; counter < datasetSpec->rank; counter++)
     {
       index[counter] = 0L;
     }


    /* prepare to read in the field values from file */
    if ( !SetCurField(FIELD1, dataset) ) 
      {
        Free(index);
        Free(bounds);
        return NULL;
      }

    /* adjust the index and read in the corresponding field value */ 
    {
      real curData;
      char dataStr[32];
      int n; /* counter */
      
      counter = 0;

      if ( datasetSpec->fieldType == STF_SCALAR ||
                   datasetSpec->permutation == STF_NONINTERLACED )
        {
         for ( n = 0; n < nComponents; n++ )
           {
             counter = 0;
             while(getStr(inFile, dataStr)) 
             {
                int first = 1;
                int j;

                if ( !ParseReal( &curData, dataStr ) )
                  {
                    FileFormatError("ReadSTFFile", "Invalid numeric Value");
                    /* clean up and return an NULL ptr */
                    Free(index);
                    Free(bounds);
                    return NULL;
                  }

                PutFieldComponent(FIELD1, n, index, (real) curData);

                for (j = 0; first && j < datasetSpec->rank; j++)
                   first = index[j] == 0;
                if (first)
                  {
                    bounds[2 * n] = curData;
                    bounds[2 * n + 1] = curData;
                  }
                else
                  {
                    /* looking for bounds */
                    if ( curData < bounds[2*n] )
                      {
                        bounds[2*n] = curData;
                      }
                    else if ( curData > bounds[2*n+1] )
                      {
                        bounds[2*n+1] = curData;
                      }
                  }

                /* adjust the index according to column or row major */
                if( datasetSpec->order == STF_COLUMNMAJOR)
                  {
                    for(counter = 0; counter < datasetSpec->rank; counter++)
                      {
                        index[counter]++;
                        if (index[counter] == datasetSpec->dimensions[counter])
                          {
                           index[counter] = 0;
                           continue;
                          }
                        break;
                      }
                    if (counter == datasetSpec->rank)
                      {
                        /* This signifies the nth component data is all read */
                        break;
                      }
                  }
                else
                  {
                    for(counter = datasetSpec->rank - 1; counter > -1 ; counter--)
                      {
                        index[counter]++;
                        if (index[counter] == datasetSpec->dimensions[counter])
                          {
                           index[counter] = 0;
                           continue;
                          }
                        break;
                      }
                    if (counter == -1)
                      {
                        /* This signifies the nth component value is all read */
                        break;
                      }
                  }

             } /* while */

           }  /* for ( n = 0; n < nComponents; n++ ) */

         /* this indicates that the data number is less than dimensions */
         if ( ( datasetSpec->fieldType == STF_VECTOR && n != nComponents) ||
              ( datasetSpec->order == STF_COLUMNMAJOR && counter != datasetSpec->rank )||
              ( datasetSpec->order == STF_ROWMAJOR && counter != -1 ) )
           {
            FileFormatError("ReadSTFFile", "Inconsistent Data File, number of data values is less than that specified by dimensions");
            Free(index);
            Free(bounds);
            return NULL;
           }

        } /* if fieldType == STF_SCALAR || permutation == STF_NONINTERLACED */
      else /* fieldType is VECTOR and permutation is INTERLACED */ 
        {
           n = 0;
           counter = 0;
           while(getStr(inFile, dataStr))
           {
             int first = 1;
             int j;

              if ( !ParseReal( &curData, dataStr ) )
                {
                  FileFormatError("ReadSTFFile", "Invalid numeric Value");
                  /* clean up and return NULL */
                  Free(index);
                  Free(bounds);
                  return NULL;
                }

              PutFieldComponent(FIELD1, n, index, (real) curData);

              for (j = 0; first && j < datasetSpec->rank; j++)
                 first = index[j] == 0;
              if (first)
                {
                  bounds[2 * n] = curData;
                  bounds[2 * n + 1] = curData;
                }
              else
                {
                  /* looking for bounds */
                  if ( curData < bounds[2*n] )
                    {
                      bounds[2*n] = curData;
                    }
                  else if ( curData > bounds[2*n+1] )
                    {
                      bounds[2*n+1] = curData;
                    }
                 }

              n++;

              if ( n == nComponents ) /* time to update index */
                {
                 /* adjust the index according to column or row major */
                 if( datasetSpec->order == STF_COLUMNMAJOR)
                   {
                     for(counter = 0; counter < datasetSpec->rank; counter++)
                       {
                         index[counter]++;
                         if (index[counter] == datasetSpec->dimensions[counter])
                           {
                            index[counter] = 0;
                            continue;
                           }
                         break;
                       }
                     if (counter == datasetSpec->rank)
                       {
                         /* This signifies the field data is all read */
                         break;
                       }
                   }
                 else
                   {
                     for(counter = datasetSpec->rank - 1; counter > -1 ; counter--)
                       {
                         index[counter]++;
                         if (index[counter] == datasetSpec->dimensions[counter])
                           {
                            index[counter] = 0;
                            continue;
                           }
                         break;
                       }
                     if (counter == -1)
                       {
                         /* This signifies the field data is all read */
                         break;
                       }
                   }
 
                 n = 0;

                } /* if n == nComponents */

             } /* for fscanf */

         /* this indicates that the data number is less than dimensions */
         if ( ( datasetSpec->order == STF_COLUMNMAJOR && counter != datasetSpec->rank )||
              ( datasetSpec->order == STF_ROWMAJOR && counter != -1 ) )
           {
            FileFormatError("ReadSTFFile", "Inconsistent Data File, number of data values is less than that specified by dimensions");
            Free(index);
            Free(bounds);
            return NULL;
           }

        } /* else where fieldType is VECTOR and INTERLACED */


    }

    Free(index);
    return bounds;

}/* readValue */

/* for use with ftp reader only to accomodate time dependent grids */
#ifdef HAVE_PROTOTYPES
static ObjPtr RegisterNewTimeDataset(ObjPtr dataset, float time)
#else
static ObjPtr RegisterNewTimeDataset(dataset, time)
ObjPtr dataset;
float time;
#endif
{
  ObjPtr form, formDataset, formDatasetData, data;
  ObjPtr formTimeSteps, dataTimeSteps, formTimeData, dataTimeData;
  ObjPtr var, namevar;
  real *elements, *bounds;
  char *name;
  int rank, i, j;
  long *dims;
  int nComponents;

  ObjPtr priorDataset, priorDatasetData;
  ObjPtr priorForm, priorFormDataset, priorFormDatasetData;
  real *priorBounds, *priorDimElements;
  int priorRank, priorNComponents;

  /* get the DATA object of the dataset */
  MakeVar(dataset, DATA);
  data = GetVar(dataset, DATA);

  /* get the bounds of the form */
  MakeVar(dataset, DATAFORM);
  form = GetVar(dataset, DATAFORM);
  MakeVar(form, BOUNDS);
  var = GetArrayVar("RegisterNewTimeDataset", form, BOUNDS);
  rank = (DIMS(var))[0];
  if (!(bounds = (real *)Alloc(sizeof(real) * rank)))
    {
      return ObjFalse;
    }
  elements = ELEMENTS(var);
  for (i = 0; i < rank; i++)
     {
       bounds[i] = elements[i];
     }
  rank /= 2;

  /* get dimensions of the formDataset */
  if (!(dims = (long *) Alloc(sizeof(long) * rank)))
    {
      /* HWUHERE */
      return ObjFalse;
    }
  MakeVar(form, DIMENSIONS);
  var = GetArrayVar("RegisterNewTimeDataset", form, DIMENSIONS);
  elements = ELEMENTS(var);
  for (i = 0; i < rank; i++)
     {
       dims[i] = (long) elements[i];
     }

  MakeVar(form, NAME);
  var = GetVar(form, NAME);
  name = GetString(var);

  /* get the dataset used to make the dataform */
  MakeVar(form, DATA);
  formDataset = GetVar(form, DATA);
  if (NULLOBJ != formDataset)
    {
      /* get the DATA part of the formDataset */
      MakeVar(formDataset, DATA);
      formDatasetData = GetVar(formDataset, DATA);

      /* get the NCOMPONENTS of the formDataset */
      MakeVar(formDataset, NCOMPONENTS);
      var = GetIntVar("RegisterNewTimeDataset", formDataset, NCOMPONENTS);
      nComponents = GetInt(var);

      /* Create a new form for this to avoid change the formDataset of the 
	 original dataset in case the caller still depends on it create another 
	 dataset */
      formDataset = NewDataset(name, rank, dims, nComponents);
      SetVar(formDataset, DATA, formDatasetData);
      form = NewCurvilinearDataForm(name, rank, dims, bounds, formDataset);
      SetDatasetForm(dataset, form);
    }
  else
    {
      /* NULLOBJ == formDataset */
      /* must be a regular grid, construct a separable grid for it */
      real **scales;

      MakeVar(form, NAME);
      var = GetVar(form, NAME);
      name = GetString(var);

      /* build the scales for new separable dataform */
      if (!(scales = (real **) Alloc(sizeof(real *) * rank)))
        {
	  /* HWUHERE */
          return ObjFalse;
        }
      for (i = 0; i < rank; i++)
         {
           if (!(scales[i] = (real *) Alloc(sizeof(real) * dims[i])))
	     {
	       for (j = i - 1; j >= 0; j--)
		  {
		    Free(scales[j]);
		  }
               Free(scales);

	       /* HWUHERE */
	       return ObjFalse;
             }
         }
      for (i = 0; i < rank; i++)
	 {
	   real interval;
	   real min, max;

           min = bounds[2 * i];
	   max = bounds[2 * i + 1];
	   if (dims[i] > 1)
	     {
	       interval = (max - min) / ((float) dims[i] - 1.0);
	     }
           for (j = 0; j < dims[i] - 1; j++)
	      {
	        scales[i][j] = min + (real) j * interval;
	      }
           scales[i][j] = max;
	 }
      form = NewSeparableDataForm(name, rank, dims, scales);
      MakeVar(form, DATA);
      formDataset = GetVar(form, DATA);
      MakeVar(formDataset, DATA);
      formDatasetData = GetVar(formDataset, DATA);

      /* get the NCOMPONENTS of the formDataset */
      MakeVar(formDataset, NCOMPONENTS);
      var = GetIntVar("RegisterNewTimeDataset", formDataset, NCOMPONENTS);
      nComponents = GetInt(var);

      /* set this as the new form */
      SetDatasetForm(dataset, form);

      /* dispose of the memory */
      for (i = 0; i < rank; i++)
         {
           Free(scales[i]);
         }
      Free(scales);
    }

  /*Find out if there already exists a timed dataset like this*/
  MakeVar(dataset, NAME);
  namevar = GetVar(dataset, NAME);
  name = GetString(namevar);
  i = 1;
  while (NULLOBJ != (priorDataset = FindDatasetByName(GetString(namevar))))
       {
	 char nameString[TEMPSTRSIZE + 1];

         /* see if this is a timed dataset, if not then looking for name (2), name (3), ... 
	    until a timed dataset is found or NULLOBJ is reached */
         MakeVar(priorDataset, DATA);
         priorDatasetData = GetVar(priorDataset, DATA);
         if (!IsTimedObj(priorDatasetData))
           {
	     goto reset;
           }

         /* see if the RANK, NCOMPONENTS and DIMENSIONS of the priorDataset matches
	    the form of the dataset to be registered.  The least I can check! */
         MakeVar(priorDataset, DATAFORM);
         priorForm = GetVar(priorDataset, DATAFORM);

	 /* get the rank from the dim[0] of the bounds */
         MakeVar(priorForm, BOUNDS);
         var = GetArrayVar("RegisterNewTimeDataset", priorForm, BOUNDS);
         priorBounds = ELEMENTS(var);
         priorRank = DIMS(var)[0];
         priorRank /= 2;
         /* check to see if its rank matches the rank of the dataset to be registered */
         if (rank != priorRank)
           {
	     goto reset;
           }

          /* get dimensions of the priorFormDataset */
          MakeVar(priorForm, DIMENSIONS);
          var = GetArrayVar("RegisterNewTimeDataset", priorForm, DIMENSIONS);
          priorDimElements = ELEMENTS(var);
	  for (j = 0; j < priorRank; j++)
	     {
	       if (dims[j] != (long) priorDimElements[j])
		 {
		   goto reset;
		 }
	     }

         /* check spatial dimensionality now, i.e. ncomponents */
         MakeVar(priorForm, DATA);
         if (NULLOBJ == (priorFormDataset = GetVar(priorForm, DATA)))
	   {
	     priorNComponents = priorRank;
	   }
         else
	   {
	     MakeVar(priorFormDataset, NCOMPONENTS);
	     var = GetIntVar("RegisterNewTimeDataset", priorFormDataset, NCOMPONENTS);
	     priorNComponents = GetInt(var);
	   }
         if (nComponents != priorNComponents)
	   {
	     goto reset;
	   }
         
	 /* If everything matches so far, then assume the dataset found is the one
	    we are going to append the new time dataset to. */
         break;

reset:   i++;
         sprintf(nameString, "%s (%d)", name, i);
         namevar = NewString(nameString);
       }

  if (priorDataset)
    {
      /* check the DATAFORM of the matched time dependent dataset, to make sure
	 other file readers does the necessaries to set up the timedobj of the 
	 as the DATA of the dataset composing the dataform.  
       SIDE EFFECT: 
	 This will change the DATA of DATAFORM of the matched dataset to a 
	 time dependent dataset if it's not.
      */

      if (NULLOBJ != priorFormDataset)
	{
          MakeVar(priorFormDataset, DATA);
          priorFormDatasetData = GetVar(priorFormDataset, DATA);
	}
      else
      /* NULLOBJ == priorFormDataset */
	{
          real **scales;
          long *priorDims;

	  /* The prior form must be a regular grid, build it as though it's a 
	     separable grid, and build up a timesteps for it */
          /* and point the priorFormData to it */
	  if (NULL == (priorDims = (long *) Alloc(sizeof(long) * priorRank)))
	    {
	      /* HWUHERE */
	      return ObjFalse;
	    }
          for (i = 0; i < priorRank; i++)
	     {
	       priorDims[i] = (long) priorDimElements[i];
	     }

          MakeVar(priorForm, NAME);
          var = GetVar(priorForm, NAME);
          name = GetString(var);

          if (!(scales = (real **) Alloc(sizeof(real *) * priorRank)))
            {
              return ObjFalse;
            }
          for (i = 0; i < priorRank; i++)
             {
               if (!(scales[i] = (real *) Alloc(sizeof(real) * priorDims[i])))
	         {
	           for (j = i - 1; j >= 0; j--)
		      {
		        Free(scales[j]);
		      }
                   Free(scales);
	           return ObjFalse;
	         }
             }
          for (i = 0; i < priorRank; i++)
	     {
	       real interval;
	       real min, max;

               min = priorBounds[2 * i];
	       max = priorBounds[2 * i + 1];
	       if (priorDimElements[i] > 1.0)
	         {
	           interval = (max - min) / ((float) priorDims[i] - 1.0);
	         }
               for (j = 0; j < priorDims[i] - 1; j++)
	          {
	            scales[i][j] = min + (real) j * interval;
	          }
               scales[i][j] = max;
	     }

          priorForm = NewSeparableDataForm(name, priorRank, priorDims, scales);
          MakeVar(priorForm, DATA);
          priorFormDataset = GetVar(priorForm, DATA);
	  MakeVar(priorFormDataset, DATA);
	  priorFormDatasetData = GetVar(priorFormDataset, DATA);

	  /* set this as the new form */
	  SetDatasetForm(priorDataset, priorForm);

	  for (i = 0; i < priorRank; i++)
	     {
	       Free(scales[i]);
	     }
          Free(scales);
	  Free(priorDims);
	}

      /* built the DATA part of the formDataset as a timed object */
      if (!IsTimedObj(priorFormDatasetData))
	{
	  int nSteps, n;
	  real *dataTimes, *dims;
	  real *elements;

          MakeVar(priorDatasetData, TIMESTEPS);
          dataTimeSteps = GetArrayVar("InsertTimeSlice", var, TIMESTEPS);
	  dataTimes = ELEMENTS(dataTimeSteps);
	  nSteps = (int) (DIMS(dataTimeSteps))[0];

          formTimeData = NewList();
	  formTimeSteps = NewList();

	  for (n = 0; n < nSteps; n++)
	     {
               PrefixList(formTimeSteps, NewReal(dataTimes[n]));
	       /* repeatingly use the same form data set */
               PrefixList(formTimeData, priorFormDatasetData);
	     }
          SetVar(priorFormDataset, DATA, NewTimedObject(formTimeSteps, formTimeData));
	}


      /* reset the bounds */
      for (i = 0; i < rank; i++)
	 {
	   if (bounds[2 * i] < priorBounds[2 * i])
	     {
	       priorBounds[2 * i] = bounds[2 * i];
	     }
           if (bounds[2 * i + 1] > priorBounds[2 * i + 1])
	     {
	       priorBounds[2 * i + 1] = bounds[2 * i + 1];
	     }
	 }

      /* insert the dataset and formDataset to their list respectively */
      InsertDatasetTimeStep(priorFormDataset, formDataset, time); 
      InsertDatasetTimeStep(priorDataset, dataset, time);
    }
  else
    {
      /*Make a new dataset*/
      dataTimeSteps = NewList();
      dataTimeData = NewList();
      formTimeSteps = NewList();
      formTimeData = NewList();

      PrefixList(dataTimeSteps, NewReal(time));
      PrefixList(formTimeSteps, NewReal(time));
      PrefixList(dataTimeData, data);
      PrefixList(formTimeData, formDatasetData);
         
      SetVar(dataset, DATA, NewTimedObject(dataTimeSteps, dataTimeData));
      SetVar(formDataset, DATA, NewTimedObject(formTimeSteps, formTimeData));
      RegisterNewDataset(dataset);
    }

  Free(bounds);
  Free(dims);
  return ObjTrue;
}

static ObjPtr ReadSTFFile(reader, name)
ObjPtr reader;
char *name;
{
    int datasetType; /* DATA, GRID, END or ERROR */
    STF_spec datasetSpec;
    int specFlag = 0;

    FILE *inFile;

    real *gridBounds = NULL;
    real *dataMinmaxes = NULL;

    int isGrid = 0, isData = 0;
    int gridComponents; /* the number of components for the grid if any */
    int counter;
    int more = 1;

    ObjPtr dataForm;
    ObjPtr gridDataset;
    ObjPtr retVal;

    inFile = fopen(name, "r");
    if ( !inFile )
      {
        Error("ReadSTFFile", OPENFILEERROR, name);
        return NULLOBJ;
      }
     
    /* set default dataset specs */
    datasetSpec.rank = 0;
    datasetSpec.nComponents = 0;
    datasetSpec.order = STF_COLUMNMAJOR;
    datasetSpec.permutation = STF_INTERLACED;
    strcpy(datasetSpec.datasetName, name);

    while( more ) 
      {
        /* parse header lines */
        datasetType = parseHeader(inFile, &specFlag, &datasetSpec);

        switch ( datasetType )
          {
            case STF_ERROR: 
               FileFormatError("Error in dataset", datasetSpec.datasetName);
               more = 0;
               break;
            case STF_GRID:
               /* if isGrid then error, double grid */
               if ( isGrid )
                 {
                   FileFormatError("ReadSTFFile", "Multiple Grid Spec");
                   FileFormatError("Error in dataset", datasetSpec.datasetName);
                   more = 0;
                   break;
                 }

               gridComponents = datasetSpec.nComponents;

               /* create the grid dataset */
               gridDataset = NewStructuredDataset(datasetSpec.datasetName, datasetSpec.rank, datasetSpec.dimensions, gridComponents);

               /* something wrong with grid value format */
               if ( !(gridBounds = readValue( inFile, &datasetSpec, gridDataset)))
                 {
                   FileFormatError("Error in dataset", datasetSpec.datasetName);
                   more = 0;
                   break;
                 }

               /* reset specFlag */
               specFlag &= (~STF_ISTYPE & ~STF_ISPERMUTATION & ~STF_ISORDER);
               isGrid = 1;
               break;
            case STF_DATA:
               /* if isData then error, double data */
               if ( isData )
                 {
                   FileFormatError("ReadSTFFile", "Multiple Data Spec");
                   FileFormatError("Error in dataset", datasetSpec.datasetName);
                   more = 0;
                   break;
                 }

               /* create the dataset */
               retVal = NewStructuredDataset(datasetSpec.datasetName, datasetSpec.rank, datasetSpec.dimensions, datasetSpec.nComponents);

               /* something wrong with data value format */
               if ( !(dataMinmaxes = readValue( inFile, &datasetSpec, retVal)) )
                 {
                   FileFormatError("Error in dataset", datasetSpec.datasetName);
                   more = 0;
                   break;
                 }
               else
                 {
                   /* has no use with dataMinmaxes */
                   Free(dataMinmaxes);
                 }

               /* reset specFlag */
               specFlag &= (~STF_ISTYPE & ~STF_ISPERMUTATION & ~STF_ISORDER);
               isData = 1;
               break;
            case STF_END:
               /* if isData && isGrid then use curvilinear dataForm */
               /* if isData && !isGrid then regular dataForm */
               /* for this case, the bounds are required, error if no */
               /* if !isData && isGrid the error */
               /* if !isData && !isGrid then means end of file */
                       
               if ( !isData && !isGrid )
                 /* end of dataset file reached */
                 {
                   more = 0;
                   break;
                 }

               if ( !isData && isGrid )
                 {
                   FileFormatError("ReadSTFFile", "NO DATA");
                   FileFormatError("Error in dataset", datasetSpec.datasetName);
                   more = 0;
                   break;
                 }

               if ( specFlag & STF_ISBOUNDS )
                 {
                   int nBounds; /* number of bounds */

                   /* get rid of calculated grid bounds if any and 
                      use the bounds provided by the user */
                   /* determine the number of bounds */
                   if ( isGrid )
                     {
                       Free(gridBounds);
                       gridBounds = NULL;
                       nBounds = 2 * gridComponents;
                     }
                   else
                     {
                       nBounds = 2 * datasetSpec.rank;
                     }

                   /* test if the user provide correct bounds */
                   if ( nBounds != datasetSpec.nBounds )
                     {
                       FileFormatError("ReadSTFFile", "BOUNDS format error");
                       FileFormatError("Error in dataset", datasetSpec.datasetName);
                       more = 0;
                       break;
                     }

                   gridBounds = datasetSpec.bounds;
                 }


               if ( isData && isGrid )
                 /* for this case, bounds are optional */
                 {
                   dataForm = NewCurvilinearDataForm(datasetSpec.datasetName, datasetSpec.rank, datasetSpec.dimensions, gridBounds, gridDataset);
                 }
               else
                 {
                   /* isData && !isGrid */      
                   /* for this case bounds are required */
                   if ( !(specFlag & STF_ISBOUNDS) )
                     {
                       FileFormatError("ReadSTFFile", "BOUNDS required for dataset with no grid");
                       FileFormatError("Error in dataset", datasetSpec.datasetName);
                       more = 0;
                       break;
                     }
                   dataForm = NewRegularDataForm(datasetSpec.datasetName, datasetSpec.rank, datasetSpec.dimensions, gridBounds);
                 }

               SetDatasetForm(retVal, dataForm);
               if ( specFlag & STF_ISTIME )
                 {
		   RegisterNewTimeDataset(retVal, datasetSpec.time);
#if 0
                   RegisterTimeDataset(retVal, datasetSpec.time);
#endif
                 }
               else
                 {
                   RegisterDataset(retVal);
                 }

               isGrid = 0;
               isData = 0;
               specFlag = 0;
               Free(datasetSpec.dimensions);
               Free(gridBounds);
               gridBounds = NULL;

               /* set default dataset specs */
               datasetSpec.rank = 0;
               datasetSpec.nComponents = 0;
               datasetSpec.order = STF_COLUMNMAJOR;
	       datasetSpec.permutation = STF_INTERLACED;

               break;
            default:     
               /* dummy case */
               break;
          } /* switch */
      } /* while */

   if ( specFlag & STF_ISDIMS )
     {
       Free(datasetSpec.dimensions);
     }

   if ( specFlag & STF_ISBOUNDS )
     {
       Free(datasetSpec.bounds);
     }
   if ( isGrid )
     {
       Free(gridBounds);
     }
   fclose(inFile);
   return NULLOBJ;
}

/*******************************************************************************/
/********************************P3D READER ************************************/
/*******************************************************************************/
/* Tzong-Yow Hwu */
/* 10/26/1992 */
/* plot3d file reader */
/* 9/7/1993, added SGI binary format */

/* file types */
#define P3D_GRID        0
#define P3D_FUNC        1
#define P3D_Q           2

/* file format types, default on IRIS 4D is unformatted and on UNIX is binary */
/* for both grid and Q files */
#define P3D_DAT         0  /* default */ 
#define P3D_BIN         1 
#define P3D_FMT         2

/* data permutation method, default is the whole */
/* PLANES has no effect on reading other than 3D dataset */
#define P3D_WHOLE       0   /* default */   /* the whole dataset at once */
#define P3D_PLANE       1   /* one plane at a time */

/* Control over default checking of Q data, default is to check */
/* for Q files only */
#define P3D_CHECK       0  /* default */
#define P3D_NOCHECK     1


/* Jacobian control */
/* for Q files only */
#define P3D_NOJACOBIAN  0 /* default */
#define P3D_JACOBIAN    1 

/* IBLANK control */
/* for grid files only */
#define P3D_NOBLANK     0
#define P3D_BLANK       1

#define P3D_MAXNAMELEN 256 

typedef struct P3D_grid
{
  char name[P3D_MAXNAMELEN];
  int fileType;
  int permutation;
  int blank;
} p3d_grid;

typedef struct P3D_q
{
  char name[P3D_MAXNAMELEN];
  int isMDataset;
  int fileType;
  int permutation;
  int check;
  int jacobian;
} p3d_q;

typedef struct P3D_func
{
  char name[P3D_MAXNAMELEN];
  int fileType;
  int permutation;
} p3d_func;

/* a p3d set consists of a grid spec and several Q and/or func spec's */
typedef struct P3D_set
   {
     int ndim;
     int isMGrid;    /* 0 for no and 1 for yes */
     int isGrid;     /* is the grid present */
     p3d_grid *grid; /* pointer to the grid spec */
     int nQ;         /* number of Q files for this p3d set */
     p3d_q *qlist;   /* pointer to an array of Q spec */
     int nFunc;      /* number of func files for this p3d set */
     p3d_func *funclist; /* pointer to an array of func spec */
   } p3d_set;

/*
   1. getToken read to next occurrence of slash '/' or eof, and pass
      the string back, return 0 if end of file is encountered, -1 
      if error format, otherwise the number of chars in token.
   2. only process READ commands 
   3. newline without preceding '-' is treated as command terminator
*/
#ifdef HAVE_PROTOTYPES
static int getToken(FILE *inFile, char *token)
#else
static int getToken(inFile, token) 
FILE *inFile;
char *token;
#endif
{
  int i;
  int c;
  int inQuotes = 0;

  for ( i = 0; ( c = fgetc(inFile) ) != EOF ; i++ ) 
     {
       if (c == ' ' || c == '\t')
         {
          i--;
          continue;
         }
       else if (c == '\n')
         {
           if (i > 0 && token[i - 1] == '-')
             /* if continuation on next line */
             {
               i -= 2;
               continue;
             }
           else if (inQuotes)
             {
               fprintf(stderr, "File format error, unbalanced quote\n");
               return -1;
             }
           else if (i == 0)
             /* skip blank lines */
             {
               i--;
               continue;
             }
           else
             {
               break;
             }
         }
       else if (c == '\"')
         {
           /* flip inQuotes flag */
           inQuotes = !inQuotes;
           /* discard the '"' char */
           i--;
           continue;
         }
       else if (c == '/')
         {
           if (!inQuotes)
             {
               break;
             }
           /* otherwise, fall through to get the char */
         }
       token[i] = c;
     }
  
  if ( i == 0 )
    if ( c == EOF )
      {
        /* end of file reached */
        return 0;
      }
    else
      {
        /* error file formate */
        /* call ERROR routine */
        fprintf(stderr, "File format error\n");
        return -1;
      }

  token[i] = '\0';
  return i;
} /* getToken */ 

/* getvalue return the number of values repeated in the form
   of num*val and returns 0 if error or EOF encountered */
#ifdef HAVE_PROTOTYPES
static int getValue(FILE *inFile, double *val)
#else
static int getValue(inFile, val)
FILE *inFile;
double *val;
#endif
{
  int c;
  int count;
  char s[180];
  char *end, *newstr;

  /* skip blanks and commas */
  while(isspace(c = fgetc(inFile)) || c == ',')
     ;
  if (c != EOF)
    ungetc(c, inFile);
  else
    return EOF;

  /* get the numerical string */
  if (fscanf(inFile, "%s", s) != 1)
    {
      return 0;
    }
  *val = strtod(s, &end);
  if (end == s)
    {
      /* nothing is converted */
      return 0;
    }
  if (end[0] == '\0' || (end[0] == ',' && end[1] == '\0'))
    {
      return 1;
    }
  else if (*end == '*')
    {
      count = (int) *val;
      newstr = end + 1;
      *val = strtod(newstr, &end);
      if (end == newstr || !(*end == '\0' || (*end = ',' && *(end + 1) == '\0')))
        {
          /* nothing is converted || format error */
          return 0;
        }
      return count;
    }
  else
    {
      return 0;
    }
} /* getValue */

/* strncmp2 compares two string for n chars, case insensitive */
#ifdef HAVE_PROTOTYPES
static int strncmp2(char *s, char *t, int n)
#else
static int strncmp2(s, t, n)
char *s;
char *t;
int n;
#endif
{
  int i;

  for (i = 0; i < n && toupper(s[i]) == toupper(t[i]); i++)
     {
       if (s[i] == '\0')
         return 0;
     }

  if (i == n)
    return 0;
  else
    return (s[i] - t[i]);
}

#ifdef HAVE_PROTOTYPES
static void cleanup(int nSets, p3d_set *sets)
#else
static void cleanup(nSets, sets) 
int nSets;
p3d_set *sets;
#endif
/* cleanup free the memory alloced in the program */
{
  int i;

  if (!sets)
    return;

  for (i = 0; i < nSets; i++)
     {
       if ((sets + i)->nQ)
          {
            Free((sets + i)->qlist);
          }
       if ((sets + i)->nFunc)
          {
            Free((sets + i)->funclist);
          }
       if ((sets + i)->isGrid)
         {
            Free((sets + i)->grid);
         }
     }

  Free(sets);

  return;
} /* cleanup */ 

#ifdef HAVE_PROTOTYPES
static void compress(int *nSets, p3d_set **sets)
#else
static void compress(nSets, sets) 
int *nSets;
p3d_set **sets;
#endif
{
  int i, j;

  /* now, sort the sets according to their grid filename first */
  for (i = 1; i < *nSets; i++)
    {
     p3d_set temp;

     temp = (*sets)[i];
     /* filename is case sensitive, so use strcmp not strcmp2 */
     for (j = i - 1; j >= 0 && strcmp((*sets)[j].grid->name,
                                       temp.grid->name) > 0; j--)
        {
          (*sets)[j + 1] = (*sets)[j];
        }
     (*sets)[j + 1] = temp;
    }

  /* second compress the sets, so that if any of these sets use the
     same grid file, making them as the same set.
  */
  {
    p3d_set *last = *sets;
    p3d_set *temp;

    for (i = 1, j = 1; i < *nSets; i++)
      {
        p3d_set *curset = *sets + i;

        if (!strcmp(curset->grid->name, last->grid->name))
          {
            /* combine this set into last set */
            if (curset->nQ > 0)
              {
                int k;

                if (!(last->qlist = (p3d_q *) Realloc(last->qlist, 
                                  (last->nQ + curset->nQ) * sizeof(p3d_q))))
                  {
                    fprintf(stderr, "Unable to alloc memory, aborted\n");
                    return;
                  }
                for (k = 0; k < curset->nQ; k++)
                  {
                    last->qlist[last->nQ + k] = curset->qlist[k];
                  }
                Free(curset->qlist);
                curset->qlist = NULL;
                last->nQ += curset->nQ;
              }

            if (curset->nFunc > 0)
              {
                int k;

                if (!(last->funclist = (p3d_func *) Realloc(last->funclist, 
                            (last->nFunc + curset->nFunc) * sizeof(p3d_func))))
                  {
                    fprintf(stderr, "Unable to alloc memory, aborted\n");
                    return;
                  }
                for (k = 0; k < curset->nFunc; k++)
                  {
                    last->funclist[last->nFunc + k] = curset->funclist[k];
                  }
                Free(curset->funclist);
                curset->funclist = NULL;
                last->nFunc += curset->nFunc;
              }
            
            /* mark this set has been combined */
            Free(curset->grid);
            curset->grid = NULL;
          }
        else
          {
            last[j++] = *curset; 
          }
      } /* for */
      

    /* reset sets and nSets */
    /* cannot use realloc for *sets, since if it fails, then the q and func lists 
       is unaccessible and becomes garbages occupying the memory */
    if (!(temp = (p3d_set *) Alloc(j * sizeof(p3d_set))))
      {
        fprintf(stderr, "Unable to alloc memory, aborted\n");
        return;
      }
    for (i = 0; i < j; i++)
      {
        temp[i] = (*sets)[i];
      }
    Free(*sets); 
    *sets = temp;
    *nSets = j;
  }

#ifdef DEBUG

  for (i = 0; i < *nSets; i++)
    {

      p3d_set *curset = *sets + i;
      int isGrid = curset->isGrid;
      p3d_grid *grid = curset->grid;
      int nQ = curset->nQ;
      p3d_q *qlist = curset->qlist;
      int nFunc = curset->nFunc;
      p3d_func *funclist = curset->funclist;
      int j;

      printf("READ");
      /* ndim */
      if (curset->ndim == 1)
        {
          printf("/1D"); 
        }
      else if (curset->ndim == 2)
        {
          printf("/2D");
        }
      else if (curset->ndim == 3)
        {
          printf("/3D");
        }
    
      /* MGRID */
      if (curset->isMGrid)
        {
          printf("/MGRID");
        }

      /* process the XYZ file */
      if (isGrid)
        {
          /* file format */
          if (grid->fileType == P3D_DAT)
            {
              printf("/UNFORMATTED");
            }
          else if (grid->fileType == P3D_BIN)
            {
              printf("/BINARY");
            }
          else if (grid->fileType == P3D_FMT)
            {
              printf("/FORMATTED");
            }

          /* permutation */
          if (grid->permutation == P3D_WHOLE)
            {
              printf("/WHOLE");
            }
          else if (grid->permutation == P3D_PLANE)
            {
              printf("/PLANES");
            }

          /* BLANK */
          if (grid->permutation == P3D_BLANK)
            {
              printf("/BLANK");
            }
          else if (grid->permutation == P3D_NOBLANK)
            {
              printf("/NOBLANK");
            }

          printf("/XYZ=%s", grid->name);
        }
      else
         /* should not be the case */
        {
          fprintf(stderr, "XYZ file missing\n");
          continue;  /* do nothing for this set */
        }

      /* process the Q files */
      for (j = 0; j < nQ; j++)
        {
          /* MDATASET */
          if (qlist[j].isMDataset)
            {
              printf("/MDATASET");
            }

          /* file format */
          if (qlist[j].fileType == P3D_DAT)
            {
              printf("/UNFORMATTED");
            }
          else if (qlist[j].fileType == P3D_BIN)
            {
              printf("/BINARY");
            }
          else if (qlist[j].fileType == P3D_FMT)
            {
              printf("/FORMATTED");
            }

          /* permutation */
          if (qlist[j].permutation == P3D_WHOLE)
            {
              printf("/WHOLE");
            }
          else if (qlist[j].permutation == P3D_PLANE)
            {
              printf("/PLANES");
            }

          /* CHECK */
          if (qlist[j].check)
            {
              printf("/CHECK");
            }
          else
            {
              printf("/NOCHECK");
            }

          /* JACOBIAN */
          if (qlist[j].jacobian)
            {
              printf("/JACOBIAN");
              printf("\nsorry, jacobian is not implemented\n");
            }
          else
            {
              printf("/NOJACOBIAN");
            }

          printf("/Q=%s", qlist[j].name);
        }

      /* process the FUNCTION files */
      for (j = 0; j < nFunc; j++)
        {
          /* file format */
          if (funclist[j].fileType == P3D_DAT)
            {
              printf("/UNFORMATTED");
            }
          else if (funclist[j].fileType == P3D_BIN)
            {
              printf("/BINARY");
            }
          else if (funclist[j].fileType == P3D_FMT)
            {
              printf("/FORMATTED");
            }

          /* permutation */
          if (funclist[j].permutation == P3D_WHOLE)
            {
              printf("/WHOLE");
            }
          else if (funclist[j].permutation == P3D_PLANE)
            {
              printf("/PLANES");
            }

          printf("/FUNCTION=%s", funclist[j].name);
        }

      printf("\n");
    } /* outter for */

#endif   /* for ifdef DEBUG */

  return;
} /* compress */ 


/* return 0 on failure */
#ifdef HAVE_PROTOTYPES
static int parseP3DMeta(char *p3dname, int *nSets, p3d_set **sets)
#else
static int parseP3DMeta(p3dname, nSets, sets)
char *p3dname;
int *nSets;
p3d_set **sets;
#endif
{

  FILE *inFile;
  char token[P3D_MAXNAMELEN];
  int isNdim = 0, ndim = 3;  /* default is 3D */
  int isMDataset = 0;  /* time dependent datasets */
  int isMGrid = 0;    /* multiple grid files */  /* How do deal with it in Scian? */
  int isFileType = 0, fileType = P3D_DAT; /* file type, default is UNFORMATTED */
  int isPermutation = 0, permutation = P3D_WHOLE; /* data permutation way */
  int isCheck = 0, check = P3D_NOCHECK;  /* check control for Q file */
  int isJacobian = 0, jacobian = P3D_NOJACOBIAN;  /* jacobian control for grid file*/
  int isBlank = 0, blank = P3D_NOBLANK;  /* blank control for grid file*/
  int status;
  int valid = 1; /* if this file format is valid */


  if ( !(inFile = fopen(p3dname, "r")) )
    {
      fprintf(stderr, "ReadP3DFile: %s\n", p3dname);
      return 0;
    }

  /* the first thing the .p3d file should be READ */
  if ((status = getToken(inFile, token)) <= 0 || strcmp2(token, "READ"))
    { 
      fprintf(stderr, "ReadP3DFile: %s, file format error\n", p3dname);
      fclose(inFile);
      return 0;
    }

  /* now since we know that this is plot 3d file, lets process it */
  /* alloc memory for the first p3d set */
  if (!(*sets = (p3d_set *) Alloc(sizeof(p3d_set))))
    {
      fprintf(stderr, "Unable to alloc memory, aborted\n");
      fclose(inFile);
      return 0;
    }
  (*sets)->ndim = 3;
  (*sets)->isMGrid = 0;
  (*sets)->isGrid = 0;
  (*sets)->grid = NULL;
  (*sets)->nQ = 0;
  (*sets)->qlist = NULL;
  (*sets)->nFunc = 0;
  (*sets)->funclist = NULL;
  *nSets = 1;

  while ( (status = getToken(inFile, token)) > 0 )
    {
       /* for dimensions */
       if (!strcmp2(token, "1D") || !strcmp2(token, "1") ||
           !strcmp2(token, "2D") || !strcmp2(token, "2") ||
           !strcmp2(token, "3D") || !strcmp2(token, "3"))
         {
           p3d_set *curset = *sets + *nSets - 1;

           /* all number of dimensions for the same set must be consistent */
           if (isNdim && ndim != token[0] - '0')
             {
               fprintf(stderr, "ReadP3DFile: %s, inconsistent dim spec.\n", p3dname);
               valid = 0;
               break; /* get out of the while loop */ 
             }

           isNdim = 1;
           ndim = token[0] - '0';
           continue;
         }

       /* for MDATASET */
       if (!strcmp2(token, "MD") || !strcmp2(token, "MDA") ||
           !strcmp2(token, "MDAT") || !strcmp2(token, "MDATA") ||
           !strcmp2(token, "MDATAS") || !strcmp2(token, "MDATASE") ||
           !strcmp2(token, "MDATASET"))
         {
           isMDataset = 1;
           continue;
         }

       /* for MGRID */
       if (!strcmp2(token, "MG") || !strcmp2(token, "MGR") ||
           !strcmp2(token, "MGRI") || !strcmp2(token, "MGRID"))
         {
           /* all files for this set must be specified as mgrid */
           isMGrid = 1;
           continue;
         }

       /* for file format */
       if (!strcmp2(token, "FO") || !strcmp2(token, "FOR") ||
           !strcmp2(token, "FORM") || !strcmp2(token, "FORMA") || 
           !strcmp2(token, "FORMAT") || !strcmp2(token, "FORMATT") || 
           !strcmp2(token, "FORMATTE") || !strcmp2(token, "FORMATTED"))

         {
           if (isFileType && fileType != P3D_FMT)
             {
                fprintf(stderr, "ReadP3DFile: %s, double file format spec.\n", p3dname);
                valid = 0;
                break; /* get out of the while loop */ 
             }
           isFileType = 1;
           fileType = P3D_FMT;
           continue;
         }
       if (!strcmp2(token, "U") || !strcmp2(token, "UN") ||
           !strcmp2(token, "UNF") || !strcmp2(token, "UNFO") ||
           !strcmp2(token, "UNFOR") || !strcmp2(token, "UNFORM") ||
           !strcmp2(token, "UNFORMA") || !strcmp2(token, "UNFORMAT") ||
           !strcmp2(token, "UNFORMATT") || !strcmp2(token, "UNFORMATTE") ||
           !strcmp2(token, "UNFORMATTED"))
         {
           if (isFileType && fileType != P3D_DAT)
             {
                fprintf(stderr, "ReadP3DFile: %s, double file format spec.\n", p3dname);
                valid = 0;
                break; /* get out of the while loop */ 
             }
           isFileType = 1;
           fileType = P3D_DAT;
           continue;
         }
       if (!strcmp2(token, "BI") || !strcmp2(token, "BIN") ||
           !strcmp2(token, "BINA") || !strcmp2(token, "BINAR") ||
           !strcmp2(token, "BINARY"))
         {
           if (isFileType && fileType != P3D_BIN)
             {
                fprintf(stderr, "ReadP3DFile: %s, double file format spec.\n", p3dname);
                valid = 0;
                break; /* get out of the while loop */ 
             }
           isFileType = 1;
           fileType = P3D_BIN;
           continue;
         }


       /* for data permutation */
       if (!strcmp2(token, "P") || !strcmp2(token, "PL") ||
           !strcmp2(token, "PLA") || !strcmp2(token, "PLAN") ||
           !strcmp2(token, "PLANE") || !strcmp2(token, "PLANES"))
         {
           if (isPermutation && permutation != P3D_PLANE)
             {
                fprintf(stderr, "ReadP3DFile: %s, double permutation spec.\n", p3dname);
                valid = 0;
                break; /* get out of the while loop */ 
             }
           isPermutation = 1;
           permutation = P3D_PLANE;
           continue;
         }
       if (!strcmp2(token, "W") || !strcmp2(token, "WH") ||
           !strcmp2(token, "WHO") || !strcmp2(token, "WHOL") ||
           !strcmp2(token, "WHOLE"))
         {
           if (isPermutation && permutation != P3D_WHOLE)
             {
                fprintf(stderr, "ReadP3DFile: %s, double permutation spec.\n", p3dname);
                valid = 0;
                break; /* get out of the while loop */ 
             }
           isPermutation = 1;
           permutation = P3D_WHOLE;
           continue;
         }

       /* for CHECK or not */
       if (!strcmp2(token, "C") || !strcmp2(token, "CH") ||
           !strcmp2(token, "CHE") || !strcmp2(token, "CHEC") ||
           !strcmp2(token, "CHECK"))
         {
           if (isCheck && check != P3D_CHECK)
             {
                fprintf(stderr, "ReadP3DFile: %s, double check spec.\n", p3dname);
                valid = 0;
                break; /* get out of the while loop */ 
             }
           isCheck = 1;
           check = P3D_CHECK;
           continue;
         }
       if (!strcmp2(token, "NOC") || !strcmp2(token, "NOCH") ||
           !strcmp2(token, "NOCHE") || !strcmp2(token, "NOCHEC") ||
           !strcmp2(token, "NOCHECK"))
         {
           if (isCheck && check != P3D_NOCHECK)
             {
                fprintf(stderr, "ReadP3DFile: %s, double check spec.\n", p3dname);
                valid = 0;
                break; /* get out of the while loop */ 
             }
           isCheck = 1;
           check = P3D_NOCHECK;
           continue;
         }

       /* for Jacobian spec */
           /* jacobian is not implemented */
       if (!strcmp2(token, "J") || !strcmp2(token, "JA") ||
           !strcmp2(token, "JAC") || !strcmp2(token, "JACO") ||
           !strcmp2(token, "JACOB") || !strcmp2(token, "JACOBI") ||
           !strcmp2(token, "JACOBIA") || !strcmp2(token, "JACOBIAN"))
         {
           /* 
           isJacobian = 1;
           check = P3D_JACOBIAN;
           continue;
           */
           fprintf(stderr, "ReadP3DFile: %s, Jacobian is unimplemented\n", p3dname);
           valid = 0;
           break;
         }
       if (!strcmp2(token, "NOJ") || !strcmp2(token, "NOJA") ||
           !strcmp2(token, "NOJAC") || !strcmp2(token, "NOJACO") ||
           !strcmp2(token, "NOJACOB") || !strcmp2(token, "NOJACOBI") ||
           !strcmp2(token, "NOJACOBIA") || !strcmp2(token, "NOJACOBIAN"))
         {
           /*
           if (isJacobian)
             {
                fprintf(stderr, "ReadP3DFile: %s, double Jacobian spec.\n", p3dname);
                valid = 0;
                break; 
             }
           isJacobian = 1;
           check = P3D_NOJACOBIAN;
           */
           /* this is default */
           continue;
         }

       /* for blank spec */
       if (!strcmp2(token, "BL") || !strcmp2(token, "BLA") ||
           !strcmp2(token, "BLAN") || !strcmp2(token, "BLANK"))
         {
           if (isBlank && blank != P3D_BLANK)
             {
                fprintf(stderr, "ReadP3DFile: %s, double Blank spec.\n", p3dname);
                valid = 0;
                break; /* get out of the while loop */ 
             }
           isBlank = 1;
           blank = P3D_BLANK;
           continue;
         }
       if (!strcmp2(token, "NOB") || !strcmp2(token, "NOBL") ||
           !strcmp2(token, "NOBLA") || !strcmp2(token, "NOBLAN") ||
           !strcmp2(token, "NOBLANK"))
         {
           if (isBlank && blank != P3D_NOBLANK)
             {
                fprintf(stderr, "ReadP3DFile: %s, double Blank spec.\n", p3dname);
                valid = 0;
                break; /* get out of the while loop */ 
             }
           isBlank = 1;
           blank = P3D_NOBLANK;
           continue;
         }

       /* for Q file */
       if ( !strncmp2(token, "Q=", 2) )
         {
           p3d_set *curset = *sets + *nSets - 1; /* points to the current set */
           int nQ = curset->nQ; /* the number of Q files so far */
           p3d_q **qlistptr = &(curset->qlist); /* points to the qlist */

           /* if there is a filename */
           if (!token[2])
             {
                fprintf(stderr, "ReadP3DFile: %s, missing Q file name.\n", p3dname);
                valid = 0;
                break; /* get out of the while loop */ 
             }
           
           /* alloc a storage for qlist and update the values */
           /* reset the array of qlist */
           if (nQ == 0)
             {
               *qlistptr = (p3d_q *) Alloc(sizeof(p3d_q));
             }
           else
             {
               *qlistptr = (p3d_q *) Realloc(*qlistptr, (nQ + 1) * sizeof(p3d_q));
             }

           /* if memory alloc is OK, then */
           if (!qlistptr)
             {
               fprintf(stderr, "unable to alloc memory\n");
               valid = 0;
               break;
             }

           /* use the old nQ value for array index */
           /* process the Q file spec */
           strcpy((*qlistptr)[nQ].name, token+2);
           (*qlistptr)[nQ].isMDataset = isMDataset;
           (*qlistptr)[nQ].fileType = fileType;
           (*qlistptr)[nQ].permutation = permutation; 
           (*qlistptr)[nQ].check = check;
           (*qlistptr)[nQ].jacobian = jacobian;

           /* increment the number of q file for the current p3d set */
           (curset->nQ)++;

           /* reset the flags exclusive to Q files */
           isFileType = 0;
           isPermutation = 0;
           isCheck = 0;
           isJacobian = 0;

           continue;
         }

       /* for FUNCTION file */
       if (!strncmp2(token, "F=", 2) || !strncmp2(token, "FU=", 3) ||
           !strncmp2(token, "FUN=", 4) || !strncmp2(token, "FUNC=", 5) ||
           !strncmp2(token, "FUNCT=", 6) || !strncmp2(token, "FUNCTI=", 7) ||
           !strncmp2(token, "FUNCTIO=", 8) || !strncmp2(token, "FUNCTION=", 9))
         {
           p3d_set *curset = *sets + *nSets - 1; /* points to the current set */
           int nFunc = curset->nFunc; /* the number of func files so far */
           p3d_func **funclistptr = &(curset->funclist); /* points to the funclist */
           int i = 0;

           /* search for the index of the file name */
           while(token[i++] != '=')
               ;
           /* if there is a filename */
           if (!token[i])
             {
                fprintf(stderr, "ReadP3DFile: %s, missing FUNCTION file name.\n",
                         p3dname);
                valid = 0;
                break; /* get out of the while loop */ 
             }
           
           /* alloc a storage for qlist and update the values */
           /* reset the array of qlist */
           if (nFunc == 0)
             {
               *funclistptr = (p3d_func *) Alloc(sizeof(p3d_func));
             }
           else
             {
               *funclistptr = (p3d_func *) Realloc(*funclistptr,
                                      (nFunc + 1) * sizeof(p3d_func));
             }

           /* if memory alloc is OK, then */
           if (!funclistptr)
             {
               fprintf(stderr, "unable to alloc memory\n");
               valid = 0;
               break;
             }

           /* use the old nFunc value for array index */
           /* process the function file spec */
           strcpy((*funclistptr)[nFunc].name, token+i);
           (*funclistptr)[nFunc].fileType = fileType;
           (*funclistptr)[nFunc].permutation = permutation; 

           /* increment the number of function file for the current p3d set */
           (curset->nFunc)++;

           /* reset default values */
           isFileType = 0;
           isPermutation = 0;

           continue;
         }


       /* for GRID file */
       if (!strncmp2(token, "X=", 2) || !strncmp2(token, "XY=", 3) ||
           !strncmp2(token, "XYZ=", 4))
         {
           p3d_set *curset = *sets + *nSets - 1; /* points to the current set */
           p3d_grid **gridptr = &(curset->grid); /* points to the grid spec */
           int i = 0;

           while(token[i++] != '=')
               ;

           /* if there is a filename */
           if (!token[i])
             {
               fprintf(stderr, "ReadP3DFile: %s, missing XYZ file name.\n", p3dname);
               valid = 0;
               break; /* get out of the while loop */ 
             }

           /* it's an error to specify double XYZ for a single set */
           if (curset->isGrid)
             {
               fprintf(stderr, "ReadP3DFile: %s, double XYZ file spec.\n", p3dname);
               valid = 0;
               break; /* get out of the while loop */ 
             }

           /* alloc memory for grid file spec */
           *gridptr = (p3d_grid *) Alloc(sizeof(p3d_grid));

           /* if memory alloc is OK. then turn on the isGrid */
           if (!gridptr)
             {
               fprintf(stderr, "unable to alloc memory\n");
               valid = 0;
               break;
             }

           curset->isGrid = 1;

           /* process the grid file spec */
           strcpy((*gridptr)->name, token + i);
           (*gridptr)->fileType = fileType;
           (*gridptr)->permutation = permutation; 
           (*gridptr)->blank = blank;

           /* reset default values */
           isFileType = 0;
           isPermutation = 0;
           isBlank = 0;

           continue;
         }

       if (!strcmp2(token, "READ"))
         {
           p3d_set *curset = *sets + *nSets - 1;
           p3d_set *temp;
           int k;

           /* check the set processed so far */
           if (!curset->isGrid)
             {
               fprintf(stderr, "ReadP3DFile: in %s XYZ file missing.\n", p3dname);
               valid = 0;
               break;
             }
           else if(!curset->nQ && !curset->nFunc)
             {
               fprintf(stderr, "ReadP3DFile: %s, no Q or function file.\n", p3dname);
               valid = 0;
               break;
             }

           curset->ndim = ndim;
           curset->isMGrid = isMGrid;

           /* another set begins here */
           if (!(temp = (p3d_set *) Alloc(sizeof(p3d_set) * (*nSets + 1))))
             {
               fprintf(stderr, "Unable to alloc memory, aborted\n");
               valid = 0;
               break; /* get out of the while loop */ 
             }
           for (k = 0; k < *nSets; k++)
             {
               temp[k] = (*sets)[k];
             }
           Free(*sets); 
           *sets = temp;
           (*nSets)++;
           curset = *sets + *nSets - 1;

           /* initialize some values */
           curset->ndim = 3;
           curset->isMGrid = 0;
           curset->isGrid = 0;
           curset->grid = NULL;
           curset->nQ = 0;
           curset->qlist = NULL;
           curset->nFunc = 0;
           curset->funclist = NULL;

           /* reset default values */
           isNdim = 0, ndim = 3;
           isMDataset = 0;
           isMGrid = 0;
           isFileType = 0, fileType = P3D_DAT;
           isPermutation = 0, permutation = P3D_WHOLE;
           isCheck = 0, check = P3D_NOCHECK;
           isJacobian = 0, jacobian = P3D_NOJACOBIAN;
           isBlank = 0, blank = P3D_NOBLANK;

           continue;
         }
       
       fprintf(stderr, "ReadP3DFile: %s, unknown spec %s.\n", p3dname, token);
       valid = 0;
       break; /* get out of the while loop */ 
    } /* while */

   fclose(inFile);

   if (!valid || status == -1)
     {
       fprintf(stderr, "ReadP3DFile: %s, file format error.\n", p3dname);
       return 0;
     }

   /* also, check current set, about nQ, nFunc, and isGrid, to see
      if there is any error */
   if ((*sets + *nSets - 1)->isGrid == 0)
     {
       fprintf(stderr, "ReadP3DFile: %s, XYZ missing in number %d read command.\n", p3dname, *nSets);
       return 0;
     }
   else if((*sets + *nSets - 1)->nQ == 0 && (*sets + *nSets - 1)->nFunc == 0)
     {
       fprintf(stderr, "ReadP3DFile: %s, no Q or function file.\n", p3dname);
       return 0;
     }
   else
     {
       (*sets + *nSets - 1)->ndim = ndim;
       (*sets + *nSets - 1)->isMGrid = isMGrid;
     }

   return 1;
} /* parseP3DMeta */

/* using row major order arrangement */
#ifdef HAVE_PROTOTYPES
static void putIblank(int iblank, int **ipp, long *indices, long *dims,
        int ndim, int igrid) 
#else
static void putIblank(iblank, ipp, indices, dims, ndim, igrid) 
int iblank;
int **ipp;
long *indices;
long *dims;
int ndim;  /* rank */
int igrid;
#endif
/* return 0 when failed */
{
  int i;
  long offset;

  offset = indices[ndim - 1];
  for (i = ndim - 1; i > 0; i--)
     {
       offset = offset * dims[i - 1] + indices[i - 1];
     }
  ipp[igrid][offset] = iblank;

  return;
} /* putIblank */

#ifdef HAVE_PROTOTYPES
static void getIblank(int *iblank, int **ipp, long *indices,
                                  long *dims, int ndim, int igrid) 
#else
static void getIblank(iblank, ipp, indices, dims, ndim, igrid)
int *iblank;
int **ipp;
long *indices;
long *dims;
int ndim;
int igrid;
#endif
{
  int i;
  long offset;

  offset = indices[ndim - 1];
  for (i = ndim - 1; i > 0; i--)
     {
       offset = offset * dims[i - 1] + indices[i - 1];
     }
  *iblank = ipp[igrid][offset];

  return;
} /* getIblank */

#ifdef HAVE_PROTOTYPES
static void freeiblank(int ngrid, int **ipp)
#else
static void freeiblank(ngrid, ipp)
int ngrid;
int **ipp;
#endif
{
  int i;

  for (i = 0; i < ngrid; i++)
     Free(ipp[i]);
  Free(ipp);

  return;
}

#ifdef HAVE_PROTOTYPES
static void freedims(int ngrid, long **dimsptr)
#else
static void freedims(ngrid, dimsptr)
int ngrid;
long **dimsptr;
#endif
{
  int i;

  for (i = 0; i < ngrid; i++)
     Free(dimsptr[i]);
  Free(dimsptr);

  return;
}

/*********************formatted p3d reading routines, start *********************/

#ifdef HAVE_PROTOTYPES
static int RdFmtXYZ(int ndim, int isMGrid, long ***dimsptr,
                               p3d_grid *gridptr, ObjPtr **formptr, int ***ipp)
#else
static int RdFmtXYZ(ndim, isMGrid, dimsptr, gridptr, formptr, ipp)
int ndim;
int isMGrid;
long ***dimsptr;
p3d_grid *gridptr;
ObjPtr **formptr;
int ***ipp;
#endif
{
  FILE *inFile;
  int ngrid;     /* number of grids */
  long *indices;
  real *bounds;
  int i, j;
  int n;
  char *name;  /* dataset name starting point */
  int namelen;
  int repeat;
  double dval;

  if (!(inFile = fopen(gridptr->name, "r")))
    {
      fprintf(stderr, "Unable to open %s, aborted\n", gridptr->name);
      return 0;
    }

  /* determine the dataset name now */
  for (i = (int) strlen(gridptr->name) - 1; i >= 0 && gridptr->name[i] != '/'; i--)
     ;
  if (i == -1 || gridptr->name[i] == '/')
    i++;
  name = gridptr->name + i;
  for (namelen = 0; name[namelen] != '\0' && name[namelen] != '.'; namelen++)
     ;

  if (isMGrid)
    {
      SkipBlanksAndCommas(inFile);
      if (fscanf(inFile, "%d", &ngrid) != 1 || ngrid <= 0)
        {
          fprintf(stderr, "%s format error, invalid NGRID value\n", gridptr->name);
          fclose(inFile);
          return 0;
        } 
    }
  else
    {
      ngrid = 1;
    }
  
  /* alloc dimsptr and indices storage */
  if (!(indices = (long *) Alloc(ndim * sizeof(long))))
    {
      fprintf(stderr, "Unable to alloc memory\n");
      fclose(inFile);
      return 0;
    }
  if (!(*dimsptr = (long **) Alloc(ngrid * sizeof(long *))))
    {
      fprintf(stderr, "Unable to alloc memory\n");
      Free(indices);
      fclose(inFile);
      return 0;
    }
  repeat = 0;
  for (i = 0; i < ngrid; i++)
    {
      if (!((*dimsptr)[i] = (long *) Alloc(ndim * sizeof(long))))
        {
          fprintf(stderr, "Unable to alloc memory\n");
          Free(indices);
          freedims(i, *dimsptr);
          fclose(inFile);
          return 0;
        }

      for (j = 0; j < ndim; j++)
         {
           if (!repeat)
             {
                if ((repeat = getValue(inFile, &dval)) <= 0 || (int) dval <= 0)
                  {
                    fprintf(stderr, "Invalid dimension value in %s\n",
                            gridptr->name);
                    Free(indices);
                    freedims(i, *dimsptr);
                    fclose(inFile);
                    return 0;
                  }
             }
           repeat--;
           *((*dimsptr)[i] + j) = (long) dval; 
         }
    }

  /* alloc strorage for pointer to bounds */
  if (!(bounds = (real *) Alloc(ndim * 2 * sizeof(real))))
    {
      fprintf(stderr, "Unable to alloc memory\n");
      Free(indices);
      freedims(ngrid, *dimsptr);
      fclose(inFile);
      return 0;
    }
  /* alloc storage for pointers to dataforms */
  if (!(*formptr = (ObjPtr *) Alloc(ngrid * sizeof(ObjPtr))))
    {
      fprintf(stderr, "Unable to alloc memory\n");
      Free(indices);
      Free(bounds);
      freedims(ngrid, *dimsptr);
      fclose(inFile);
      return 0;
    }

  /* alloc iblankptr if IBLANK */
  if (gridptr->blank == P3D_BLANK)
    {
      if (!(*ipp = (int **) Alloc(ngrid * sizeof(int *))))
        {
          fprintf(stderr, "Unable to alloc memory\n");
          Free(indices);
          Free(bounds);
          Free(*formptr);
          freedims(ngrid, *dimsptr);
          fclose(inFile);
          return 0;
        }

    }

  repeat = 0;
  for (n = 0; n < ngrid; n++)
    {
      long size, count; /* number of grid values promised */
      ObjPtr dataset;
      char datasetname[P3D_MAXNAMELEN];
      int whichcomponent;

      strncpy(datasetname, name, namelen);
      datasetname[namelen] = '\0';
      if (ngrid > 1)
        {
          sprintf(datasetname + namelen, "%d", n + 1);
        }

      /* setup a dataset for the grid data */
      dataset = NewStructuredDataset(datasetname, ndim, (*dimsptr)[n], ndim);
      if (!SetCurField(FIELD1, dataset))
        {
          fprintf(stderr, "Unable to SetCurField\n");
          Free(indices);
          Free(bounds);
          Free(*formptr);
          freedims(ngrid, *dimsptr);
          if (gridptr->blank == P3D_BLANK)
            {
              freeiblank(n-1, *ipp);
            }
          fclose(inFile);
          return 0;
        }

      size = 1;
      for (i = 0; i < ndim; i++)
         {
           size *= (*dimsptr)[n][i];
         }
      size *= (long)(ndim + (gridptr->blank == P3D_BLANK));

      if (gridptr->blank == P3D_BLANK)
        {
          /* alloc storage for the iblank arrays */
          if (!((*ipp)[n] = (int *) Alloc(size * sizeof(int))))
            {
              fprintf(stderr, "Unable to alloc memory\n");
              freeiblank(n, *ipp);
              Free(indices);
              Free(bounds);
              Free(*formptr);
              freedims(ngrid, *dimsptr);
              fclose(inFile);
              return 0;
            }
        }

      for (i = 0; i < ndim; i++)
         {
           indices[i] = 0;
         }
      whichcomponent = 0;
      /* read in all the grid values */
      for (count = 0; count < size; count++)
        {
          int first = 1;

          if (!repeat)
            {
              if ((repeat = getValue(inFile, &dval)) <= 0)
                {
                  fprintf(stderr, "Error reading %s, invalid numeric value\n",
                           gridptr->name);
                  Free(indices);
                  Free(bounds);
                  Free(*formptr);
                  freedims(ngrid, *dimsptr);
                  if (gridptr->blank == P3D_BLANK)
                    {
                      freeiblank(n + 1, *ipp);
                    }
                  fclose(inFile);
                  return 0;
                }
            }

          repeat--;
          if (whichcomponent < ndim) /* read grid values */
            {
              PutFieldComponent(FIELD1, whichcomponent, indices, (real) dval);  
              for (j = 0; first && j < ndim; j++)
                 first = indices[j] == 0;

              if (first)
                {
                  bounds[whichcomponent * 2] = (real) dval;
                  bounds[whichcomponent * 2 + 1] = (real) dval;
                }
              else
                {
                  if (dval < bounds[whichcomponent * 2])
                    bounds[whichcomponent * 2] = (real) dval;
                  else if (dval > bounds[whichcomponent * 2 + 1])
                    bounds[whichcomponent * 2 + 1] = (real) dval;
                }
            }
          else /* read iblank values */
            {
              putIblank((int) dval, *ipp, indices, (*dimsptr)[n], ndim, n);
            }
          /* adjust the indices and whichcomponent */
          if (gridptr->permutation == P3D_WHOLE)
            {
              for (j = 0; j < ndim; j++)
                 {
                   indices[j]++;
                   if (indices[j] == (*dimsptr)[n][j])
                     {
                       indices[j] = 0;
                       continue;
                     }
                   break;
                 }
              if (j == ndim) /* all indices are set to zero now */ 
                {
                  whichcomponent++;
                }
            }
          else  /* P3D_PLANE */
            {
              for (j = 0; j < ndim - 1; j++)
                 {
                   indices[j]++;
                   if (indices[j] == (*dimsptr)[n][j])
                     {
                       indices[j] = 0;
                       continue;
                     }
                   break;
                 }
              if (j == ndim - 1)  
                {
                  whichcomponent++;
                  if ((gridptr->blank == P3D_BLANK) ? 
                     (whichcomponent == ndim + 1) : (whichcomponent == ndim))
                  /* leave a space for reading iblanks */
                    {
                      whichcomponent = 0;
                      indices[j]++;
                    }
                }
            }
        } /* for read grid values */

     (*formptr)[n] = NewCurvilinearDataForm(datasetname, ndim, (*dimsptr)[n],
                    bounds, dataset);

    } /* outer for */

  /* indices and bounds are no longer of any use */
  Free(indices);
  Free(bounds);
  fclose(inFile);
  return n;
} /* RdFmtXYZ */

#ifdef HAVE_PROTOTYPES
static int RdFmtSol(ObjPtr reader, int ndim, int isMGrid, int ngrid, int blank,
                 long **dimsptr, p3d_q *qptr, ObjPtr *formptr, int **ipp)
#else
static int RdFmtSol(reader, ndim, isMGrid, ngrid, blank, dimsptr, qptr, formptr, ipp)
ObjPtr reader;
int ndim;
int isMGrid;
int ngrid;
int blank;
long **dimsptr;
p3d_q *qptr;
ObjPtr *formptr;
int **ipp;
#endif
{
  FILE *inFile;
  int i, j, n;
  char *name;
  long *indices;
  int numsets; /* number of datasets read so far, for multiple datasets in a file*/
  int namelen; /* the length of dataset name going to be used */
  int repeat;  /* to deal with <num>*<value> format */
  double dval;

  if (!(inFile = fopen(qptr->name, "r")))
    {
      fprintf(stderr, "Unable to open %s, aborted\n", qptr->name);
      return 0;
    }

  if (!(indices = (long *) Alloc(ndim * sizeof(long))))
    {
      fprintf(stderr, "Unable to alloc memory, aborted\n");
      fclose(inFile);
      return 0;
    }

  /* determine the dataset name from the file name */
  for (i = (int) strlen(qptr->name) - 1; i >= 0 && qptr->name[i] != '/'; i--)
     ;
  if (i == -1 || qptr->name[i] == '/')
    i++;
  name = qptr->name + i;
  for (namelen = 0; name[namelen] != '\0' && name[namelen] != '.'; namelen++)
     ;

  numsets = 0;

  do
    {
      repeat = 0;
      /* get ngrid value if it is a MGRID q file */
      if (isMGrid)
        {
           repeat = getValue(inFile, &dval);
           if (repeat == EOF)
             {
              if (numsets == 0)
                 {
                  fprintf(stderr, "Error in %s: empty file\n", qptr->name);
                 }
              fclose(inFile);
              Free(indices);
              return (numsets ? ngrid : 0);
             }
          else if (repeat == 0 || (n = (int) dval) <= 0)
             {
              fprintf(stderr, "Error in %s: Invalid NGRID spec.\n", qptr->name);
              fclose(inFile);
              Free(indices);
              return 0;
             }
          repeat--;
        }
      else
        {
          n = 1;
        }

      /* if ngrid != this file's ngrid then error */
      if (ngrid != n)
        {
          fprintf(stderr, "Error in %s: NGRID value does not match its grid\n",
                                                             qptr->name);
          fclose(inFile);
          Free(indices);
          return 0;
          
        }
      /* if dims does not match the dims read in XYZ file then error */
      for (n = 0; n < ngrid; n++)
        {
          for (i = 0; i < ndim; i++)
            {
              if (!repeat)
                {
                  if ((repeat = getValue(inFile, &dval)) <= 0)
                    {
                      /* if this is not multiple grid, then the first
                      dimension value is used to indicate whether there is
                      more datasets or not */
                      if (!isMGrid && i == 0 && repeat == EOF)
                        {
                         if (numsets == 0)
                           {
                            fprintf(stderr, "Error in %s: empty file\n", qptr->name);
                           }
                          fclose(inFile);
                          Free(indices);
                          return (numsets ? ngrid : 0);
                        }
                      fprintf(stderr, "Invalid dimension value in %s\n", qptr->name);
                      fclose(inFile);
                      Free(indices);
                      return 0;
                    }
                }
              repeat--;
              if ((long) dval != dimsptr[n][i])
                {
                  fprintf(stderr,
                  "Dimensions mismatch between %s and its XYZ file\n", qptr->name);
                  fclose(inFile);
                  Free(indices);
                  return 0;
                }
            }
        }
    
      for (n = 0; n < ngrid; n++)
        {
          long size, count; /* number of data values promised */
          ObjPtr density, energy, momentum;
          char momentumName[P3D_MAXNAMELEN];
          char densityName[P3D_MAXNAMELEN];
          char energyName[P3D_MAXNAMELEN];
          char nameBase[P3D_MAXNAMELEN];
          double time;
          int whichcomponent;
          ObjPtr whichdataset;
    
          /* process the FAMACH..TIME line */
          for (i = 1; i <= 4; i++)
            {
              if (!repeat)
                {
                  if ((repeat = getValue(inFile, &dval)) <= 0)
                    {
                     fprintf(stderr,
                     "Error in %s: invalid FSMACH,ALPHA,RE,TIME line\n", qptr->name);
                     Free(indices);
                     fclose(inFile);
                     return 0;
                    }
                }
              repeat--;
              if (i == 4)
                {
                  time = dval;
                }
            }
    
	  strncpy(nameBase, name, namelen);
	  nameBase[namelen] = '\0';
          /* setup datasets */
          if (ngrid > 1)
            {
              sprintf(densityName, "%s %s %d", nameBase, "Density", n + 1);
              sprintf(energyName, "%s %s %d", nameBase, "Stagnation Energy per Unit Volume", n + 1);
              sprintf(momentumName, "%s %s %d", nameBase, "Momentum", n + 1);
            }
          else
	    {
              sprintf(densityName, "%s %s", nameBase, "Density");
              sprintf(energyName, "%s %s", nameBase, "Stagnation Energy per Unit Volume");
              sprintf(momentumName, "%s %s", nameBase, "Momentum");
	    }
          density = NewStructuredDataset(densityName, ndim, dimsptr[n], 0);
          energy = NewStructuredDataset(energyName, ndim, dimsptr[n], 0);
          momentum = NewStructuredDataset(momentumName, ndim, dimsptr[n], ndim);

          size = 1;
          for (i = 0; i < ndim; i++)
             size *= dimsptr[n][i];
          size *= ndim + 2;  /* includes density and pressure */
          for (i = 0; i < ndim; i++)
             {
               indices[i] = 0;
             }
          whichcomponent = 0;
          whichdataset = density;
    
          /* read in all the q values */
          for (count = 0; count < size; count++)
            {
              int iblank = 1;
    
              if (!repeat)
                {
                  if ((repeat = getValue(inFile, &dval)) <= 0)
                    {
                      fprintf(stderr, "Error reading %s, invalid numeric value\n",
                                  qptr->name);
                      Free(indices);
                      fclose(inFile);
                      return 0;
                    }
                }
              repeat--;
              /* do check of density and pressure, and iblanks here */
              if (qptr->check == P3D_CHECK && (whichdataset == energy ||
                  whichdataset == density))
                {
                  if (dval < 0.0)
		    {
		      fprintf(stderr, "Negative %s value %f found in Q file %s and is converted to %f\n", whichdataset == energy ? "pressure" : "density", dval, qptr->name, 0.0);
                      dval = 0.0;
                    }
                }
              if (blank == P3D_BLANK)
                {
                   getIblank(&iblank, ipp, indices, dimsptr[n], ndim, n);
                }
    
              if (!SetCurField(FIELD1, whichdataset))
                {
                  fprintf(stderr, "Unable to SetCurField\n");
                  Free(indices);
                  fclose(inFile);
                  return 0;
                }
    
              PutFieldComponent(FIELD1, whichcomponent, indices, 
                          iblank == 0 ? (real) missingData : (real) dval);  
    
              /* adjust the indices, whichdataset, and whichcomponent */
              if (qptr->permutation == P3D_WHOLE)
                {
                  for (j = 0; j < ndim; j++)
                     {
                       indices[j]++;
                       if (indices[j] == dimsptr[n][j])
                         {
                           indices[j] = 0;
                           continue;
                         }
                       break;
                     }
                  if (j == ndim) /* all indices are set to zero now */ 
                    {
		      if (whichdataset == density)
			{
			  whichdataset = momentum;
                        }
                      else if (++whichcomponent == ndim)
			{
			  whichdataset = energy;
			  whichcomponent = 0;
                        }
                    }
                }
              else  /* P3D_PLANE */
                {
                  for (j = 0; j < ndim - 1; j++)
                     {
                       indices[j]++;
                       if (indices[j] == dimsptr[n][j])
                         {
                           indices[j] = 0;
                           continue;
                         }
                       break;
                     }
                  if (j == ndim - 1)  
                    {
                      if (whichdataset == density)
                        whichdataset = momentum;
                      else if (whichdataset = momentum)
                        {
                          if (++whichcomponent == ndim)
                            {
                              whichdataset = energy;
                              whichcomponent = 0;
                            }
                        }
                      else /* reading energy */
			{
			  whichdataset = density;
                          indices[j]++;
			}
                    }
                }
            } 
    
          SetDatasetForm(density, formptr[n]);
          SetDatasetForm(energy, formptr[n]);
          SetDatasetForm(momentum, formptr[n]);
          if (qptr->isMDataset)
            {
              RegisterNewTimeDataset(density, (real) time);
              RegisterNewTimeDataset(momentum, (real) time);
              RegisterNewTimeDataset(energy, (real) time);
#if 0
              RegisterTimeDataset(density, (real) time);
              RegisterTimeDataset(momentum, (real) time);
              RegisterTimeDataset(energy, (real) time);
#endif
            }
          else
            {
              RegisterDataset(density);
              RegisterDataset(momentum);
              RegisterDataset(energy);
            }
        } /* outer for */
      numsets++;
    } while(qptr->isMDataset);

  Free(indices);
  fclose(inFile);
  return ngrid;
} /* RdFmtSol */

#ifdef HAVE_PROTOTYPES
static int RdFmtFun(int ndim, int isMGrid, int ngrid, int blank,
                long **dimsptr, p3d_func *funcptr, ObjPtr *formptr, int **ipp)
#else
static int RdFmtFun(ndim, isMGrid, ngrid, blank, 
                              dimsptr, funcptr, formptr, ipp)
int ndim;
int isMGrid;
int ngrid;
int blank;
long **dimsptr;
p3d_func *funcptr;
ObjPtr *formptr;
int **ipp;
#endif
{
  FILE *inFile;/* the current function file being processed */
  int i, j, n;
  char *name; /* points to the portion of function file name after the last / */
  long *indices;
  int *nvar;  /* number of variables of this function */
  int namelen;
  int repeat;
  double dval;

  if (!(inFile = fopen(funcptr->name, "r")))
    {
      fprintf(stderr, "Unable to open %s, aborted\n", funcptr->name);
      return 0;
    }
  if (isMGrid)
    {
      SkipBlanksAndCommas(inFile);
      if (fscanf(inFile, "%d", &n) != 1 || n <= 0)
        {
          fprintf(stderr, "Error in %s: invalid NGRID value\n", funcptr->name);
          fclose(inFile);
          return 0;
        }
    }
  else
    {
      n = 1;
    }
  /* if ngrid != this file's ngrid then error */
  if (ngrid != n)
    {
      fprintf(stderr, "Error in %s: NGRID value does not match its grid\n",
                                                                 funcptr->name);
      fclose(inFile);
      return 0;
    }
  if (!(indices = (long *) Alloc(ndim * sizeof(long))))
    {
      fprintf(stderr, "Unable to alloc memory, aborted\n");
      fclose(inFile);
      return 0;
    }
  if (!(nvar = (int *) Alloc(ngrid * sizeof(int))))
    {
      fprintf(stderr, "Unable to alloc memory, aborted\n");
      Free(indices);
      fclose(inFile);
      return 0;
    }

  /* if dims does not match the dims read in XYZ file then error */
  repeat = 0;
  for (n = 0; n < ngrid; n++)
    {
      for (i = 0; i < ndim + 1; i++) /* nvar follows dimension spec's */
        {
          if (!repeat)
            {
              if ((repeat = getValue(inFile, &dval)) <= 0)
                {
                  fprintf(stderr, "Invalid dim/nvar value in %s\n", funcptr->name);
                  Free(indices);
                  Free(nvar);
                  fclose(inFile);
                  return 0;
                }
            }
          repeat--;
          if (i < ndim)  /* the value is dimension */
            {
              if ((long) dval != dimsptr[n][i])
                {
                  fprintf(stderr, "Dimensions mismatch between %s"
                                      " and its XYZ file\n", funcptr->name);
                  Free(indices);
                  Free(nvar);
                  fclose(inFile);
                  return 0;
                }
            }
          else /* i = ndim, the value is of nvar */
            {
              nvar[n] = (int) dval;
            }
        }
    }

  /* determine the dataset name */
  for (i = (int) strlen(funcptr->name) - 1; i >= 0 && funcptr->name[i] != '/'; i--)
     ;
  if (i == -1 || funcptr->name[i] == '/')
    i++;
  name = funcptr->name + i;
  for (namelen = 0; name[namelen] != '\0' && name[namelen] != '.'; namelen++)
     ;

  repeat = 0;
  for (n = 0; n < ngrid; n++)
    {
      long size, count; /* number of data values promised */
      ObjPtr dataset;
      char datasetname[P3D_MAXNAMELEN];
      int whichcomponent;

      strncpy(datasetname, name, namelen);
      datasetname[namelen] = '\0';
      /* number the datasets only if its number of grids is than 1 */
      if (ngrid > 1)
        {
	  /* ordering starts from 1 */
          sprintf(datasetname + namelen, "%d", n + 1);
        }

      size = 1;
      for (i = 0; i < ndim; i++)
         size *= dimsptr[n][i];
      size *= nvar[n];
      for (i = 0; i < ndim; i++)
         {
           indices[i] = 0;
         }

      /* setup datasets */
      dataset = NewStructuredDataset(datasetname, ndim, dimsptr[n],
                                              nvar[n] == 1 ? 0 : nvar[n]);
      if (!SetCurField(FIELD1, dataset))
        {
          fprintf(stderr, "Unable to SetCurField\n");
          Free(indices);
          Free(nvar);
          fclose(inFile);
          return 0;
        }
      whichcomponent = 0;

      /* read in all the func values */
      for (count = 0; count < size; count++)
        {
          int iblank = 1;
     
          if (!repeat)
            {
              if ((repeat = getValue(inFile, &dval)) <= 0)
                {
                  fprintf(stderr, "Error reading %s, invalid numeric value\n",
                                funcptr->name);
                  Free(indices);
                  Free(nvar);
                  fclose(inFile);
                  return 0;
                }
            }
          repeat--;
          if (blank == P3D_BLANK)
            {
               getIblank(&iblank, ipp, indices, dimsptr[n], ndim, n);
            }
    
          PutFieldComponent(FIELD1, whichcomponent, indices, 
                  iblank == 0 ? (real) missingData : (real) dval);  

          /* adjust the indices and whichcomponent */
          if (funcptr->permutation == P3D_WHOLE)
            {
              for (j = 0; j < ndim; j++)
                {
                  indices[j]++;
                  if (indices[j] == dimsptr[n][j])
                    {
                      indices[j] = 0;
                      continue;
                    }
                  break;
                }
              if (j == ndim) /* all indices are set to zero now */
                {
                  whichcomponent++;
                }
            }
          else  /* P3D_PLANE */
            {
              for (j = 0; j < ndim - 1; j++)
                {
                  indices[j]++;
                  if (indices[j] == dimsptr[n][j])
                    {
                      indices[j] = 0;
                      continue;
                    }
                  break;
                }
              if (j == ndim - 1 && ++whichcomponent == nvar[n])
                {
                  indices[j]++;
                  whichcomponent = 0;
                }
            }
        }     /* for size */

      SetDatasetForm(dataset, formptr[n]);
      RegisterDataset(dataset);
    } /* outer for */

  Free(indices);
  Free(nvar);
  fclose(inFile);
  return ngrid;
} /* RdFmtFun */

/*********************formatted p3d reading routines, end *********************/

/*********************unformatted p3d reading routines, start *********************/
/* only available when there's fortran compiler */

#ifdef FORTRAN

#ifdef HAVE_PROTOTYPES
static int RdUnfmtXYZ(int ndim, int isMGrid, long ***dimsptr,
                      p3d_grid *gridptr, ObjPtr **formptr, int ***ipp);
static int RdUnfmtSol(ObjPtr reader, int ndim, int isMGrid, int ngrid, int blank, 
		      long **dimsptr, p3d_q *qptr, ObjPtr *formptr, int **ipp);
static int RdUnfmtFun(int ndim, int isMGrid, int ngrid, int blank, long **dimsptr,
		      p3d_func *funcptr, ObjPtr *formptr, int **ipp);
#ifdef FORTRAN_
/* fortran subroutines */
extern void foropn_(int *lun, int intname[], int *length, int *opstat);
extern void forcls_(int *lun);
extern void rngrid_(int *lun, int *ngrid, int *opstat);
extern void rddims_(int *lun, int *ndim, int *ngrid, int *indims, int *opstat);
extern void rdgrid_(int *lun, int *ndim, int *indims, int *isiblk, int *iblank,
		    int *iperm, float *gdvals, long *size, int *opstat); 
extern void rdsolu_(int *lun, int *ndim, int *indims, int *iperm, 
		     float *slvals, long *size, int *opstat);
extern void rdtime_(int *lun, float *time, int *opstat);
extern void rdfunc_(int *lun, int *ndim, int *indims, int *nvar, int *iperm,
		    float *fnvals, long *size, int *opstat); 
#else
/* fortran subroutines */
extern void foropn(int *lun, int iname[], int *length, int *opstat);
extern void forcls(int *lun);
extern void rngrid(int *lun, int *ngrid, int *opstat);
extern void rddims(int *lun, int *ndim, int *ngrid, int *indims, int *opstat);
extern void rdgrid(int *lun, int *ndim, int *indims, int *isiblk, int *iblank,
		   int *iperm, float *gdvals, long *size, int *opstat); 
extern void rdsolu(int *lun, int *ndim, int *indims, int *iperm, 
		    float *slvals, long *size, int *opstat);
extern void rdtime(int *lun, float *time, int *opstat);
extern void rdfunc(int *lun, int *ndim, int *indims, int *nvar, int *iperm,
		   float *fnvals, long *size, int *opstat); 
#endif

#else

static int RdUnfmtXYZ();
static int RdUnfmtSol();
static int RdUnfmtFun();
#ifdef FORTRAN_
/* fortran subroutines */
extern void foropn_();
extern void forcls_();
extern void rngrid_();
extern void rddims_();
extern void rdgrid_();
extern void rdsolu_();
extern void rdtime_();
extern void rdfunc_();
#else
/* fortran subroutines */
extern void foropn();
extern void forcls();
extern void rngrid();
extern void rddims();
extern void rdgrid();
extern void rdsolu();
extern void rdtime();
extern void rdfunc();
#endif
#endif


/* for reading fortran unformatted grid file */
#ifdef HAVE_PROTOTYPES
static int RdUnfmtXYZ(int ndim, int isMGrid, long ***dimsptr,
                               p3d_grid *gridptr, ObjPtr **formptr, int ***ipp)
#else
static int RdUnfmtXYZ(ndim, isMGrid, dimsptr, gridptr, formptr, ipp)
int ndim;
int isMGrid;
long ***dimsptr;
p3d_grid *gridptr;
ObjPtr **formptr;
int ***ipp;
#endif
{
  int ins[P3D_MAXNAMELEN]; /* ascii coding of file name */
  int namelen;         /* file name length */
  int lun = 2;         /* logical unit number of fortran opened file */
  int opstat=0;        /* operation status of fortran file op subroutine */
  int *indims, ngrid;
  int i, j;
  
/* try to open the xyz file, ascii coding of the xyz file name */
  namelen = (int) strlen(gridptr->name);
  for (i = 0; i <= namelen; i++)
     {
       ins[i] = gridptr->name[i];
     }
/* open file with fortran subroutine */
#ifdef FORTRAN_
  foropn_(&lun, ins, &namelen, &opstat);
#else
  foropn(&lun, ins, &namelen, &opstat);
#endif

  if (!opstat)
    {
      fprintf(stderr, "Unable to open file %s, aborted\n", gridptr->name);
      return 0;
    }

/* reading ngrid */
  if (isMGrid)
    {
      opstat = 0;
#ifdef FORTRAN_
      rngrid_(&lun, &ngrid, &opstat);
#else
      rngrid(&lun, &ngrid, &opstat);
#endif
      if (opstat <= 0 || ngrid <=0)
        {
	  fprintf(stderr, "%s format error, invalid NGRID value\n",
	  gridptr->name);
#ifdef FORTRAN_
          forcls_(&lun);
#else
          forcls(&lun);
#endif
	  return 0;
        }
    }
  else
    {
      ngrid = 1;
    }

/* reading dimensions */
  indims = (int *) 0;
  *dimsptr = (long **) 0;
  if (!(indims = (int *) Alloc(ndim*ngrid*sizeof(int))) ||
      !(*dimsptr = (long **) Alloc(ndim * sizeof(long *))))
    {
      fprintf(stderr, "Unable to alloc mem, aborted\n");
#ifdef FORTRAN_
          forcls_(&lun);
#else
          forcls(&lun);
#endif
      Free(indims);
      return 0;
    }
  opstat=0;
#ifdef FORTRAN_
  rddims_(&lun, &ndim, &ngrid, indims, &opstat);
#else
  rddims(&lun, &ndim, &ngrid, indims, &opstat);
#endif
  if (opstat > 0)
    {
      for (i = 0; i < ngrid; i++)
	 {
           if (!((*dimsptr)[i] = (long *) Alloc(ndim * sizeof(long))))
	     {
	       fprintf(stderr, "Unable to alloc memory, aborted\n");
	       i--;
	       opstat = 0;  /* turn off opstat to indicate failure */
	       break;
             }
	   for (j = 0; j < ndim; j++)
             {
	       *((*dimsptr)[i]+j)= (long) indims[i*ndim + j];
             }
         }
    }
  if (opstat <= 0)
    {
      fprintf(stderr, "Unable to Read Dimension Values in %s\n",
	      gridptr->name);
#ifdef FORTRAN_
          forcls_(&lun);
#else
          forcls(&lun);
#endif
      freedims(i, *dimsptr);
      Free(indims);
      return 0;
    }

/* reading grid values*/
    {
      long *indices;
      real *bounds;
      char *name;   /* dataset name, extracted from file name */

      /* extracting the dataset name from file name */
      for(i = namelen -1;
	  i >= 0 && gridptr->name[i] != '/';
	  i--) ;
      if (i == -1 || gridptr->name[i] == '/')
	i++;
      name = gridptr->name + i;
      for (namelen = 0; 
	   name[namelen] != '\0' && name[namelen] != '.';
	   namelen++) ;
         
      /* memory allocation for indices, bounds, formptr, ipp  */
      indices = (long *) 0;
      bounds = (real *) 0;
      *formptr = (ObjPtr *) 0;
      *ipp = (int **) 0;
      if ( !(indices = (long *) Alloc(ndim * sizeof(long))) ||
	   !(bounds = (real *) Alloc(ndim * 2 * sizeof(real))) ||
	   !(*formptr = (ObjPtr *) Alloc(ngrid * sizeof(ObjPtr))) ||
	   !(*ipp = (int **) Alloc(ngrid *sizeof(int *))) )
           /* alloc *ipp regardless of the value of gridptr->blank, but
	      remember to free it if it is of no use */
        {
          fprintf(stderr, "Unable to alloc mem, aborted\n");
#ifdef FORTRAN_
          forcls_(&lun);
#else
          forcls(&lun);
#endif
          Free(indims);
	  freedims(ngrid, *dimsptr);
	  Free(*formptr);
	  Free(*ipp);
	  Free(indices);
	  Free(bounds);
          return 0;
	}

      for(i = 0; i < ngrid; i++)
	{
          long size, count;
          float *gdvals;
	  ObjPtr dataset;
	  char datasetname[P3D_MAXNAMELEN];

	  /* create a dataset for the grid to be read in */
          strncpy(datasetname, name, namelen);
	  datasetname[namelen] = '\0';
	  if (ngrid > 1)
	    sprintf(datasetname + namelen, "%d", i + 1);
          dataset=NewStructuredDataset(datasetname, ndim, (*dimsptr)[i], ndim);
	  if (!SetCurField(FIELD1, dataset))
	    {
                 fprintf(stderr, "Unable to SETCURFIELD, aborted\n");
#ifdef FORTRAN_
                 forcls_(&lun);
#else
                 forcls(&lun);
#endif
                 Free(indims);
	         freedims(ngrid, *dimsptr);
		 if (gridptr->blank == P3D_BLANK)
		   {
	             freeiblank(i-1, *ipp);
		   }
                 else
		   {
		     Free(*ipp);
		   }
	         Free(*formptr);
	         Free(indices);
	         Free(bounds);
                 return 0;
	    }

          /* calculate size of one component */
          size=1;
          for(j = 0; j < ndim; j++)
           {
             size *= (*dimsptr)[i][j];
           }
	  (*ipp)[i] = (int *) 0;
          if (gridptr->blank == P3D_BLANK)
	    {
	      if (!((*ipp)[i] = (int *) Alloc(size*sizeof(int))))
		{
                 fprintf(stderr, "Unable to alloc mem, aborted\n");
#ifdef FORTRAN_
                 forcls_(&lun);
#else
                 forcls(&lun);
#endif
                 Free(indims);
	         freedims(ngrid, *dimsptr);
	         freeiblank(i, *ipp);
	         Free(indices);
	         Free(bounds);
	         Free(*formptr);
                 return 0;
                }
	    }

	  /* read in grid values */
          if (!(gdvals = (float *) Alloc(ndim*size*sizeof(float))))
	    {
                 fprintf(stderr, "Unable to alloc mem, aborted\n");
#ifdef FORTRAN_
                 forcls_(&lun);
#else
                 forcls(&lun);
#endif
                 Free(indims);
	         freedims(ngrid, *dimsptr);
		 if (gridptr->blank == P3D_BLANK)
		   {
	             freeiblank(i, *ipp);
		   }
                 else
		   {
		     Free(*ipp);
		   }
	         Free(indices);
	         Free(bounds);
	         Free(*formptr);
                 return 0;
	    }
	  opstat=0;
#ifdef FORTRAN_
          rdgrid_(&lun, &ndim, indims+i*ndim, &gridptr->blank, (*ipp)[i],
		  &(gridptr->permutation), gdvals, &size, &opstat);
#else
          rdgrid(&lun, &ndim, indims+i*ndim, &gridptr->blank, (*ipp)[i],
		  &(gridptr->permutation), gdvals, &size, &opstat);
#endif
          if (!opstat)
	    {
	      fprintf(stderr, "Read error on %d grid of file %s\n", 
		       i, gridptr->name);
#ifdef FORTRAN_
              forcls_(&lun);
#else
              forcls(&lun);
#endif
              Free(indims);
	      freedims(ngrid, *dimsptr);
	      if (gridptr->blank == P3D_BLANK)
	        {
	          freeiblank(i, *ipp);
		}
              else
	        {
		  Free(*ipp);
		}
	      Free(indices);
	      Free(bounds);
	      Free(*formptr);
	      Free(gdvals);
              return 0;
	    }


	  /* assign gdvals to dataset and calculate bounds here */
          for(j = 0; j < ndim; j++)
	    {
	      indices[j] = 0;
	      /* give bounds initial values */
	      bounds[j*2] = bounds[j*2 + 1] = gdvals[j*size];
            }
          for(count = 0; count < size; count++)
	    {
	      int k, whichcomponent;

	      /* assign grid values */
	      for(whichcomponent = 0; whichcomponent < ndim; whichcomponent++)
		{
		  real rval;

		  rval = (real) gdvals[whichcomponent*size + count];
		  PutFieldComponent(FIELD1, whichcomponent, indices, rval);
                  if (rval < bounds[whichcomponent*2])
		    bounds[whichcomponent*2] = rval;
                  else if (rval > bounds[whichcomponent*2 + 1])
		    bounds[whichcomponent*2 + 1] = rval;
		}
	      /* update indices */
	      /* row major */
              for(k = ndim - 1; k >= 0; k--)
		{
		  indices[k]++;
		  if (indices[k] == (*dimsptr)[i][k])
		    {
		      indices[k] = 0;
		      continue;
		    }
		  break;
		}
	    }

          Free(gdvals);
	  (*formptr)[i]=NewCurvilinearDataForm(datasetname, ndim,
	                           (*dimsptr)[i], bounds, dataset);
	}
      Free(indices);
      Free(bounds);
    }

/* close file with fortran subroutine */
#ifdef FORTRAN_
  forcls_(&lun);
#else
  forcls(&lun);
#endif

  Free(indims); 
  if (gridptr->blank != P3D_BLANK)
    Free(*ipp);
  return i;
} /* RdUnfmtXYZ */


/* for reading fortran unformatted solution file */
#ifdef HAVE_PROTOTYPES
static int RdUnfmtSol(ObjPtr reader, int ndim, int isMGrid, int ngrid, int blank, 
                          long **dimsptr, p3d_q *qptr, ObjPtr *formptr, int **ipp)
#else
static int RdUnfmtSol(reader, ndim, isMGrid, ngrid, blank, dimsptr, qptr, formptr, ipp)
ObjPtr reader;
int ndim;
int isMGrid;
int ngrid;
int blank;
long **dimsptr;
p3d_q *qptr;
ObjPtr *formptr;
int **ipp;
#endif
{
  int ins[P3D_MAXNAMELEN]; /* ascii coding of file name */
  int namelen;         /* file name length */
  int lun = 2;         /* logical unit number of fortran opened file */
  int opstat=0;        /* operation status of fortran file op subroutine */
  int *indims;
  int numsets;         /* number of datasets read so far */
  int i, j;
  int myngrid;         /* the ngrid read from the solution file */
  long *indices;
  char *name;   /* dataset name, extracted from file name */


/* try to open the solution file, ascii coding of the solution file name */
  namelen = (int) strlen(qptr->name);
  for (i = 0; i <= namelen; i++)
     {
       ins[i] = qptr->name[i];
     }
/* open file with fortran subroutine */
#ifdef FORTRAN_
  foropn_(&lun, ins, &namelen, &opstat);
#else
  foropn(&lun, ins, &namelen, &opstat);
#endif
  if (!opstat)
    {
      fprintf(stderr, "Unable to open file %s, aborted\n", qptr->name);
      return 0;
    }

  /* extracting the dataset name from file name */
  for(i = namelen -1; i >= 0 && qptr->name[i] != '/'; i--)
    ;
  if (i == -1 || qptr->name[i] == '/')
    i++;
  name = qptr->name + i;

  for (namelen = 0; 
       name[namelen] != '\0' && name[namelen] != '.';
       namelen++) ;

  /* allocate storage for indices */
  if (!(indices = (long *) Alloc(ndim*sizeof(long))))
    {
      fprintf(stderr, "Unable to alloc memory, aborted\n");
#ifdef FORTRAN_
      forcls_(&lun);
#else
      forcls(&lun);
#endif
      return 0;
    }

  /* 0 dataset read so far */
  numsets = 0;

/* do multiple datasets starting here */
  do
    {
      /* reading myngrid */
      if (isMGrid)
        {
          opstat = 0;
#ifdef FORTRAN_
          rngrid_(&lun, &myngrid, &opstat);
#else
          rngrid(&lun, &myngrid, &opstat);
#endif
          if (opstat == EOF)
    	{
    	  if (numsets == 0)
    	    {
    	      fprintf(stderr, "Error in %s: Empty file\n", qptr->name);
    	    }
#ifdef FORTRAN_
              forcls_(&lun);
#else
              forcls(&lun);
#endif
              Free(indices);
    	  /* if numsets >= 1 then this is an EOF, else error */
    	  return (numsets ? ngrid : 0);
    	}
          else if (!opstat || myngrid <=0)
            {
    	  fprintf(stderr, "%s format error, invalid NGRID value\n", qptr->name);
#ifdef FORTRAN_
              forcls_(&lun);
#else
              forcls(&lun);
#endif
              Free(indices);
    	  return 0;
            }
        }
      else
        {
          myngrid = 1;
        }
    
      /* if myngrid is not the same as ngrid then error */
      if (myngrid != ngrid)
        {
          fprintf(stderr, "Error in %s: NGRID value does not match that"
    	      "of its XYZ file\n", qptr->name);
#ifdef FORTRAN_
          forcls_(&lun);
#else
          forcls(&lun);
#endif
          Free(indices);
          return 0;
        }
    
       /* reading dimensions */
      indims = (int *) 0;
      if (!(indims = (int *) Alloc(ndim*myngrid*sizeof(int))))
        {
          fprintf(stderr, "Unable to alloc mem, aborted\n");
#ifdef FORTRAN_
          forcls_(&lun);
#else
          forcls(&lun);
#endif
          Free(indices);
          return 0;
        }
      opstat=0;
#ifdef FORTRAN_
      rddims_(&lun, &ndim, &ngrid, indims, &opstat);
#else
      rddims(&lun, &ndim, &ngrid, indims, &opstat);
#endif
      if (opstat == EOF && !isMGrid)
        {
          if (numsets == 0)
    	{
    	  fprintf(stderr, "Error in %s: Empty file\n", qptr->name);
    	}
#ifdef FORTRAN_
          forcls_(&lun);
#else
          forcls(&lun);
#endif
          Free(indices);
          Free(indims);
          return (numsets ? ngrid : 0);
        }
      else if (opstat > 0)
        {
          /* Check dimension consistency with dims read from XYZ file */
          for (i = 0; i < ngrid; i++)
    	 {
    	   for (j = 0; j < ndim; j++)
                 {
    	       if (dimsptr[i][j] != (long) indims[i*ndim + j])
    		 {
    		   fprintf(stderr, "Dimensions mismatch between %s "
    		   "and its XYZ file\n", qptr->name);
#ifdef FORTRAN_
                       forcls_(&lun);
#else
                       forcls(&lun);
#endif
                       Free(indices);
                       Free(indims);
                       return 0;
                     }
                 }
             }
        }
      else
      /* opstat is 0 or (opstat is EOF and isMGrid is true) */
        {
          fprintf(stderr, "Unable to Read Dimension Values in %s\n", qptr->name);
#ifdef FORTRAN_
              forcls_(&lun);
#else
              forcls(&lun);
#endif
          Free(indices);
          Free(indims);
          return 0;
        }

      for(i = 0; i < ngrid; i++)
	{
          long size, count;
          float *slvals;
	  ObjPtr density, energy, momentum;
	  char densityName[P3D_MAXNAMELEN];
	  char energyName[P3D_MAXNAMELEN];
	  char momentumName[P3D_MAXNAMELEN];
	  char nameBase[P3D_MAXNAMELEN];
	  float time;
	  int whichcomponent;

	  /* get the time data of the current grid solutions */
          opstat = 0;
#ifdef FORTRAN_
          rdtime_(&lun, &time, &opstat);
#else
          rdtime(&lun, &time, &opstat);
#endif
          if (!opstat)
	    {
	      fprintf(stderr, "Error in %s, invalid FAMACH...TIME value\n",
		      qptr->name);
#ifdef FORTRAN_
              forcls_(&lun);
#else
              forcls(&lun);
#endif
              Free(indices);
              Free(indims);
              return 0;
	    }

          /* calculate size of one component */
          size=1;
          for(j = 0; j < ndim; j++)
           {
             size *= dimsptr[i][j];
           }
	  /* read in solution, density, and pressure values */
          if (!(slvals = (float *) Alloc((ndim+2)*size*sizeof(float))))
	    {
                 fprintf(stderr, "Unable to alloc mem, aborted\n");
#ifdef FORTRAN_
                 forcls_(&lun);
#else
                 forcls(&lun);
#endif
                 Free(indims);
	         Free(indices);
                 return 0;
	    }
	  opstat=0;
#ifdef FORTRAN_
          rdsolu_(&lun, &ndim, indims+i*ndim, &(qptr->permutation), 
		  slvals, &size, &opstat);
#else
          rdsolu(&lun, &ndim, indims+i*ndim, &(qptr->permutation), 
		  slvals, &size, &opstat);
#endif
          if (!opstat)
	    {
	      fprintf(stderr, "Read error on file %s\n", qptr->name);
#ifdef FORTRAN_
              forcls_(&lun);
#else
              forcls(&lun);
#endif
              Free(indims);
	      Free(indices);
	      Free(slvals);
              return 0;
	    }

          strncpy(nameBase, name, namelen);
	  nameBase[namelen] = '\0';
	  /* create datasets */
          if (ngrid > 1)
            {
              sprintf(densityName, "%s %s %d", nameBase, "Density", i + 1);
              sprintf(energyName, "%s %s %d", nameBase, "Stagnation Energy per Unit Volume", i + 1);
              sprintf(momentumName, "%s %s %d", nameBase, "Momentum", i + 1);
            }
          else
	    {
              sprintf(densityName, "%s %s", nameBase, "Density");
              sprintf(energyName, "%s %s", nameBase, "Stagnation Energy per Unit Volume");
              sprintf(momentumName, "%s %s", nameBase, "Momentum");
	    }
          density = NewStructuredDataset(densityName, ndim, dimsptr[i], 0);
          energy = NewStructuredDataset(energyName, ndim, dimsptr[i], 0);
          momentum = NewStructuredDataset(momentumName, ndim, dimsptr[i], ndim);

	  /* assign slvals to densiity, momentum, and energy */
          for(whichcomponent = 0; whichcomponent < ndim + 2; whichcomponent++)
	    {
	      ObjPtr whichdataset;

	      if (whichcomponent == 0)
		/* for density */
		{
		  whichdataset = density;
	        }
              else if (whichcomponent <= 3)
		/* for momentum */
		{
		  whichdataset = momentum;
		}
              else
		/* for stagnation energy */
		{
		  whichdataset = energy;
		}

	      if (!SetCurField(FIELD1, whichdataset))
	        {
                  fprintf(stderr, "Unable to SETCURFIELD, aborted\n");
#ifdef FORTRAN_
                  forcls_(&lun);
#else
                  forcls(&lun);
#endif
                  Free(indims);
	          Free(indices);
		  Free(slvals);
                  return 0;
	        }
              for(j = 0; j < ndim; j++)
	        {
	          indices[j] = 0;
                }
              for(count = 0; count < size; count++)
	        {
	          int k, iblank = 1;
		  real rval;

	          rval = (real) slvals[whichcomponent*size + count];
	          /* do iblank if needed */
	          if (blank == P3D_BLANK)
	            {
	              getIblank(&iblank, ipp, indices, dimsptr[i], ndim, i);
		    }
	          /* do check of negative values of density or pressure 
		     if required */
		  if ((whichcomponent == 0 || whichcomponent == 4) && 
		      qptr->check == P3D_CHECK && rval < (real) 0.0)
		    {
		      fprintf(stderr, "Negative %s value %f found in Q file %s and is converted to %f\n", whichdataset == energy ? "pressure" : "density", rval, qptr->name, 0.0);
		      rval = (real) 0.0;
		    }
		  PutFieldComponent(FIELD1, 
		     (whichcomponent == 0 || whichcomponent == 4) ? 0 : whichcomponent - 1, 
		      indices, iblank == 0 ? (real) missingData : rval);
	          /* update indices */
		  /* row major way */
                  for(k = ndim - 1; k >= 0; k--)
		    {
		      indices[k]++;
		      if (indices[k] == dimsptr[i][k])
		        {
		          indices[k] = 0;
		          continue;
		        }
		      break;
		    }
		} /* for count = 0 ~ size */
	    } /* for whichcomponent = 0 ~ ndim+2 */

          Free(slvals);
	  SetDatasetForm(density, formptr[i]);
	  SetDatasetForm(momentum, formptr[i]);
	  SetDatasetForm(energy, formptr[i]);
	  if ( qptr->isMDataset )
	    {
	      RegisterNewTimeDataset(density, (real) time);
	      RegisterNewTimeDataset(momentum, (real) time);
	      RegisterNewTimeDataset(energy, (real) time);
#if 0
	      RegisterTimeDataset(density, (real) time);
	      RegisterTimeDataset(momentum, (real) time);
	      RegisterTimeDataset(energy, (real) time);
#endif
	    }
          else
	    {
	      RegisterDataset(density);
	      RegisterDataset(momentum);
	      RegisterDataset(energy);
	    }
	}/* for loop for i = 0, ngrid - 1 */
      numsets++;
    } while(qptr->isMDataset);

/* close file with fortran subroutine */
#ifdef FORTRAN_
  forcls_(&lun);
#else
  forcls(&lun);
#endif

  Free(indims); 
  Free(indices);
  return ngrid;
} /* RdUnfmtSol */

/* for reading fortran unformatted function file */
#ifdef HAVE_PROTOTYPES
static int RdUnfmtFun(int ndim, int isMGrid, int ngrid, int blank, long **dimsptr,
		      p3d_func *funcptr, ObjPtr *formptr, int **ipp)
#else
static int RdUnfmtFun(ndim, isMGrid, ngrid, blank, dimsptr, funcptr, formptr, ipp)
int ndim;
int isMGrid;
int ngrid;
int blank;
long **dimsptr;
p3d_func *funcptr;
ObjPtr *formptr;
int **ipp;
#endif
{
  int ins[P3D_MAXNAMELEN]; /* ascii coding of file name */
  int namelen;         /* file name length */
  int lun = 2;         /* logical unit number of fortran opened file */
  int opstat=0;        /* operation status of fortran file op subroutine */
  int *indims;
  int i, j;
  int myngrid;         /* the ngrid read from the solution file */
  int nvardim;         /* number of dimensions + 1 */
  long *indices;
  char *name;   /* dataset name, extracted from file name */


/* try to open the solution file, ascii coding of the solution file name */
  namelen = (int) strlen(funcptr->name);
  for (i = 0; i <= namelen; i++)
     {
       ins[i] = funcptr->name[i];
     }
/* open file with fortran subroutine */
#ifdef FORTRAN_
  foropn_(&lun, ins, &namelen, &opstat);
#else
  foropn(&lun, ins, &namelen, &opstat);
#endif
  if (!opstat)
    {
      fprintf(stderr, "Unable to open file %s, aborted\n", funcptr->name);
      return 0;
    }

  /* read myngrid */
  if (isMGrid)
    {
      opstat = 0;
#ifdef FORTRAN_
      rngrid_(&lun, &myngrid, &opstat);
#else
      rngrid(&lun, &myngrid, &opstat);
#endif
      if (opstat <= 0 || myngrid <=0)
        {
	  fprintf(stderr, "%s format error, invalid NGRID value\n", funcptr->name);
#ifdef FORTRAN_
          forcls_(&lun);
#else
          forcls(&lun);
#endif
	  return 0;
        }
    }
  else
    {
      myngrid = 1;
    }

  /* if myngrid is not the same as ngrid then error */
  if (myngrid != ngrid)
    {
      fprintf(stderr, "Error in %s: NGRID value does not match that"
	      "of its XYZ file\n", funcptr->name);
#ifdef FORTRAN_
      forcls_(&lun);
#else
      forcls(&lun);
#endif
      return 0;
    }

  /* extracting the dataset name from file name */
  for(i = namelen -1; i >= 0 && funcptr->name[i] != '/'; i--)
    ;
  if (i == -1 || funcptr->name[i] == '/')
    i++;
  name = funcptr->name + i;

  for (namelen = 0; 
       name[namelen] != '\0' && name[namelen] != '.';
       namelen++) ;

  /* allocate storage for indices */
  if (!(indices = (long *) Alloc(ndim*sizeof(long))))
    {
      fprintf(stderr, "Unable to alloc memory, aborted\n");
#ifdef FORTRAN_
      forcls_(&lun);
#else
      forcls(&lun);
#endif
      return 0;
    }

  /* reading dimensions */
  indims = (int *) 0;
  nvardim = ndim + 1;
  if (!(indims = (int *) Alloc(nvardim*myngrid*sizeof(int))))
    {
      fprintf(stderr, "Unable to alloc mem, aborted\n");
#ifdef FORTRAN_
      forcls_(&lun);
#else
      forcls(&lun);
#endif
      Free(indices);
      return 0;
    }
  opstat=0;
#ifdef FORTRAN_
  rddims_(&lun, &nvardim, &ngrid, indims, &opstat);
#else
  rddims(&lun, &nvardim, &ngrid, indims, &opstat);
#endif
  if (opstat > 0)
    {
      /* Check dimension consistency with dims read from XYZ file */
      for (i = 0; i < ngrid; i++)
	 {
	   for (j = 0; j < ndim; j++)
             {
	       if (dimsptr[i][j] != (long) indims[i*nvardim + j])
		 {
		   fprintf(stderr, "Dimensions mismatch between %s "
		   "and its XYZ file\n", funcptr->name);
#ifdef FORTRAN_
                   forcls_(&lun);
#else
                   forcls(&lun);
#endif
                   Free(indices);
                   Free(indims);
                   return 0;
                 }
             }
         }
    }
  else
  /* opstat <= 0 */
    {
      fprintf(stderr, "Unable to Read Dimension Values in %s\n", funcptr->name);
#ifdef FORTRAN_
          forcls_(&lun);
#else
          forcls(&lun);
#endif
      Free(indices);
      Free(indims);
      return 0;
    }

      for(i = 0; i < ngrid; i++)
	{
          long size, count;
          float *fnvals;
	  ObjPtr dataset;
	  char datasetname[P3D_MAXNAMELEN];

          /* calculate size of one component */
          size=1;
          for(j = 0; j < ndim; j++)
           {
             size *= dimsptr[i][j];
           }
	  /* read in function values */
          if (!(fnvals = (float *) Alloc(nvardim*size*sizeof(float))))
	    {
                 fprintf(stderr, "Unable to alloc mem, aborted\n");
#ifdef FORTRAN_
                 forcls_(&lun);
#else
                 forcls(&lun);
#endif
                 Free(indims);
	         Free(indices);
                 return 0;
	    }
	  opstat=0;
#ifdef FORTRAN_
          rdfunc_(&lun, &ndim, indims+i*nvardim, indims+i*nvardim+ndim, 
	           &(funcptr->permutation), fnvals, &size, &opstat);
#else
          rdfunc(&lun, &ndim, indims+i*nvardim, indims+i*nvardim+ndim, 
	           &(funcptr->permutation), fnvals, &size, &opstat);
#endif
          if (!opstat)
	    {
	      fprintf(stderr, "Read error on file %s\n", funcptr->name);
#ifdef FORTRAN_
              forcls_(&lun);
#else
              forcls(&lun);
#endif
              Free(indims);
	      Free(indices);
	      Free(fnvals);
              return 0;
	    }

	  /* create datasets for function */
          strncpy(datasetname, name, namelen);
	  datasetname[namelen] = '\0';
	  if (ngrid > 1)
	    sprintf(datasetname + namelen, "%d", i + 1);
          dataset=NewStructuredDataset(datasetname, ndim, dimsptr[i], 
		   indims[i*nvardim+ndim] == 1 ? 0 : indims[i*nvardim+ndim]);

	  /* assign fnvals to dataset */
          for(j = 0; j < ndim; j++)
	    {
	      indices[j] = 0;
            }
          if (!SetCurField(FIELD1, dataset))
            {
              fprintf(stderr, "Unable to SETCURFIELD, aborted\n");
#ifdef FORTRAN_
              forcls_(&lun);
#else
              forcls(&lun);
#endif
              Free(indims);
              Free(indices);
              Free(fnvals);
              return 0;
	    }
          for(count = 0; count < size; count++)
	    {
	      int whichcomponent, k;
	      int iblank = 1;

	      /* do iblank if needed */
	      if (blank == P3D_BLANK)
	        {
	          getIblank(&iblank, ipp, indices, dimsptr[i], ndim, i);
		}
	      /* assign solution values */
	      for(whichcomponent = 0; whichcomponent < indims[i*nvardim+ndim];
		  whichcomponent++)
		{
		  real rval;

		  rval = (real) fnvals[whichcomponent*size + count];

		  PutFieldComponent(FIELD1, whichcomponent, indices, 
		  iblank == 0 ? (real) missingData : rval);
		}
	      /* update indices */
	      /* row major way */
              for(k = ndim - 1; k >= 0; k--)
		{
		  indices[k]++;
		  if (indices[k] == dimsptr[i][k])
		    {
		      indices[k] = 0;
		      continue;
		    }
		  break;
		}
            }
          Free(fnvals);
          SetDatasetForm(dataset, formptr[i]);
       	  RegisterDataset(dataset);
	}/* for loop for i = 0, ngrid - 1 */

/* close file with fortran subroutine */
#ifdef FORTRAN_
  forcls_(&lun);
#else
  forcls(&lun);
#endif

  Free(indims); 
  Free(indices);
  return ngrid;
} /* RdUnfmtFun */

#endif /* for ifdef FORTRAN */

/*********************unformatted p3d reading routines, end ***********************/

/*********************binary p3d reading routines, start *********************/

#ifdef HAVE_PROTOTYPES
static int RdBinXYZ(int ndim, int isMGrid, long ***dimsptr,
                               p3d_grid *gridptr, ObjPtr **formptr, int ***ipp)
#else
static int RdBinXYZ(ndim, isMGrid, dimsptr, gridptr, formptr, ipp)
int ndim;
int isMGrid;
long ***dimsptr;
p3d_grid *gridptr;
ObjPtr **formptr;
int ***ipp;
#endif
{
  FILE *inFile;
  int ngrid;            /* number of grids */

  int i, j;
  int n;

  char *name;           /* dataset name starting point */
  int namelen;

  int ival;             /* an integer value */
  float fval;           /* a float value */

  long *indices;
  real *bounds;

  if (!(inFile = fopen(gridptr->name, "r")))
    {
      fprintf(stderr, "Unable to open %s, aborted\n", gridptr->name);
      return 0;
    }

  /* determine the dataset name now */
  for (i = (int) strlen(gridptr->name) - 1; i >= 0 && gridptr->name[i] != '/'; i--)
     ;
  if (i == -1 || gridptr->name[i] == '/')
    i++;
  name = gridptr->name + i;
  for (namelen = 0; name[namelen] != '\0' && name[namelen] != '.'; namelen++)
     ;

  if (isMGrid)
    {
      if (1 != fread((void *) &ival, sizeof(int), 1, inFile))
        {
          fprintf(stderr, "Error reading NGRID value in %s\n", gridptr->name);
          fclose(inFile);
          return 0;
        } 
      else if (ival <= 0)
        {
          fprintf(stderr, "%s format error, invalid NGRID value\n", gridptr->name);
          fclose(inFile);
          return 0;
        } 
      else
	{
	  ngrid = ival;
	}
    }
  else
    {
      ngrid = 1;
    }
  
  /* alloc dimsptr and indices storage */
  if (!(indices = (long *) Alloc(ndim * sizeof(long))))
    {
      fprintf(stderr, "Unable to alloc memory\n");
      fclose(inFile);
      return 0;
    }
  if (!(*dimsptr = (long **) Alloc(ngrid * sizeof(long *))))
    {
      fprintf(stderr, "Unable to alloc memory\n");
      Free(indices);
      fclose(inFile);
      return 0;
    }

  /* read in all the dimensions */
  for (i = 0; i < ngrid; i++)
    {
      if (!((*dimsptr)[i] = (long *) Alloc(ndim * sizeof(long))))
        {
          fprintf(stderr, "Unable to alloc memory\n");
          Free(indices);
          freedims(i, *dimsptr);
          fclose(inFile);
          return 0;
        }

      for (j = 0; j < ndim; j++)
         {
           if (1 != fread((void *) &ival, sizeof(int), 1, inFile) || ival <= 0)
             {
               fprintf(stderr, "Invalid dimension value in %s\n", gridptr->name);
               Free(indices);
               freedims(i, *dimsptr);
               fclose(inFile);
               return 0;
             }
           *((*dimsptr)[i] + j) = (long) ival; 
         }
    }

  /* alloc strorage for pointer to bounds */
  if (!(bounds = (real *) Alloc(ndim * 2 * sizeof(real))))
    {
      fprintf(stderr, "Unable to alloc memory\n");
      Free(indices);
      freedims(ngrid, *dimsptr);
      fclose(inFile);
      return 0;
    }
  /* alloc storage for pointers to dataforms */
  if (!(*formptr = (ObjPtr *) Alloc(ngrid * sizeof(ObjPtr))))
    {
      fprintf(stderr, "Unable to alloc memory\n");
      Free(indices);
      Free(bounds);
      freedims(ngrid, *dimsptr);
      fclose(inFile);
      return 0;
    }

  /* alloc iblankptr if IBLANK */
  if (gridptr->blank == P3D_BLANK)
    {
      if (!(*ipp = (int **) Alloc(ngrid * sizeof(int *))))
        {
          fprintf(stderr, "Unable to alloc memory\n");
          Free(indices);
          Free(bounds);
          Free(*formptr);
          freedims(ngrid, *dimsptr);
          fclose(inFile);
          return 0;
        }
    }

  for (n = 0; n < ngrid; n++)
    {
      long size, count; /* number of grid values promised */
      ObjPtr dataset;
      char datasetname[P3D_MAXNAMELEN];
      int whichcomponent;

      strncpy(datasetname, name, namelen);
      datasetname[namelen] = '\0';
      if (ngrid > 1)
        {
          sprintf(datasetname + namelen, "%d", n + 1);
        }

      /* setup a dataset for the grid data */
      dataset = NewStructuredDataset(datasetname, ndim, (*dimsptr)[n], ndim);
      if (!SetCurField(FIELD1, dataset))
        {
          fprintf(stderr, "Unable to SetCurField\n");
          Free(indices);
          Free(bounds);
          Free(*formptr);
          freedims(ngrid, *dimsptr);
          if (gridptr->blank == P3D_BLANK)
            {
              freeiblank(n-1, *ipp);
            }
          fclose(inFile);
          return 0;
        }

      size = 1;
      for (i = 0; i < ndim; i++)
         {
           size *= (*dimsptr)[n][i];
         }
      size *= (long)(ndim + (gridptr->blank == P3D_BLANK));

      if (gridptr->blank == P3D_BLANK)
        {
          /* alloc storage for the iblank arrays */
          if (!((*ipp)[n] = (int *) Alloc(size * sizeof(int))))
            {
              fprintf(stderr, "Unable to alloc memory\n");
              freeiblank(n, *ipp);
              Free(indices);
              Free(bounds);
              Free(*formptr);
              freedims(ngrid, *dimsptr);
              fclose(inFile);
              return 0;
            }
        }

      for (i = 0; i < ndim; i++)
         {
           indices[i] = 0;
         }
      whichcomponent = 0;
      /* read in all the grid values */
      for (count = 0; count < size; count++)
        {
          int first = 1;

          if (whichcomponent < ndim) /* read grid values */
            {
              if (1 != fread((void *) &fval, sizeof(float), 1, inFile))
                {
	          if (feof(inFile))
		    {
                      fprintf(stderr, "Short of data in %s\n", gridptr->name);
		    }
                  else
		    {
                      fprintf(stderr, "Error reading %s\n", gridptr->name);
		    }
                  Free(indices);
                  Free(bounds);
                  Free(*formptr);
                  freedims(ngrid, *dimsptr);
                  if (gridptr->blank == P3D_BLANK)
                    {
                      freeiblank(n + 1, *ipp);
                    }
                  fclose(inFile);
                  return 0;
                }
              PutFieldComponent(FIELD1, whichcomponent, indices, (real) fval);  
              for (j = 0; first && j < ndim; j++)
                 first = indices[j] == 0;

              if (first)
                {
                  bounds[whichcomponent * 2] = (real) fval;
                  bounds[whichcomponent * 2 + 1] = (real) fval;
                }
              else
                {
                  if (fval < bounds[whichcomponent * 2])
                    bounds[whichcomponent * 2] = (real) fval;
                  else if (fval > bounds[whichcomponent * 2 + 1])
                    bounds[whichcomponent * 2 + 1] = (real) fval;
                }
            }
          else /* read iblank values */
            {
              if (1 != fread((void *) &ival, sizeof(int), 1, inFile))
                {
	          if (feof(inFile))
		    {
                      fprintf(stderr, "Short of data in %s\n", gridptr->name);
		    }
                  else
		    {
                      fprintf(stderr, "Error reading %s\n", gridptr->name);
		    }
                  Free(indices);
                  Free(bounds);
                  Free(*formptr);
                  freedims(ngrid, *dimsptr);
                  if (gridptr->blank == P3D_BLANK)
                    {
                      freeiblank(n + 1, *ipp);
                    }
                  fclose(inFile);
                  return 0;
                }
              putIblank(ival, *ipp, indices, (*dimsptr)[n], ndim, n);
            }

          /* adjust the indices and whichcomponent */
          if (gridptr->permutation == P3D_WHOLE)
            {
              for (j = 0; j < ndim; j++)
                 {
                   indices[j]++;
                   if (indices[j] == (*dimsptr)[n][j])
                     {
                       indices[j] = 0;
                       continue;
                     }
                   break;
                 }
              if (j == ndim) /* all indices are set to zero now */ 
                {
                  whichcomponent++;
                }
            }
          else  /* P3D_PLANE */
            {
              for (j = 0; j < ndim - 1; j++)
                 {
                   indices[j]++;
                   if (indices[j] == (*dimsptr)[n][j])
                     {
                       indices[j] = 0;
                       continue;
                     }
                   break;
                 }
              if (j == ndim - 1)  
                {
                  whichcomponent++;
                  if ((gridptr->blank == P3D_BLANK) ? 
                     (whichcomponent == ndim + 1) : (whichcomponent == ndim))
                  /* leave a space for reading iblanks */
                    {
                      whichcomponent = 0;
                      indices[j]++;
                    }
                }
            }
        } /* for read grid values */

     (*formptr)[n] = NewCurvilinearDataForm(datasetname, ndim, (*dimsptr)[n],
                    bounds, dataset);

    } /* outer for */

  /* indices and bounds are no longer of any use */
  Free(indices);
  Free(bounds);
  fclose(inFile);
  return n;
} /* RdBinXYZ */

#ifdef HAVE_PROTOTYPES
static int RdBinSol(ObjPtr reader, int ndim, int isMGrid, int ngrid, int blank,
                 long **dimsptr, p3d_q *qptr, ObjPtr *formptr, int **ipp)
#else
static int RdBinSol(reader, ndim, isMGrid, ngrid, blank, dimsptr, qptr, formptr, ipp)
ObjPtr reader;
int ndim;
int isMGrid;
int ngrid;
int blank;
long **dimsptr;
p3d_q *qptr;
ObjPtr *formptr;
int **ipp;
#endif
{
  FILE *inFile;
  int i, j, n;
  char *name;
  long *indices;
  int numsets; /* number of datasets read so far, for multiple datasets in a file*/
  int namelen; /* the length of dataset name going to be used */
  int ival;
  float fval;

  if (!(inFile = fopen(qptr->name, "r")))
    {
      fprintf(stderr, "Unable to open %s, aborted\n", qptr->name);
      return 0;
    }

  if (!(indices = (long *) Alloc(ndim * sizeof(long))))
    {
      fprintf(stderr, "Unable to alloc memory, aborted\n");
      fclose(inFile);
      return 0;
    }

  /* determine the dataset name from the file name */
  for (i = (int) strlen(qptr->name) - 1; i >= 0 && qptr->name[i] != '/'; i--)
     ;
  if (i == -1 || qptr->name[i] == '/')
    i++;
  name = qptr->name + i;
  for (namelen = 0; name[namelen] != '\0' && name[namelen] != '.'; namelen++)
     ;

  numsets = 0;

  /* We have to consider multiple datasets here. */
  do
    {
      if (isMGrid)
	{
          if (1 != fread((void *) &ival, sizeof(int), 1, inFile))
            {
              if (feof(inFile))
                {
                  if (numsets == 0)
                    {
                      fprintf(stderr, "Error in %s: empty file\n", qptr->name);
                    } /* else it reads to the end of a Mdata file already */
                  fclose(inFile);
                  Free(indices);
                  return (numsets ? ngrid : 0);
                }
              else 
                {
                  fprintf(stderr, "Error reading %s.\n", qptr->name);
                  fclose(inFile);
                  Free(indices);
                  return 0;
                }
            }
          else if (ival <= 0)
	    {
              fprintf(stderr, "Error in %s: Invalid NGRID spec.\n", qptr->name);
              fclose(inFile);
              Free(indices);
              return 0;
	    }
          else
	    {
	      n = ival;
	    }
        }
      else
        {
          n = 1;
        }

      /* if ngrid != this file's ngrid then error */
      if (ngrid != n)
        {
          fprintf(stderr, "Error in %s: NGRID value does not match its grid\n",
                                                             qptr->name);
          fclose(inFile);
          Free(indices);
          return 0;
        }

      /* reading dimensions */
      for (n = 0; n < ngrid; n++)
        {
          for (i = 0; i < ndim; i++)
            {
              if (1 != fread((void *) &ival, sizeof(int), 1, inFile))
                {
                  /* if this is not multiple grid, then the first dimension value 
		     is used to indicate whether there is more datasets or not */
                  if (!isMGrid && i == 0 && feof(inFile))
                    {
                      if (numsets == 0)
                        {
                          fprintf(stderr, "Error in %s: empty file\n", qptr->name);
                        }
                      fclose(inFile);
                      Free(indices);
                      return (numsets ? ngrid : 0);
                    }
                  fprintf(stderr, "Error reading dimensions in %s\n", qptr->name);
                  fclose(inFile);
                  Free(indices);
                  return 0;
                }
              /* if dims does not match the dims read in XYZ file then error */
              if ((long) ival != dimsptr[n][i])
                {
                  fprintf(stderr,
                  "Dimensions mismatch between %s and its XYZ file\n", qptr->name);
                  fclose(inFile);
                  Free(indices);
                  return 0;
                }
            }
        }

      for (n = 0; n < ngrid; n++)
        {
          long size, count; /* number of data values promised */
          ObjPtr density, energy, momentum;
          char momentumName[P3D_MAXNAMELEN];
          char densityName[P3D_MAXNAMELEN];
          char energyName[P3D_MAXNAMELEN];
          char nameBase[P3D_MAXNAMELEN];
          double time;
          int whichcomponent;
          ObjPtr whichdataset;
    
          /* process the FAMACH..TIME line */
          for (i = 1; i <= 4; i++)
            {
              if (1 != fread((void *) &fval, sizeof(float), 1, inFile))
                {
                  fprintf(stderr,
                     "Error in %s: invalid FSMACH,ALPHA,RE,TIME line\n", qptr->name);
                  Free(indices);
                  fclose(inFile);
                  return 0;
                }
              if (i == 4)
                {
                  time = (double) fval;
                }
            }
    
	  strncpy(nameBase, name, namelen);
	  nameBase[namelen] = '\0';
	  /* create datasets */
          if (ngrid > 1)
            {
              sprintf(densityName, "%s %s %d", nameBase, "Density", n + 1);
              sprintf(energyName, "%s %s %d", nameBase, "Stagnation Energy per Unit Volume", n + 1);
              sprintf(momentumName, "%s %s %d", nameBase, "Momentum", n + 1);
            }
          else
	    {
              sprintf(densityName, "%s %s", nameBase, "Density");
              sprintf(energyName, "%s %s", nameBase, "Stagnation Energy per Unit Volume");
              sprintf(momentumName, "%s %s", nameBase, "Momentum");
	    }
          density = NewStructuredDataset(densityName, ndim, dimsptr[n], 0);
          energy = NewStructuredDataset(energyName, ndim, dimsptr[n], 0);
          momentum = NewStructuredDataset(momentumName, ndim, dimsptr[n], ndim);

          size = 1;
          for (i = 0; i < ndim; i++)
             size *= dimsptr[n][i];
          size *= ndim + 2;  /* includes density and energy */
          for (i = 0; i < ndim; i++)
             {
               indices[i] = 0;
             }
          whichcomponent = 0;
          whichdataset = density;
    
          /* read in all the q values */
          for (count = 0; count < size; count++)
            {
              int iblank = 1;
    
              if (1 != fread((void *) &fval, sizeof(float), 1, inFile))
                {
		  if (feof(inFile))
		    {
                      fprintf(stderr, "Error reading %s, short of data\n", 
			      qptr->name);
		    }
                  else
		    {
                      fprintf(stderr, "Error reading %s\n", qptr->name);
		    }
                  Free(indices);
                  fclose(inFile);
                  return 0;
                }
              /* do check of density and energy, and iblanks here */
              if (qptr->check == P3D_CHECK && (whichdataset == energy ||
                  whichdataset == density))
                {
                  if (fval < 0.0)
		    {
		      fprintf(stderr, "Negative %s value %f found in Q file %s and is converted to %f\n", whichdataset == energy ? "pressure" : "density", fval, qptr->name, 0.0);
                      fval = 0.0;
                    }
                }
              if (blank == P3D_BLANK)
                {
                   getIblank(&iblank, ipp, indices, dimsptr[n], ndim, n);
                }
    
              if (!SetCurField(FIELD1, whichdataset))
                {
                  fprintf(stderr, "Unable to SetCurField\n");
                  Free(indices);
                  fclose(inFile);
                  return 0;
                }
    
              PutFieldComponent(FIELD1, whichcomponent, indices, 
                          iblank == 0 ? (real) missingData : (real) fval);  
    
              /* adjust the indices, whichdataset, and whichcomponent */
              if (qptr->permutation == P3D_WHOLE)
                {
                  for (j = 0; j < ndim; j++)
                     {
                       indices[j]++;
                       if (indices[j] == dimsptr[n][j])
                         {
                           indices[j] = 0;
                           continue;
                         }
                       break;
                     }
                  if (j == ndim) /* all indices are set to zero now */ 
                    {
		      if (whichdataset == density)
			{
			  whichdataset = momentum;
                        }
                      else if (++whichcomponent == ndim)
			{
			  whichdataset = energy;
			  whichcomponent = 0;
                        }
                    }
                }
              else  /* P3D_PLANE */
                {
                  for (j = 0; j < ndim - 1; j++)
                     {
                       indices[j]++;
                       if (indices[j] == dimsptr[n][j])
                         {
                           indices[j] = 0;
                           continue;
                         }
                       break;
                     }
                  if (j == ndim - 1)  
                    {
                      if (whichdataset == density)
			{
                          whichdataset = momentum;
			}
                      else if (whichdataset == momentum)
			{
			  if (++whichcomponent == ndim)
			    {
			      whichdataset = energy;
			      whichcomponent = 0;
			    }
			}
                      else /* whichdataset is energy now */
                        {
                          whichdataset = density;
                          indices[j]++;
                        }
                    }
                }
            } 
    
          SetDatasetForm(density, formptr[n]);
          SetDatasetForm(momentum, formptr[n]);
          SetDatasetForm(energy, formptr[n]);
          if (qptr->isMDataset)
            {
              RegisterNewTimeDataset(density, (real) time);
              RegisterNewTimeDataset(momentum, (real) time);
              RegisterNewTimeDataset(energy, (real) time);
#if 0
              RegisterTimeDataset(density, (real) time);
              RegisterTimeDataset(momentum, (real) time);
              RegisterTimeDataset(energy, (real) time);
#endif
            }
          else
            {
              RegisterDataset(density);
              RegisterDataset(momentum);
              RegisterDataset(energy);
            }
        } /* outer for */
      numsets++;
    } while(qptr->isMDataset);

  Free(indices);
  fclose(inFile);
  return ngrid;
} /* RdBinSol */

#ifdef HAVE_PROTOTYPES
static int RdBinFun(int ndim, int isMGrid, int ngrid, int blank,
                long **dimsptr, p3d_func *funcptr, ObjPtr *formptr, int **ipp)
#else
static int RdBinFun(ndim, isMGrid, ngrid, blank, 
                              dimsptr, funcptr, formptr, ipp)
int ndim;
int isMGrid;
int ngrid;
int blank;
long **dimsptr;
p3d_func *funcptr;
ObjPtr *formptr;
int **ipp;
#endif
{
  FILE *inFile;/* the current function file being processed */
  int i, j, n;
  char *name; /* points to the portion of function file name after the last / */
  long *indices;
  int *nvar;  /* number of variables of this function */
  int namelen;
  int ival;
  float fval;

  if (!(inFile = fopen(funcptr->name, "r")))
    {
      fprintf(stderr, "Unable to open %s, aborted\n", funcptr->name);
      return 0;
    }
  if (isMGrid)
    {
      if (1 != fread((void *) &ival, sizeof(int), 1, inFile))
        {
          fprintf(stderr, "Error reading NGRID value in %s\n", funcptr->name);
          fclose(inFile);
          return 0;
        } 
      else if (ival <= 0)
        {
          fprintf(stderr, "%s format error, invalid NGRID value\n", funcptr->name);
          fclose(inFile);
          return 0;
        } 
      else
	{
	  n = ival;
	}
    }
  else
    {
      n = 1;
    }
  /* if ngrid != this file's ngrid then error */
  if (ngrid != n)
    {
      fprintf(stderr, "Error in %s: NGRID value does not match its grid\n",
                                                                 funcptr->name);
      fclose(inFile);
      return 0;
    }
  if (!(indices = (long *) Alloc(ndim * sizeof(long))))
    {
      fprintf(stderr, "Unable to alloc memory, aborted\n");
      fclose(inFile);
      return 0;
    }
  if (!(nvar = (int *) Alloc(ngrid * sizeof(int))))
    {
      fprintf(stderr, "Unable to alloc memory, aborted\n");
      Free(indices);
      fclose(inFile);
      return 0;
    }

  /* if dims does not match the dims read in XYZ file then error */
  for (n = 0; n < ngrid; n++)
    {
      for (i = 0; i < ndim + 1; i++) /* nvar follows dimension spec's */
        {
          if (1 != fread((void *) &ival, sizeof(int), 1, inFile))
            {
              fprintf(stderr, "Error reading dim/nvar value in %s\n", funcptr->name);
              Free(indices);
              Free(nvar);
              fclose(inFile);
              return 0;
            }
          if (i < ndim)  /* the value is dimension */
            {
              if ((long) ival != dimsptr[n][i])
                {
                  fprintf(stderr, "Dimensions mismatch between %s"
                                      " and its XYZ file\n", funcptr->name);
                  Free(indices);
                  Free(nvar);
                  fclose(inFile);
                  return 0;
                }
            }
          else /* i = ndim, the value is of nvar */
            {
              nvar[n] = ival;
            }
        }
    }

  /* determine the dataset name */
  for (i = (int) strlen(funcptr->name) - 1; i >= 0 && funcptr->name[i] != '/'; i--)
     ;
  if (i == -1 || funcptr->name[i] == '/')
    i++;
  name = funcptr->name + i;
  for (namelen = 0; name[namelen] != '\0' && name[namelen] != '.'; namelen++)
     ;

  for (n = 0; n < ngrid; n++)
    {
      long size, count; /* number of data values promised */
      ObjPtr dataset;
      char datasetname[P3D_MAXNAMELEN];
      int whichcomponent;

      strncpy(datasetname, name, namelen);
      datasetname[namelen] = '\0';
      /* number the datasets only if its number of grids is than 1 */
      if (ngrid > 1)
        {
	  /* ordering starts from 1 */
          sprintf(datasetname + namelen, "%d", n + 1);
        }

      size = 1;
      for (i = 0; i < ndim; i++)
         size *= dimsptr[n][i];
      size *= nvar[n];
      for (i = 0; i < ndim; i++)
         {
           indices[i] = 0;
         }

      /* setup datasets */
      dataset = NewStructuredDataset(datasetname, ndim, dimsptr[n],
                                              nvar[n] == 1 ? 0 : nvar[n]);
      if (!SetCurField(FIELD1, dataset))
        {
          fprintf(stderr, "Unable to SetCurField\n");
          Free(indices);
          Free(nvar);
          fclose(inFile);
          return 0;
        }
      whichcomponent = 0;

      /* read in all the func values */
      for (count = 0; count < size; count++)
        {
          int iblank = 1;
     
          if (1 != fread((void *) &fval, sizeof(float), 1, inFile))
            {
	      if (feof(inFile))
		{
                  fprintf(stderr, "Error reading %s, short of data\n", 
			  funcptr->name);
		}
              else
		{
                  fprintf(stderr, "Error reading %s\n", funcptr->name);
                }
              Free(indices);
              Free(nvar);
              fclose(inFile);
              return 0;
            }
          if (blank == P3D_BLANK)
            {
               getIblank(&iblank, ipp, indices, dimsptr[n], ndim, n);
            }
    
          PutFieldComponent(FIELD1, whichcomponent, indices, 
                  iblank == 0 ? (real) missingData : (real) fval);  

          /* adjust the indices and whichcomponent */
          if (funcptr->permutation == P3D_WHOLE)
            {
              for (j = 0; j < ndim; j++)
                {
                  indices[j]++;
                  if (indices[j] == dimsptr[n][j])
                    {
                      indices[j] = 0;
                      continue;
                    }
                  break;
                }
              if (j == ndim) /* all indices are set to zero now */
                {
                  whichcomponent++;
                }
            }
          else  /* P3D_PLANE */
            {
              for (j = 0; j < ndim - 1; j++)
                {
                  indices[j]++;
                  if (indices[j] == dimsptr[n][j])
                    {
                      indices[j] = 0;
                      continue;
                    }
                  break;
                }
              if (j == ndim - 1 && ++whichcomponent == nvar[n])
                {
                  indices[j]++;
                  whichcomponent = 0;
                }
            }
        }     /* for size */

      SetDatasetForm(dataset, formptr[n]);
      RegisterDataset(dataset);
    } /* outer for */

  Free(indices);
  Free(nvar);
  fclose(inFile);
  return ngrid;
} /* RdBinFun */

/*********************binary p3d reading routines, end *********************/

#ifdef HAVE_PROTOTYPES
static void parseDataFiles(ObjPtr reader, int nSets, p3d_set *sets)
#else
static void parseDataFiles(reader, nSets, sets)
ObjPtr reader;
int nSets;
p3d_set *sets;
#endif
{
  int i, j;

/* If no fortran compiler, then unable to handle unformatted p3d data files. */
  for (i = 0; i < nSets; i++)
     {
       p3d_set *curset = sets + i;
       ObjPtr *formptr;           /* an array of ngrid dataforms */
       int **ipp;                 /* points to ngrid arrays of iblank ints */ 
       long **dimsptr;            /* points to ngrid arrays of dims */
       int ngrid;

       if (curset->grid->fileType == P3D_FMT)
	 {
           ngrid = RdFmtXYZ(curset->ndim, curset->isMGrid, &dimsptr, curset->grid, 
			    &formptr, &ipp);
	 }
       else if (curset->grid->fileType == P3D_BIN)
	 {
           ngrid = RdBinXYZ(curset->ndim, curset->isMGrid, &dimsptr, curset->grid, 
			    &formptr, &ipp);
	 }
       else
	 {
#ifdef FORTRAN
           ngrid = RdUnfmtXYZ(curset->ndim, curset->isMGrid, &dimsptr, curset->grid, 
			    &formptr, &ipp);
#else
           fprintf(stderr, "Error in %s, %s data reader is not available\n",
		   curset->grid->name, "Unformatted");
           return;
#endif
         }
       if (!ngrid)
         {
           FileFormatError("parseDataFiles", curset->grid->name);
           return;  /* something wrong with the grid file */
         }

       /* read all the q files of this set and construct datasets for them */
       for (j = 0; j < curset->nQ; j++)
          {
	    int status;

            if ((curset->qlist + j)->fileType == P3D_FMT)
	      {
                status = RdFmtSol(reader, curset->ndim, curset->isMGrid, ngrid,
                  curset->grid->blank, dimsptr, curset->qlist + j, formptr, ipp);
	      }
            else if ((curset->qlist + j)->fileType == P3D_BIN)
	      {
                status = RdBinSol(reader, curset->ndim, curset->isMGrid, ngrid,
                  curset->grid->blank, dimsptr, curset->qlist + j, formptr, ipp);
	      }
            else
	      {
#ifdef FORTRAN
                status = RdUnfmtSol(reader, curset->ndim, curset->isMGrid, ngrid,
                      curset->grid->blank, dimsptr, curset->qlist + j, formptr, ipp);
#else
                fprintf(stderr, "Error in %s, %s data reader is not available\n",
		(curset->qlist + j)->name, "Unformatted");
                return;
#endif
              }
            /* return a value of 0 in status if failed */
            if (!status)
              {
                FileFormatError("Q file", curset->qlist[j].name);
                Free(formptr);
                freedims(ngrid, dimsptr);
                if (curset->grid->blank == P3D_BLANK)
                  freeiblank(ngrid, ipp);
                return;
              }
          }

       for (j = 0; j < curset->nFunc; j++)
          {
	    int status;

            if ((curset->funclist + j)->fileType == P3D_FMT)
	      {
                status = RdFmtFun(curset->ndim, curset->isMGrid, ngrid,
                  curset->grid->blank, dimsptr, curset->funclist + j, formptr, ipp);
	      }
            else if ((curset->funclist + j)->fileType == P3D_BIN)
	      {
                status = RdBinFun(curset->ndim, curset->isMGrid, ngrid,
                  curset->grid->blank, dimsptr, curset->funclist + j, formptr, ipp);
	      }
            else
	      {
#ifdef FORTRAN
                status = RdUnfmtFun(curset->ndim, curset->isMGrid, ngrid,
                  curset->grid->blank, dimsptr, curset->funclist + j, formptr, ipp);
#else
                fprintf(stderr, "Error in %s, %s data reader is not available\n",
		(curset->funclist + j)->name, "Unformatted");
                return;
#endif
              }
            /* return a value of 0 in status if failed */
            if (!status)
              {
                FileFormatError("Function file", curset->funclist[j].name);
                Free(formptr);
                freedims(ngrid, dimsptr);
                if (curset->grid->blank == P3D_BLANK)
                  freeiblank(ngrid, ipp);
                return;
              }
          }

       /* clean up */
       if (curset->grid->blank == P3D_BLANK)
         {
           freeiblank(ngrid, ipp);
           ipp = NULL;
         }
       freedims(ngrid, dimsptr);
       Free(formptr);
       dimsptr = NULL;
       formptr = NULL;
    } /* outer for */
  return;
} /* parseDataFiles */

static ObjPtr ReadP3DFile(reader, p3dname)
ObjPtr reader;
char *p3dname;
{
  int nSets = 0; /* number of p3d sets */
  p3d_set *sets = NULL; /* pointer to an array of p3d sets */ 

  /* parse p3d meta file */
  if (!parseP3DMeta(p3dname, &nSets, &sets))
    {
      FileFormatError("ReadP3DFile", p3dname);
    }
  else
    {
      /* find out if any of the sets using the same XYZ file */
      compress(&nSets, &sets);
      parseDataFiles(reader, nSets, sets);
    }

  cleanup(nSets, sets);
  return NULLOBJ;
} /* ReadP3DFile */

/*******************************************************************************/
/*******************************************************************************/
/*******************************************************************************/

#ifndef RELEASE
/*******************************************************************************/
/********************************AVS READER ************************************/
/*******************************************************************************/
/* AVS field format input file reader */
/* Tzong-Yow Hwu */
/* 11/15/1993 */
/* file extension name .fld */
/* Time dependent file name format is implemented.*/
/* 
   AVS field format input file header format:
   1. The first line must begin with "# AVS".
   2. ndim=<int>
      dim1=<int>
      dim2=<int>
      .
      .
      .
      nspace=<int>
      veclen=<int>
      data={byte, integer, float, double, xdr_integer, xdr_float, xdr_double}
      field={uniform, rectilinear, irregular}
      [min_ext=<float> <float> ...]                 # <nspace> values
      [max_ext=<float> <float> ...]                 # <nspace> values
      label=<string> <string> ...                   # <veclen> values
      min_val=<val> <val> ...                       # <veclen> number of values
      max_val=<val> <val> ...                       # <veclen> number of values
   3. Data are arranged in interlaced, column-major order followed by coordinates in
      noninterlaced, cloumn-major order, with nspace components for irregular grid, 
      in x, y, ... scales format for rectilinear grid, and in minx, maxx, miny, maxy, 
      ... format for uniform grid.  All coordinate are of single precision floating value.
   4. For uniform grid, if the image is cropped, downsized, or interpolated, then a 
      different set of min_ext and max_ext is specified in the header for the uniform dataform.
   5. For uniform or rectilinear grid, <ndim> must be the same as <nspace>.
   6. Number of labels, if specified, must be the same as the value of veclen.
   7. Treat the data as scalar when ndim=1.
*/

/* Required AVS header specifications */
#define AVS_ISNDIM          1     /* for ndim= */
#define AVS_ISDIMS          2     /* for dim?= */
#define AVS_ISNSPACE        4     /* for nspace= */
#define AVS_ISDATA          8     /* for data=   */
#define AVS_ISFIELD        16     /* for field=  */
#define AVS_ISVECLEN       32     /* for veclen = */
/* Optional ones */
#define AVS_ISMINEXT       64     /* for min_ext=  */
#define AVS_ISMAXEXT      128     /* for max_ext=  */
#define AVS_ISLABEL       256     /* for label= */
#define AVS_ISMINVAL      512     /* for min_val= */
#define AVS_ISMAXVAL     1024     /* for max_val= */

/* DATA format */
#define AVS_BYTE          1
#define AVS_INTEGER       2
#define AVS_FLOAT         3
#define AVS_DOUBLE        4
#define AVS_XDR_INTEGER   5
#define AVS_XDR_FLOAT     6
#define AVS_XDR_DOUBLE    7

/* FIELD type */
#define AVS_UNIFORM       1
#define AVS_RECTILINEAR   2
#define AVS_IRREGULAR     3

/* Miscellaneous */
#define AVS_MAXLINELEN    256 
#define AVS_MAXNAMELEN    256 

#define AVS_SKIPBLANKS(s)    while(*(s) == ' ' || *(s) == '\t') (s)++; 

typedef struct avs_spec 
  {
     int mask;           /* indication of what's here and what's not */
     int fieldType;      /* field type  */
     int dataFormat;     /* data format */
     int  rank;          /* for computational dimensions*/
     long *dimensions; 
     int  nSpace;        /* for spatial dimensions */
     real *bounds;  
     int  nComponents;   /* for data components */
     real *minmax;
     char **labels;      /* for labels of each component of data sample */ 
     char name[AVS_MAXNAMELEN];     /* for file name */
  } AVS_Spec;

#ifdef HAVE_PROTOTYPES
static int avs_parseHeader(FILE *inFile, AVS_Spec *pHeader);
static ObjPtr avs_parseData(FILE *inFile, AVS_Spec *pHeader);
static ObjPtr avs_parseCoordinates(FILE *inFile, AVS_Spec *pHeader);
static void avs_cleanUp(AVS_Spec *pHeader);
static real **avs_createScales(int rank, long *dims);
static void avs_freeScales(real **scales, int rank);
static char *avs_fgets(char *s, int n, FILE *stream);
#else
static int avs_parseHeader();
static ObjPtr avs_parseData();
static ObjPtr avs_parseCoordinates();
static void avs_cleanUp();
static real **avs_createScales();
static void avs_freeScales();
static char *avs_fgets();
#endif

#ifdef HAVE_PROTOTYPES
static char *avs_fgets(char *s, int n, FILE *stream)
#else
static char *avs_fgets(s, n, stream)
char *s;
int n;
FILE *stream;
#endif
{
  /* used to read header lines */
  /* special case with form feed \f char */
  int c;

  if ('\f' == (c = fgetc(stream)))
    {
      *s = (char) c;
      *(s + 1) = fgetc(stream);
      *(s + 2) = '\0';
      return s;
    }
  else
    {
      ungetc(c, stream);
      return fgets(s, n, stream);
    }
}

/* return 0 on error */
#ifdef HAVE_PROTOTYPES
static int avs_parseHeader(FILE *inFile, AVS_Spec *pHeader)
#else
static int avs_parseHeader(inFile, pHeader)
FILE *inFile;
AVS_Spec *pHeader;
#endif
{
  int mask;
  char aLine[AVS_MAXLINELEN];
  char trail[AVS_MAXLINELEN];
  char *pc;
  int isError;

  int data, field;
  int nspace, veclen, ndim;
  int whichDim;
  char **labels;
  int nLabels;

  int maxDim;                    /* the max dimensionality of the dimensions specified */
  int nDimLines;                 /* number of dimensions specified */
  unsigned long dimMask;         /* which of the dimensions are specified */
  long *dims;                    /* value of dimensions */

  float *minExts, *maxExts;      /* bounds of space coordinates */
  int nMinExts, nMaxExts;
  float *minVals, *maxVals;      /* bounds of data samples */
  int nMinVals, nMaxVals;
  int whichBound;                /* min_ext, max_ext, min_val or max_val */
  int i;


  while(1)
  {
    /* The first non-blank line of the header must begin with # AVS */
    if (NULL == avs_fgets(aLine, (int) AVS_MAXLINELEN, inFile))
      {
        FileFormatError("avs_parseHeader", "Empty file.");
        return 0;
      }
    pc = aLine;
    AVS_SKIPBLANKS(pc);
    if (*pc == '\n')
      {
	/* blank line */
	continue;
      }
    if (strncmp2("# AVS", pc, 5))
      {
        FileFormatError("avs_parseHeader", "The first non-blank line of AVS file must begin with # AVS");
        return 0;
      }
    break;
  }

  /* initializations */
  mask = 0;
  isError = 0;
  nLabels = 0;
  nDimLines = 0;
  dimMask = 0;
  nMinExts = 0;
  nMaxExts = 0;
  nMinVals = 0;
  nMaxVals = 0;
  pHeader->nSpace = 0;
  pHeader->nComponents = 0;

  while(avs_fgets(aLine, AVS_MAXLINELEN, inFile))
  {
    pc = aLine;
    AVS_SKIPBLANKS(pc);
    if (pc[0] == '#' || pc[0] == '\n')
      {
	/* It's a comment line or a blank line. */
        continue;
      }

    trail[0] = '\0';  /* see if anything trailing after a regular header spec. */

    if ( (0 == strncmp2(pc, "ndim", 4)   && (whichDim = AVS_ISNDIM))   ||
         (0 == strncmp2(pc, "nspace", 6) && (whichDim = AVS_ISNSPACE)) ||
	 (0 == strncmp2(pc, "veclen", 6) && (whichDim = AVS_ISVECLEN)) )
      { 
	int toSkip;
	int iVal;

	/* double specifications */
        if ( ((whichDim == AVS_ISNDIM)   && (mask & AVS_ISNDIM))   ||
             ((whichDim == AVS_ISNSPACE) && (mask & AVS_ISNSPACE)) ||
             ((whichDim == AVS_ISVECLEN) && (mask & AVS_ISVECLEN)) ) 
          {
            FileFormatError("avs_parseHeader, double specification of", aLine); 
	    isError = 1;
            break;
          }

        if (whichDim == AVS_ISNDIM)
	  {
	    toSkip = 4;
	  }
        else if (whichDim == AVS_ISNSPACE)
	  {
	    toSkip = 6;
	  }
        else if (whichDim == AVS_ISVECLEN)
	  {
	    toSkip = 6;
	  }

	/* find the '=' char and read value*/
        pc += toSkip;
	AVS_SKIPBLANKS(pc);
        if (*pc != '=' || 0 >= sscanf(pc + 1, "%d%s", &iVal, trail))
          {
            FileFormatError("avs_parseHeader, Invalid header", aLine);
	    isError = 1;
            break;
          }
        
	/* check the value of iVal */
        if (iVal < 1)
          {
            FileFormatError("avs_parseHeader, Invalid header", aLine);
	    isError = 1;
            break;
          }

        if (whichDim == AVS_ISNDIM)
	  {
	    ndim = iVal;
	  }
        else if (whichDim == AVS_ISNSPACE)
	  {
	    nspace = iVal;
	  }
        else if (whichDim == AVS_ISVECLEN)
	  {
	    veclen = iVal;
	  }
	mask |= whichDim;
      }
    else if (0 == strncmp2(pc, "data", 4))
      {
	char temp[AVS_MAXLINELEN];

        if (mask & AVS_ISDATA)
          {
            FileFormatError("avs_parseHeader, double specification", aLine); 
	    isError = 1;
            break;
          }

	/* find the '=' char and read value*/
        pc += 4;
	AVS_SKIPBLANKS(pc);
        if (*pc != '=' || 0 >= sscanf(pc + 1, "%s%s", temp, trail))
          {
            FileFormatError("avs_parseHeader, Invalid header", aLine);
	    isError = 1;
            break;
          }
        
        if (0 == strncmp2(temp, "byte", 4) || 0 == strncmp2(temp, "xdr_byte", 8))
	  {
	    data = AVS_BYTE;
	  }
        else if (0 == strncmp2(temp, "integer", 7))
	  {
	    data = AVS_INTEGER;
	  }
        else if (0 == strncmp2(temp, "xdr_integer", 11))
	  {
	    data = AVS_XDR_INTEGER;
	  }
        else if (0 == strncmp2(temp, "float", 5))
	  {
	    data = AVS_FLOAT;
	  }
        else if (0 == strncmp2(temp, "xdr_float", 9))
	  {
	    data = AVS_XDR_FLOAT;
	  }
        else if (0 == strncmp2(temp, "double", 6))
	  {
	    data = AVS_DOUBLE;
	  }
        else if (0 == strncmp2(temp, "xdr_double", 10))
	  {
	    data = AVS_XDR_DOUBLE;
	  }
        else
          {
            FileFormatError("avs_parseHeader, Invalid header", aLine);
	    isError = 1;
            break;
          }
	mask |= AVS_ISDATA;
      }
    else if (0 == strncmp2(pc, "field", 5))
      {
	char temp[AVS_MAXLINELEN];

        if (mask & AVS_ISFIELD)
          {
            FileFormatError("avs_parseHeader, double specification", aLine); 
	    isError = 1;
            break;
          }

	/* find the '=' char and read value*/
        pc += 5;
	AVS_SKIPBLANKS(pc);
        if (*pc != '=' || 0 >= sscanf(pc + 1, "%s%s", temp, trail))
          {
            FileFormatError("avs_parseHeader, Invalid header", aLine);
	    isError = 1;
            break;
          }
        
        if (0 == strncmp2(temp, "uniform", 7))
	  {
	    field = AVS_UNIFORM;
	  }
        else if (0 == strncmp2(temp, "rectilinear", 11))
	  {
	    field = AVS_RECTILINEAR;
	  }
        else if (0 == strncmp2(temp, "irregular", 9))
	  {
	    field = AVS_IRREGULAR;
	  }
        else
          {
            FileFormatError("avs_parseHeader, Invalid header", aLine);
	    isError = 1;
            break;
          }
	mask |= AVS_ISFIELD;
      }
    else if (0 == strncmp2(pc, "dim", 3))
      {
	int curDim;
	 
	pc += 3;
        if (!isdigit(*pc) || *pc == '0')
	  {
	    FileFormatError("avs_parseHeader, Invalid header", aLine);
	    isError = 1;
	    break;
	  }
        curDim = *pc - '0';

	/* if there is already a dim?= line */
        if (mask & AVS_ISDIMS)
          {
	    if (dimMask & (unsigned long)(1 << curDim))
	      {
                FileFormatError("avs_parseHeader, double specification", aLine);
		isError = 1;
		break;
	      }
	    dimMask |= (unsigned long)(1 << curDim);
	    if (maxDim < curDim)
	      {
	    	 maxDim = curDim;
		 if (!(dims = (long *) Realloc((void *) dims, sizeof(long) * maxDim)))
		   {
		     OMErr();
		     isError = 1;
		     nDimLines = 0;
		     break;
		   }
	      }
            nDimLines++;
          }
        else
	  {
            if (!(dims = (long *) Alloc(sizeof(long) * curDim)))
	      {
	        OMErr();
	        isError = 1;
	        break;
	      }
	    maxDim = curDim;
	    mask |= AVS_ISDIMS;
	    dimMask |= (unsigned long)(1 << curDim);
	    nDimLines = 1;
	  }

	/* find the '=' char and read value*/
        pc ++;
	AVS_SKIPBLANKS(pc);
        if (*pc != '=' || 0 >= sscanf(pc + 1, "%ld%s", dims + curDim - 1, trail))
          {
            FileFormatError("avs_parseHeader, Invalid header\n", aLine);
	    isError = 1;
            break;
          }
	/* check its value */
        if ( dims[curDim - 1] < 1)
          {
            FileFormatError("avs_parseHeader, Invalid header", aLine);
	    isError = 1;
            break;
          }
      }
    else if ((0 == strncmp2(pc, "min_ext", 7) && (whichBound = AVS_ISMINEXT)) ||
	     (0 == strncmp2(pc, "max_ext", 7) && (whichBound = AVS_ISMAXEXT)) ||
             (0 == strncmp2(pc, "min_val", 7) && (whichBound = AVS_ISMINVAL)) ||
	     (0 == strncmp2(pc, "max_val", 7) && (whichBound = AVS_ISMAXVAL)))
      {
	float fVal, *vals;
	char valStr[AVS_MAXLINELEN];
	int nVals;

        if ((whichBound == AVS_ISMINEXT && mask & AVS_ISMINEXT)  ||
	    (whichBound == AVS_ISMAXEXT && mask & AVS_ISMAXEXT)  ||
            (whichBound == AVS_ISMINVAL && mask & AVS_ISMINVAL)  ||
	    (whichBound == AVS_ISMAXVAL && mask & AVS_ISMAXVAL))
          {
            FileFormatError("avs_parseHeader, double specification", aLine); 
	    isError = 1;
            break;
          }

	/* find the '=' char and read value*/
        pc += 7;
	AVS_SKIPBLANKS(pc);
        if (*pc != '=')
          {
            FileFormatError("avs_parseHeader, Invalid header", aLine);
	    isError = 1;
            break;
          }
	/* skip over the '=' character */
	pc++;

	/* read in values from aLine */
	AVS_SKIPBLANKS(pc);
			/* if there's something immediately trailing a number, then done */
	for (nVals = 0; trail[0] == '\0' && 0 < sscanf(pc, "%s", valStr); nVals++)
	   {
	     if (0 >= sscanf(valStr, "%f%s", &fVal, trail))
	       {
		 strcpy(trail, valStr); /* since it's not a number, keep a copy of it */
		 break;
	       }
             pc += strlen(valStr);
	     AVS_SKIPBLANKS(pc);

	     if (nVals)
	       {
		 if (!(vals = (float *) Realloc(vals, sizeof(float) * (nVals + 1))))
		   {
		     OMErr();
		     isError = 1;
		     break;
		   }
	       }
             else
	       {
		 if (!(vals = (float *) Alloc(sizeof(float))))
		   {
		     OMErr();
		     isError = 1;
		     break;
		   }
	       }
             vals[nVals] = fVal;
	   }
	if (isError)
	  {
	    /* This is an error caused by out of memory. */
	    break;
	  }
        else if (nVals == 0)
	  {
            FileFormatError("avs_parseHeader, Invalid header", aLine);
	    isError = 1;
            break;
	  }

        if (whichBound == AVS_ISMINEXT)
	  {
	    minExts = vals;
	    nMinExts = nVals;
	  }
        else if (whichBound == AVS_ISMAXEXT)
	  {
	    maxExts = vals;
	    nMaxExts = nVals;
	  }
        else if (whichBound == AVS_ISMINVAL)
	  {
	    minVals = vals;
	    nMinVals = nVals;
	  }
        else if (whichBound == AVS_ISMAXVAL)
	  {
	    maxVals = vals;
	    nMaxVals = nVals;
	  }
	mask |= whichBound;
      }
    else if (0 == strncmp2(pc, "label", 5))
      {
	char **ptemps;

        if (mask & AVS_ISLABEL)
          {
            FileFormatError("avs_parseHeader, double specification", aLine); 
	    isError = 1;
            break;
          }

	/* find the '=' char and read value*/
        pc += 5;
	AVS_SKIPBLANKS(pc);
        if (*pc != '=')
          {
            FileFormatError("avs_parseHeader, Invalid header", aLine);
	    isError = 1;
            break;
          }
	/* skip over the '=' character */
	pc++;

	/* read in labels from aLine */
	AVS_SKIPBLANKS(pc);
	for (nLabels = 0; 0 < sscanf(pc, "%s", trail) && trail[0] != '#'; nLabels++, trail[0] = '\0')
	   {
	     int length;

	     length = strlen(trail);
             pc += length;
	     AVS_SKIPBLANKS(pc);
 
	     /* Can not use realloc since if it fails we lose the pointer to ptemps contains. */
	     if (!(ptemps = (char **) Alloc(sizeof(char *) * (nLabels + 1))) ||
		 !(ptemps[nLabels] = (char *) Alloc(sizeof(char) * (length + 1))) )
	       {
	         for (i = 0; i < nLabels; i++)
		    {
		      Free(labels[i]);
		    }
                 if (nLabels) 
		   { 
		     Free(labels); 
		     SAFEFREE(ptemps); /* If allocation of ptemps failed, it will be NULL anyway */
		   }  
		 OMErr();
	         isError = 1;
		 nLabels = 0;
		 break;
	       }
             else
	       {
		 strcpy(ptemps[nLabels], trail); 
	         for (i = 0; i < nLabels; i++)
		    {
		      ptemps[i] = labels[i];
		    }
                 if (nLabels) Free(labels);
                 labels = ptemps;
	       }
	   }
	if (isError)
	  {
	    /* This is an error caused by out of memory. */
	    break;
	  }
        else if (nLabels == 0)
	  {
            FileFormatError("avs_parseHeader, value missing", aLine);
	    isError = 1;
            break;
	  }
	mask |= AVS_ISLABEL;
      }
    else if (0 == strncmp(pc, "unit", 4))
      {
	continue;
      }
    else if (0 == strncmp(pc, "\f\f", 2))
      {
	/* header section terminator */
	break;
      }
    else
      {
        FileFormatError("avs_parseHeader, Invalid header", aLine);
	isError = 1;
        break;
      }

    /* check if anything remaining in the line */
    if (trail[0] == '\0' || trail[0] == '#') 
      {
        /* a comment */
	continue;
      }
    else
      {
        FileFormatError("avs_parseHeader, invalid header", aLine);
        isError = 1;
        break;
      }

  } /* while */

  /* check if all the required headers are present */
  if (!isError)
    {
      if (!(mask & AVS_ISNDIM))
	{
          FileFormatError("avs_parseHeader", "header of ndim= missing");
          isError = 1;
	}
      if (!(mask & AVS_ISDIMS))
	{
          FileFormatError("avs_parseHeader", "headers of dim?= missing");
          isError = 1;
	}
      if (!(mask & AVS_ISNSPACE))
	{
          FileFormatError("avs_parseHeader", "headers of nspace= missing");
          isError = 1;
	}
      if (!(mask & AVS_ISDATA))
	{
          FileFormatError("avs_parseHeader", "headers of data= missing");
          isError = 1;
	}
      if (!(mask & AVS_ISFIELD))
	{
          FileFormatError("avs_parseHeader", "headers of field= missing");
          isError = 1;
	}
      if ( (!(mask & AVS_ISMINEXT) &&  (mask & AVS_ISMAXEXT)) ||
           ( (mask & AVS_ISMINEXT) && !(mask & AVS_ISMAXEXT)) ) 
	{
          FileFormatError("avs_parseHeader", "Either both of min_ext and max_ext or neither must be present");
          isError = 1;
	}
      if ( (!(mask & AVS_ISMINVAL) &&  (mask & AVS_ISMAXVAL)) ||
           ( (mask & AVS_ISMINVAL) && !(mask & AVS_ISMAXVAL)) ) 
	{
          FileFormatError("avs_parseHeader", "Either both of min_val and max_val or neither must be present");
          isError = 1;
	}
    }

  /* check the consistency */
  if (!isError)
    {
      if (ndim != nDimLines ||   /* makes sure there are adequate dim?= lines */
	  ndim != maxDim)   /* makes sure there is no spec of higher dimenionality than necessary */
        {
          FileFormatError("avs_parseHeader", "Number of dim header lines is not the same as the value of ndim");
          isError = 1;
        }
      if ((field == AVS_UNIFORM || field == AVS_RECTILINEAR) && (ndim != nspace))
        {
          FileFormatError("avs_parseHeader", "the value of ndim should be the same as that of nspace for uniform or rectilinear field");
          isError = 1;
        }
      if (mask & AVS_ISMINEXT) 
	{
          if (nMaxExts != nMinExts)
	    {
              FileFormatError("avs_parseHeader", "the number of values of min_ext and max_ext are different");
              isError = 1;
	    }
          else if (nMinExts != nspace)
	    {
              FileFormatError("avs_parseHeader", "the number of values of min_ext and max_ext are different from nspace");
              isError = 1;
	    }
	}
      if (mask & AVS_ISMINVAL) 
	{
          if (nMaxVals != nMinVals)
	    {
              FileFormatError("avs_parseHeader", "the number of values of min_val and max_val are different");
              isError = 1;
	    }
          else if (nMinVals != veclen)
	    {
              FileFormatError("avs_parseHeader", "the number of values of min_val and max_val are different from veclen");
              isError = 1;
	    }
	}
      if (mask & AVS_ISLABEL)
	{
	  if (nLabels != veclen)
	    {
              FileFormatError("avs_parseHeader", "number of labels is different from the value of veclen");
              isError = 1;
	    }
	}
    } 

  /* reassign bounds value */
  if (!isError && mask & AVS_ISMINEXT)
    {
      if (!(pHeader->bounds = (real *) Alloc(sizeof(real) * nspace *2)))
        {
          OMErr();
          isError = 1;
	}
      else
        {
          for (i = 0; i < nspace; i++)
             {
	       pHeader->bounds[2 * i] = (real) minExts[i];
	       pHeader->bounds[2 * i + 1] = (real) maxExts[i];
             }
          Free(minExts);
          Free(maxExts);
	  pHeader->nSpace = nspace;
        }
    }
  if (!isError && mask & AVS_ISMINVAL)
    {
      if (!(pHeader->minmax = (real *) Alloc(sizeof(real) * veclen * 2)))
        {
          OMErr();
          isError = 1;
	}
      else
        {
          for (i = 0; i < veclen; i++)
             {
	       pHeader->minmax[2 * i] = (real) minVals[i];
	       pHeader->minmax[2 * i + 1] = (real) maxVals[i];
             }
          Free(minVals);
          Free(maxVals);
          pHeader->nComponents = veclen;
        }
    }

  if (isError)
    {
      /* something wrong about the file headers */
      if (nDimLines) Free(dims);
      if (nMinExts) Free(minExts);
      if (nMaxExts) Free(maxExts);
      if (nMinVals) Free(minVals);
      if (nMaxVals) Free(maxVals);
      if (nLabels)
	{
	  for (i = 0; i < nLabels; i++)
	     {
	       Free(labels[i]);
	     }
	  Free(labels);
	}
      if (pHeader->nComponents) Free(pHeader->minmax);
      if (pHeader->nSpace) Free(pHeader->bounds);
    }
  else
    {
      /* if everything goes alright */
      pHeader->mask = mask;
      pHeader->rank = ndim;
      pHeader->dimensions = dims;
      if (!(mask & AVS_ISMINVAL)) pHeader->nComponents = veclen;
      if (!(mask & AVS_ISMINEXT)) pHeader->nSpace = nspace;
      pHeader->dataFormat = data;
      pHeader->fieldType = field;
      pHeader->labels = labels;
    }
      
  return !isError;
} /* avs_parseHeader */

#ifdef HAVE_PROTOTYPES
static void avs_cleanUp(AVS_Spec *pHeader)
#else
static void avs_cleanUp(pHeader)
AVS_Spec *pHeader;
#endif
{
  /* free some header storage */
  if (pHeader->mask & AVS_ISDIMS)    Free(pHeader->dimensions);
  if (pHeader->mask & AVS_ISMINEXT)  Free(pHeader->bounds);
  if (pHeader->mask & AVS_ISMINVAL)  Free(pHeader->minmax); 
  if (pHeader->mask & AVS_ISLABEL)
    {
      int i;

      for (i = 0; i < pHeader->nComponents; i++)
	 {
	   Free(pHeader->labels[i]);
	 }
      Free(pHeader->labels);
    }
  return;
}

/* return a dataset */
#ifdef HAVE_PROTOTYPES
static ObjPtr avs_parseData(FILE *inFile, AVS_Spec *pHeader)
#else
static ObjPtr avs_parseData(inFile, pHeader)
FILE *inFile;
AVS_Spec *pHeader;
#endif
{
  ObjPtr retVal;
  int whichComponent;
  long *indices;
  int i;
  size_t dataSize;
  size_t nData;
  void *data;

  /* create a dataset as the retVal */
  if (NULLOBJ == (retVal = NewStructuredDataset( pHeader->name, pHeader->rank, pHeader->dimensions, 
      (pHeader->nComponents == 1 ? 0 : pHeader->nComponents) )))
    {
      return NULLOBJ;
    }
  if (!SetCurField(FIELD1, retVal)) 
    {
      return NULLOBJ;
    }

  /* for traversing the data values */
  if (!(indices = (long *) Alloc(sizeof(long) * pHeader->rank)))
    {
      OMErr();
      return NULLOBJ;
    }
  for (i = 0; i < pHeader->rank; i++)
     {
       indices[i] = 0L;
     }

  /* what is the size per data */
  switch (pHeader->dataFormat) 
     {
       case AVS_BYTE:
	   dataSize = sizeof(char);
	   break;
       case AVS_INTEGER:  case AVS_XDR_INTEGER:
	   dataSize = sizeof(int);
	   break;
       case AVS_FLOAT:    case AVS_XDR_FLOAT:
	   dataSize = sizeof(float);
	   break;
       case AVS_DOUBLE:   case AVS_XDR_DOUBLE:
	   dataSize = sizeof(double);
	   break;
       default:
	   break;
     }

  nData = 1;
  for (i = 0; i < pHeader->rank; i++)
     {
       nData *= pHeader->dimensions[i]; 
     }
  nData *= pHeader->nComponents;

  /* if there's still enough memory then read all the data at one time */
  if (data = (void *) Alloc(dataSize * nData))
    {
      long whichData;

      if (nData != fread(data, dataSize, nData, inFile))
        {
          FileFormatError("avs_parseData", "Error reading binary data values");
          Free(indices);
          Free(data);
          return NULLOBJ;
        }

      whichData = 0;
      while(1)
	{
          for (whichComponent = 0; whichComponent < pHeader->nComponents; whichComponent++)
	     {
               switch (pHeader->dataFormat) 
                  {
                    case AVS_BYTE:
                        PutFieldComponent(FIELD1, whichComponent, indices, 
					  (real) *((char *) data + whichData));
	                break;
                    case AVS_INTEGER:  case AVS_XDR_INTEGER:
                        PutFieldComponent(FIELD1, whichComponent, indices, 
					  (real) *((int *) data + whichData));
	                break;
                    case AVS_FLOAT:    case AVS_XDR_FLOAT:
                        PutFieldComponent(FIELD1, whichComponent, indices, 
					  (real) *((float *) data + whichData));
	                break;
                    case AVS_DOUBLE:   case AVS_XDR_DOUBLE:
                        PutFieldComponent(FIELD1, whichComponent, indices, 
					  (real) *((double *) data + whichData));
	                break;
                    default:
	                break;
                  }
               whichData++;
	     }
	  for (i = 0; i < pHeader->rank; i++)
	     {
	       indices[i]++;
	       if (indices[i] == pHeader->dimensions[i])
		 {
		   indices[i] = 0;
		   continue;
		 }
               break;
	     }
          if (i == pHeader->rank)
	    {
	      /* all indices are traversed */
	      break;
	    }
	}  /* while loop */
    } /* if */
  else if (data = (void *) Alloc(dataSize))
      /* read in one data value at a time */
    {
      /* assign the data to the dataset while reading one data value at a time */
      while(1)
	{
          for (whichComponent = 0; whichComponent < pHeader->nComponents; whichComponent++)
	     {
               if (1 != fread(data, dataSize, (size_t) 1, inFile))
                 {
                   FileFormatError("avs_parseData", "Error reading binary data values");
                   Free(indices);
                   Free(data);
                   return NULLOBJ;
                 }
               switch (pHeader->dataFormat) 
                  {
                    case AVS_BYTE:
                        PutFieldComponent(FIELD1, whichComponent, indices, 
					  (real) *((char *) data));
	                break;
                    case AVS_INTEGER:  case AVS_XDR_INTEGER:
                        PutFieldComponent(FIELD1, whichComponent, indices, 
					  (real) *((int *) data));
	                break;
                    case AVS_FLOAT:    case AVS_XDR_FLOAT:
                        PutFieldComponent(FIELD1, whichComponent, indices, 
					  (real) *((float *) data));
	                break;
                    case AVS_DOUBLE:   case AVS_XDR_DOUBLE:
                        PutFieldComponent(FIELD1, whichComponent, indices, 
					  (real) *((double *) data));
	                break;
                    default:
	                break;
                  }
	     }
	  for (i = 0; i < pHeader->rank; i++)
	     {
	       indices[i]++;
	       if (indices[i] == pHeader->dimensions[i])
		 {
		   indices[i] = 0;
		   continue;
		 }
               break;
	     }
          if (i == pHeader->rank)
	    {
	      /* all indices are traversed */
	      break;
	    }
	}  /* while loop */
    }
  /* out of memory */
  else
    {
      OMErr();
      Free(indices);
      return NULLOBJ;
    }

  Free(data);
  Free(indices);

  return retVal; 
}/* avs_parseData */

/* return a dataForm ObjPtr*/
#ifdef HAVE_PROTOTYPES
static ObjPtr avs_parseCoordinates(FILE *inFile, AVS_Spec *pHeader)
#else
static ObjPtr avs_parseCoordinates(inFile, pHeader)
FILE *inFile;
AVS_Spec *pHeader;
#endif
{
  ObjPtr retVal;
  int i, j;

  switch(pHeader->fieldType)
    {
      case AVS_UNIFORM:
	  {
	    real *realBounds;
	    float *floatBounds;
            int nBounds;
	   
	    /* for UNIFORM field the pHeader->nSpace equals to pHeader->rank */
	    if (!(realBounds = (real *) Alloc(sizeof(real) * pHeader->rank * 2)))
	      {
		OMErr();
		return NULLOBJ;
	      }

	    /* default bounds, HWU */
            for(i = 0; i < pHeader->rank; i++)
	       {
		 realBounds[2 * i] = 0.0;
		 realBounds[2 * i + 1] = 1.0;
	       }

	    /* read in all the bounds at a time */
            if (floatBounds = (float *) Alloc(sizeof(float) * pHeader->rank * 2))
	      {
		/*
		if (pHeader->rank * 2 != fread(floatBounds, sizeof(float), (size_t) pHeader->rank * 2, inFile))
                */
		nBounds = fread(floatBounds, sizeof(float), (size_t) pHeader->rank * 2, inFile);
		if (nBounds && pHeader->rank * 2 != nBounds)
		  {
                     FileFormatError("avs_parseCoordinates", "Error: number of bounds provided is inconsistent with the rank.");
		     Free(realBounds);
		     Free(floatBounds);
                     return NULLOBJ;
		  }
                for (i = 0; i < 2 * pHeader->rank; i++)
	           {
                     realBounds[i] = (real) floatBounds[i];
	           }
	      }
	    /* read in bounds one at a time */
            else if (floatBounds = (float *) Alloc(sizeof(float)))
	      {
                for (nBounds = 0; nBounds < 2 * pHeader->rank; nBounds++)
	           {
                     if (1 != fread(floatBounds, sizeof(float), (size_t) 1, inFile))
	               {
			 if (nBounds)
			   {
                             FileFormatError("avs_parseCoordinates", "Error reading coordinate values");
		             Free(realBounds);
			     Free(floatBounds);
                             return NULLOBJ;
			   }
                         else
			   {
			     break;
                           }
	               }
                     realBounds[i] = (real) *floatBounds;
	           }
              }
            else
	      {
		OMErr();
		Free(realBounds);
		return NULLOBJ;
	      }

	    /* see if the bounds is the same as the ones specified in the header if any */
	    if ((pHeader->mask & AVS_ISMINEXT) && nBounds)
	      {
		int equalBounds = 1;

		for (i = 0; equalBounds && i < 2 * pHeader->rank; i++)
		   {
		     equalBounds = (realBounds[i] == pHeader->bounds[i]);
		   }
                if (equalBounds)
		  {
		    retVal = NewRegularDataForm(pHeader->name, pHeader->rank, pHeader->dimensions, realBounds);
		  } /* if equalBounds */
                else
		  {
		    /* use scales calculated from original bounds */
		    real **scales;
		    real *elements;
		    ObjPtr boundObj;
                   
                    if (NULLOBJ == (boundObj = NewRealArray(1, (long) pHeader->rank * 2L)))
		      {
			Free(realBounds);
			return NULLOBJ;
                      }
		    if (!(scales = (real **) avs_createScales(pHeader->rank, pHeader->dimensions)))
		      {
			Free(realBounds);
			return NULLOBJ;
		      }
                    for (i = 0; i < pHeader->rank; i++)
		       {
		         real interval;
			 real min, max;

                         min = pHeader->bounds[2 * i];
			 max = pHeader->bounds[2 * i + 1];
			 if (pHeader->dimensions[i] > 1)
			   {
			     interval = (max - min) / (real) (pHeader->dimensions[i] - 1);
			   }
                         for (j = 0; j < pHeader->dimensions[i] - 1; j++)
			    {
			      scales[i][j] = min + (real) j * interval;
			    }
                         scales[i][j] = max;
		       }
                    retVal = NewSeparableDataForm(pHeader->name, pHeader->rank, pHeader->dimensions, scales);
		    /* reset the dataForm's bounds */
		    elements = ELEMENTS(boundObj);
		    for (i = 0; i < pHeader->rank * 2; i++)
		       {
			 elements[i] = realBounds[i];
		       }
	            SetVar(retVal, BOUNDS, boundObj);
		    avs_freeScales(scales, pHeader->rank);
		  } /* else of if equalBounds */
	      } /* if AVS_ISMINEXT */
            else if (nBounds)
	      {
		retVal = NewRegularDataForm(pHeader->name, pHeader->rank, pHeader->dimensions, realBounds);
	      } /* if nBounds */
            else if (pHeader->mask & AVS_ISMINEXT)
	      {
		retVal = NewRegularDataForm(pHeader->name, pHeader->rank, pHeader->dimensions, pHeader->bounds);
	      }
            else
	      {
		/* use 0, 1, 0, 1, ... as bounds */
	        /* default bounds, HWU */
		retVal = NewRegularDataForm(pHeader->name, pHeader->rank, pHeader->dimensions, realBounds);
		/* or treat as error */
		/* 
                FileFormatError("avs_parseCoordinates", "Coordinate values missing");
                Free(realBounds);
		return NULLOBJ;
		*/
	      }
            Free(realBounds);
            Free(floatBounds);
	  }
	  break;
      case AVS_RECTILINEAR:
	  {
	    real **scales;
	    float *coords;
	    long nCoords;
                  
	    /* for rectilinear dataform, pHeader->rank equals to pHeader->nSpace */
	    if (!(scales = (real **) avs_createScales(pHeader->rank, pHeader->dimensions)))
	      {
		return NULLOBJ;
	      }
	    nCoords = 0L;
            for (i = 0; i < pHeader->rank; i++)
	       {
		 nCoords += pHeader->dimensions[i];
	       }

	    /* read in all scale values at one time */
            if (coords = (float *) Alloc(sizeof(float) * nCoords))
	      {
		long sofar;

		if (nCoords != fread(coords, sizeof(float), (size_t) nCoords, inFile))
	          {
                    FileFormatError("avs_parseCoordinates", "Error reading coordinate values");
		    avs_freeScales(scales, pHeader->rank);
		    return NULLOBJ;
		  }
                sofar = 0;
                for (i = 0; i < pHeader->rank; i++)
	           {
                     for (j = 0; j < pHeader->dimensions[i]; j++)
		        {
                          scales[i][j] = coords[j + sofar];
		        }
                     sofar += pHeader->dimensions[i];
	           }
	      }
	    /* read in scale values one at a time */
            else if (coords = (float *) Alloc(sizeof(float)))
	      {
                for (i = 0; i < pHeader->rank; i++)
	           {
                     for (j = 0; j < pHeader->dimensions[i]; j++)
		        {
		          if (1 != fread(coords, sizeof(float), (size_t) 1, inFile))
			    {
                              FileFormatError("avs_parseCoordinates", "Error reading coordinate values");
			      avs_freeScales(scales, pHeader->rank);
		              return NULLOBJ;
			    }
                          scales[i][j] = (real) *coords;
		        }
	           }
	      }
            else  /* out of memory */
	      {
                avs_freeScales(scales, pHeader->rank);
		OMErr();
		return NULLOBJ;
	      }
            retVal = NewSeparableDataForm(pHeader->name, pHeader->rank, pHeader->dimensions, scales);
	    Free(coords);
	    avs_freeScales(scales, pHeader->rank);
	  }
	  break;
      case AVS_IRREGULAR:
	  {
	    ObjPtr formDataset;
	    real *bounds;
            long nCoords;
	    int whichCoord;
            float *coords;
            long *indices;
            int i;

	    if (NULLOBJ == (formDataset = NewStructuredDataset(pHeader->name, pHeader->rank, pHeader->dimensions, pHeader->nSpace)))
	      {
		return NULLOBJ;
	      }
            if (!SetCurField(FIELD1, formDataset)) 
              {
                return NULLOBJ;
              }

            /* for traversing the data values */
            if (!(indices = (long *) Alloc(sizeof(long) * pHeader->rank)))
              {
                OMErr();
                return NULLOBJ;
              }
            for (i = 0; i < pHeader->rank; i++)
               {
                 indices[i] = 0L;
               }

	    /* see if we need to calculate bounds */
            if (!(pHeader->mask & AVS_ISMINEXT))
	      {
		if (!(bounds = (real *)Alloc(sizeof(real) * pHeader->nSpace * 2)))
		  {
                    OMErr();
		    Free(indices);
                    return NULLOBJ;
		  }
	      }
            else
	      {
		bounds = pHeader->bounds;
	      }

            nCoords = 1;
            for (i = 0; i < pHeader->rank; i++)
               {
                 nCoords *= pHeader->dimensions[i]; 
               }
            nCoords *= pHeader->nComponents;

            /* if there's still enough memory then read all the coordinates at one time */
            if (coords = (float *) Alloc(sizeof(float) * nCoords))
              {
                long sofar;

                if (nCoords != fread(coords, sizeof(float), nCoords, inFile))
                  {
                    FileFormatError("avs_parseCoordinates", "Error reading coordinates");
		    if (!(pHeader->mask & AVS_ISMINEXT))  Free(bounds);
                    Free(indices);
                    Free(coords);
                    return NULLOBJ;
                  }
                sofar = 0;
                for (whichCoord = 0; whichCoord < pHeader->nSpace; whichCoord++)
	           {
                     /* need to calculate grid bounds if no bounds is provided */
		     if (!(pHeader->mask & AVS_ISMINEXT))
		       {
		         bounds[2 * whichCoord] = bounds[2 * whichCoord + 1] = (real) coords[sofar];
		       }

		     while(1)
		       {
                         /* need to calculate grid bounds if no bounds is provided */
		         if (!(pHeader->mask & AVS_ISMINEXT))
		           {
		             if ((real) coords[sofar] > bounds[2 * whichCoord + 1])
		               {
		                 bounds[2 * whichCoord + 1] = (real) coords[sofar];
		               }
                             else if ((real) coords[sofar] < bounds[2 * whichCoord])
		               {
		                 bounds[2 * whichCoord] = (real) coords[sofar];
		               }
		           }
                         PutFieldComponent(FIELD1, whichCoord, indices, (real) *(coords + sofar));
		         sofar++;
	                 for (i = 0; i < pHeader->rank; i++)
	                    {
	                      indices[i]++;
	                      if (indices[i] == pHeader->dimensions[i])
		                {
		                  indices[i] = 0;
		                  continue;
		                }
                              break;
	                    }
                         if (i == pHeader->rank)
	                   {
	                     /* all indices are traversed */
	                     break;
	                   }
		       } /* while */
                   }
              } /* if */
            else if (coords = (float *) Alloc(sizeof(float)))
            /* read in one data value at a time */
              {
		int curCoord= -1;

                for (whichCoord = 0; whichCoord < pHeader->nSpace; whichCoord++)
	           {
		     while(1)
		       {
                         if (1 != fread(coords, sizeof(float), (size_t) 1, inFile))
                           {
                             FileFormatError("avs_parseCoordinates", "Error reading coordinates");
		             if (!(pHeader->mask & AVS_ISMINEXT))  Free(bounds);
                             Free(indices);
                             Free(coords);
                             return NULLOBJ;
                           }
                         PutFieldComponent(FIELD1, whichCoord, indices, (real) *coords);

                         /* need to calculate grid bounds if no bounds is provided */
			 if (!(pHeader->mask & AVS_ISMINEXT))
			   {
			     if (curCoord == whichCoord)
			       {
			         if ((real) *coords > bounds[2 * whichCoord + 1])
			           {
				     bounds[2 * whichCoord + 1] = (real ) *coords;
			           }
                                 else if ((real) *coords < bounds[2 * whichCoord])
			           {
				     bounds[2 * whichCoord] = (real ) *coords;
			           }
			       }
                             else
			       {
			         bounds[2 * whichCoord] = bounds[2 * whichCoord + 1] = (real) *coords;
				 curCoord = whichCoord;
			       }
			   }
	                 for (i = 0; i < pHeader->rank; i++)
	                    {
	                      indices[i]++;
	                      if (indices[i] == pHeader->dimensions[i])
		                {
		                  indices[i] = 0;
		                  continue;
		                }
                              break;
	                    }
                         if (i == pHeader->rank)
	                   {
	                     /* all indices are traversed */
	                     break;
	                   }
		       } /* while */
                   }
              }
            /* out of memory */
            else
              {
                OMErr();
		if (!(pHeader->mask & AVS_ISMINEXT))  Free(bounds);
                Free(indices);
                return NULLOBJ;
              }

	    retVal = NewCurvilinearDataForm(pHeader->name, pHeader->rank, pHeader->dimensions, bounds, formDataset);
            Free(coords);
            Free(indices);
	    if (!(pHeader->mask & AVS_ISMINEXT)) Free(bounds);
	  }
	  break;
      default:
	  break;
    }
  return retVal; 
}/* avs_parseCoordinates */

#ifdef HAVE_PROTOTYPES
static real **avs_createScales(int rank, long *dims)
#else
static real **avs_createScales(rank, dims)
int rank;
long *dims;
#endif
{
  int i;
  real **retVal;

  if (!(retVal = (real **) Alloc(sizeof(real *) * rank)))
    {
      return (real **) NULL;
    }
  for (i = 0; i < rank; i++)
     {
       if (!(retVal[i] = (real *) Alloc(sizeof(real) * dims[i])))
	 {
	   avs_freeScales(retVal, i);
	   return (real **) NULL;
	 }
     }
  return retVal;
}

#ifdef HAVE_PROTOTYPES
static void avs_freeScales(real **scales, int rank)
#else
static void avs_freeScales(scales, rank);
real **scales;
int rank;
#endif
{
  int i;

  for (i = 0; i < rank; i++)
     {
       Free(scales[i]);
     }
  Free(scales);
}

static ObjPtr ReadAVSFile(reader, avsname)
ObjPtr reader;
char *avsname;
{
  struct avs_spec header;
  FILE *inFile;
  int i;
  ObjPtr dataset, dataForm;

  if (!(inFile = fopen(avsname, "r")))
    {
      Error("ReadAVSFile", OPENFILEERROR, avsname);
      return NULLOBJ;
    }

  /* get dataset name */
  {
    ObjPtr datasetName;

    if (NULLOBJ == (datasetName = MakeDatasetName(avsname)))
      {
	return NULLOBJ;
      }
    strcpy(header.name, GetString(datasetName));
  }

  if (avs_parseHeader(inFile, &header))
    {
      if (NULLOBJ != (dataset = avs_parseData(inFile, &header)) &&
          NULLOBJ != (dataForm = avs_parseCoordinates(inFile, &header)) )
        {
          /* HWU? Do we always need to set the minmax if it's available */
          if (header.mask & AVS_ISMINVAL)
	    {
	      ObjPtr minmax;
    
	      /* set the dataset's MINMAX if it is provided */
	      if (NULLOBJ == (minmax = NewRealArray(1, (long) header.nComponents * 2L)))
	        {
	          return NULLOBJ;
	        }
	      CArray2Array(minmax, header.minmax); 
	      SetVar(dataset, MINMAX, minmax);
	    }
          SetDatasetForm(dataset, dataForm);
          RegisterDataset(dataset);
        }
      else
        {
          fprintf(stderr, "Error reading %s\n", avsname);
        }
      avs_cleanUp(&header);
    }
  else
    {
      fprintf(stderr, "Error reading header of %s\n", avsname);
    }
  fclose(inFile);
  return NULLOBJ;
} /* ReadAVSFile */
/*******************************************************************************/
/*******************************************************************************/
/*******************************************************************************/
#endif

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

  /* PLOT3D reader */
  fileReader = NewFileReader("P3D");
  SetVar(fileReader, EXTENSION, NewString("p3d"));
  SetVar(fileReader, TIMEDDATASETS, NewInt(0));
  SetVar(fileReader, HELPSTRING, 
	 NewString("This file reader reads scientific datasets using \
the PLOT3D data format.  No calculation is performed on the data stored in\
solution files.  Consequently, only the raw data stored in solution files, \
i.e, density, momemtum vector and stagnation energy per unit volume, are visualized.\
Fortunatelly, the reader is capable of reading PLOT3D function file.\
To visualize other quantities derived from the data stored in solution files, you \
only need to pre-compute the data and store the result in a funtion file of \
PLOT3D format.")); 
  SetMethod(fileReader, READALL, ReadP3DFile);
  ApplySavedSettings(fileReader);

  /* STF reader */
  fileReader = NewFileReader("STF");
  SetVar(fileReader, EXTENSION, NewString("stf"));
  SetVar(fileReader, TIMEDDATASETS, NewInt(0));
  SetVar(fileReader, HELPSTRING, 
	 NewString("This file reader reads scientific datasets using \
the Simple Text Format created by Eric Pepke at Supercomputer Computation Research \
Institute of Florida State University.  Refer to SciAn manuals for detailed \
information on the STF format."));
  SetMethod(fileReader, READALL, ReadSTFFile);
  ApplySavedSettings(fileReader);

#ifndef RELEASE
  /* AVS reader */
  fileReader = NewFileReader("AVS");
  SetVar(fileReader, EXTENSION, NewString("fld"));
  SetVar(fileReader, HELPSTRING, 
	 NewString("This file reader reads scientific datasets using \
the AVS read field data format."));
  SetMethod(fileReader, READALL, ReadAVSFile);
  ApplySavedSettings(fileReader);
#endif
}

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