/* cdfrayrv -- work on radial velocities along ray. gcc -c -O cdfrayrv.c && gcc -O -o cdfrayrv cdfrayrv.o -lm -lcd */ #include #include #include "cdfhdr.h" #define GATES 512 /*Number of gates in ray*/ #define MAXCHAR 100 /*Max characters on the calling line*/ #define D2R 0.01745329278 /* pi/180 conversion factor (degrees to radians)*/ /* assumed names of radial velocity and reflectivity fields */ #define VR "vr" #define DBZ "dbz" /* assumed names defining positioning */ #define RANGE "range" /* range */ #define ALT "alt" /* altitude */ #define EL "el" /* elevation */ #define AZ "az" /* Azimuth */ /* assumed names defining velocities */ #define U_AC "ua" /* ac-velocity x-dir */ #define V_AC "va" /* ac-velocity y-dir */ #define W_AC "wa" /* air-velocity z-dir */ #define U_AIR "u" /* air-velocity x-dir */ #define V_AIR "v" /* air-velocity y-dir */ #define W_AIR "w" /* air-velocity z-dir */ /* assumed name for Nyquist Velocity */ /* #define NYQ_V "nyq_vel" */ #define NYQ_V "nyqvel" float nyqvel_0 = 0; /*critical t value for alpha =0.xxx from abramowitz & Stegun for convinience, tc_xxx was set to -1 so the index corresponds to degrees of freedom from 1 to 30 */ static float tc_001[31] = {-1,636.619,31.598,12.924,8.610,6.869,5.959,5.408,5.041,4.781,4.587,4.437,4.318,4.221,4.140,4.073,4.015,3.965,3.922,3.883,3.850,3.819,3.792,3.768,3.745,3.725,3.707,3.690,3.674,3.659,3.646}; static float tc_001_lnu[4]={3.551,3.460,3.373,3.291}; static float tc_005[31] = {-1,127.321,14.089,7.453,5.598,4.773,4.317,4.029,3.833,3.690,3.581,3.497,3.428,3.372,3.326,3.286,3.252,3.223,3.197,3.174,3.153,3.135,3.119,3.104,3.090,3.078,3.067,3.057,3.047,3.038,3.030}; static float tc_005_lnu[4]={2.971,2.915,2.860,2.807}; static float tc_010[31] = {-1,63.657,9.925,5.841,4.604,4.032,3.707,3.499,3.355,3.250,3.169,3.106,3.055,3.012,2.977,2.947,2.921,2.898,2.878,2.861,2.845,2.831,2.819,2.807,2.797,2.787,2.779,2.771,2.763,2.756,2.750}; static float tc_010_lnu[4]={2.704,2.660,2.617,2.576}; static float tc_050[31] = {-1,12.706,4.303,3.182,2.776,2.571,2.447,2.365,2.306,2.262,2.228,2.201,2.179,2.160,2.145,2.131,2.120,2.110,2.101,2.093,2.086,2.080,2.074,2.069,2.064,2.060,2.056,2.052,2.048,2.045,2.042}; static float tc_050_lnu[4]={2.021,2.000,1.980,1.960}; /* candis stuff */ static char hbuff[HBMAX][LINE]; static char cline[LINE]; static float bad,badlim; static float *sbuff,*vbuff; static long nsbuff,nvbuff; static struct field *fp; /* names for pointers to input fields */ static float *nyqv; /* --> NYQ_V */ static float *range; /* --> RANGE (index field) */ static float *dbz; /* --> DBZ */ static float *vr; /* --> VR */ static float *el; /* --> EL */ static float *az; /* --> AZ */ static float *alt; /* --> ALT */ static float *uac; /* --> U_AC */ static float *vac; /* --> V_AC */ static float *wac; /* --> W_AC */ static float *uair; /* --> U_AIR */ static float *vair; /* --> V_AIR */ static float *wair; /* --> W_AIR */ static float in_vr[GATES]; /* not required radial velocity --> VR */ /* functions */ void fun_make_header(); void fun_alloc(); void fun_make_s_slice(); void fun_stop(char *s); void fun_usage(char *s); void fun_cgeta(int argc,char **argv); int fun_find_good_segments(); void fun_dbzthresh(); void fun_unfoldit(); void fun_unfoldit_lnyqi(); void fun_acorr(float *v, float *cr); void fun_vel_thresh(); void fun_basic_statistics(float *x,float *mx, float *stdx, float *nx); /* line arguments -- treated as global for simplicity*/ static char argu[MAXCHAR]; static float dbz_thresh; /* reflectivity threshold (dbz) */ static float vel_thresh; /* reflectivity threshold (dbz) */ static float min_range; /* minimum range for which we accept gates */ static float cl_h; /*clutter layer height (km) */ static int flag_r = 0; static int flag_u = 0; static int flag_a = 0; static int flag_t = 0; /*other globals */ static int good1v[GATES]; static int good2v[GATES]; static int disposses = 0; static long int unfo_tot_cnt; static long int unfo_bad_cnt; static long acor_bad_cnt; static long acor_tot_cnt; float *tc_t,*tc_t_lnu; int conf_level = 999; main(int argc,char **argv) { float refl; /* reflectivity factor */ float ac; /* auto correlation coeficient*/ float stat; /*aux to hold percent of ray rejections*/ int i; if ( conf_level == 999 ) { tc_t = tc_001; tc_t_lnu = tc_001_lnu; } if ( conf_level == 995 ) { tc_t = tc_005; tc_t_lnu = tc_005_lnu; } if ( conf_level == 990 ) { tc_t = tc_010; tc_t_lnu = tc_010_lnu; } if ( conf_level == 950 ) { tc_t = tc_050; tc_t_lnu = tc_050_lnu; } /* for (i=0; i< 31; i++){ fprintf(stderr,"DEBUG_MAIN:tc_t=%f\n",tc_t[i]);} for (i=0; i< 4; i++) { fprintf(stderr,"DEBUG_MAIN:tc_t_lnu=%f\n",tc_t_lnu[i]);} exit(1); */ /* check and get the command line Parameters */ fun_cgeta(argc,argv); /* Candis header */ fun_make_header(); /* allocate memory for static & variable fields */ fun_alloc(); /* Candis static slice */ fun_make_s_slice(); /* loop through input variable slices */ while (getslice(stdin,nvbuff,vbuff) != EOF) { if (*nyqv != nyqvel_0) { fprintf(stderr,"NYQUISTVEL CHANGE: %f --> %f\n",nyqvel_0,*nyqv); nyqvel_0 = *nyqv; } if (flag_r) fun_dbzthresh(); if (flag_u) fun_unfoldit_lnyqi(); /* fun_unfoldit(); fun_unfoldit_lnyqi(); */ if (flag_a) fun_acorr(vr, &ac); if (flag_t) fun_vel_thresh(); /*must be done after unfolding*/ putslice(stdout,nvbuff,vbuff); } /* read_slice_while */ if (flag_u && unfo_bad_cnt > 0){ /*UNFO rejection statistics*/ stat = 100.0 * ((float)unfo_bad_cnt/(float)unfo_tot_cnt); fprintf(stderr,"UNFO_STAT:%d\t%d\t\t%f\n",unfo_tot_cnt,unfo_bad_cnt,stat); /* use this when looking for stats by the ray unfo_tot_cnt = unfo_bad_cnt=0; */ } if (flag_a && acor_bad_cnt > 0 ) { /*AUTO_CORR rejection Statistics*/ stat = 100.0 * ((float)acor_bad_cnt/(float)acor_tot_cnt); fprintf(stderr,"ACOR_STAT:%d\t%d\t%f\n",acor_bad_cnt,acor_tot_cnt,stat); } exit(0); } /* EOmain */ /*--------------------------------------------------------------------*/ void fun_acorr(float *v, float *cr){ /* * This function calculate the autocorrelation (lag =1 )of the given * vector to the confidence level (or betterfor N > 32) specified in * main. The input vector is set to bad when the confidence criteria * is not met. ---------------------------------------------------------------------- Bad Data along the ray are treated in a way that is equivalent to "single inputation with mean value", where the mean calculated ignoring bad data. The input vector (of original size N) may contain M bad data values in it in any order. In this case, the autocorrelation is calculated by skiping terms that involve bad data, which is equivalent to single fill of missing data with the mean value of the input vector. This procedure does not affect the value of the autocorrelation coeficient, which can be seen from the following argument: In general the autocorrelation coefficient is given by (see for instance http://www.asp.ucar.edu/colloquium/1992/notes/part1/course.html by William A. Cooper and Thomas W. Schlatter ) sum_{i=1}^{i=N-k} (v_{i} - OVM)(v_{i+k} - OVM) r= --------------------------------------------------- sum_{i=1}^{i=N} (v_i - OVM)^2 where OVM is the mean of the original vector. The mean is calculted using only good values, which is equivalent to set the bad data value to OVM. The effect of setting bad values of v to OVM in the expresion for r is simply to skip terms that contain bad data. Therefore, we may simply skip those terms regardless the value assigned to bad with out modifying r. ---------------------------------------------------------------------- Confidence in the autocorrelation: Let t_c the critical value of the t-student distribution, such that A(nu,t_c) = 1 - alpha/2 where alpha is the significance level for a two tailed test of the no-correlation hypotesis. In general t is related to r as follows: t = r*sqrt((k-2)/(1-r^2)) where k = nu+2 is the number of data in the sample -- good data that is. For a given alpha, t_c is found from tables, abramowitz&stegun in our case. Since our tables only offer values for t_c for nu in [1,30] and {40,60,120,inf}, we will use the actual critical value from the table to accept or reject the null hypothesis of no-correlation for rays with less than 33 data. Rays with more data, are subjet to a more stringet criteria(see coding below) and their confidence level exeeds the criteria. In doing this we might loose some good rays however, is probably best some good ones than to put trash in the analysis. So in setting the confidence level, one should be based on the principle that accepting the null hypothesis means to bring garbache into the analysis. However you are fre to experiment. -------------------------------------------------------------------- */ /*bad, badlim, and GATES are globals and so are tc_t and tc_t_lnu*/ int i,nu; float nx; float xi,xipl; float mx,varx,stdx,maxdev; float acc = 0.0; float t_d = 0.0; float tx = 0.0; float snx = 0.0; float outl_cnt = 0.0; float SET_RAY_TO_BAD = 0.0; /*get mean, standard deviation, and number of good data from input vector*/ fun_basic_statistics(v,&mx,&stdx,&nx); if (nx > 2 ) { /*minumun amount of good data required */ maxdev=2.5*stdx; /*maximun deviation allowed from the mean of v*/ for (i = 0; i < GATES; i++ ) { /*check for outliers*/ if (fabs(v[i] - mx) > maxdev) { v[i] = bad; outl_cnt++;} } if (outl_cnt > 0) { /*recompute basic statistics*/ fun_basic_statistics(v,&mx,&stdx,&nx); } /* estimate acc*/ if ( stdx > 1e-3 ) { snx = 0.0; for (i = 0; i < (GATES - 1); i++ ) { xi = v[i]; xipl = v[i+1]; if (fabs(xi) < badlim && fabs(xipl) < badlim) { snx += (xi - mx)*(xipl - mx); } } varx = (nx -1.0)* stdx*stdx; acc = snx/varx; } /*else deside what to do about acc -- currently zero at this branch */ /*check for BUGS in the code that my lead to this*/ if (fabs(acc) > 1) {fprintf(stderr,"acc=%f > 1\n",acc); exit(1); } /*---check statistics --------------------------------------- Null_Hypothesis of No-correlation @ the chosen confidence level -- probability of accepting thrash is set to that level --*/ nu = ceil(nx - 2.0); /*get degrees of freedom from number of good data*/ t_d = acc*sqrt((nx-2.0)/(1.0-acc*acc)); /* t from data*/ if (nu <= 30) tx = tc_t[nu]; else if (nu < 40 ) tx = tc_t[30]; else if (nu < 60 ) tx = tc_t_lnu[0]; /*tc_[40]*/ else if (nu < 120 ) tx = tc_t_lnu[1]; /*tc_[60]*/ else tx = tc_t_lnu[2]; /*tc_[120]*/ if ( fabs(t_d) < tx) {SET_RAY_TO_BAD = 1.0;}/*couldn't reject NullHypotesis*/ /* fprintf(stderr,"%f\t%i\t%f\t%f\t%f\t%f\n",nx,nu,t_d,tx,acc,SET_RAY_TO_BAD); */ } else {SET_RAY_TO_BAD = 1.0; /*Not enough good data*/} /*----Discard Ray if necessary------------------------------------------*/ if (SET_RAY_TO_BAD > 0.0) { /*seting whole ray to bad data*/ if (nx > 0) acor_bad_cnt++; for (i =0; i < GATES; i++){ v[i] = bad; } } acor_tot_cnt++; } /*EOFU*/ /*--------------------------------------------------------------------*/ void fun_basic_statistics(float *x,float *mx, float *stdx, float *nx) { /* This function calculates the mean(mx), standard deviation (stdx) * and number of good data(nx) in the input vector(x) */ int i; float xi,sx,sxx,lx; *mx=*stdx=bad; *nx=0.0; lx=sx=sxx=0.0; for (i = 0; i < GATES; i++ ) { xi = x[i]; if (fabs(xi) < badlim ) { sx += xi; sxx += (xi*xi); lx++; } } *nx = lx; /* number of good gates in ray*/ if (lx > 1.0 ) { *mx=sx/lx; /* mean */ *stdx=sqrt((sxx - lx * (*mx) * (*mx))/(lx -1.0)); /*standard deviation*/ } if (lx > 0.0 ) *mx=sx/lx; /* mean */ /* fprintf(stderr,"%f\t%f\t%f\n",mx,stdx,nx); */ } /*--------------------------------------------------------------------*/ void fun_dbzthresh() { /*only allow radial velocities on gates that have good reflectivities*/ int i; for (i=0; i < GATES; i++) { if ( fabs(dbz[i]) > badlim ) vr[i] = bad; } } /*--------------------------------------------------------------------*/ void fun_vel_thresh() { /*after unfolding, allow only radial velocities smaller than vel_thresh*/ int i; for (i=0; i < GATES; i++) { if ( fabs(vr[i]) > vel_thresh ) vr[i] = bad; /*oviously, bad is larger than vel_thresh*/ } } /*--------------------------------------------------------------------*/ void fun_unfoldit_lnyqi() { /* This function is to do the unfolding when the NYQI (nyquist interval) * is so large that we don't have to worry about true radial velocities * being outside the NYQI. Although the retrived radial velocities * are relative to the plane, they are contaminated with the aircraft * velocity most of the time. Hence, they are observed as larger than * the Nyquist velocity and consequently artificially folded into the * NYQI. However when we are sure that the true radial velocity for * the targets are always inside the NYQ interval, we can unfold the * measured velocities to their right place without using the * continuity-along-the-ray approach. * Note that continuity can be done afterwards, by a different function, * when true radial velocities are expected to be larger than the * Nyquist velocity. * * Rationale -- * In theory, the velocity measured by the radar (v_m) is the target velocity relative to the aircraft folded into the nyquist velocity interval. Let, v_t, v_ac, and v_rel respectibely be the target velocity, aircraft velocity, and relative velocity of the target with respect to the aircraft. The vector relationship among this velocities is v_rel = (v_t - v_ac) The dot product of this equation with a unit vector along the ray, n, is v_rel . n = v_t . n - v_ac . n Since v_rel is in the ray direction, calling v_t_r = v_t . n ,and v_ac_r = v_ac . n the radial velocities of the target and aircraft, respectively, results in v_rel = v_t_r - v_ac_r what gets measure is v_rel folded in to the Nyquist interval: v_m = v_t_r - v_ac_r +- k * nyq_inter where k is an unknown integer. The radial velocity of the target can therefore be obtained as folows for some n: v_t_r = v_m + v_ac_r -+ k * nyq_inter The Algorith to find k in this function has basically two steps: 1) add the projection of the aircraft velocity onto the ray to the measured velocity -- v_m + v_ac_r; 2) add or substract an interger times of the NYQI to until the result is inside the NYQI -- fabs(v_m + v_ac_r -+ k * nyq_inter) <= NYQI. * */ /* Globals: int GATES,unfo_tot_cnt,unfo_bad_cnt; float *nyqv,*el,*az,uac,vac,wac,bad,badlim;*/ int g; float radel,radaz; float nx,ny,nz; float v_ac_r; float half_nyq_inter; /* half of the unfolding interval */ float nyq_inter; /* unfolding interval */ float nyqi_left,nyqi_right; float vr_curr_gate; float vr_prev_gate; if (fabs(*el) < badlim && fabs(*az) < badlim && fabs(*uac) < badlim && fabs(*vac) < badlim && fabs(*wac) < badlim ) { /*check if basic info is available*/ /* compute aircraft speed along the ray */ radel = (*el)*D2R; radaz = (*az)*D2R; nx = cos(radel)*sin(radaz); ny = cos(radel)*cos(radaz); nz = sin(radel); v_ac_r = (*uac)*nx + (*vac)*ny + (*wac)*nz; /*get size of nyquist interval*/ half_nyq_inter = (*nyqv); nyq_inter = 2.0 * half_nyq_inter; /*extremes of the NYQI centered in zero*/ nyqi_right = 0.0 + half_nyq_inter; nyqi_left = 0.0 - half_nyq_inter; for (g=0; g < GATES; g++) { /*cycle over ray*/ vr_curr_gate = vr[g]; if ( vr_curr_gate < badlim ) { /*only unfold good radial velocities*/ /* add aircraft velocity before unfolding*/ vr_curr_gate += v_ac_r; /* do unfolding */ while (vr_curr_gate > nyqi_right) { vr_curr_gate -= nyq_inter; } while (vr_curr_gate < nyqi_left) { vr_curr_gate += nyq_inter; } vr[g] = vr_curr_gate; /*Update Value of radial velocity*/ } } /*finish with ray*/ } /* not enough info to unfold velocities */ else { unfo_tot_cnt++; for (g=0; g < GATES; g++) vr[g] = bad; } unfo_tot_cnt++; /*count number of rays proccessed*/ } /*END_OF_fun_unfoldit_lnyqi*/ /*--------------------------------------------------------------------*/ void fun_unfoldit() { int good,i,g,start,end; float radel,radaz; float nx,ny,nz; float v_ac_r,v_inswind_r; float vcorre,vmean,vmean_prev_block; float half_nyq_inter; /* half of the unfolding interval */ float nyq_inter; /* unfolding interval */ float nyqi_left,nyqi_right; float vr_curr_gate; float vr_prev_gate; float vr_prev_block; /* compute components of aircraft and insitu-wind vel along ray */ if (fabs(*el) < badlim && fabs(*az) < badlim && fabs(*uac) < badlim && fabs(*vac) < badlim && fabs(*wac) < badlim && fabs(*uair) < badlim && fabs(*vair) < badlim && fabs(*wair) < badlim ) { radel = (*el)*D2R; radaz = (*az)*D2R; nx = cos(radel)*sin(radaz); ny = cos(radel)*cos(radaz); nz = sin(radel); v_ac_r = (*uac)*nx + (*vac)*ny + (*wac)*nz; v_inswind_r = (*uair)*nx + (*vair)*ny + (*wair)*nz; half_nyq_inter = (*nyqv); nyq_inter = 2.0 * half_nyq_inter; /* fprintf(stderr,"UNFO:%f\t%f\t%f\t%f\n",v_inswind_r,v_ac_r,half_nyq_inter,nyq_inter); */ good = fun_find_good_segments(); /*aplly unfolding so consecutive gates in each good block(i.e. contiguos gates by definition of good block)have radial velocities in the same Nyquist interval. In theory, the velocity measured by the radar (v_m) is the target velocity relative to the aircraft folded into the nyquist velocity interval. if v_t -- target velocity v_ac -- aircraft velocity v_rel -- relative velocity of the target with respect to the ac. v_rel = (v_t - v_ac) if n is a unit vector in the ray direction, since v_rel is in the ray direction, we have v_rel . n = v_t . n - v_ac . n calling v_t_r = v_t . n ,and v_ac_r = v_ac . n the radial velocities of the target and aircraft, respectively, v_rel = v_t_r - v_ac_r Now v_rel gets folded into the nyquist interval, that is v_m = v_t_r - v_ac_r +- n * nyq_inter where n is an unknown integer. so, the radial velocity of the target can be obtained as folows: v_t_r = v_m + v_ac_r -+ n * nyq_inter the way we find n is by continuity along the ray assuming that velocities of contiguos gates must be inside the same nyquist interval we further assume that gates closer to the aircraft must have a their velocity in the same nyquist interval as the wind velocity measured by the aircraft insitu instrumentation. */ /*Initialize vr_prev_block to the radial-ac-wind speed*/ vr_prev_block = v_inswind_r; for (g=0; g <= good; g++) { /*cycle over good block*/ /*get index of the interval extremes*/ start = good1v[g]; end=good2v[g]; /*Initialize vr_prev_gate to velocity value in previous block*/ vr_prev_gate = vr_prev_block; for (i = start ; i <= end; i++) { /*Do block */ /*find extremes of the NYQI centered in previous gate*/ nyqi_right = vr_prev_gate + half_nyq_inter; nyqi_left = vr_prev_gate - half_nyq_inter; /* add aircraft velocity before unfolding*/ vr_curr_gate = vr[i] + v_ac_r; /* do unfolding */ while (vr_curr_gate > nyqi_right) { vr_curr_gate -= nyq_inter; } while (vr_curr_gate < nyqi_left) { vr_curr_gate += nyq_inter; } vr[i] = vr_prev_gate = vr_curr_gate; /*Update Values*/ } /*finish with block*/ /*get vel at the end of this block to update vr_prev_block -- this culd be made more robust by taking some weighted average of velocities in the block */ vr_prev_block = vr[end]; } /*finish with all blocks in this ray*/ /* first compute mean in good block */ /* vmean = 0.; for (i = start; i <= end; i++) vmean += vr[i]; vmean /= (float)(end - start + 1); */ } /*not enough info to unfold it; as much as I hate it, I've get rid of it*/ else { for (i=0; i < GATES; i++) vr[i] = bad; disposses++; } /* if ( unfo_bad_cnt > 0 ){ fprintf(stderr,"DEBUG_MAIN:%d\t%d\t\t%f\n",unfo_tot_cnt,unfo_bad_cnt,100.0 *((float)unfo_bad_cnt/(float)unfo_tot_cnt)); } */ /* unfo_tot_cnt = unfo_bad_cnt=0; */ } /*--------------------------------------------------------------------*/ int fun_find_good_segments() { /*return Num. of good segments in ray. The vectors goodv and goodv2 contain the positions for starting and ending points for all good segments found in the ray */ int i = 0; int good = 0; while (i < GATES) { good1v[good] = good2v[good] = GATES; /* fprintf(stderr,"GOOD:%i\t%i\t%i\n",good,good1v[good],good2v[good]); */ while ( (i < GATES) && (fabs(vr[i]) > badlim) ) { i++;} /*skip bads*/ if ( i < GATES) { good1v[good] = i;} while ( (i < GATES) && (fabs(vr[i]) < badlim) ) { i++;} if ( i < GATES) { good2v[good] = i-1;} else {good--;} /*fprintf(stderr,"DEBUG_good_seg:%i\t%i\t%i\n",good,good1v[good],good2v[good]);*/ good++; } good -= 1; return(good); } /*--------------------------------------------------------------------*/ void fun_cgeta(int argc,char **argv) { /* fun_cgeta stuff */ char *STDMESSAGE=" "; int i,argtc; /* total */ int maxarg = 5; /* maximum number of arguments */ /* getopt stuff */ extern char *optarg; /* extern int optind = 1; */ int c; int j,k,l; if (argc > maxarg+1) {fprintf(stderr,"argc = %i\n",argc); exit(1);} /* Process Arguments */ while ((c = getopt(argc,argv,"ruat:h")) != -1) { switch (c) { case 'r' : flag_r = 1.0; /*activate fun_dbzthresh*/ argtc+=1; break; case 'u' : flag_u = 1.0; /*activate unfolding*/ argtc+=1; break; case 'a' : flag_a = 1.0; /*activate fun_acorr -- autocorrelation*/ argtc+=1; break; case 't' : flag_t = 1.0; /*activate fun_vel_thresh */ if (sscanf(optarg,"%f",&vel_thresh) != 1) fun_usage("cdfrayrv: bad -t option"); fprintf(stderr,"vel_thresh=%f\n", vel_thresh); argtc+=2; break; case 'h' : fun_usage(STDMESSAGE); argtc+=1; break; default : fun_usage(STDMESSAGE); break; } } /*get string with all called arguments in it*/ l = 0; for ( k=0 ; k < argc; k++) { j=0; while ( (c = argv[k][j++]) != '\0' ) { if (l < MAXCHAR) argu[l++] = c; else fun_stop("Not enough room to hold data"); } argu[l++] = ' '; } argu[l] = '\n'; } /* EOfun_cgeta */ /*--------------------------------------------------------------------*/ void fun_usage(char *s) { char *USAGE = "Usage:cdfrayrv [-r][-u][-a][-h]\n"; if (s != " ") fprintf(stderr,"%s\n",s); fprintf(stderr,"%s",USAGE); exit(1); } /* EOfun_usage */ /*--------------------------------------------------------------------*/ void fun_stop(char *s) { fprintf(stderr,"cdfrayrv@fun_%s",s); exit(1); } /*--------------------------------------------------------------------*/ void fun_make_header() { /* get header of input file */ gethdr(stdin,hbuff,HBMAX); if (getfmt(hbuff,HBMAX) != 'f') { fun_stop("make_header:format expected"); } /*get bad and badlim */ if (seekpar(hbuff,HBMAX,"bad",cline) == OK) bad = atof(cline); else {fprintf(stderr,"Could't find bad value\n"); exit(1);} if (seekpar(hbuff,HBMAX,"badlim",cline) == OK) badlim = atof(cline); else {fprintf(stderr,"Could't find badlim value\n"); exit(1);} /* make header of output file */ /* comments */ addcline(hbuff,HBMAX,argu); /* write output header */ puthdr(stdout,hbuff,HBMAX); } /* EOfun_make_header */ /*--------------------------------------------------------------------*/ void fun_alloc() { /* allocate buffers */ nsbuff = elemcnt(hbuff,HBMAX,'s'); nvbuff = elemcnt(hbuff,HBMAX,'v'); sbuff = getbuff(nsbuff); vbuff = getbuff(nvbuff); /* get pointers for index fields */ range = getptr2(hbuff,HBMAX,sbuff,'s','d',RANGE,&fp); if (fp->dim != 1 || strcmp(fp->fname,fp->dname1) != 0) { fprintf(stderr,"@alloc: %s is not an index field\n",fp->fname); exit(1); } if ( fp->dsize1 != GATES ) { fprintf(stderr,"#of gates is not 512!! \n"); exit(1); } /* get pointers for static fields */ /* nyqv = getptr(hbuff,HBMAX,sbuff,'s','d',NYQ_V); */ nyqv = getptr(hbuff,HBMAX,vbuff,'v','d',NYQ_V); /* get pointers for var fields */ dbz = getptr(hbuff,HBMAX,vbuff,'v','d',DBZ); vr = getptr(hbuff,HBMAX,vbuff,'v','d',VR); /* get pointer elevation */ el = getptr(hbuff,HBMAX,vbuff,'v','d',EL); az = getptr(hbuff,HBMAX,vbuff,'v','d',AZ); alt = getptr(hbuff,HBMAX,vbuff,'v','d',ALT); uac = getptr(hbuff,HBMAX,vbuff,'v','d',U_AC); vac = getptr(hbuff,HBMAX,vbuff,'v','d',V_AC); wac = getptr(hbuff,HBMAX,vbuff,'v','d',W_AC); uair = getptr(hbuff,HBMAX,vbuff,'v','d',U_AIR); vair = getptr(hbuff,HBMAX,vbuff,'v','d',V_AIR); wair = getptr(hbuff,HBMAX,vbuff,'v','d',W_AIR); } /* EOfun_alloc */ /*--------------------------------------------------------------------*/ void fun_make_s_slice() { /* get input static slice */ getslice(stdin,nsbuff,sbuff); /* write output static slice */ putslice(stdout,nsbuff,sbuff); } /* EOfun_make_s_slice */ /*--------------------------------------------------------------------*/