/* cdfray_ct -- commonn thresholds for reflectivity and radial along the ray. gcc -c -O cdfray_ct.c && gcc -O -o cdfray_ct cdfray_ct.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 */ /* 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 *range; /* --> RANGE (index field) */ static float *dbz; /* --> DBZ */ static float *vr; /* --> VR */ static float *el; /* --> EL */ static float *alt; /* --> ALT */ static float in_vr[512]; /* 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); void fun_dbzthresh(); void fun_minrange(); void fun_seaclutter(); void fun_cut_surf_layer(); /* line arguments -- treated as global for simplicity*/ static char argu[MAXCHAR]; static float min_range; /* minimum range for which we accept gates */ static float cl_h; /*clutter layer height (km) */ static int flag_c = 0; /* Flags for different arguments options */ static int flag_r = 0; static int flag_t = 0; static int flag_s = 0; main(int argc,char **argv) { float refl; /* reflectivity factor */ /* 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 (flag_r) fun_minrange(); if (flag_s) fun_seaclutter(); if (flag_c) fun_cut_surf_layer(); /*has to be done after*/ putslice(stdout,nvbuff,vbuff); } /* read_slice_while */ fprintf(stderr,"\nDONE IT\n"); exit(0); } /* EOmain */ /*--------------------------------------------------------------------*/ void fun_minrange(){ /*set dbz and vr to bad for getes within the min_rage*/ float aux,faux; int iaux,loop; /*find last gate to be supressed*/ for (loop = 0; loop < GATES; loop++) { /* loop on range */ if (range[loop] > min_range) { iaux = loop; break; } } /*supress dbz and vr values from first to iaux gate.*/ for (loop = 0; loop < iaux; loop++) { dbz[loop] = vr[loop] = bad; } /* range_loop */ } /*--------------------------------------------------------------------*/ void fun_cut_surf_layer(){ /* Cleans-up everything from surface to the altitude cl_h (specified in the command line) cl_h must be lower than the flight level (i.e. only downward pointing rays are examined.) NOTE: It has to be done after the sea clutter is removed, from reflectivity data collected above flight level. */ int loop; int n_c = -1; for (loop = 0; loop < GATES; loop++) { if ( (*el < 0) && (*alt - (range[loop] * sin((fabs(*el) + 1.0) * D2R )) < cl_h) ) { n_c = loop - 1; break; } } if (n_c > 0 ) { /*clean up everything beyound n_c gate.*/ for (loop = n_c; loop < GATES; loop++) { dbz[loop] = vr[loop] = bad; } } } /*--------------------------------------------------------------------*/ /* Attempt to despecle a ray. Not done yet!! This maybe relevant only for unfolding velocities. I will implement this in cdfrayvr where, tsegments of good data is already done. Because, the only thing required is to get rid of good length-one segments. */ /*--------------------------------------------------------------------*/ void fun_seaclutter(){ /* elimination of specular reflection by discarding gates whose range value would put them beyoud the surface regardles where the ray is pointing to. This is an attempt of implementing the thechnique described by Testud ... complete the reference*/ #define N_t 5 #define N_d 10 #define DBZ_INIT -100 /*initial value in the search of max in dbz*/ int loop,n1_d,n2_d,n_max; float rs,aux,faux,dbz_max; /* get nominal ray-distance to surf based on fabs(el)*/ if (fabs(*alt) < badlim && fabs(*el) < badlim && (*el) != 0 ) { rs = *alt/fabs(sin((*el) * D2R)); if ( (rs > range[N_d]) && (rs < range[GATES-N_d -1]) ) { /* get interval centered in rs. This interval might be contaminated by sea clutter regardless to where the ray is pointing*/ for (loop = N_d ; loop < GATES; loop++ ) { if ( range[loop] > rs) { n1_d = loop - N_d; n2_d = loop + N_d; if ( n2_d > (GATES - 1)) n2_d = GATES - 1 ; break; } } /* find max in dbz within the interval centered in rs*/ dbz_max = DBZ_INIT; n_max = -1; for (loop = n1_d; loop <= n2_d; loop++ ) { aux = dbz[loop]; faux = fabs(aux); if ( faux < badlim && aux > dbz_max) { dbz_max = aux; n_max = loop; n1_d = n_max - N_t; n2_d = n_max + N_t; } } if (n_max > 0) { /* discard everything between n1_d and n2_d regardless where the ray is pointing to*/ for (loop = n1_d; loop <= n2_d ; loop++ ) { dbz[loop] = vr[loop] = bad;} if ( *el < 0 ){ /*Further discard everything beyound gate n2_d*/ for (loop = n2_d; loop < GATES; loop++ ) { dbz[loop] = vr[loop] = bad; } } } } } } /*--------------------------------------------------------------------*/ void fun_cgeta(int argc,char **argv) { /* fun_cgeta stuff */ char *STDMESSAGE=" "; int i,argtc; /* total */ int fixarg = 1; /* fixed arguments */ 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,"c:r:sh")) != -1) { switch (c) { case 'c' : flag_c = 1.0; if (sscanf(optarg," %f",&cl_h) != 1) fun_usage(": incorrect format for clutter_layer_height"); argtc+=2; break; case 'r' : flag_r = 1.0; if (sscanf(optarg," %f",&min_range) != 1) fun_usage(": incorrect format for minimum_range"); argtc+=2; break; case 's' : flag_s = 1.0; /*specular reflection*/ argtc+=1; break; case 'h' : fun_usage(STDMESSAGE); 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 input parameters data"); } argu[l++] = ' '; } argu[l] = '\n'; } /* EOfun_cgeta */ /*--------------------------------------------------------------------*/ void fun_usage(char *s) { char *USAGE = "Usage:cdfray_ct [-c h(km)][-r min_range][-s][-h]\n"; if (s != " ") fprintf(stderr,"%s\n",s); fprintf(stderr,"%s",USAGE); exit(1); } /* EOfun_usage */ /*--------------------------------------------------------------------*/ void fun_stop(char *s) { fprintf(stderr,"@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 */ /*sprintf(cline,": %s\n",argu); */ 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 fields depending on range */ 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); alt = getptr(hbuff,HBMAX,vbuff,'v','d',ALT); } /* 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 */ /*--------------------------------------------------------------------*/