/* radcart6 -- same as radcart5 except that I've included * the interpolation of the field "rsgp" (reflectivity seen grid point) * to the cartesian grid. * This field is range dependent and is build from the reflectivity * field before any postprocessing threshold is applied to the data. * It is assigned the value of 1, wherever the reflectivity * has a good value, and it is set to bad wherever the reflectivity * has a bad value. Since its values are constant, the result of its * interpolation will assign 1 to any grid point that has relfectivities * measurements made close to the grid point . If no measurements were * made in that vecinity, rsgp will have a bad value for that grid * point. The idea is to mark all grid points * for which the radar istrumentation has assigned a good value of * reflectivity in their neiborhoood. This will give an indication * of regions that were actually seen by the radar but does not * show up on the reflectivity map because of thresholds applied * in postprocessing. In this way we can distinguish * regions of radar malfunctinig from regions that did not meet * our threshold criteria and also give us the oportunity to see * how much are we loosing due to a very stringent threshold. * * I also have include the field "rara" for the rain rate. Normally * rain rate is calculated from reflectivity after the interpolation is * carried out. "rara", however, will be calculated on a gate-by-gate * basis and then interpolated to the cartesian grid. The downside of * this is that whenever a more accurate ZvsR relationship is developed, * we have to redo the whole interpolation proccess. * "rara" is created from the reflectivity input. The function fun_rara * is in charge of this calculation. * * April-21-2002(under test): since Z-R relation dependes basically on two parameters * I'm going to include two input parameters to feed this relationship. * care should be made to easily extend this for more parameters in case * a different Z-R or several Z-R are used. * * Oct 9,2002: The elimination of radial velocities based on * reflectivity threshold ( if dbz is bad then vr is bad) is * been eliminated. This if required has to be done beforehand * using, for instance, the cdfraydbz filter. This will allowd * to analyse radial velocities in cases where the reflectivity field * is not good. This is Ok since radial velocities and reflectivities * are independen measures. * * Oct 11, 2002: I've included the field "vsgp" (velocity seen grid point) * which has the same purpose of rsgp but for the radial * velocities. * * radcart5 is based on Carlos Lopez-Carrillo's modifications * to Lucio Lopez's Cartesianization program, which is based on the * Cartesianization program of Sharon Lewis, which is based on the * original program by Dave Raymond. It does only dual Doppler. * Modifications made by Dave Raymond. * * The output consists of the following variables: * x,y,z - static fields holding the grid * time - overall time for data collections * avgtime - time for data collection @ each point * dbz - reflectivity factor from the weighted average of * the interpolated reflectivities @ each point * ux0,uy0 - Dual doppler velocities assuming zero vertical velocity @ * each point * gx,gy - from ux=ux0-gx*vz and uy=uy0-gy*vz, in dual analysis * for postprocessing after solving the continuity equation * eur,evr - standard errors squared in dual doppler velocities assuming * zero error in precip vertical velocity * eup,evp - standard errors squared in dual doppler velocities assuming * zero error in radial velocities * pcount - number of data used for least square technique @ each point * * rara - rain rate in mm/hr * * rsgp - reflectivity seen grid point * * vsgp - velocity seen grid point * * */ #include #include #include "cdfhdr.h" #define D2R 0.017453 /* degrees to radians conversion */ #define TRUNC(a0,da,aloc,aa) ((int)(aa = (aloc - a0)/da + 1000.) - 1000) #define FRACT(ia,aa,ba) (aa -= ia + 1000, ba = 1. - aa) /* assumed names of radial velocity and reflectivity fields */ #define VR "vr" #define DBZ "dbz" #define RR "rara" #define RSGP "rsgp" #define VSGP "vsgp" /* other assumed names defining positioning (meaning self evident) */ #define RANGE "range" /* range */ #define LON "lon" /* longitude */ #define LAT "lat" /* latitude */ #define ALT "alt" /* altitude */ #define TIME "monotime" /* monotonically guaranteed time */ #define ALTTIME "time" /* alternate time if monotime doesn't exist */ #define AZ "az" /* azimuth */ #define EL "el" /* elevation */ /* candis stuff */ static char hbin[HBMAX][LINE]; static char hbout[HBMAX][LINE]; static char cline[LINE],pline[LINE]; static float bad,badlim; static float *sbin,*vbin,*sbout,*vbout; static long nsbin,nvbin,nsbout,nvbout; static struct field *fp; /* names for pointers to input fields */ static float *range; /* --> RANGE (index field) */ static float *time; /* --> TIME */ static float *lon; /* --> LON */ static float *lat; /* --> LAT */ static float *alt; /* --> ALT */ static float *az; /* --> AZ */ static float *el; /* --> EL */ /* fields to be cartesianized */ static float *in_dbz; /* raflectivity factor --> DBZ */ static float *in_vr; /* radial velocity --> VR */ static float *in_rsgp; /* */ static float *in_vsgp; /* */ /* fields having to do with the grid */ static float *xg,*yg,*zg; /* grid point coordinates (index fields) */ static float *initime; /* initial time @ each point */ static float *avgtime; /* time interval to get all data @ each pt. */ static float *dbz; /* cartesianized output-fields */ static float *rara; /* cartesianized output-fields */ static float *rsgp; /* cartesianized output-fields */ static float *vsgp; /* cartesianized output-fields */ static float *nfac_dbz; /* normalization factor for dbz */ static float *nfac_rsgp; /* normalization factor for rsgp */ static float *nfac_vsgp; /* normalization factor for vsgp */ static float *pcount; /* # of vr good instances @ each point */ /* dual doppler stuff @ each grid point -- all them are output fields */ static float *ux0,*uy0; /* dual doppler velocities assuming wp = 0 */ static float *gx,*gy; /* velocity correction factors for wp != 0 */ static float *eur,*evr,*eup,*evp; /*errors in dual doppler */ /* other global variables */ static float *ctime; /* pendiente?????? */ static float orax,oray,oraz; /* orientation ray angles.(among ray & axis) */ /* internal fields having to do with the doppler analysis @ each point */ static float *a00,*a01,*a02,*a11,*a12; /*A-elements eq.(pend) */ static float *w00,*w01,*w11; /*G-elements eq.(pend)*/ static float *b0,*b1; /*have to do with right hand side of eq.(pend) */ /* Parameters for the Z-R relationship*/ static float alpha,beta; /*parameters for the Z-R */ /* functions needed */ void fun_make_header(float *args, float x0,float y0,float z0, int nx,int ny,int nz); void fun_alloc(int *nrange); void fun_fill_index(float x0,float dx,int nx,float *xg, float y0,float dy,int ny,float *yg, float z0,float dz,int nz,float *zg); void fun_make_s_slice(); void fun_get_buffer(long num); void fun_ini_gridpoints(long num, float *rara, float *rsgp,float *nfac_rsgp, float *vsgp,float *nfac_vsgp, float *dbz,float *nfac_dbz, float *ux0,float *uy0,float *gx, float *gy,float *eur,float *evr,float *eup, float *evp,float *pcount, float *initime,float *avgtime,float *a00, float *a01,float *a02,float *a11,float *a12, float *b0,float *b1,float *w00,float *w01, float *w11); /* special functions */ void fun_stop(char *s); void fun_usage(char *s); void fun_assign(int iloc,float weight,float timenow, float in_vr,float refl,float in_rsgp,float in_vsgp); void fun_weight(int nx,int ny,int nz,int ix,int iy, int iz, float cx, float cy, float cz,int *np,float *we); void fun_ploc(float range, float xl, float yl, float zl, float deltax, float deltay, float *xloc, float *yloc,float *zloc); void fun_rloc(float lon0,float lat0,float alt0,float *xl,float *yl,float *zl); void fun_ora(float *orax,float *oray,float *oraz); void fun_filldbz_thresh(int flag_p,float dbz0,float r0, float *range,int maxrange, float *dbz_thresh); void fun_cgeta(int argc,char **argv,float *args); void fun_dualdoppa(float a00,float a01,float a02, float a11,float a12, float b0, float b1,float w00, float w01, float w11, float pcount,long loop); float fun_rara(float Z); main(int argc,char **argv) { int nx,ny,nz; /* grid dimensions */ long num; /* total ponints on the grid */ float x0,y0,z0; /* cartesian Origin */ float xl,yl,zl; /* radar location from Origin */ float xloc,yloc,zloc; /* point location from Origin */ int ix,iy,iz; /* grid location */ float cx,cy,cz; /* parameters to define weigth function */ float refl; /* reflectivity factor */ int tind; /* if set then time goes in average for volume */ float timenow; /* slice's time */ float deltatime; /* difference between slice's time and reftime */ float deltax,deltay; /* drift due to storm velocity */ int np[8]; /* neighborhood points */ float we[8]; /* weights for neighborhood points*/ int corner; /* loop variable over neighborhood points */ float *dbz_thresh; /* dbz threshold as function of range */ long cnt; /* # of slices that have at least 1 data in the grid */ long lloop; /* loop var for long loops */ int iloop; /* loop var for int loops */ int nrange; /* maxrange -- no more than 512 gates are expected */ /* line arguments */ float lon0,lat0,alt0; /* grid origin */ float dx,dy,dz; /* step size */ float xkm,ykm,zkm; /* domain size */ float svx,svy,reftime; /* storm velocities and reference time */ float hsc_thresh; /* ht(km) for surface clutter thresh at 0 range */ float hsc_coef; /* surface clutter thresh range coefficient */ float mingatez; /* minimum elevation to accept gates for surface clutter suppression -- function of range */ float r0; /* reference point to define power threshold */ float dbz0; /* minimum value for dbz at r0 */ float el_thresh; /* threshold for elevation angles (in degrees) */ float min_range; /* minimum range for which we accept gates */ int flag_p,flag_c,flag_s,flag_e,flag_r; /* Flags for different arguments options */ float args[26]; /* store parameters after check */ /* check and get the command line Parameters */ fun_cgeta( argc,argv,args); /* INITIALIZATION */ lon0 = args[1]; dx = args[2]; xkm = args[3]; lat0 = args[4]; dy = args[5]; ykm = args[6]; alt0 = args[7]; dz = args[8]; zkm = args[9]; flag_p = args[10]; r0 = args[11]; dbz0 = args[12]; flag_c = args[13]; hsc_thresh = args[14]; hsc_coef = args[15]; flag_s = args[16]; svx = args[17]; svy = args[18]; reftime = args[19]; flag_e = args[20]; el_thresh = args[21]; flag_r = args[22]; min_range = args[23]; if (args[24] > 0) { alpha = args[24]; beta = args[25]; } else { alpha = 153.0; beta = 1.53; } /* make x0,y0,z0 equal to zero, then the origin of the grid is at lon0,lat0,alt0. */ x0 = y0 = z0 = 0.0 ; /* the +1 is to compensate for truncation */ nx = xkm/dx + 1.5; ny = ykm/dy + 1.5; nz = zkm/dz + 1.5; /* read inheader and make outheader */ fun_make_header(args,x0,y0,z0,nx,ny,nz); /* allocate memory (get buffer) for static & variable fields */ fun_alloc(&nrange); /* fill in index fields *xg,*yg,*zg */ fun_fill_index(x0,dx,nx,xg,y0,dy,ny,yg,z0,dz,nz,zg); /* get s slice and make the newone */ fun_make_s_slice(); /* dbz thresholds: allocate space & fill dbz_thresh-vector dbz_thresh = getbuff(nrange); fun_filldbz_thresh(flag_p,dbz0,r0,range,nrange,dbz_thresh); */ /* calculate How many points does the grid have */ num = nx*ny*nz; /* GET BUFFERS FOR INTERNAL point-dependent field buffers for all output and input fields were gotten @ fun_alloc */ fun_get_buffer(num); /* INITIALIZE EVERY point-dependent field */ fun_ini_gridpoints(num,rara,rsgp,nfac_rsgp,vsgp,nfac_vsgp,dbz,nfac_dbz, ux0,uy0,gx,gy,eur,evr,eup,evp, pcount,initime,avgtime, a00,a01,a02,a11,a12,b0,b1,w00,w01,w11); /* loop through input variable slices */ deltatime = 0.; cnt = 0; *ctime = 0.; while (getslice(stdin,nvbin,vbin) != EOF) { /* initialize tind for this ray (default: ray outside grid space) */ tind = 0; /* get location of radar respect to grid origin -- these distances are approximated by the arc lengths */ fun_rloc(lon0,lat0,alt0,&xl,&yl,&zl); /* made elevation thresholding if ( (flag_e == 0) || fabs(*el) <= el_thresh ) { */ /* get the orientation angules for this ray */ fun_ora(&orax,&oray,&oraz); /* time is assumed to be monotonic -- hrdrad insures this if monotime was available */ timenow = *time; /* Get the time difference between this data slice and the reference time */ deltatime = timenow - reftime; /* get the correction due to storm velocity */ deltax = svx*deltatime; deltay = svy*deltatime; /* loop on range */ for (iloop = 0; iloop < nrange; iloop++) { /*check if there is somthing to interpolate*/ if (fabs(in_vr[iloop]) < badlim || fabs(in_dbz[iloop]) < badlim || fabs(in_rsgp[iloop]) < badlim || fabs(in_vsgp[iloop]) < badlim) { /* convert reflectivity from log (dbz) to linear (refl) scale*/ if (fabs(in_dbz[iloop]) < badlim) { refl = pow(10.,(in_dbz[iloop]/10.)); } else refl = bad; /* Get location of point respect to grid origin */ fun_ploc(range[iloop],xl,yl,zl,deltax,deltay,&xloc,&yloc,&zloc); /* find grid location */ ix = TRUNC(x0,dx,xloc,cx); iy = TRUNC(y0,dy,yloc,cy); iz = TRUNC(z0,dz,zloc,cz); /* check that the data point is in the grid space */ if (ix >= 0 && ix < nx - 1 && iy >= 0 && iy < ny - 1 && iz >= 0 && iz < nz - 1) { tind = 1; /* define neighbourhood points ,np[8], & their weights, we[8] */ fun_weight(nx,ny,nz,ix,iy,iz,cx,cy,cz,np,we); /* fill the eight neighbourhood points around */ for(corner = 0; corner < 8 ; corner++ ) { fun_assign(np[corner],we[corner],timenow, in_vr[iloop],refl,in_rsgp[iloop],in_vsgp[iloop]); } } /*if point inside*/ } /*Nothing good to interpolate */ /* if tind is set then account for time */ if (tind) { /* means at least 1 point of this ray was interpolated */ *ctime += timenow; /* increase recollection-data time */ ++cnt; /* account for another slice--remember 1ray/slice */ } /* else this slice was skipped */ } /* range_loop */ /*} slice skipedd due to elevation angle thresholding*/ } /* read_slice_while */ /* POST PROCESSING */ if (cnt > 0) { /* check if some grid point was hit */ *ctime /= (float)cnt; /* normalize ctime */ /* filling up output-grid fields */ for (lloop = 0; lloop < num; lloop++) { /* loop over grid points */ /* normalize dbz, rara & restore reflectivity factor */ if (nfac_dbz[lloop] > 0.0) { dbz[lloop] /= nfac_dbz[lloop]; rara[lloop] /= nfac_dbz[lloop]; dbz[lloop] = 10.*log10(dbz[lloop]); } else { dbz[lloop] = rara[lloop] = bad; } /* normalize rsgp */ if (nfac_rsgp[lloop] > 0.0) rsgp[lloop] /= nfac_rsgp[lloop]; else rsgp[lloop] = bad; /* normalize vsgp */ if (nfac_vsgp[lloop] > 0.0) vsgp[lloop] /= nfac_vsgp[lloop]; else vsgp[lloop] = bad; /* check if velocities were interpolated to this point*/ if (avgtime[lloop] < badlim) { /* get the time interval in secs for this point */ avgtime[lloop] = (avgtime[lloop] - initime[lloop])*1000.; fun_dualdoppa(a00[lloop],a01[lloop],a02[lloop],a11[lloop], a12[lloop],b0[lloop],b1[lloop],w00[lloop], w01[lloop],w11[lloop],pcount[lloop],lloop); } } /* loop over grid points */ /* write variable slice */ putslice(stdout,nvbout,vbout); exit(0); } else fun_stop("radcart6: No data for this grid"); } /* EOmain */ /*--------------------------------------------------------------------*/ void fun_cgeta(int argc,char **argv,float *args) { /* fun_cgeta stuff */ char *STDMESSAGE=" "; int i,argtc; /* total */ int fixarg = 9; /* fixed arguments */ int maxarg = 24; /* maximum number of arguments */ /* getopt stuff */ extern char *optarg; extern int optind; int c; if (argc == 1) fun_usage(STDMESSAGE); if (argc < fixarg+1) fun_usage("radcart_x: too few arguments"); if (argc > maxarg+1) fun_usage("radcart_x: too many arguments"); /* Process Arguments */ /* initialization */ for (i = 0; i <= maxarg; i++) args[i] = 0.0; /* set default value for optional arguments(if != 0.0) */ /* Process Optional Arguments (if any) */ optind = 1; while ((c = getopt(argc,argv,"p:c:s:e:r:z:")) != -1) { switch (c) { case 'p' : args[10] = 1.0; if (sscanf(optarg,"%f:%f",&args[11],&args[12]) != 2) fun_usage("radcart_x: incorrect format for power_thresh"); if (args[11] <= 0.0) fun_usage("radcart_x: r0 must be positive"); argtc += 2; break; case 'c' : args[13] = 1.0; if (sscanf(optarg," %f:%f",&args[14],&args[15]) != 2) fun_usage("radcart_x: incorrect format for clutter_thresh"); argtc+=2; break; case 's' : args[16] = 1.0; if (sscanf(optarg,"%f:%f:%f",&args[17],&args[18],&args[19]) != 3) fun_usage("radcart_x: incorrect format for storm_movement"); argtc+=2; break; case 'e' : args[20] = 1.0; if (sscanf(optarg," %f",&args[21]) != 1) fun_usage("radcart_x: incorrect format for elevation_thresh"); argtc+=2; break; case 'r' : args[22] = 1.0; if (sscanf(optarg," %f",&args[23]) != 1) fun_usage("radcart_x: incorrect format for minimum_range"); argtc+=2; break; case 'z' : if (sscanf(optarg," %f:%f",&args[24],&args[25]) != 2) fun_usage("radcart_x: incorrect format for Z-R parameters"); argtc+=2; break; default : fun_usage(STDMESSAGE); break; } } /* set fixed arguments */ if (argc != optind + fixarg) fun_usage("radcart_x: incorrect number of args"); for (i = 1; i <= fixarg; i++) { args[i] = atof(argv[i + optind - 1]); } } /* EOfun_cgeta */ /*--------------------------------------------------------------------*/ void fun_usage(char *s) { char *USAGE = "Usage: radcart_x [-p r0:dbz0] [-c zmin0:zmincoef] [-s svx:svy:reftime] [-e el] [-r min_range] [-z alpha:beta] -- lon0 dx xkm lat0 dy ykm alt0 dz zkm\n"; if (s != " ") fprintf(stderr,"%s\n",s); fprintf(stderr,"%s",USAGE); exit(1); } /* EOfun_usage */ /*--------------------------------------------------------------------*/ void fun_stop(char *s) { fprintf(stderr,"radcart_x@fun_%s",s); exit(1); } /*--------------------------------------------------------------------*/ void fun_make_header(float *args,float x0,float y0,float z0, int nx,int ny,int nz) { /* get header of input file */ gethdr(stdin,hbin,HBMAX); if (getfmt(hbin,HBMAX) != 'f') { fun_stop("make_header:format expected"); } /* make header of output file */ /* comments */ nullhdr(hbout,HBMAX,"float"); copycmt(hbout,HBMAX,hbin,HBMAX); sprintf(cline,"radcart_x:\n % 5.2f % 5.2f % 5.2f\n",args[1],args[2],args[3]); addcline(hbout,HBMAX,cline); sprintf(cline," % 5.2f % 5.2f % 5.2f\n",args[4],args[5],args[6]); addcline(hbout,HBMAX,cline); sprintf(cline," % 5.2f % 5.2f % 5.2f\n",args[7],args[8],args[9]); addcline(hbout,HBMAX,cline); if (args[10] >= 1.) { sprintf(cline," -p % 5.2f % 5.2f\n",args[11],args[12]); addcline(hbout,HBMAX,cline); } if (args[13] >= 1.) { sprintf(cline," -c % 5.2f % 5.3f\n",args[14],args[15]); addcline(hbout,HBMAX,cline); } if (args[16] >= 1.) { sprintf(cline," -s % 5.2f % 5.2f % 5.2f\n",args[17],args[18],args[19]); addcline(hbout,HBMAX,cline); } if (args[20] >= 1.) { sprintf(cline," -e % 5.2f\n",args[21]); addcline(hbout,HBMAX,cline); } if (args[22] >= 1.) { sprintf(cline," -r % 5.1f\n",args[23]); addcline(hbout,HBMAX,cline); } if (args[24] >= 0.) { sprintf(cline,"alpha % 6.3f\n",args[24]); addcline(hbout,HBMAX,cline); sprintf(cline,"beta % 5.3f\n",args[25]); addcline(hbout,HBMAX,cline); } /* parameters */ copypar(hbout,HBMAX,hbin,HBMAX); /* bad data values */ if (seekpar(hbin,HBMAX,"bad",pline) == OK) bad = atof(pline); else { bad = BAD; sprintf(pline,"%g",bad); addpar(hbout,HBMAX,"bad",pline); } if (seekpar(hbin,HBMAX,"badlim",pline) == OK) badlim = atof(pline); else { badlim = BADLIM; sprintf(pline,"%g",badlim); addpar(hbout,HBMAX,"badlim",pline); } /* index parameters */ sprintf(pline,"%g",args[1]); addpar(hbout,HBMAX,"lon0",pline); sprintf(pline,"%g",x0); addpar(hbout,HBMAX,"x0",pline); sprintf(pline,"%g",args[2]); addpar(hbout,HBMAX,"dx",pline); sprintf(pline,"%g",args[4]); addpar(hbout,HBMAX,"lat0",pline); sprintf(pline,"%g",y0); addpar(hbout,HBMAX,"y0",pline); sprintf(pline,"%g",args[5]); addpar(hbout,HBMAX,"dy",pline); sprintf(pline,"%g",args[7]); addpar(hbout,HBMAX,"alt0",pline); sprintf(pline,"%g",z0); addpar(hbout,HBMAX,"z0",pline); sprintf(pline,"%g",args[8]); addpar(hbout,HBMAX,"dz",pline); /* add index fields */ addfld(hbout,HBMAX,'s',"x",1000.,0.,'l',1,"x",nx); addfld(hbout,HBMAX,'s',"y",1000.,0.,'l',1,"y",ny); addfld(hbout,HBMAX,'s',"z",1000.,0.,'l',1,"z",nz); /* add variable fields with range dimension to be cartesianized other fields are dropped */ /* Add radial reflectivity */ if ( (fp = seekfld(hbin,HBMAX,'v',DBZ)) != NULL) if (fp->dim == 1 && strcmp(fp->dname1,RANGE) == 0) addfld(hbout,HBMAX,'v',fp->fname,fp->smul,fp->sadd,'c',3, "x",nx,"y",ny,"z",nz); else fun_stop("radcart_x: DBZ is not 1-dim or does not depend on range"); else fun_stop("radcart_x: can't find DBZ"); /* Add rsgp */ if ( (fp = seekfld(hbin,HBMAX,'v',RSGP)) != NULL) if (fp->dim == 1 && strcmp(fp->dname1,RANGE) == 0) addfld(hbout,HBMAX,'v',fp->fname,fp->smul,fp->sadd,'c',3, "x",nx,"y",ny,"z",nz); else fun_stop("radcart_x: RSGP is not 1-dim or does not depend on range"); else fun_stop("radcart_x: can't find RSGP"); /* Add vsgp */ if ( (fp = seekfld(hbin,HBMAX,'v',VSGP)) != NULL) if (fp->dim == 1 && strcmp(fp->dname1,RANGE) == 0) addfld(hbout,HBMAX,'v',fp->fname,fp->smul,fp->sadd,'c',3, "x",nx,"y",ny,"z",nz); else fun_stop("radcart_x: VSGP is not 1-dim or does not depend on range"); else fun_stop("radcart_x: can't find VSGP"); /* add new fields */ addfld(hbout,HBMAX,'v',RR,10.,0.,'s',3,"x",nx,"y",ny,"z",nz); addfld(hbout,HBMAX,'v',"ux0",10.,0.,'s',3,"x",nx,"y",ny,"z",nz); addfld(hbout,HBMAX,'v',"uy0",10.,0.,'s',3,"x",nx,"y",ny,"z",nz); addfld(hbout,HBMAX,'v',"gx",10.,0.,'s',3,"x",nx,"y",ny,"z",nz); addfld(hbout,HBMAX,'v',"gy",10.,0.,'s',3,"x",nx,"y",ny,"z",nz); addfld(hbout,HBMAX,'v',"eur",10.,0.,'s',3,"x",nx,"y",ny,"z",nz); addfld(hbout,HBMAX,'v',"evr",10.,0.,'s',3,"x",nx,"y",ny,"z",nz); addfld(hbout,HBMAX,'v',"eup",10.,0.,'s',3,"x",nx,"y",ny,"z",nz); addfld(hbout,HBMAX,'v',"evp",10.,0.,'s',3,"x",nx,"y",ny,"z",nz); addfld(hbout,HBMAX,'v',"avgtime",10.,0.,'s',3,"x",nx,"y",ny,"z",nz); addfld(hbout,HBMAX,'v',"pcount",10.,0.,'s',3,"x",nx,"y",ny,"z",nz); addfld(hbout,HBMAX,'v',"time",1000.,0.,'l',0); /* write output header */ puthdr(stdout,hbout,HBMAX); } /* EOfun_make_header */ /*--------------------------------------------------------------------*/ void fun_alloc(int *nrange) { /* allocate buffers */ nsbin = elemcnt(hbin,HBMAX,'s'); nsbout = elemcnt(hbout,HBMAX,'s'); nvbin = elemcnt(hbin,HBMAX,'v'); nvbout = elemcnt(hbout,HBMAX,'v'); sbin = getbuff(nsbin); sbout = getbuff(nsbout); vbin = getbuff(nvbin); vbout = getbuff(nvbout); /* get pointers for index fields */ range = getptr2(hbin,HBMAX,sbin,'s','d',RANGE,&fp); if (fp->dim != 1 || strcmp(fp->fname,fp->dname1) != 0) { fprintf(stderr,"radcart_x@alloc: %s is not an index field\n",fp->fname); exit(1); } *nrange = fp->dsize1; xg = getptr(hbout,HBMAX,sbout,'s','d',"x"); yg = getptr(hbout,HBMAX,sbout,'s','d',"y"); zg = getptr(hbout,HBMAX,sbout,'s','d',"z"); /* fill in index fields */ /* get pointers for fields to be cartesianized */ in_dbz = getptr(hbin,HBMAX,vbin,'v','d',DBZ); in_vr = getptr(hbin,HBMAX,vbin,'v','d',VR); in_rsgp= getptr(hbin,HBMAX,vbin,'v','d',RSGP); in_vsgp= getptr(hbin,HBMAX,vbin,'v','d',VSGP); dbz = getptr(hbout,HBMAX,vbout,'v','d',DBZ); rara = getptr(hbout,HBMAX,vbout,'v','d',RR); rsgp = getptr(hbout,HBMAX,vbout,'v','d',RSGP); vsgp = getptr(hbout,HBMAX,vbout,'v','d',VSGP); /* get pointers for new fields */ ux0 = getptr(hbout,HBMAX,vbout,'v','d',"ux0"); uy0 = getptr(hbout,HBMAX,vbout,'v','d',"uy0"); gx = getptr(hbout,HBMAX,vbout,'v','d',"gx"); gy = getptr(hbout,HBMAX,vbout,'v','d',"gy"); eur = getptr(hbout,HBMAX,vbout,'v','d',"eur"); evr = getptr(hbout,HBMAX,vbout,'v','d',"evr"); eup = getptr(hbout,HBMAX,vbout,'v','d',"eup"); evp = getptr(hbout,HBMAX,vbout,'v','d',"evp"); pcount = getptr(hbout,HBMAX,vbout,'v','d',"pcount"); avgtime = getptr(hbout,HBMAX,vbout,'v','d',"avgtime"); /* get pointers for radar location */ lon = getptr(hbin,HBMAX,vbin,'v','d',LON); lat = getptr(hbin,HBMAX,vbin,'v','d',LAT); alt = getptr(hbin,HBMAX,vbin,'v','d',ALT); /* get pointers for time */ if((time = getptr(hbin,HBMAX,vbin,'v','r',TIME)) == NULL) { fprintf(stderr,"radcart_x: warning: %s not available, %s used\n", TIME,ALTTIME); time = getptr(hbin,HBMAX,vbin,'v','r',ALTTIME); } ctime = getptr(hbout,HBMAX,vbout,'v','d',"time"); /* get pointers for horizontal azimuth and elevation */ az = getptr(hbin,HBMAX,vbin,'v','d',AZ); el = getptr(hbin,HBMAX,vbin,'v','d',EL); } /* EOfun_alloc */ /*--------------------------------------------------------------------*/ void fun_fill_index(float x0,float dx,int nx,float *xg, float y0,float dy,int ny,float *yg, float z0,float dz,int nz,float *zg) { float auxvf; int iloop; auxvf = xg[0] = x0; for (iloop = 1; iloop < nx; iloop++) { auxvf += dx; xg[iloop] = auxvf; } auxvf = yg[0] = y0; for (iloop = 1; iloop < ny; iloop++) { auxvf += dy; yg[iloop] = auxvf; } auxvf = zg[0] = z0; for (iloop = 1; iloop < nz; iloop++){ auxvf += dz; zg[iloop] = auxvf; } } /* EOfun_fill_index */ /*--------------------------------------------------------------------*/ void fun_make_s_slice() { /* get input static slice */ getslice(stdin,nsbin,sbin); /* write output static slice */ putslice(stdout,nsbout,sbout); } /* EOfun_make_s_slice */ /*--------------------------------------------------------------------*/ void fun_get_buffer(long num) { /* Get buffers for internal point-dependent fields. The buffers for all output and input fields were gotten @ fun_alloc */ /* get buffer for normalization factors to each ofield */ nfac_dbz = getbuff(num); nfac_rsgp = getbuff(num); nfac_vsgp = getbuff(num); /* get buffer space for matriz A eq.(2.12) --for each point */ a00 = getbuff(num); a01 = getbuff(num); a02 = getbuff(num); a11 = getbuff(num); a12 = getbuff(num); /* get buffer space for vector B eq.(A.12) --for each point */ b0 = getbuff(num); b1 = getbuff(num); /* get buffer space for matriz G eq.(A.12) --for each point */ w00 = getbuff(num); w01 = getbuff(num); w11 = getbuff(num); /* get buffer space for initial time at each point */ initime = getbuff(num); } /* fun_get_buffer */ /*--------------------------------------------------------------------*/ void fun_ini_gridpoints(long num,float *rara, float *rsgp,float *nfac_rsgp, float *vsgp,float *nfac_vsgp, float *dbz,float *nfac_dbz, float *ux0,float *uy0,float *gx, float *gy,float *eur,float *evr, float *eup,float *evp, float *pcount,float *initime, float *avgtime,float *a00,float *a01,float *a02, float *a11,float *a12, float *b0,float *b1,float *w00,float *w01, float *w11) { long lloop; /* INITIALIZE EVERY point-dependent field */ for(lloop = 0; lloop < num; lloop++) { /*loop over grid points */ /* zero a bunch of stuff */ b0[lloop] = b1[lloop] = 0.; a00[lloop] = a01[lloop] = a02[lloop] = a11[lloop] = a12[lloop] = 0.; w00[lloop] = w01[lloop] = w11[lloop] = 0.; pcount[lloop] = 0.; /* number of terms in A-elements */ nfac_dbz[lloop] = 0.; nfac_rsgp[lloop] = 0.; nfac_vsgp[lloop] = 0.; dbz[lloop] = 0.0; rara[lloop] = 0.0; rsgp[lloop] = 0.0; vsgp[lloop] = 0.0; /* make bad all the fields for the dual doppler analysis */ ux0[lloop] = uy0[lloop] = bad; gx[lloop] = gy[lloop] = bad; eur[lloop] = evr[lloop] = bad; eup[lloop] = evp[lloop] = bad; initime[lloop] = bad; /* initial time -- for each point */ avgtime[lloop] = bad; /* time interval-- for each point */ } } /* fun_ini_gridpoints */ /*--------------------------------------------------------------------*/ void fun_assign(int iloc,float weight,float timenow, float in_vr,float refl,float in_rsgp,float in_vsgp) { /* fun_assign -- do the assignment of data values to a grid point */ /* global: rsgp,nfac_rsgp,vsgp,nfac_vsgp,dbz,nfac_dbz,a00,a01,a02,a11,a12,b0,b1, w00,w01,w11,pcount,initime,avgtime */ float fxfx,fxfy,fxfz,fyfy,fyfz,fzfz; float rayvr; /* update dbz and rara -- */ if (fabs(refl) < badlim) { dbz[iloc] += weight*refl; rara[iloc] += weight*fun_rara(refl); /* update the normalization factors */ nfac_dbz[iloc] += weight; } /* update rsgp -- */ if (fabs(in_rsgp) < badlim ) { rsgp[iloc] += weight*in_rsgp; nfac_rsgp[iloc] += weight; } /* update vsgp -- */ if (fabs(in_vsgp) < badlim ) { vsgp[iloc] += weight*in_vsgp; nfac_vsgp[iloc] += weight; } /* check for a good radial velocity */ if (fabs(in_vr) < badlim) { rayvr=weight*in_vr; /* update stuff for dual doppler analysis */ fxfx =weight*orax*orax; fxfy =weight*orax*oray; fxfz =weight*orax*oraz; fyfy =weight*oray*oray; fyfz =weight*oray*oraz; /* update sums on left hand side */ a00[iloc] += fxfx; a01[iloc] += fxfy; a02[iloc] += fxfz; a11[iloc] += fyfy; a12[iloc] += fyfz; /* update sums on the right hand side */ b0[iloc] += orax*rayvr; b1[iloc] += oray*rayvr; /* update error terms */ w00[iloc] += weight*fxfx; w01[iloc] += weight*fxfy; w11[iloc] += weight*fyfy; /* update counter of terms in vr average */ pcount[iloc] += 1.0; /* do time stuff */ if (fabs(initime[iloc]) < badlim) avgtime[iloc] = timenow; /* update avgtime */ else /* means 1st time @ iloc, so set avgtime & initime */ initime[iloc] = avgtime[iloc] = timenow; } /* else Skip this input radial velocity from the vr-average */ } /* EOfun_assign */ /*--------------------------------------------------------------------*/ #define EARTH_R 6370.0 void fun_rloc(float lon0,float lat0,float alt0,float *xl,float *yl,float *zl) { /* Radar location form the origin of the coordinate system */ /*global: *lat,*lon,*alt */ *xl = (EARTH_R)*cos(*lat*D2R)*((*lon - lon0)*D2R); *yl = (EARTH_R)*((*lat - lat0)*D2R); *zl = (*alt - alt0); } /* EOfun_radloc_g */ /*--------------------------------------------------------------------*/ void fun_ora(float *orax,float *oray,float *oraz) { /* angular factors from the radar */ /*global: *el,*az */ *orax = cos(*el*D2R); *oray = cos(*az*D2R)*(*orax); *orax = sin(*az*D2R)*(*orax); *oraz = sin(*el*D2R); } /* EOfun_ora */ /*--------------------------------------------------------------------*/ #define EARTH_D 12740.0 void fun_ploc(float range,float xl, float yl, float zl, float deltax, float deltay, float *xloc,float *yloc,float *zloc) { float delx,dely,delz; /* get locations of range point relative to radar */ delx = range*orax; dely = range*oray; delz = range*oraz; /* do corrections for earth's curvature */ delz += (delx*delx + dely*dely)/EARTH_D; /* get location relative to coordinate system origin */ *xloc = xl + delx; *yloc = yl + dely; *zloc = zl + delz; /* translate all points to their location at the common time of reftime */ *xloc = *xloc - (deltax); *yloc = *yloc - (deltay); } /* EOfun_ploc */ /*--------------------------------------------------------------------*/ void fun_weight(int nx,int ny,int nz,int ix,int iy, int iz, float cx, float cy, float cz,int *np,float *we) { int jx,jy,jz; float bx,by,bz; FRACT(ix,cx,bx); FRACT(iy,cy,by); FRACT(iz,cz,bz); jx = ix + 1; jy = iy + 1; jz = iz + 1; np[0] = iz + nz*iy + nz*ny*ix; we[0] = bx*by*bz; np[1] = iz + nz*iy + nz*ny*jx; we[1] = cx*by*bz; np[2] = iz + nz*jy + nz*ny*ix; we[2] = bx*cy*bz; np[3] = jz + nz*iy + nz*ny*ix; we[3] = bx*by*cz; np[4] = iz + nz*jy + nz*ny*jx; we[4] = cx*cy*bz; np[5] = jz + nz*iy + nz*ny*jx; we[5] = cx*by*cz; np[6] = jz + nz*jy + nz*ny*ix; we[6] = bx*cy*cz; np[7] = jz + nz*jy + nz*ny*jx; we[7] = cx*cy*cz; } /* EOfun_weight */ /*--------------------------------------------------------------------*/ #define DBZMIN -100.0 #define LASTBADGATE 0 void fun_filldbz_thresh(int flag_p,float dbz0,float r0, float *range,int maxrange, float *dbz_thresh) { /* set up dbz thresholds */ float auxdbm; int gate; int jrange; /*I expect no more than 512 gates=maxrange */ /* bad value is assumed for power at range = 0 */ dbz_thresh[0] = bad; for (gate = 1; gate <= LASTBADGATE; gate++) dbz_thresh[gate] = bad; if (flag_p == 1) { auxdbm = dbz0 - 20.0*log10(r0); for (jrange = 1; jrange < maxrange; jrange++) { dbz_thresh[jrange] = auxdbm + 20.0*log10(range[jrange]); } } else /* no power thresholding */ for (jrange = 1; jrange < maxrange; jrange++) dbz_thresh[jrange] = DBZMIN; } /* EOfun_filldbz_thresh */ /*--------------------------------------------------------------------*/ #define MINDET 1.e-10 void fun_dualdoppa(float a00,float a01, float a02, float a11, float a12, float b0, float b1,float w00, float w01, float w11, float pcount,long loop) { /* globals: ux0,uy0,gx,gy,eur,evr,eup,evp,icount */ float det1,aux; /* if there are no data at the point return to main */ if (pcount <= 0.) return; /* compute the determinant of the system */ det1 = a00*a11 - a01*a01; /* if det1 is ok make the analysis */ if (fabs(det1) > MINDET) { /* calculate the variances */ /* first, assuming error in wp is zero */ eur[loop] = (a11*a11*w00 - 2.*a11*a01*w01 + a01*a01*w11)/(det1*det1); evr[loop] = (a01*a01*w00 - 2.*a01*a00*w01 + a00*a00*w11)/(det1*det1); /* next, assuming error in vr is zero */ aux = (a02*a11 - a12*a01)/det1; eup[loop] = aux*aux; aux = (a12*a00 - a02*a01)/det1; evp[loop] = aux*aux; /* this test is here because ``bad'' values of these quantities sometimes occur when ``good'' values of ux0, uy0, etc. are calculated */ if (fabs(eur[loop]) < badlim && fabs(evr[loop]) < badlim && fabs(eup[loop]) < badlim && fabs(evp[loop]) < badlim) { /* Dual Doppler velocities */ ux0[loop] = (b0*a11 - a01*b1)/det1; uy0[loop] = (a00*b1 - b0*a01)/det1; /* now calculate the gx and gy */ gx[loop] = (a01*a12 - a11*a02)/det1; gy[loop] = (a01*a02 - a00*a12)/det1; } } } /* EOfun_dualdoppa */ /*----------------------------------------------------------------*/ float fun_rara(float Z) { /* typically there is a ZvsR relationship of the form Z=alpha*R^beta with alpha and beta reflecting a particular climatology of a region as well as the precipitation regime -- namely convective or stratiform. At this point I don't see how to tell apart convective from stratiform maybe looking at the reflectivity gradient along the ray will help in this regard. We must keep in mind the attenuation!! Attenuation oly reduces reflectivity. The Following alpha& beta are for Convective Toga-Coare conditions ALPHA = 139.0; BETA = 1.43; Mean values for Epic provided by Darrel Baumgardner Centro de Ciencias de la Atmosfera, Universidad Nacional Autonoma de Mexico ALPHA = 153; BETA = 1.53; */ float ALPHA,BETA,rr; rr = pow(Z/alpha,1.0/beta); return(rr); }