/***********************************************************************/
/* This program demonstrates how to read a field from a DAO HDF-EOS    */
/* product using the HDF library (HDF-EOS not required).  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, 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) */

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

/* Find the index and ID of the SDS for the given variable name. */

sds_index = SDnametoindex (sd_id, vname);
if (sds_index < 0) {
   printf ("Could not find %s\n",vname);
   exit(-1);
}

sds_id = SDselect (sd_id,sds_index);

/* 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

*/

/* Read the data into data_array */

status = SDreaddata (sds_id, start, stride, edges, (VOIDP) 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 = SDendaccess (sds_id);
status = SDend (sd_id);
       
}