X-Git-Url: http://pileus.org/git/?p=~andy%2Frsl;a=blobdiff_plain;f=volume.c;h=a30368b4796258f72416061a22602b32a54f791d;hp=2ff6aa9650a06e1e38a6be432d9f725e0ff58a01;hb=d08a7f8a699a044bc4ac5f93917aa7f6c463923b;hpb=0a59ff8412de3e114b7b5fd081a5fd39102542c1 diff --git a/volume.c b/volume.c index 2ff6aa9..a30368b 100644 --- a/volume.c +++ b/volume.c @@ -114,7 +114,7 @@ extern int radar_verbose_flag; float DZ_F(Range x) { if (x >= F_OFFSET) /* This test works when Range is unsigned. */ - return (((float)x-F_OFFSET)/F_FACTOR - F_DZ_RANGE_OFFSET); /* Default wsr88d. */ + return (((float)x-F_OFFSET)/F_FACTOR - F_DZ_RANGE_OFFSET); /* Default wsr88d. */ if (x == 0) return BADVAL; if (x == 1) return RFVAL; if (x == 2) return APFLAG; @@ -125,9 +125,9 @@ float DZ_F(Range x) { float VR_F(Range x) { float val; if (x >= F_OFFSET) { /* This test works when Range is unsigned. */ - val = (((float)x-F_OFFSET)/F_FACTOR - F_VR_OFFSET); /* Default wsr88d coding. */ - /* fprintf(stderr, "x=%d, val=%f\n", x, val); */ - return val; + val = (((float)x-F_OFFSET)/F_FACTOR - F_VR_OFFSET); /* Default wsr88d coding. */ + /* fprintf(stderr, "x=%d, val=%f\n", x, val); */ + return val; } if (x == 0) return BADVAL; if (x == 1) return RFVAL; @@ -139,8 +139,8 @@ float VR_F(Range x) { float DR_F(Range x) { /* Differential reflectivity */ float val; if (x >= F_OFFSET) { /* This test works when Range is unsigned. */ - val = (((float)x-F_OFFSET)/F_DR_FACTOR - F_DR_OFFSET); - return val; + val = (((float)x-F_OFFSET)/F_DR_FACTOR - F_DR_OFFSET); + return val; } if (x == 0) return BADVAL; if (x == 1) return RFVAL; @@ -151,7 +151,7 @@ float DR_F(Range x) { /* Differential reflectivity */ float LR_F(Range x) {/* From MCTEX */ if (x >= F_OFFSET) /* This test works when Range is unsigned. */ - return (float) (x - 250.)/6.; + return (float) (x - 250.)/6.; if (x == 0) return BADVAL; if (x == 1) return RFVAL; if (x == 2) return APFLAG; @@ -228,13 +228,23 @@ float KD_F(Range x) return (x-32768.)/100.; } -float NP_F(Range x) { /* Normalized Coherent Power (DORADE) */ +/* Normalized Coherent Power (DORADE) */ +float NP_F(Range x) +{ + if (x == 0) return BADVAL; + return (float)(x - 1) / 100.; +} + +/* Standard Deviation (for Dual-pole QC testing.) */ +float SD_F(Range x) +{ if (x == 0) return BADVAL; return (float)x / 100.; } /* Signal Quality Index */ -float SQ_F(Range x) { +float SQ_F(Range x) +{ if (x == 0) return BADVAL; return (float)(x-1) / 65533.; } @@ -360,28 +370,34 @@ Range PH_INVF(float x) { /* KD_INVF for 1 or 2 byte data. */ Range KD_INVF(float x) { if (x == BADVAL) return (Range)0; -/****** Commented-out code for 1-byte Sigmet native data format. + return (Range)(x * 100. + 32768. + 0.5); +/****** Old code for 1-byte Sigmet native data format commented-out: if (x == RFVAL) return (Range)1; if (x == APFLAG) return (Range)2; if (x == NOECHO) return (Range)3; if (rsl_kdp_wavelen == 0.0) return (Range)0; if (x < 0) { - x = 127 - - 126 * (log((double)-x) - log((double)(0.25/rsl_kdp_wavelen))) / - log((double)600.0) + - 0.5; + x = 127 - + 126 * (log((double)-x) - log((double)(0.25/rsl_kdp_wavelen))) / + log((double)600.0) + + 0.5; } else if (x > 0) { - x = 129 + - 126 * (log((double)x) - log((double)0.25/rsl_kdp_wavelen)) / - log((double)600.0) + - 0.5; + x = 129 + + 126 * (log((double)x) - log((double)0.25/rsl_kdp_wavelen)) / + log((double)600.0) + + 0.5; } else { - x = 128; + x = 128; } x += F_OFFSET; ******/ - return (Range)(x * 100. + 32768. + 0.5); - +} + +/* Standard Deviation (for Dual-pole QC testing.) */ +Range SD_INVF(float x) +{ + if (x == BADVAL) return (Range)0; + return (Range)(x * 100.); } /* Signal Quality Index */ @@ -391,11 +407,11 @@ Range SQ_INVF(float x) return (Range)(x * 65533. + 1. +.5); } - -Range NP_INVF(float x) /* Normalized Coherent Power (DORADE) */ +/* Normalized Coherent Power (DORADE) */ +Range NP_INVF(float x) { if (x == BADVAL) return (0); - return (Range)(x * 100.); + return (Range)(x * 100. + 1.); } Range TI_INVF(float x) /* MCTEX */ @@ -497,9 +513,9 @@ typedef struct { #define STATIC STATIC int RSL_max_sweeps = 0; /* Initial allocation for sweep_list. - * RSL_new_sweep will allocate the space first - * time around. - */ + * RSL_new_sweep will allocate the space first + * time around. + */ STATIC int RSL_nsweep_addr = 0; /* A count of sweeps in the table. */ STATIC Sweep_list *RSL_sweep_list = NULL; STATIC int RSL_nextents = 0; @@ -516,7 +532,7 @@ void FREE_HASH_TABLE(Hash_table *table) int i; if (table == NULL) return; for (i=0; inindexes; i++) - FREE_HASH_NODE(table->indexes[i]); /* A possible linked list of Rays. */ + FREE_HASH_NODE(table->indexes[i]); /* A possible linked list of Rays. */ free(table->indexes); free(table); } @@ -527,7 +543,7 @@ void REMOVE_SWEEP(Sweep *s) int j; /* Find where it goes, split the list and slide the tail down one. */ for (i=0; i= RSL_max_sweeps) { /* Current list is too small. */ - RSL_nextents++; - new_list = (Sweep_list *) calloc(100*RSL_nextents, sizeof(Sweep_list)); - if (new_list == NULL) { - perror("INSERT_SWEEP"); - exit(2); - } + RSL_nextents++; + new_list = (Sweep_list *) calloc(100*RSL_nextents, sizeof(Sweep_list)); + if (new_list == NULL) { + perror("INSERT_SWEEP"); + exit(2); + } /* Copy the old list to the new one. */ - for (i=0; ii; j--) - RSL_sweep_list[j] = RSL_sweep_list[j-1]; + RSL_sweep_list[j] = RSL_sweep_list[j-1]; RSL_sweep_list[i].s_addr = s; RSL_sweep_list[i].hash = NULL; @@ -581,7 +597,7 @@ int SWEEP_INDEX(Sweep *s) /* Simple linear search; but this will be a binary search. */ int i; for (i=0; iray = (Ray **) calloc(max_rays, sizeof(Ray*)); if (s->ray == NULL) perror("RSL_new_sweep, Ray*"); s->h.nrays = max_rays; /* A default setting. */ + s->h.elev = -999.; + s->h.azimuth = -999.; return s; } @@ -633,7 +651,7 @@ Sweep *RSL_clear_sweep(Sweep *s) int i; if (s == NULL) return s; for (i=0; ih.nrays; i++) { - RSL_clear_ray(s->ray[i]); + RSL_clear_ray(s->ray[i]); } return s; } @@ -642,7 +660,7 @@ Volume *RSL_clear_volume(Volume *v) int i; if (v == NULL) return v; for (i=0; ih.nsweeps; i++) { - RSL_clear_sweep(v->sweep[i]); + RSL_clear_sweep(v->sweep[i]); } return v; } @@ -664,7 +682,7 @@ void RSL_free_sweep(Sweep *s) int i; if (s == NULL) return; for (i=0; ih.nrays; i++) { - RSL_free_ray(s->ray[i]); + RSL_free_ray(s->ray[i]); } if (s->ray) free(s->ray); REMOVE_SWEEP(s); /* Remove from internal Sweep list. */ @@ -676,9 +694,9 @@ void RSL_free_volume(Volume *v) if (v == NULL) return; for (i=0; ih.nsweeps; i++) - { - RSL_free_sweep(v->sweep[i]); - } + { + RSL_free_sweep(v->sweep[i]); + } if (v->sweep) free(v->sweep); free(v); } @@ -711,7 +729,7 @@ Sweep *RSL_copy_sweep(Sweep *s) n_sweep->h = s->h; for (i=0; ih.nrays; i++) { - n_sweep->ray[i] = RSL_copy_ray(s->ray[i]); + n_sweep->ray[i] = RSL_copy_ray(s->ray[i]); } return n_sweep; } @@ -728,7 +746,7 @@ Volume *RSL_copy_volume(Volume *v) new_vol->h = v->h; for (i=0; ih.nsweeps; i++) { - new_vol->sweep[i] = RSL_copy_sweep(v->sweep[i]); + new_vol->sweep[i] = RSL_copy_sweep(v->sweep[i]); } return new_vol; } @@ -817,9 +835,9 @@ Ray *RSL_get_next_ccwise_ray(Sweep *s, Ray *ray) double cwise_angle_diff(float x,float y) { /* Returns the clockwise angle difference of x to y. - * If x = 345 and y = 355 return 10. - * If x = 345 and y = 335 return 350 - */ + * If x = 345 and y = 355 return 10. + * If x = 345 and y = 335 return 350 + */ double d; d = (double)(y - x); @@ -836,9 +854,9 @@ double cwise_angle_diff(float x,float y) double ccwise_angle_diff(float x,float y) { /* Returns the counterclockwise angle differnce of x to y. - * If x = 345 and y = 355 return 350. - * If x = 345 and y = 335 return 10 - */ + * If x = 345 and y = 355 return 350. + * If x = 345 and y = 335 return 10 + */ double d; d = (double)(x - y); @@ -862,48 +880,48 @@ Azimuth_hash *the_closest_hash(Azimuth_hash *hash, float ray_angle) if (hash == NULL) return NULL; /* Set low pointer to hash index with ray angle just below - * requested angle and high pointer to just above requested - * angle. - */ + * requested angle and high pointer to just above requested + * angle. + */ /* set low and high pointers to initial search locations*/ low = hash; high = hash->ray_high; /* Search until clockwise angle to high is less then clockwise - * angle to low. - */ + * angle to low. + */ clow = cwise_angle_diff(ray_angle,low->ray->h.azimuth); chigh = cwise_angle_diff(ray_angle,high->ray->h.azimuth); cclow = ccwise_angle_diff(ray_angle,low->ray->h.azimuth); while((chigh > clow) && (clow != 0)) - { - if (clow < cclow) - { - low = low->ray_low; - high = low->ray_high; /* Not the same low as line before ! */ - } - else - { - low = low->ray_high; - high = low->ray_high; /* Not the same low as line before ! */ - } - - clow = cwise_angle_diff(ray_angle,low->ray->h.azimuth); - chigh = cwise_angle_diff(ray_angle,high->ray->h.azimuth); - cclow = ccwise_angle_diff(ray_angle,low->ray->h.azimuth); - } + { + if (clow < cclow) + { + low = low->ray_low; + high = low->ray_high; /* Not the same low as line before ! */ + } + else + { + low = low->ray_high; + high = low->ray_high; /* Not the same low as line before ! */ + } + + clow = cwise_angle_diff(ray_angle,low->ray->h.azimuth); + chigh = cwise_angle_diff(ray_angle,high->ray->h.azimuth); + cclow = ccwise_angle_diff(ray_angle,low->ray->h.azimuth); + } if(chigh <= cclow) - { - return high; - } + { + return high; + } else - { - return low; - } + { + return low; + } } @@ -1025,8 +1043,8 @@ int hash_bin(Hash_table *table,float angle) * why bother? */ while(table->indexes[hash] == NULL) { - hash++; - if(hash >= table->nindexes) hash = 0; + hash++; + if(hash >= table->nindexes) hash = 0; } return hash; @@ -1038,13 +1056,13 @@ Hash_table *hash_table_for_sweep(Sweep *s) i = SWEEP_INDEX(s); if (i==-1) { /* Obviously, an unregistered sweep. Most likely the - * result of pointer assignments. - */ - i = INSERT_SWEEP(s); + * result of pointer assignments. + */ + i = INSERT_SWEEP(s); } if (RSL_sweep_list[i].hash == NULL) { /* First time. Construct the table. */ - RSL_sweep_list[i].hash = construct_sweep_hash_table(s); + RSL_sweep_list[i].hash = construct_sweep_hash_table(s); } return RSL_sweep_list[i].hash; @@ -1061,7 +1079,7 @@ Ray *RSL_get_closest_ray_from_sweep(Sweep *s,float ray_angle, float limit) /* * Return closest Ray in Sweep within limit (angle) specified * in parameter list. Assume PPI mode. - */ + */ int hindex; Hash_table *hash_table; Azimuth_hash *closest; @@ -1078,8 +1096,8 @@ Ray *RSL_get_closest_ray_from_sweep(Sweep *s,float ray_angle, float limit) closest = the_closest_hash(hash_table->indexes[hindex],ray_angle); /* Is closest ray within limit parameter ? If - * so return ray, else return NULL. - */ + * so return ray, else return NULL. + */ close_diff = angle_diff(ray_angle,closest->ray->h.azimuth); @@ -1134,13 +1152,13 @@ float RSL_get_value_from_ray(Ray *ray, float r) if (ray == NULL) return BADVAL; if(ray->h.gate_size == 0) - { - if(radar_verbose_flag) - { - fprintf(stderr,"RSL_get_value_from_ray: ray->h.gate_size == 0\n"); - } - return BADVAL; - } + { + if(radar_verbose_flag) + { + fprintf(stderr,"RSL_get_value_from_ray: ray->h.gate_size == 0\n"); + } + return BADVAL; + } /* range_bin1 is range to center of first bin */ bin_index = (int)(((rm - ray->h.range_bin1)/ray->h.gate_size) + 0.5); @@ -1254,10 +1272,10 @@ Ray *RSL_get_ray_above(Volume *v, Ray *current_ray) i++; while( i < v->h.nsweeps) - { - if(v->sweep[i] != NULL) break; - i++; - } + { + if(v->sweep[i] != NULL) break; + i++; + } if(i >= v->h.nsweeps) return NULL; @@ -1283,10 +1301,10 @@ Ray *RSL_get_ray_below(Volume *v, Ray *current_ray) i--; while( i >= 0) - { - if(v->sweep[i] != NULL) break; - i--; - } + { + if(v->sweep[i] != NULL) break; + i--; + } if(i < 0) return NULL; @@ -1332,13 +1350,13 @@ Ray *RSL_get_first_ray_of_sweep(Sweep *s) smallest_ray_num = 9999999; if (s == NULL) return r; for (i=0; ih.nrays; i++) - if (s->ray[i]) { - if (s->ray[i]->h.ray_num <= 1) return s->ray[i]; - if (s->ray[i]->h.ray_num < smallest_ray_num) { - r = s->ray[i]; - smallest_ray_num = r->h.ray_num; - } - } + if (s->ray[i]) { + if (s->ray[i]->h.ray_num <= 1) return s->ray[i]; + if (s->ray[i]->h.ray_num < smallest_ray_num) { + r = s->ray[i]; + smallest_ray_num = r->h.ray_num; + } + } return r; } @@ -1361,7 +1379,7 @@ Sweep *RSL_get_first_sweep_of_volume(Volume *v) int i; if (v == NULL) return NULL; for (i=0; ih.nsweeps; i++) - if (RSL_get_first_ray_of_sweep(v->sweep[i])) return v->sweep[i]; + if (RSL_get_first_ray_of_sweep(v->sweep[i])) return v->sweep[i]; return NULL; } @@ -1384,6 +1402,8 @@ int rsl_qfield[MAX_RADAR_VOLUMES] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1 }; @@ -1438,26 +1458,26 @@ void RSL_select_fields(char *field_type, ...) if (radar_verbose_flag) fprintf(stderr,"Selected fields for ingest:"); while (c_field) { - /* CHECK EACH FIELD. This is a fancier case statement than C provides. */ - if (radar_verbose_flag) fprintf(stderr," %s", c_field); - if (strcasecmp(c_field, "all") == 0) { - for (i=0; i> specified.\n", c_field); - } - } - c_field = va_arg(ap, char *); + /* CHECK EACH FIELD. This is a fancier case statement than C provides. */ + if (radar_verbose_flag) fprintf(stderr," %s", c_field); + if (strcasecmp(c_field, "all") == 0) { + for (i=0; i> specified.\n", c_field); + } + } + c_field = va_arg(ap, char *); } if (radar_verbose_flag) fprintf(stderr,"\n"); @@ -1497,15 +1517,15 @@ int rsl_query_field(char *c_field) RSL_invf_list[0] = RSL_invf_list[0]; for (i=0; i> specified.\n", c_field); + * properly written ingest code. + */ + fprintf(stderr, "rsl_query_field: Invalid field name <<%s>> specified.\n", c_field); } /* 'i' is the index. Is it set? */ @@ -1515,8 +1535,8 @@ int rsl_query_field(char *c_field) /* Could be static and force use of 'rsl_query_sweep' */ int *rsl_qsweep = NULL; /* If NULL, then read all sweeps. Otherwise, - * read what is on the list. - */ + * read what is on the list. + */ #define RSL_MAX_QSWEEP 500 /* It'll be rediculious to have more. :-) */ int rsl_qsweep_max = RSL_MAX_QSWEEP; @@ -1555,35 +1575,35 @@ void RSL_read_these_sweeps(char *csweep, ...) if (radar_verbose_flag) fprintf(stderr,"Selected sweeps for ingest:"); - for (;c_sweep; c_sweep = va_arg(ap, char *)) - { - /* CHECK EACH FIELD. This is a fancier case statement than C provides. */ - if (radar_verbose_flag) fprintf(stderr," %s", c_sweep); - if (strcasecmp(c_sweep, "all") == 0) { - for (i=0; i RSL_MAX_QSWEEP) { - if (radar_verbose_flag) fprintf(stderr,"\nRSL_read_these_sweeps: parameter %s not in [0,%d). Ignoring.\n", c_sweep, RSL_MAX_QSWEEP); - continue; - } - - if (isweep > rsl_qsweep_max) rsl_qsweep_max = isweep; - rsl_qsweep[isweep] = 1; - } + for (;c_sweep; c_sweep = va_arg(ap, char *)) + { + /* CHECK EACH FIELD. This is a fancier case statement than C provides. */ + if (radar_verbose_flag) fprintf(stderr," %s", c_sweep); + if (strcasecmp(c_sweep, "all") == 0) { + for (i=0; i RSL_MAX_QSWEEP) { + if (radar_verbose_flag) fprintf(stderr,"\nRSL_read_these_sweeps: parameter %s not in [0,%d). Ignoring.\n", c_sweep, RSL_MAX_QSWEEP); + continue; + } + + if (isweep > rsl_qsweep_max) rsl_qsweep_max = isweep; + rsl_qsweep[isweep] = 1; + } } if (radar_verbose_flag) fprintf(stderr,"\n"); @@ -1639,9 +1659,9 @@ void RSL_add_dbz_offset_to_ray(Ray *r, float dbz_offset) if (r == NULL) return; for (ibin=0; ibinh.nbins; ibin++) { - val = r->h.f(r->range[ibin]); - if ( val >= (float)NOECHO ) continue; /* Invalid value */ - r->range[ibin] = r->h.invf(val + dbz_offset); + val = r->h.f(r->range[ibin]); + if ( val >= (float)NOECHO ) continue; /* Invalid value */ + r->range[ibin] = r->h.invf(val + dbz_offset); } } @@ -1655,7 +1675,7 @@ void RSL_add_dbz_offset_to_sweep(Sweep *s, float dbz_offset) int iray; if (s == NULL) return; for (iray=0; irayh.nrays; iray++) - RSL_add_dbz_offset_to_ray(s->ray[iray], dbz_offset); + RSL_add_dbz_offset_to_ray(s->ray[iray], dbz_offset); } /*********************************************************************/ @@ -1668,5 +1688,5 @@ void RSL_add_dbz_offset_to_volume(Volume *v, float dbz_offset) int isweep; if (v == NULL) return; for (isweep=0; isweeph.nsweeps; isweep++) - RSL_add_dbz_offset_to_sweep(v->sweep[isweep], dbz_offset); + RSL_add_dbz_offset_to_sweep(v->sweep[isweep], dbz_offset); }