/*Name: cdffmean 
 *Status: Exp
 *Date:June-15-1999
 *
 *Description: This program is a modification to cdfmean so that, the mean value
 * of uniform sampled  functions can be calculated in the following way:
 *
 * If F(z_i) is a function  sampled at n equally spaced points,  
 *
 * Then
 *
 *        mean[F(z_i)] = int(F(z_i))/int(dz) 
 *
 * instead of
 *
 *   mean[F(z_i)] = sum(F(z_i))/n
 *
 * The Integral of the function and increments are aproximated using 
 * The Composite Trapazoide Rule using the function ctrap_rg.
 *
 * Compilation line:
 *
 * gcc -c -O cdffmean.c && gcc -o cdffmean cdffmean.o ctrap_rg.o -lm -lcdf
 *
 */

#include <stdio.h>
#include <math.h>
#include "cdfhdr.h"

static char hbuff[HBMAX][LINE];              /* header buffer */
static char cline[3*LINE];                   /* line buffer */
static char pline[LINE];
static float *sbuff,*vbuff;                  /* slice buffers */
static long nsbin,nsbout,nvbin,nvbout;       /* slice buffer lengths */
static struct field *fp;                     /* field pointer temporary */
static float bad,badlim;                     /* bad data values */
static char *dimension,*suffix,*vars[HBMAX]; /* input strings */
static long loop,poop;                       /* looping indices */
static long nvars;                           /* number of variables affected */
static long indexsize;                       /* size of index field */
static long bingo;                           /* flag indicating dim exists */
static char ndname[4][LINE],*odname[4];      /* new and old dimension names */
static long ndsize[4],odsize[4];             /* new and old dimension sizes */
static long dimpos;                          /* counter dimension position */
static float *index_var,*field,*mfield;      /* subjects of operation */
static long start,incr,size,index,instance;  /* for subspace1 call */
static double fmean;                         /* for averaging */
static float goodelems;                      /* for averaging */
static char buffer[3*LINE];                  /* name buffer for new name */
static char *newname();                      /* cat suffix to name */
static float element,last_good_element;      /* temporaries */
static long start0;
static int badval;
static long first_index_val;
static float deltadim;
static float integ_of_field,integ_of_index;

float ctrap_rg(float *field, long start, long incr, long size, 
	                        float delta_dim, float bad, float badlim);

