#include <stdio.h>
#include <avs/avs.h>
#include <avs/port.h>
#include <avs/field.h>
/* John's Includes */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
/* Local Includes */
#include <IO.hh>
#include <IEEEIO.hh>
#ifdef WITH_HDF
#include <HDFIO.hh>
#endif
// #include <Reader.hh> // an even higher level interface, but must skip for now

extern "C" {
int  readIEEE(AVSfield **out,char *Filename,int RecordNumber,
	   int SetExtents,char *SetCoordinates);
}

int readIEEE(AVSfield **out,char *Filename,int RecordNumber,
	     int SetExtents,char *SetCoordinates){
  static IObase *infile=0;
  // Handle case where Filename has changed
  if(AVSparameter_changed("Filename") || !infile){
    // Filename has changed, so close any currently open data files
    // (this is done by delete'ing here.)
    if(infile) delete infile;
    infile=0;
    // Just in case the new file is bogus, hide the dial widgets so 
    // that they can't be messed with if the data isn't valid.
    AVSparameter_visible("RecordNumber",0);
#ifdef WITH_HDF
    if(strstr(strchr(Filename,','),"hdf"))
      infile = new HDFIO(Filename,IObase::Read);
    else
#endif
      infile = new IEEEIO(Filename,IObase::Read); // open the datafile

    if(!infile->isValid()){ // if the file is not IEEEIO, exit
      delete infile;
      infile=0;
      fprintf(stderr,"%s is not an IEEEIO file",Filename);
      return 0;
    }
    if(infile->nDatasets()>1){// don't set widget dial for min=max=0
      AVSmodify_parameter("RecordNumber",AVS_VALUE|AVS_MAXVAL,
			  0,0,infile->nDatasets()-1);
      AVSparameter_visible("RecordNumber",1);
    } // so dial widget remains hidden if dataset contains 1 or fewer datasets
    RecordNumber=0;
  }
  // Failsafe sanity check (things are really hosed if this check fails)
  if(!infile) {
    puts("Not an IEEEIO file");
    AVSparameter_visible("RecordNumber",0);
    return 0;
  }
  
  // now do memory allocation for AVS
  char AllocationDescription[64];
  char *typname;
  IObase::DataType datatype;
  int rank,dims[5];
  int nelements;

  infile->seek(RecordNumber); // seek to desired record
  infile->readInfo(datatype,rank,dims); // get info about the dataset
  // we are building the AVS memory description text string to do the allocation.
  switch(datatype){
  case IObase::Float32:
    puts("Data Type float");
    typname="float";
    break;
  case IObase::Float64:
    puts("Data Type double");
    typname="double";
    break;
  case IObase::Int32:
    typname="integer";
    break;
  case IObase::uChar:
  case IObase::Int8:
    typname="byte";
    break;
  }

  if(!strcmp("Set Edges",SetCoordinates))
    sprintf(AllocationDescription,"field %uD %u-space scalar rectilinear %s",
	    rank,rank,typname); // put it all together into the allocation string
  else
    sprintf(AllocationDescription,"field %uD %u-space scalar uniform %s",
	    rank,rank,typname); // put it all together into the allocation string

  //printf("AVS allocating %s\n",AllocationDescription);
  //printf("With dims[0]=%d dims[1]=%d dims[2]=%d\n",dims[0],dims[1],dims[2]);
  if(*out) AVSfield_free(*out); // free old memory if its still around
  *out = (AVSfield *)AVSdata_alloc(AllocationDescription,dims); // and allocate!
  if(!*out){
    fprintf(stderr,"Allocation of %s failed\n",
	    AllocationDescription);
    return 0; // allocation failed
  }

  /*---
    Note: So if this read doesn't work, you'll need to use a switch statement
    with exact casts to the proper datatypes. 
    I skipped using the switch statement because it looks really stupid
    but its necessary for some machines which align vectors
    differently depending on the size of the type the vector
    points to. (yes folks... this is what the "align" keyword is about)
    Hopefully we won't use such poorly designed machines.

    This need not be done for the parameters though since that field
    of the avsfld is always float.
  ---*/
  AVSfield_float *fld=(AVSfield_float*)(*out);
  puts("Preparing to read data");
  infile->read(fld->data); // read the data
  puts("read data");
  int minindex,maxindex;
  /*---
    In this section we are just setting the extents and edge points information
    for AVS (the physical size of the dataset in floating point coordinates).
    There are just two Attribute pairs that it will recognize to do this.

    1) min_ext,max_ext
       These are vectors of size <rank of dataset> which contain the
       coordinates of the minimum and maximum corners of the bounding box.
    2) origin,delta
       The GR people seem to like this.  Its the origin (same as the min_ext)
       of the dataset and then the delta (or "dx") between grid points in the
       dataset.

    It is now able to handle these parameters as Float32 or Float64 data.
    Integer is ignored.  The dataset will simply be ranged from 0-dim[x]
    if the above Attributes are not found.
  ---*/
  {
    float min_ext[5],max_ext[5];

    if((minindex=infile->readAttributeInfo("min_ext",datatype,nelements))>0 &&
       (maxindex=infile->readAttributeInfo("max_ext",datatype,nelements))>0 &&
       (datatype==IObase::Float32 || datatype==IObase::Float64)) {
      if(datatype==IObase::Float32){
	infile->readAttribute(maxindex,min_ext);
	infile->readAttribute(minindex,max_ext);
      }
      else if(datatype==IObase::Float64){
	double maxext[5],minext[5];
	infile->readAttribute(maxindex,maxext);
	infile->readAttribute(minindex,minext);
	for(int i=0;i<rank;i++){
	  min_ext[i]=(float)(minext[i]);
	  max_ext[i]=(float)(maxext[i]);
	}
      }
      
    } //***** OK, if we cant find extents, then try origin & delta
    else if((minindex=infile->readAttributeInfo("origin",datatype,nelements))>0 &&
	    (maxindex=infile->readAttributeInfo("delta",datatype,nelements))>0 &&
	    (datatype==IObase::Float32 || datatype==IObase::Float64))  {
      if(datatype==IObase::Float32){
	infile->readAttribute(maxindex,fld->min_extent);
	infile->readAttribute(minindex,fld->max_extent);
	for(int i=0;i<rank;i++){
	  max_ext[i]*=(float)dims[i];
	  max_ext[i]+=min_ext[i];
	}
      }
      else if(datatype==IObase::Float64){
	double maxext[5],minext[5];
	infile->readAttribute(maxindex,maxext);
	infile->readAttribute(minindex,minext);
	for(int i=0;i<rank;i++){
	  maxext[i]*=(double)dims[i];
	  maxext[i]+=minext[i];
	  min_ext[i]=(float)(minext[i]);
	  max_ext[i]=(float)(maxext[i]);
	}
      } 
    }
    else { // No information about edges.
      for(int i=0;i<rank;i++){
	min_ext[i]=0.0;
	max_ext[i]=(float)(dims[i]-1);
      }
    }
    if(SetExtents){
      for(int i=0;i<rank;i++){
	fld->min_extent[i]=min_ext[i];
	fld->max_extent[i]=max_ext[i];
      }
    }
    if(!strcmp("Set Min/Max",SetCoordinates)){
      for(int i=0;i<rank;i++){
	fld->points[i*2]=min_ext[i];
	fld->points[i*2 + 1]=max_ext[i];
      }
    }
    else if(!strcmp("Set Edges",SetCoordinates)){
      for(int index=0,i=0;i<rank;i++){
	double dx = (max_ext[i]-min_ext[i])/((float)(fld->dimensions[i]));
	printf("SetEdges[%u]: dx=%lf\n",i,dx);
	for(int j=0;j<fld->dimensions[i];j++,index++){
	  fld->points[index] = min_ext[i] + ((double)j)*dx;
	}
      }
    }
  }
  return 1; // for now
}