/*Name: ctrap-rg.c -- composite-trapaziode rule for regular grids. Status: Exp Date:June-22-1999 Author: Carlos Lopez-Carrillo Description: This function is ment to be used as a core routine in candis programs that perform integrals of functions sampled on regular grids. The order of integration preceeds in the order of the index variable -- direct order. Calling statement: integral = ctrap_rg(float *iv, long start, long incr, long size, float delta_dim, float bad, float badlim) parameter description: *iv -- pointer to the vector that contains the sample values for the variable being integrated. This is normally obtained by the user with getptr2 on the calling program. start -- the starting point for integration var on the input vector. incr -- step between values of the integration var on the input vector. size -- number of vlues of the integration var on the input vector. start, incr and size are accses information required to locate, on the input vector, the values of the variable being integrated. Normally the user can get this info using the function subspace1 on the calling program. If the integral is not to be done over the entire grid, then the appropiate value of size has to be given in order to stop the integration on the desired point. This would be the right thing to do for example in cdfinteg. If only an integration on some subinterval of the domain is required, then start as well as the size values have to be entered approperly in order to reflect the rigth subinterval. delta_dim -- grid interval. This is normally gather from parameters using the function seekparameter on the calling program. bad -- bad data value. badlim -- badlim value. Pocedure: A loop over the input vector is runed to find continuos segments of good values. Every time a segment is find the integral over that segment is calculated as follows. If the number of good values in the segment is bigger than two, then the composite trapazoide rule is applied to calculate the integral from the first to the last good value. If the number of good values in the segment is equal to two then trapazoide rule is applied to calculate the integral. If the number of good values is equal to one then it is skipped. This is the same as saying that spikes do not contribute to the final integral. if there are not good values on the input vector, then the integral is set to bad. Every time a integral over a segment is calculated it is added to the final value of the integral over the the entire domain of the index field. MOD-MAY-12-200: If the final value of the integral exeeds badlim then a warning message is issued. This may occur when integrating over lage domains and/or fields with large values. The normal solution is to scale down the field being integrated. Compilation line: gcc -c -O ctrap_rg.c */ #include #include #include /*------------------------------------------------------------------------------------- F U N C T I O N H E A D E R S -------------------------------------------------------------------------------------*/ 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); /*-------------------------------------------------------------------------------------*/ /* */ /* F U N C T I O N I M P L E M E N T A T I O S */ /*-------------------------------------------------------------------------------------*/ float ctrap_rg(float *field, long start, long incr, long size, float delta_dim, float bad, float badlim) { float integral; /*final value of the integral of the variable over the specified grid*/ float sub_integral; int done,good_intervals; long start0,starti,loop,first_good,last_good,goodelements,index_f,index_l; float aux_integral; aux_integral=0.0; integral=0.0; starti=0; done=0; good_intervals=0; while ( done == 0 ) { done = good_interval(field,starti,start,incr,size,badlim,&first_good,&last_good,&goodelements); starti=last_good+1; switch ( goodelements ) { case 0 : if (good_intervals == 0 ) integral = bad; if (done != 1 ) { /*leaved for testing purposes only*/ fprintf(stderr,"SOMETHIG_WRONG: inegral was set to bad but it is not done yet\n"); fprintf(stderr,"start=%i<--",start); exit(1); } break; case 1 : /*spikes do not contribute to the integration*/ integral = 0.0; good_intervals++; break; default : /*Note that if goodelements=2 then the for loop is not executed so this reduces to a simple trapezoide rule for two points if goodelements > 2 then composite trapezoide rule is performed. */ sub_integral=0.0; index_f =start + incr*first_good; index_l =start + incr*last_good; sub_integral = (field[index_f] + field[index_l]) * 0.5 ; start0=first_good+1; for (loop = start0; loop < last_good ; loop++ ) { sub_integral += field[start + incr*loop]; } sub_integral *= delta_dim; integral += sub_integral; /*integrals over good intervals add up to the final value of the integral*/ /* MOD-MAY-12-2000 Even though he integral is comming from additions of good values of the field over the grid, it can very well exeed the balim value at some point. In order be aware of this event, we use the auxiliar variable aux_integral. */ aux_integral = integral; good_intervals++; } } /*MOD-MAY-12-2000*/ if ( aux_integral >= badlim ) { /*This test can not be done using the variable "integral" as it could had been set to bad in cases where there are no good points. What we want to check here is the event of the integral becoming larger of badlim as a result of adding good values of the field. This is why we need to use the auxiliar variable*/ fprintf(stderr,"\n ctrap_rg !!! WARNING !!! value of integral exeeds badlim\n"); fprintf(stderr," badlim = %f;\tintegral = %f\n",badlim,aux_integral); fprintf(stderr,"A possible solution is to scale down the field being integrated\n"); } return(integral); } /*End-Of-ctrap_rg*/ int good_interval(float *field,long starti, long start, long incr, long size, float badlim, long *first_good, long *last_good, long *goodelements) { float element; long start0,loop,index; /* find first good value */ *goodelements=0; for (start0 =starti; start0 < size; start0++) { index = start + incr*start0; element = field[index]; if (fabs(element) < badlim) { *first_good=start0; *last_good= start0; /* initialize last good to first good just in case there is only one good */ (*goodelements)++; start0++; /*In order to have the usual behavior of counter increment, we have to manually increase the counter before the break */ break; } } /*find last good of this segment*/ for (loop = start0; loop < size; loop++) { index = start + incr*loop; element = field[index]; if (fabs(element) < badlim) { *last_good=loop; (*goodelements)++; } else break; } if (loop == size) return(1); else return(0); } /*End-Of-good_interval*/