/***********************************************************************/
/* This program demonstrates how to read a field from a DAO HDF-EOS    */
/* product using the HDF-EOS library.  It will take a file name and    */
/* field name on the command line, read the first time of the given    */
/* field, calculate an average of that time and print the average.     */
/*                                                                     */
/*             usage:  avg <file name> <field name>                    */
/*                                                                     */
/*                                                     Rob Lucchesi    */
/*                                            rob.lucchesi@nasa.gov    */
/***********************************************************************/

#include "hdf.h"
#include "mfhdf.h"
#include <stdio.h>


main(int argc,char *argv[]) {
   
int32 sd_id, sds_id, status;
int32 sds_index;
int32 start[4], edges[4], stride[4];
char *fname, *vname;
float32 *data_array;
float32 avg, sum;
int32 i;
int32 file_id, gd_id;
int32 xdim, ydim, zdim, len;

if (argc != 3) {
   printf("Usage: avg <filename> <field> \n");
   exit (-1);
}

fname = argv[1];
vname = argv[2];

/* Open the file (read-only) */

file_id = GDopen (fname, DFACC_RDONLY);
if (file_id < 0) {
   printf ("Could not open %s\n",fname);
   exit(-1);
}

/* Attach to the EOS grid contained within the file. */
/* The DAO uses the generic name "EOSGRID" for the grid in all products. */

gd_id = GDattach (file_id,"EOSGRID");  
if (gd_id < 0) {
   printf ("Could not open %s\n",fname);
   exit(-1);
}

/* Set positioning arrays to read the entire field at the first time. */

start[0] = 0;
start[1] = 0;
start[2] = 0;
start[3] = 0;

stride[0] = 1;
stride[1] = 1;
stride[2] = 1;
stride[3] = 1;

edges[0] = 1;
edges[1] = zdim;
edges[2] = ydim;
edges[3] = xdim;

/*
   In this program, we read the entire field.  By manipulating the start 
   and edges arrays, it is possible to read a subset of the entire array.  
   For example, to read a 3D section defined by x=100,225; y=50,150; 
   z=15,16 you would set the start and edges arrays to the following:

   start[0] = 0;    time start location
   start[1] = 15;   z-dim start location
   start[2] = 50;   y-dim start location
   start[3] = 100;  x-dim start location

   edges[0] = 1;    time length
   edges[1] = 2;    z-dim length
   edges[2] = 100;  y-dim length
   edges[3] = 125;  x-dim length

*/

len = xdim*ydim*zdim;
data_array = (float32 *) malloc(len);

/* Read the data into data_array */

status = GDreadfield (gd_id, vname, start, stride, edges, data_array);
printf ("Read status=%d\n",status);

/* Calculate and print the average */

sum=0.0;
for (i=0; i<len; i++) sum += data_array[i];
avg = sum/(float32)len;

printf ("Average of %s in 3 dimensions is=%f\n",vname,avg);

/* Close file. */

status = GDdetach (gd_id);
status = GDclose (file_id);

}