main (argc,argv)
int argc;
char *argv[];
{

/* check args */
  if (argc < 4) {
    fprintf(stderr,"Usage: cdffmean dimension suffix var1 var2 ...\n");
    exit(1);
  }
  dimension = argv[1];
  suffix = argv[2];
  nvars = argc - 3;
  for (loop = 0; loop < nvars; loop++) {
    vars[loop] = argv[loop + 3];
  }

/* get a header from standard input */
  gethdr(stdin,hbuff,HBMAX);
  if (getfmt(hbuff,HBMAX) != 'f') {
    fprintf(stderr,"cdffmean float format expected\n");
    exit(1);
  }

/* get bad data values */
  if (seekpar(hbuff,HBMAX,"bad",cline) == OK) bad = atof(cline);
  else bad = BAD;
  if (seekpar(hbuff,HBMAX,"badlim",cline) == OK) badlim = atof(cline);
  else badlim = BADLIM;

/* check to see if dimension exists */
  if ((fp = seekfld(hbuff,HBMAX,'s',dimension)) == NULL ||
       strcmp(fp->fname,fp->dname1) != 0) {
    fprintf(stderr,"cdffmean %s doesn't exist or is not an index field\n",
         dimension);
    exit(1);
  }
  if (fp->tsize < 2 || fp->dim != 1) {
     fprintf(stderr,"cdfdefint_exp field %s too small or not one dimensional\n",dimension);
     exit(1);
  }
  indexsize = fp->dsize1;
  first_index_val = fp->ptr;
  
  /* get buffer for index_var */
  index_var = getbuff(indexsize);
  
  /* compute deltadim from index parameters if available */
  sprintf(cline,"d%s",dimension);
  deltadim=-1.0;
  if (seekpar(hbuff,HBMAX,cline,pline) == OK) deltadim = atof(pline);

/* get input buffer sizes */
  nsbin = elemcnt(hbuff,HBMAX,'s');
  nvbin = elemcnt(hbuff,HBMAX,'v');

/* check consistency for each variable -- does it have the right dimension? */
  for (loop = 0; loop < nvars; loop++) {
    if ((fp = seekfld(hbuff,HBMAX,'v',vars[loop])) == NULL) {
      fprintf(stderr,"cdffmean variable field %s doesn't exist\n",vars[loop]);
      exit(1);
    }
    bingo = 0;
    if (fp->dim > 0 && strcmp(fp->dname1,dimension) == 0) bingo = 1;
    if (fp->dim > 1 && strcmp(fp->dname2,dimension) == 0) bingo = 2;
    if (fp->dim > 2 && strcmp(fp->dname3,dimension) == 0) bingo = 3;
    if (fp->dim > 3 && strcmp(fp->dname4,dimension) == 0) bingo = 4;
    if (bingo == 0) {
      fprintf(stderr,"cdffmean variable %s doesn't have dimension %s\n",
           vars[loop],dimension);
      exit(1);
    }

/* construct new entry information for reduced variable */
    dimpos = 0;
    odname[0] = fp->dname1;
    odname[1] = fp->dname2;
    odname[2] = fp->dname3;
    odname[3] = fp->dname4;
    odsize[0] = fp->dsize1;
    odsize[1] = fp->dsize2;
    odsize[2] = fp->dsize3;
    odsize[3] = fp->dsize4;
    for (poop = 0; poop < fp->dim; poop++) {
      if (poop != bingo - 1) {
        strcpy(ndname[dimpos],odname[poop]);
        ndsize[dimpos] = odsize[poop];
        dimpos++;
      }
    }

/* add the new field */
    addfld(hbuff,HBMAX,'v',newname(vars[loop],suffix),fp->smul,fp->sadd,
         fp->bytes,fp->dim - 1,ndname[0],ndsize[0],ndname[1],ndsize[1],
         ndname[2],ndsize[2],ndname[3],ndsize[3]);
  }

/* add a comment */
  sprintf(cline,"cdffmean %s %s\n",dimension,suffix);
  addcline(hbuff,HBMAX,cline);
  loop = 3;
  while (loop < argc) {
    strcpy(cline,"");
    for (poop = 0; poop < 3; poop++) {
      if (loop >= argc) break;
      strcat(cline," ");
      strcat(cline,vars[loop - 3]);
      loop++;
    }
    strcat(cline,"\n");
    addcline(hbuff,HBMAX,cline);
  }

/* get output buffer sizes and create buffers */
  nsbout = elemcnt(hbuff,HBMAX,'s');
  nvbout = elemcnt(hbuff,HBMAX,'v');
  sbuff = getbuff(nsbout);
  vbuff = getbuff(nvbout);

/* put header and transfer static slice */
  puthdr(stdout,hbuff,HBMAX);
  getslice(stdin,nsbin,sbuff);

        /* compute deltadim from index field if it has not been set from
	 * params*/
        if (deltadim <  0.0 )
	    deltadim = sbuff[first_index_val + 1] - sbuff[first_index_val];

  putslice(stdout,nsbout,sbuff);

/* loop on variable slices */
  while (getslice(stdin,nvbin,vbuff) != EOF) {

/* loop on variables to be averaged */
    for (loop = 0; loop < nvars; loop++) {
      field = getptr2(hbuff,HBMAX,vbuff,'v','d',vars[loop],&fp);
      mfield = getptr(hbuff,HBMAX,vbuff,'v','d',newname(vars[loop],suffix));

/* loop on subspace cutting thru field in direction of averaging dimension */
      instance = 0;
      while (subspace1(fp,dimension,instance,&start,&incr,&size) != FAIL) {
	  
        /*built vector for index_var */
	for (poop = 0; poop < size; poop++) {
	    if (fabs(field[start + incr*poop]) < badlim) 
		index_var[poop] = 1.0;
	    else 
	        index_var[poop] = bad;
	}

        integ_of_field = ctrap_rg(field,start,incr,size,deltadim,bad,badlim);
	integ_of_index = ctrap_rg(index_var,0,1,indexsize,deltadim,bad,badlim);

	mfield[instance++] = integ_of_field / integ_of_index;

      }/*loop-while*/

    }/*loop-over variables to be averaged--for_loop*/

    /* write variable slice */
    putslice(stdout,nvbout,vbuff);

  }/* Loop over Var Slices*/
  exit(0);

}/*End-Of-Main*/

/* newname -- concatenate name and suffix together */
static char *newname(name,suffix)
char *name,*suffix;
{
  strcpy(buffer,name);
  strcat(buffer,suffix);
  return(buffer);
}