/* cdfdefint -- This program computes a definite integral on the
 * specified dimension of all fields with that dimension, leaving
 * them with one fewer dimensions.
 *
 * MODification:June-15-1999, This modification consist in changing the original
 * integration scheme for that of the  Composite Trapazide Rule as implemented
 * in ctrap_rg. So only regular grid are considered sofar!!
 *
 * Compilation line:
 *       gcc -c -O cdfdefint_exp.c && gcc -o cdfdefint_exp cdfdefint_exp.o ctrap_rg.o -lm -lcdf
 *
 */

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

static char hbin[HBMAX][LINE],hbout[HBMAX][LINE]; /* header buffers */
static char pline[LINE],cline[LINE];     /* line buffers */
static float *sbin,*sbout,*vbin,*vbout;  /* pointers to slice buffers */
static long nsbin,nsbout,nvbin,nvbout;   /* lengths of slice buffers */
static float bad,badlim;                 /* bad data stuff */
static int ds[4];                        /* dimension sizes */
static char dn[4][LINE];                 /* dimension names */
static long goop,poop,loop;              /* looping variables */
static long nflds;                       /* number of variable fields */
static int nlevels;                      /* number of good int levels */
static int start,badval;                 /* auxiliar at initialization step*/
static struct field *fp,*ifld;           /* field structure pointer */
static int nstat,nvar;                   /* numbers of fields */
static float *sfin,*vfin[HBMAX];         /* input field pointers */
static float *sfout,*vfout[HBMAX];       /* output field pointers */
static char smod[HBMAX],vmod[HBMAX];     /* modification flags */
static struct field ft[HBMAX];           /* table of field structures */
static long start1,incr1,size1;          /* subspace1 arguments */
static float deltadim,sum,element;       /* temporaries */
static float last_good_element;          /* temporaries */
static float *index_var;                 /* vector of ones for the index var */
static long index_size;                  /* size of index_var */
static char rmdim();                     /* dimension removing function */


int good_interval(float *filed,long starti, long start, long incr, long size,float badlim,
	          long *first_good, long *last_good, long *goodelements);
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 arguments */
  if (argc != 2) {
    fprintf(stderr,"Usage: cdfdefint_exp dimension_name\n");
    exit(1);
  }
/* get header buffer from standard input */
  gethdr(stdin,hbin,HBMAX);
/* check for float type */
  if (getfmt(hbin,HBMAX) != 'f') {
    fprintf(stderr,"cdfdefint_exp float format required\n");
    exit(1);
  }
/* allocate in buffers */
  nsbin = elemcnt(hbin,HBMAX,'s');
  sbin = getbuff(nsbin);
  nvbin = elemcnt(hbin,HBMAX,'v');
  vbin = getbuff(nvbin);

/* get static slice */
  getslice(stdin,nsbin,sbin);

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

/* compute deltadim from index parameters */
  deltadim = -1.0;
  sprintf(cline,"d%s",argv[1]);
  if (seekpar(hbin,HBMAX,cline,pline) == OK) deltadim = atof(pline);

  /*check existence, size and dimension for index filed*/
  if ((ifld = seekfld(hbin,HBMAX,'s',argv[1])) == NULL) {
     fprintf(stderr,"cdfdefint_exp dimension %s doesn't exist\n",argv[1]);
     exit(1);
  }
  if (ifld->tsize < 2 || ifld->dim != 1) {
     fprintf(stderr,"cdfdefint_exp field %s too small or not one dimensional\n",
             argv[1]);
      exit(1);
 }

 /*compute deltadim from index field if have not been set yet*/
 if (deltadim < 0 ) deltadim = sbin[ifld->ptr + 1] - sbin[ifld->ptr];

 /*get buffer for index_var to be integrated.*/
 index_size = ifld->tsize;
 index_var = getbuff(index_size);

