/* cdfinteg_exp -- Take the integral of a specified field with respect * to a specified index field. * * Modification Jun-3-1999: Feliz cumple Adri. (my sister) * I'm reviewing this so it do not reset the value of the integral * after a bad value. This is done in the effort to brig the results * of this program into agreement with those of cdfdefint at the extremes * of the interval of integration. I belive the reseting thing is not good. * By the way in doing this revision it seems clear to me that the integration * scheme used by this program is a Composite Trapazoide Rule. (Not Simpson * Rule as declared on its man page). * * Modification Jun-27-1999: This is to incorporate the function ctrap_rg * This function calculate the integral over a regular grid using * composite trapazoide rule. It is intended to be a core function for all * cdf programs performing integrations. * * Since ctrap_rg only works for regular grids, this program will only * work for those kind of grids. * * Compilation line: * gcc -c -O cdfinteg_exp.c && gcc -o cdfinteg_exp cdfinteg_exp.o -lm -lcdf * * Carlos Lopez-Carrillo * * */ #include #include #include "cdfhdr.h" static char *hbuff[HBMAX][LINE]; /* header buffer */ static char cline[LINE],pline[LINE]; /* character line buffer */ static long nsbuff,nvbin,nvbout; /* slice buffer sizes */ static float *sbuff,*vbuff; /* slice buffer pointers */ static float *ifield; /* pointer to specified index field */ static float *sfield; /* pointer to specified field */ static float *dfield; /* pointer to integral field */ static long instance; /* looping var over instances */ static long loop; /* the other looping variable */ static struct field *fp; /* field structure temporary pointer */ static struct field *fps; /* field structure pointer for specified field */ static float bad,badlim; /* bad data values */ static long start,incr,size; /* loop variables */ static int flag; /* up (+1) or down (-1) indicator */ static long index; /* field looping variable */ static double sum; /* integration temporary */ static long naux; static long start_aux,incr_aux,size_aux; static float deltadim; static int sign,factor; 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 != 5) { fprintf(stderr,"Usage: cdfinteg_exp -[tb] integrated_fld field dimension\n"); exit(1); } /* get header */ gethdr(stdin,hbuff,HBMAX); if (getfmt(hbuff,HBMAX) != 'f') { fprintf(stderr,"cdfinteg_exp float format expected\n"); exit(1); } /* get deltadim from index parameters */ deltadim = -1.0; sprintf(cline,"d%s",argv[1]); if (seekpar(hbuff,HBMAX,cline,pline) == OK) deltadim = atof(pline); /* get bad data values */ bad = (seekpar(hbuff,HBMAX,"bad",cline) == OK) ? bad = atof(cline) : BAD; badlim = (seekpar(hbuff,HBMAX,"badlim",cline) == OK) ? atof(cline):BADLIM; /* get the flag */ if (strcmp(argv[1],"-b") != 0 || strcmp(argv[1],"-t") != 0) flag = (strcmp(argv[1],"-b") == 0) ? 1 : -1; else { fprintf(stderr,"cdfinteg_exp flag type is not known\n"); exit(1); } /* compute input buffer lengths */ nsbuff = elemcnt(hbuff,HBMAX,'s'); nvbin = elemcnt(hbuff,HBMAX,'v'); /* add stuff to header, comment first */ sprintf(cline,"cdfinteg_exp %s\n",argv[1]); addcline(hbuff,HBMAX,cline); sprintf(cline," %s %s %s\n",argv[2],argv[3],argv[4]); addcline(hbuff,HBMAX,cline); /* integral field next -- same dimensions as specified field */ if ((fp = seekfld(hbuff,HBMAX,'v',argv[3])) == NULL) { fprintf(stderr,"cdfinteg_exp %s field not found in variable slice\n",argv[3]); exit(1); } addfld(hbuff,HBMAX,'v',argv[2],fp->smul,fp->sadd,fp->bytes,fp->dim, fp->dname1,fp->dsize1,fp->dname2,fp->dsize2, fp->dname3,fp->dsize3,fp->dname4,fp->dsize4); /* allocate buffers */ nvbout = elemcnt(hbuff,HBMAX,'v'); sbuff = getbuff(nsbuff); vbuff = getbuff(nvbout); /* assign field and structure pointers */ sfield = getptr2(hbuff,HBMAX,vbuff,'v','d',argv[3],&fps); dfield = getptr(hbuff,HBMAX,vbuff,'v','d',argv[2]); if ((ifield = getptr2(hbuff,HBMAX,sbuff,'s','r',argv[4],&fp)) == NULL) { fprintf(stderr,"cdfinteg_exp %s field doesn't exist\n",argv[4]); exit(1); } if ((fp->dim != 1) || (strcmp(fp->fname,fp->dname1) != 0)) { fprintf(stderr,"cdfinteg_exp %s field is not an index field\n",argv[4]); exit(1); } /* write header */ puthdr(stdout,hbuff,HBMAX); /* transfer static slice */ getslice(stdin,nsbuff,sbuff); /*compute deltadim from index field if have not been set yet*/ if (deltadim < 0 ) deltadim = sbuff[fp->ptr + 1] - sbuff[fp->ptr]; putslice(stdout,nsbuff,sbuff); /* up or down? */ if (flag == 1) { /* up! */ sign = 1.0; factor = 0.0; } else { /*down*/ sign = -1.0; factor = 1.0; } /* loop on variable slices */ while (getslice(stdin,nvbin,vbuff) != EOF) { /* loop on subspace instances */ instance = 0; while (subspace1(fps,argv[4],instance,&start,&incr,&size) != FAIL) { start_aux = start + factor * incr * (size - 1.0); incr_aux = sign*incr; /*loop over grid */ /* dfield[start_aux] = 0.0; integral got to be zero at starting point*/ for (loop = 0; loop < size; loop++ ) { index = start_aux + incr_aux*loop; /*Insure bad for bad entries other wise their values will be set to the last calculated value of the integral.*/ if (sfield[index] < badlim) { size_aux = loop + 1; dfield[index] = sign * ctrap_rg(sfield,start_aux,incr_aux,size_aux, deltadim,bad,badlim); } else dfield[index] = bad; }/*loop over grid*/ /* increment instance */ instance++; }/*loop-subspace-while*/ /* write variable slice */ putslice(stdout,nvbout,vbuff); } /*loop-slices-while*/ exit(0); }/*End-Of-main*/