/* make new header */
  nullhdr(hbout,HBMAX,"float");
  copycmt(hbout,HBMAX,hbin,HBMAX);
  sprintf(cline,"cdfdefint_exp %s\n",argv[1]);
  addcline(hbout,HBMAX,cline);
  copypar(hbout,HBMAX,hbin,HBMAX);
  nstat = 0;
  while ((fp = getfld(hbin,HBMAX,'s',nstat)) != NULL) {
    smod[nstat++] = rmdim(fp,'s',argv[1]);
  }
  nvar = 0;
  while ((fp = getfld(hbin,HBMAX,'v',nvar)) != NULL) {
    vmod[nvar++] = rmdim(fp,'v',argv[1]);
  }

  /* write new header */
  puthdr(stdout,hbout,HBMAX);

  /* out allocate buffers */
  nsbout = elemcnt(hbout,HBMAX,'s');
  sbout = getbuff(nsbout);
  nvbout = elemcnt(hbout,HBMAX,'v');
  vbout = getbuff(nvbout);

  /* transform static slice data */
  loop = 0;
  while ((fp = getfld(hbin,HBMAX,'s',loop)) != NULL) {
    strcpy(cline,fp->fname);
    sfout = getptr(hbout,HBMAX,sbout,'s','d',cline);
    sfin = getptr2(hbin,HBMAX,sbin,'s','d',cline,&fp);

    /* modified field -- integrate specified dimension */
    if (smod[loop++]) {
      poop = 0;
      while(subspace1(fp,argv[1],poop,&start1,&incr1,&size1) != FAIL) {
	   /*built index_var to be integrated */
	   for (goop = 0; goop < index_size; goop++) {
	       element = sfin[start1 + incr1*goop];
	       if (fabs(element) < badlim) index_var[goop] = 1.0;
	       else index_var[goop] = bad;
           }
	   sfout[poop++] = ctrap_rg(index_var,0,1,index_size,deltadim,bad,badlim);
      } /*loop over subspace*/
    }
    /* unmodified field -- simple transfer */
    else {
      for (poop = 0; poop < fp->tsize; poop++) sfout[poop] = sfin[poop];
    }
  }
/* put the modified static slice */
  putslice(stdout,nsbout,sbout);

/* compile information on variable slice fields */
  nflds = 0;
  while((fp = getfld(hbin,HBMAX,'v',nflds)) != NULL) {
    vfout[nflds] = getptr(hbout,HBMAX,vbout,'v','d',fp->fname);
    vfin[nflds] = getptr(hbin,HBMAX,vbin,'v','d',fp->fname);
    ft[nflds++] = *fp;
  }

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

       /* transform variable slice data */
       for (loop = 0; loop < nflds; loop++) {

           /* modified field -- integrate specified dimension */
           if (vmod[loop]) {
              poop = 0;

              while(subspace1(&ft[loop],argv[1],poop,&start1,&incr1,&size1) != FAIL) {

	         vfout[loop][poop++] = ctrap_rg(vfin[loop],start1,incr1,size1,deltadim,bad,badlim);

              } /*Loop over subspace*/
           }

           /* unmodified field -- simple transfer */
           else {
               for (poop = 0; poop < ft[loop].tsize; poop++) {
	           vfout[loop][poop] = vfin[loop][poop];
	       }
           }
       }/*loop over fields*/

       /* put the variable slice */
       putslice(stdout,nvbout,vbout);
  } /*loop over variable slices*/
  exit(0);
}

/* rmdim -- remove a dimension if it exists and add modified field
 * to output.
 */
static char rmdim(fp,type,dimname)
struct field *fp;
char type;
char *dimname;
{
  char flag;
  int dimout;
/* work through dimensions */
  dimout = 0;
  flag = 0;
  if (fp->dim > 0) {
    if (strcmp(dimname,fp->dname1) == 0) {
      flag = 1;
    }
    else {
      ds[dimout] = fp->dsize1;
      strcpy(dn[dimout++],fp->dname1);
    }
  }
  if (fp->dim > 1) {
    if (strcmp(dimname,fp->dname2) == 0) {
      flag = 1;
    }
    else {
      ds[dimout] = fp->dsize2;
      strcpy(dn[dimout++],fp->dname2);
    }
  }
  if (fp->dim > 2) {
    if (strcmp(dimname,fp->dname3) == 0) {
      flag = 1;
    }
    else {
      ds[dimout] = fp->dsize3;
      strcpy(dn[dimout++],fp->dname3);
    }
  }
  if (fp->dim > 3) {
    if (strcmp(dimname,fp->dname4) == 0) {
      flag = 1;
    }
    else {
      ds[dimout] = fp->dsize4;
      strcpy(dn[dimout++],fp->dname4);
    }
  }
/* add new field entry */
  addfld(hbout,HBMAX,type,fp->fname,fp->smul,fp->sadd,fp->bytes,dimout,
       dn[0],ds[0],dn[1],ds[1],dn[2],ds[2],dn[3],ds[3]);
/* return flag */
  return(flag);
}