3 This is the TRMM Office Radar Software Library.
4 Copyright (C) 1996, 1997
6 Space Applications Corporation
9 This library is free software; you can redistribute it and/or
10 modify it under the terms of the GNU Library General Public
11 License as published by the Free Software Foundation; either
12 version 2 of the License, or (at your option) any later version.
14 This library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 Library General Public License for more details.
19 You should have received a copy of the GNU Library General Public
20 License along with this library; if not, write to the Free
21 Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
24 * Volume functions coded in this file:
26 * Default DZ_F and DZ_INVF. For VR, SW, CZ, ZT, DR and LR too.
27 * Volume *RSL_new_volume(int max_sweeps);
28 * Sweep *RSL_new_sweep(int max_rays);
29 * Ray *RSL_new_ray(int max_bins);
30 * Ray *RSL_clear_ray(Ray *r);
31 * Sweep *RSL_clear_sweep(Sweep *s);
32 * Volume *RSL_clear_volume(Volume *v);
33 * void RSL_free_ray(Ray *r);
34 * void RSL_free_sweep(Sweep *s);
35 * void RSL_free_volume(Volume *v);
36 * Ray *RSL_copy_ray(Ray *r);
37 * Sweep *RSL_copy_sweep(Sweep *s);
38 * Volume *RSL_copy_volume(Volume *v);
39 * Ray *RSL_get_ray_from_sweep(Sweep *s, float azim);
40 * float RSL_get_value_from_sweep(Sweep *s, float azim, float r);
41 * float RSL_get_value_from_ray(Ray *ray, float r);
42 * float RSL_get_value_at_h(Volume *v, float azim, float grnd_r, float h);
43 * Sweep *RSL_get_sweep(Volume *v, float elev);
44 * Ray *RSL_get_ray(Volume *v, float elev, float azimuth);
45 * float RSL_get_value(Volume *v, float elev, float azimuth, float range);
46 * Ray *RSL_get_ray_above(Volume *v, Ray *current_ray);
47 * Ray *RSL_get_ray_below(Volume *v, Ray *current_ray);
48 * Ray *RSL_get_matching_ray(Volume *v, Ray *ray);
49 * int RSL_get_sweep_index_from_volume
51 * See image_gen.c for the Volume image generation functions.
53 * See doc/rsl_index.html for HTML formatted documentation.
55 * All routines herein coded, unless otherwise stated, by:
57 * Space Applications Corporation
60 * Dennis Flanigan, Jr.
61 * flanigan@lance.colostate.edu
68 int strcasecmp(const char *s1, const char *s2);
73 #define bin_azimuth(x, dx) (float)((float)x/dx)
74 #define bin_elevation(x, dx) (float)((float)x/dx)
75 #define bin_range(x, dx) (float)((float)x/dx)
77 extern int radar_verbose_flag;
79 /* Internal storage conversion functions. These may be any conversion and
80 * may be dynamically defined; based on the input data conversion. If you
81 * change any of the reserved values, ie. values that cannot be converted
82 * from/to internal storage (BADVAL, APFLAG, etc), be sure to change all
83 * functions, both XX_F and XX_INVF to handle it. It is best to have the
84 * reserved values stored at 0, 1, 2, on up. That way you merely need to
85 * provide an offset to the actual conversion function. See 'rsl.h'
86 * for the definition of the reserved values. Currently: BADVAL, RFVAL,
89 * The conversion functions may NOT be macros.
92 #ifdef USE_TWO_BYTE_PRECISION
93 #define F_FACTOR 100.0
94 #define F_DR_FACTOR 1000.0
95 #define F_DZ_RANGE_OFFSET 50
98 #define F_DR_FACTOR 10.0
99 #define F_DZ_RANGE_OFFSET 32
102 /* #define F_VR_OFFSET 63.5 */
103 #define F_VR_OFFSET 127.0
104 #define F_DR_OFFSET 12.0
106 /* IMPORTANT: This is the offset from reserved values. This
107 * number must be exactly (or >=) the number of
108 * reserved values in XX_F and XX_INVF.
110 * You must change nsig_to_radar.c where F_OFFSET is used for optimization.
115 float DZ_F(Range x) {
116 if (x >= F_OFFSET) /* This test works when Range is unsigned. */
117 return (((float)x-F_OFFSET)/F_FACTOR - F_DZ_RANGE_OFFSET); /* Default wsr88d. */
118 if (x == 0) return BADVAL;
119 if (x == 1) return RFVAL;
120 if (x == 2) return APFLAG;
121 if (x == 3) return NOECHO;
122 return BADVAL; /* Can't get here, but quiets the compiler. */
125 float VR_F(Range x) {
127 if (x >= F_OFFSET) { /* This test works when Range is unsigned. */
128 val = (((float)x-F_OFFSET)/F_FACTOR - F_VR_OFFSET); /* Default wsr88d coding. */
129 /* fprintf(stderr, "x=%d, val=%f\n", x, val); */
132 if (x == 0) return BADVAL;
133 if (x == 1) return RFVAL;
134 if (x == 2) return APFLAG;
135 if (x == 3) return NOECHO;
136 return BADVAL; /* Can't get here, but quiets the compiler. */
139 float DR_F(Range x) { /* Differential reflectivity */
141 if (x >= F_OFFSET) { /* This test works when Range is unsigned. */
142 val = (((float)x-F_OFFSET)/F_DR_FACTOR - F_DR_OFFSET);
145 if (x == 0) return BADVAL;
146 if (x == 1) return RFVAL;
147 if (x == 2) return APFLAG;
148 if (x == 3) return NOECHO;
149 return BADVAL; /* Can't get here, but quiets the compiler. */
152 float LR_F(Range x) {/* From MCTEX */
153 if (x >= F_OFFSET) /* This test works when Range is unsigned. */
154 return (float) (x - 250.)/6.;
155 if (x == 0) return BADVAL;
156 if (x == 1) return RFVAL;
157 if (x == 2) return APFLAG;
158 if (x == 3) return NOECHO;
162 float HC_F(Range x) { /* HydroClass (Sigmet) */
163 if (x == 0) return BADVAL;
167 /****************************
168 Sigmet RhoHV : one_byte values
169 > RohHV = sqrt((N-1)/253)
176 *******************************/
177 float RH_F(Range x) {
178 if (x == 0) return BADVAL;
179 /* return (float)(sqrt((double)((x-1.0)/253.0))); */
180 return (float)(x-1) / 65533.;
183 /*****************************
184 Sigmet PhiDP : one_byte values
185 > PhiDP (mod 180) = 180 * ((N-1)/254)
186 > The range is from 0 to 180 degrees in steps of 0.71 as follows.
192 ******************************/
193 float PH_F(Range x) {
194 if (x == 0) return BADVAL;
195 /*return (float)(180.0*((x-1.0)/254.0));*/
196 return (360.*(x-1.))/65534.;
199 /* TODO: Should this be 5. cm. instead of 0.5? Or maybe 10. cm.? */
200 float rsl_kdp_wavelen = 0.5; /* Default radar wavelen = .5 cm. See
201 * nsig_to_radar.c for the code that sets this.
203 /* KD_F for 1 or 2 byte. */
206 /****** Commented-out code for 1-byte Sigmet native data format.
209 if (rsl_kdp_wavelen == 0.0) return BADVAL;
212 -0.25 * pow((double)600.0,(double)((127-x)/126.0))
216 0.25 * pow((double)600.0,(double)((x-129)/126.0))
221 if (x == 1) return RFVAL;
222 if (x == 2) return APFLAG;
223 if (x == 3) return NOECHO;
227 if (x == 0) return BADVAL;
228 return (x-32768.)/100.;
231 /* Normalized Coherent Power (DORADE) */
234 if (x == 0) return BADVAL;
235 return (float)(x - 1) / 100.;
238 /* Standard Deviation (for Dual-pole QC testing.) */
241 if (x == 0) return BADVAL;
242 return (float)x / 100.;
245 /* Signal Quality Index */
248 if (x == 0) return BADVAL;
249 return (float)(x-1) / 65533.;
254 if (x >= F_OFFSET) return (float)(x);
255 if (x == 0) return BADVAL;
256 if (x == 1) return RFVAL;
257 if (x == 2) return APFLAG;
258 if (x == 3) return NOECHO;
262 float SW_F(Range x) { return VR_F(x); }
263 float CZ_F(Range x) { return DZ_F(x); }
264 float ZT_F(Range x) { return DZ_F(x); }
265 float ZD_F(Range x) { return DR_F(x); } /* Differential reflectivity */
266 float CD_F(Range x) { return DR_F(x); } /* Differential reflectivity */
267 float XZ_F(Range x) { return DZ_F(x); }
268 float MZ_F(Range x) { return (float)x; } /* DZ Mask */
269 float MD_F(Range x) { return MZ_F(x); } /* ZD Mask */
270 float ZE_F(Range x) { return DZ_F(x); }
271 float VE_F(Range x) { return VR_F(x); }
272 float DM_F(Range x) { return DZ_F(x); }
273 float DX_F(Range x) { return DZ_F(x); }
274 float CH_F(Range x) { return DZ_F(x); }
275 float AH_F(Range x) { return DZ_F(x); }
276 float CV_F(Range x) { return DZ_F(x); }
277 float AV_F(Range x) { return DZ_F(x); }
278 float VS_F(Range x) { return VR_F(x); }
279 float VL_F(Range x) { return VR_F(x); }
280 float VG_F(Range x) { return VR_F(x); }
281 float VT_F(Range x) { return VR_F(x); }
282 float VC_F(Range x) { return VR_F(x); }
287 /* Unfortunately, floats are stored differently than ints/shorts. So,
288 * we cannot simply set up a switch statement and we must test for
289 * all the special cases first. We must test for exactness.
291 Range DZ_INVF(float x)
293 if (x == BADVAL) return (Range)0;
294 if (x == RFVAL) return (Range)1;
295 if (x == APFLAG) return (Range)2;
296 if (x == NOECHO) return (Range)3;
297 if (x < -F_DZ_RANGE_OFFSET) return (Range)0;
298 return (Range)(F_FACTOR*(x+F_DZ_RANGE_OFFSET)+.5 + F_OFFSET); /* Default wsr88d. */
301 Range VR_INVF(float x)
303 if (x == BADVAL) return (Range)0;
304 if (x == RFVAL) return (Range)1;
305 if (x == APFLAG) return (Range)2;
306 if (x == NOECHO) return (Range)3;
307 if (x < -F_VR_OFFSET) return (Range)0;
308 return (Range)(F_FACTOR*(x+F_VR_OFFSET)+.5 + F_OFFSET); /* Default wsr88d coding. */
311 Range DR_INVF(float x) /* Differential reflectivity */
313 if (x == BADVAL) return (Range)0;
314 if (x == RFVAL) return (Range)1;
315 if (x == APFLAG) return (Range)2;
316 if (x == NOECHO) return (Range)3;
317 if (x < -F_DR_OFFSET) return (Range)0;
318 return (Range)(F_DR_FACTOR*(x + F_DR_OFFSET) + F_OFFSET + 0.5);
321 Range HC_INVF(float x) /* HydroClass (Sigmet) */
323 if (x == BADVAL) return (Range)0;
324 return (Range)(x + 0.5); /* Round */
327 Range LR_INVF(float x) /* MCTEX */
329 if (x == BADVAL) return (Range)0;
330 if (x == RFVAL) return (Range)1;
331 if (x == APFLAG) return (Range)2;
332 if (x == NOECHO) return (Range)3;
333 return (Range)((6.*x + 250) + 0.5); /* Round */
336 /**************************
337 Sigmet RhoHV : one_byte values
338 > RohHV = sqrt((N-1)/253)
345 ****************************/
346 /* RH_INVF for 1 or 2 byte data. */
347 Range RH_INVF(float x) {
348 if (x == BADVAL) return (Range)0;
349 /* return (Range)(x * x * 253.0 + 1.0 + 0.5); */
350 return (Range)(x * 65533. + 1. +.5);
353 /******************************
354 Sigmet PhiDP : one_byte values
355 > PhiDP (mod 180) = 180 * ((N-1)/254)
356 > The range is from 0 to 180 degrees in steps of 0.71 as follows.
362 *******************************/
363 Range PH_INVF(float x) {
364 if (x == BADVAL) return (Range)0;
365 /* return (Range)((x / 180.0) * 254.0 + 1.0 + 0.5); */
366 return (Range)(x*65534./360. + 1.0 + 0.5);
370 /* KD_INVF for 1 or 2 byte data. */
371 Range KD_INVF(float x) {
372 if (x == BADVAL) return (Range)0;
373 return (Range)(x * 100. + 32768. + 0.5);
374 /****** Old code for 1-byte Sigmet native data format commented-out:
375 if (x == RFVAL) return (Range)1;
376 if (x == APFLAG) return (Range)2;
377 if (x == NOECHO) return (Range)3;
378 if (rsl_kdp_wavelen == 0.0) return (Range)0;
381 126 * (log((double)-x) - log((double)(0.25/rsl_kdp_wavelen))) /
386 126 * (log((double)x) - log((double)0.25/rsl_kdp_wavelen)) /
396 /* Standard Deviation (for Dual-pole QC testing.) */
397 Range SD_INVF(float x)
399 if (x == BADVAL) return (Range)0;
400 return (Range)(x * 100.);
403 /* Signal Quality Index */
404 Range SQ_INVF(float x)
406 if (x == BADVAL) return (Range)0;
407 return (Range)(x * 65533. + 1. +.5);
410 /* Normalized Coherent Power (DORADE) */
411 Range NP_INVF(float x)
413 if (x == BADVAL) return (0);
414 return (Range)(x * 100. + 1.);
417 Range TI_INVF(float x) /* MCTEX */
419 if (x == BADVAL) return (Range)0;
420 if (x == RFVAL) return (Range)1;
421 if (x == APFLAG) return (Range)2;
422 if (x == NOECHO) return (Range)3;
427 Range SW_INVF(float x) { return VR_INVF(x); }
428 Range CZ_INVF(float x) { return DZ_INVF(x); }
429 Range ZT_INVF(float x) { return DZ_INVF(x); }
430 Range ZD_INVF(float x) { return DR_INVF(x); } /* Differential reflectivity */
431 Range CD_INVF(float x) { return DR_INVF(x); } /* Differential reflectivity */
432 Range XZ_INVF(float x) { return DZ_INVF(x); }
433 Range MZ_INVF(float x) { return (Range)x; } /* DZ Mask */
434 Range MD_INVF(float x) { return MZ_INVF(x); } /* ZD Mask */
435 Range ZE_INVF(float x) { return DZ_INVF(x); }
436 Range VE_INVF(float x) { return VR_INVF(x); }
437 Range DM_INVF(float x) { return DZ_INVF(x); }
438 Range DX_INVF(float x) { return DZ_INVF(x); }
439 Range CH_INVF(float x) { return DZ_INVF(x); }
440 Range AH_INVF(float x) { return DZ_INVF(x); }
441 Range CV_INVF(float x) { return DZ_INVF(x); }
442 Range AV_INVF(float x) { return DZ_INVF(x); }
443 Range VS_INVF(float x) { return VR_INVF(x); }
444 Range VL_INVF(float x) { return VR_INVF(x); }
445 Range VG_INVF(float x) { return VR_INVF(x); }
446 Range VT_INVF(float x) { return VR_INVF(x); }
447 Range VC_INVF(float x) { return VR_INVF(x); }
451 /**********************************************************************/
452 /* M E M O R Y M A N A G E M E N T R O U T I N E S */
453 /**********************************************************************/
454 /**********************************************************************/
460 /**********************************************************************/
461 Volume *RSL_new_volume(int max_sweeps)
464 * A volume consists of a header section and an array of sweeps.
467 v = (Volume *)calloc(1, sizeof(Volume));
468 if (v == NULL) perror("RSL_new_volume");
469 v->sweep = (Sweep **) calloc(max_sweeps, sizeof(Sweep*));
470 if (v->sweep == NULL) perror("RSL_new_volume, Sweep*");
471 v->h.nsweeps = max_sweeps; /* A default setting. */
476 * The 'Sweep_list' structure is internal to RSL. It maintains a list
477 * of sweeps allocated and it contains pointers to a hash table of Rays
478 * separately for each sweep. There is no reason to access this internal
479 * structure except when optimizing new RSL routines that access Rays.
480 * Otherwise, the RSL interfaces should suffice.
482 * The hash table is a means of finding rays, by azimuth, quickly.
483 * To find a ray is simple: use the hash function to get close
484 * to the ray, if not right on it the first time. Collisions of rays in
485 * the hash table are handled by a link list of rays from a hash entry.
486 * Typically, the first ray of the sweep is not the ray with the smallest
487 * azimuth angle. We are confident that the order of Rays in the Sweep
488 * is by azimuth angle, but that cannot be guarenteed. Therefore, this
489 * hash scheme is required.
491 * The 'Sweep_list' contains the address of the memory allocated to
492 * sweep. The list is sorted by addresses. There is no
493 * memory limit to the number of sweeps. If the number of sweeps exceeds
494 * the current allocation for the Sweep_list, then a new Sweep_list is
495 * allocated, which is bigger, and the old list copied to it.
497 * Sweep_list is at least as long as the number of sweeps allocated.
506 * By design of RSL, this should be "#define STATIC static"
508 * It is OK to "#define STATIC static", but, if you do, then
509 * the examples (run by run_tests in examples/) will fail for
510 * those programs that test these data structures. I normally,
511 * don't set this #define, for that reason.
515 STATIC int RSL_max_sweeps = 0; /* Initial allocation for sweep_list.
516 * RSL_new_sweep will allocate the space first
519 STATIC int RSL_nsweep_addr = 0; /* A count of sweeps in the table. */
520 STATIC Sweep_list *RSL_sweep_list = NULL;
521 STATIC int RSL_nextents = 0;
523 void FREE_HASH_NODE(Azimuth_hash *node)
525 if (node == NULL) return;
526 FREE_HASH_NODE(node->next); /* Tail recursive link list removal. */
530 void FREE_HASH_TABLE(Hash_table *table)
533 if (table == NULL) return;
534 for (i=0; i<table->nindexes; i++)
535 FREE_HASH_NODE(table->indexes[i]); /* A possible linked list of Rays. */
536 free(table->indexes);
540 void REMOVE_SWEEP(Sweep *s)
544 /* Find where it goes, split the list and slide the tail down one. */
545 for (i=0; i<RSL_nsweep_addr; i++)
546 if (s == RSL_sweep_list[i].s_addr) break;
548 if (i == RSL_nsweep_addr) return; /* Not found. */
549 /* This sweep is at 'i'. */
550 /* Deallocate the memory for the hash table. */
551 FREE_HASH_TABLE(RSL_sweep_list[i].hash);
554 for (j=i; j<RSL_nsweep_addr; j++)
555 RSL_sweep_list[j] = RSL_sweep_list[j+1];
557 RSL_sweep_list[RSL_nsweep_addr].s_addr = NULL;
558 RSL_sweep_list[RSL_nsweep_addr].hash = NULL;
562 int INSERT_SWEEP(Sweep *s)
564 Sweep_list *new_list;
567 if (RSL_nsweep_addr >= RSL_max_sweeps) { /* Current list is too small. */
569 new_list = (Sweep_list *) calloc(100*RSL_nextents, sizeof(Sweep_list));
570 if (new_list == NULL) {
571 perror("INSERT_SWEEP");
574 /* Copy the old list to the new one. */
575 for (i=0; i<RSL_max_sweeps; i++) new_list[i] = RSL_sweep_list[i];
576 RSL_max_sweeps = 100*RSL_nextents;
577 free(RSL_sweep_list);
578 RSL_sweep_list = new_list;
580 /* Find where it goes, split the list and slide the tail down one. */
581 for (i=0; i<RSL_nsweep_addr; i++)
582 if (s < RSL_sweep_list[i].s_addr) break;
584 /* This sweep goes at 'i'. But first we must split the list. */
585 for (j=RSL_nsweep_addr; j>i; j--)
586 RSL_sweep_list[j] = RSL_sweep_list[j-1];
588 RSL_sweep_list[i].s_addr = s;
589 RSL_sweep_list[i].hash = NULL;
594 int SWEEP_INDEX(Sweep *s)
596 /* Locate the sweep in the RSL_sweep_list. Return the index. */
597 /* Simple linear search; but this will be a binary search. */
599 for (i=0; i<RSL_nsweep_addr; i++)
600 if (s == RSL_sweep_list[i].s_addr) return i;
604 Sweep *RSL_new_sweep(int max_rays)
607 * A sweep consists of a header section and an array of rays.
610 s = (Sweep *)calloc(1, sizeof(Sweep));
611 if (s == NULL) perror("RSL_new_sweep");
613 s->ray = (Ray **) calloc(max_rays, sizeof(Ray*));
614 if (s->ray == NULL) perror("RSL_new_sweep, Ray*");
615 s->h.nrays = max_rays; /* A default setting. */
617 s->h.azimuth = -999.;
621 Ray *RSL_new_ray(int max_bins)
624 * A ray consists of a header section and an array of Range types (floats).
627 r = (Ray *)calloc(1, sizeof(Ray));
628 if (r == NULL) perror("RSL_new_ray");
629 r->range = (Range *) calloc(max_bins, sizeof(Range));
630 if (r->range == NULL) perror("RSL_new_ray, Range");
631 r->h.nbins = max_bins; /* A default setting. */
632 /* fprintf(stderr,"range[0] = %x, range[%d] = %x\n", &r->range[0], max_bins-1, &r->range[max_bins-1]);*/
636 /**********************************************************************/
642 /**********************************************************************/
643 Ray *RSL_clear_ray(Ray *r)
645 if (r == NULL) return r;
646 memset(r->range, 0, sizeof(Range)*r->h.nbins);
649 Sweep *RSL_clear_sweep(Sweep *s)
652 if (s == NULL) return s;
653 for (i=0; i<s->h.nrays; i++) {
654 RSL_clear_ray(s->ray[i]);
658 Volume *RSL_clear_volume(Volume *v)
661 if (v == NULL) return v;
662 for (i=0; i<v->h.nsweeps; i++) {
663 RSL_clear_sweep(v->sweep[i]);
667 /**********************************************************************/
673 /**********************************************************************/
674 void RSL_free_ray(Ray *r)
676 if (r == NULL) return;
677 if (r->range) free(r->range);
680 void RSL_free_sweep(Sweep *s)
683 if (s == NULL) return;
684 for (i=0; i<s->h.nrays; i++) {
685 RSL_free_ray(s->ray[i]);
687 if (s->ray) free(s->ray);
688 REMOVE_SWEEP(s); /* Remove from internal Sweep list. */
691 void RSL_free_volume(Volume *v)
694 if (v == NULL) return;
696 for (i=0; i<v->h.nsweeps; i++)
698 RSL_free_sweep(v->sweep[i]);
700 if (v->sweep) free(v->sweep);
704 /**********************************************************************/
710 /**********************************************************************/
711 Ray *RSL_copy_ray(Ray *r)
715 if (r == NULL) return NULL;
716 new_ray = RSL_new_ray(r->h.nbins);
718 memcpy(new_ray->range, r->range, r->h.nbins*sizeof(Range));
721 Sweep *RSL_copy_sweep(Sweep *s)
726 if (s == NULL) return NULL;
727 n_sweep = RSL_new_sweep(s->h.nrays);
728 if (n_sweep == NULL) return NULL;
731 for (i=0; i<s->h.nrays; i++) {
732 n_sweep->ray[i] = RSL_copy_ray(s->ray[i]);
739 Volume *RSL_copy_volume(Volume *v)
744 if (v == NULL) return NULL;
745 new_vol = RSL_new_volume(v->h.nsweeps);
748 for (i=0; i<v->h.nsweeps; i++) {
749 new_vol->sweep[i] = RSL_copy_sweep(v->sweep[i]);
755 /**********************************************************************/
756 /**********************************************************************/
757 /* G E N E R A L F U N C T I O N S */
758 /**********************************************************************/
759 /**********************************************************************/
761 double angle_diff(float x, float y)
764 d = fabs((double)(x - y));
765 if (d > 180) d = 360 - d;
769 /**********************************************************************/
771 /* RSL_get_next_cwise_ray */
772 /* Dennis Flanigan */
773 /* Mods by John Merritt 10/20/95 */
774 /**********************************************************************/
775 Ray *RSL_get_next_cwise_ray(Sweep *s, Ray *ray)
777 /* The fastest way to do this is to gain access to the hash table
778 * which maintains a linked list of sorted rays.
780 Hash_table *hash_table;
781 Azimuth_hash *closest;
785 if (s == NULL) return NULL;
786 if (ray == NULL) return NULL;
787 /* Find a non-NULL index close to hindex that we want. */
788 hash_table = hash_table_for_sweep(s);
789 if (hash_table == NULL) return NULL; /* Nada. */
790 ray_angle = ray->h.azimuth;
791 hindex = hash_bin(hash_table,ray_angle);
793 /* Find hash entry with closest Ray */
794 closest = the_closest_hash(hash_table->indexes[hindex],ray_angle);
796 return closest->ray_high->ray;
799 /**********************************************************************/
801 /* RSL_get_next_ccwise_ray */
803 /**********************************************************************/
804 Ray *RSL_get_next_ccwise_ray(Sweep *s, Ray *ray)
806 /* The fastest way to do this is to gain access to the hash table
807 * which maintains a linked list of sorted rays.
809 Hash_table *hash_table;
810 Azimuth_hash *closest;
814 if (s == NULL) return NULL;
815 if (ray == NULL) return NULL;
816 /* Find a non-NULL index close to hindex that we want. */
817 hash_table = hash_table_for_sweep(s);
818 if (hash_table == NULL) return NULL; /* Nada. */
819 ray_angle = ray->h.azimuth;
820 hindex = hash_bin(hash_table,ray_angle);
822 /* Find hash entry with closest Ray */
823 closest = the_closest_hash(hash_table->indexes[hindex],ray_angle);
825 return closest->ray_low->ray;
829 /******************************************
833 * Dennis Flanigan,Jr. 5/17/95 *
834 ******************************************/
835 double cwise_angle_diff(float x,float y)
837 /* Returns the clockwise angle difference of x to y.
838 * If x = 345 and y = 355 return 10.
839 * If x = 345 and y = 335 return 350
848 /******************************************
850 * ccwise_angle_diff *
852 * Dennis Flanigan,Jr. 5/17/95 *
853 ******************************************/
854 double ccwise_angle_diff(float x,float y)
856 /* Returns the counterclockwise angle differnce of x to y.
857 * If x = 345 and y = 355 return 350.
858 * If x = 345 and y = 335 return 10
867 /*****************************************
871 * Dennis Flanigan,Jr. 4/29/95 *
872 *****************************************/
873 Azimuth_hash *the_closest_hash(Azimuth_hash *hash, float ray_angle)
875 /* Return the hash pointer with the minimum ray angle difference. */
877 double clow,chigh,cclow;
878 Azimuth_hash *high,*low;
880 if (hash == NULL) return NULL;
882 /* Set low pointer to hash index with ray angle just below
883 * requested angle and high pointer to just above requested
887 /* set low and high pointers to initial search locations*/
889 high = hash->ray_high;
891 /* Search until clockwise angle to high is less then clockwise
895 clow = cwise_angle_diff(ray_angle,low->ray->h.azimuth);
896 chigh = cwise_angle_diff(ray_angle,high->ray->h.azimuth);
897 cclow = ccwise_angle_diff(ray_angle,low->ray->h.azimuth);
899 while((chigh > clow) && (clow != 0))
904 high = low->ray_high; /* Not the same low as line before ! */
909 high = low->ray_high; /* Not the same low as line before ! */
912 clow = cwise_angle_diff(ray_angle,low->ray->h.azimuth);
913 chigh = cwise_angle_diff(ray_angle,high->ray->h.azimuth);
914 cclow = ccwise_angle_diff(ray_angle,low->ray->h.azimuth);
928 /*******************************************************************/
930 /* get_closest_sweep_index */
932 /* Dennis Flanigan, Jr. 5/15/95 */
933 /*******************************************************************/
934 int get_closest_sweep_index(Volume *v,float sweep_angle)
938 float delta_angle = 91;
941 if(v == NULL) return -1;
945 for (i=0; i<v->h.nsweeps; i++)
948 if (s == NULL) continue;
949 check_angle = fabs((double)(s->h.elev - sweep_angle));
951 if(check_angle <= delta_angle)
953 delta_angle = check_angle;
965 /********************************************************************/
967 /* RSL_get_closest_sweep */
969 /* Dennis Flanigan, Jr. 5/15/95 */
970 /********************************************************************/
971 Sweep *RSL_get_closest_sweep(Volume *v,float sweep_angle,float limit)
973 /* Find closest sweep to requested angle. Assume PPI sweep for
974 * now. Meaning: sweep_angle represents elevation angle from
981 if (v == NULL) return NULL;
983 if((ci = get_closest_sweep_index(v,sweep_angle)) < 0)
990 delta_angle = fabs((double)(s->h.elev - sweep_angle));
992 if( delta_angle <= limit)
1002 /**********************************************************************/
1003 /* These are more specific routines to make coding hierarchical. */
1005 /* done 4/7/95 Ray *RSL_get_ray_from_sweep */
1006 /* done 3/31 float RSL_get_value_from_sweep */
1007 /* done 3/31 float RSL_get_value_from_ray */
1008 /* done 4/1 float RSL_get_value_at_h */
1010 /**********************************************************************/
1011 Ray *RSL_get_ray_from_sweep(Sweep *s, float ray_angle)
1013 /* Locate the Ray * for ray_angle in the sweep. */
1015 /* Sanity checks. */
1016 if (s == NULL) return NULL;
1017 if (ray_angle < 0) ray_angle += 360.0; /* Only positive angles. */
1018 if (ray_angle >= 360) ray_angle -= 360;
1020 return RSL_get_closest_ray_from_sweep(s,ray_angle,s->h.horz_half_bw);
1023 /**********************************************
1027 * Dennis Flanigan, Jr. 4/27/95 *
1028 **********************************************/
1029 int hash_bin(Hash_table *table,float angle)
1031 /* Internal Routine to calculate the hashing bin index
1037 res = 360.0/table->nindexes;
1038 hash = (int)(angle/res + res/2.0);/*Centered about bin.*/
1040 if(hash >= table->nindexes) hash = hash - table->nindexes;
1042 /* Could test see which direction is closer, but
1045 while(table->indexes[hash] == NULL) {
1047 if(hash >= table->nindexes) hash = 0;
1053 Hash_table *hash_table_for_sweep(Sweep *s)
1058 if (i==-1) { /* Obviously, an unregistered sweep. Most likely the
1059 * result of pointer assignments.
1061 i = INSERT_SWEEP(s);
1064 if (RSL_sweep_list[i].hash == NULL) { /* First time. Construct the table. */
1065 RSL_sweep_list[i].hash = construct_sweep_hash_table(s);
1068 return RSL_sweep_list[i].hash;
1071 /*********************************************************************/
1073 /* RSL_get_closest_ray_from_sweep */
1075 /* Dennis Flanigan 4/30/95 */
1076 /*********************************************************************/
1077 Ray *RSL_get_closest_ray_from_sweep(Sweep *s,float ray_angle, float limit)
1080 * Return closest Ray in Sweep within limit (angle) specified
1081 * in parameter list. Assume PPI mode.
1084 Hash_table *hash_table;
1085 Azimuth_hash *closest;
1088 if (s == NULL) return NULL;
1089 /* Find a non-NULL index close to hindex that we want. */
1090 hash_table = hash_table_for_sweep(s);
1091 if (hash_table == NULL) return NULL; /* Nada. */
1093 hindex = hash_bin(hash_table,ray_angle);
1095 /* Find hash entry with closest Ray */
1096 closest = the_closest_hash(hash_table->indexes[hindex],ray_angle);
1098 /* Is closest ray within limit parameter ? If
1099 * so return ray, else return NULL.
1102 close_diff = angle_diff(ray_angle,closest->ray->h.azimuth);
1104 if(close_diff <= limit) return closest->ray;
1110 /*********************************************************************/
1112 /* Rsl_get_value_from_sweep */
1114 /*********************************************************************/
1115 float RSL_get_value_from_sweep(Sweep *s, float azim, float r)
1117 /* Locate the polar point (r,azim) in the sweep. */
1119 if (s == NULL) return BADVAL;
1120 ray = RSL_get_ray_from_sweep(s, azim);
1121 if (ray == NULL) return BADVAL;
1122 return RSL_get_value_from_ray(ray, r);
1126 /*********************************************************************/
1128 /* RSL_get_range_of_range_index */
1129 /* D. Flanigan 8/18/95 */
1130 /*********************************************************************/
1131 float RSL_get_range_of_range_index(Ray *ray, int index)
1133 if (ray == NULL) return 0.0;
1134 if (index >= ray->h.nbins) return 0.0;
1135 return ray->h.range_bin1/1000.0 + index*ray->h.gate_size/1000.0;
1139 /************************************/
1140 /* RSL_get_value_from_ray */
1142 /* Updated 4/4/95 D. Flanigan */
1144 /************************************/
1145 float RSL_get_value_from_ray(Ray *ray, float r)
1152 if (ray == NULL) return BADVAL;
1154 if(ray->h.gate_size == 0)
1156 if(radar_verbose_flag)
1158 fprintf(stderr,"RSL_get_value_from_ray: ray->h.gate_size == 0\n");
1163 /* range_bin1 is range to center of first bin */
1164 bin_index = (int)(((rm - ray->h.range_bin1)/ray->h.gate_size) + 0.5);
1166 /* Bin indexes go from 0 to nbins - 1 */
1167 if (bin_index >= ray->h.nbins || bin_index < 0) return BADVAL;
1169 return ray->h.f(ray->range[bin_index]);
1173 /*********************************************************************/
1175 /* RSL_get_value_at_h */
1177 /*********************************************************************/
1178 float RSL_get_value_at_h(Volume *v, float azim, float grnd_r, float h)
1182 RSL_get_slantr_and_elev(grnd_r, h, &r, &elev);
1183 return RSL_get_value(v, elev, azim, r);
1187 /**********************************************************************/
1188 /* These take a Volume and return the appropriate structure. */
1190 /* done 4/21/95 Sweep *RSL_get_sweep */
1191 /* done 4/1 Ray *RSL_get_ray */
1192 /* done 4/1 float *RSL_get_value */
1193 /* done 5/3 Ray *RSL_get_ray_above */
1194 /* done 5/3 Ray *RSL_get_ray_below */
1195 /* done 5/12 Ray *RSL_get_ray_from_other_volume */
1197 /**********************************************************************/
1201 /*********************************************************************/
1205 /* Updated 5/15/95 Dennis Flanigan, Jr. */
1206 /*********************************************************************/
1207 Sweep *RSL_get_sweep(Volume *v, float sweep_angle)
1209 /* Return a sweep with +/- 1/2 beam_width of 'elev', if found. */
1212 if (v == NULL) return NULL;
1213 while(v->sweep[i] == NULL) i++;
1215 return RSL_get_closest_sweep(v,sweep_angle,v->sweep[i]->h.vert_half_bw);
1219 /*********************************************************************/
1223 /*********************************************************************/
1224 Ray *RSL_get_ray(Volume *v, float elev, float azimuth)
1226 /* Locate 'elev' and 'azimuth' in the Volume v by a simple +/- epsilon on
1227 * the elevation angle and azimuth angle.
1231 * 1. Locate sweep using azimuth; call RSL_get_sweep.
1232 * 2. Call RSL_get_ray_from_sweep
1235 return RSL_get_ray_from_sweep( RSL_get_sweep( v, elev ), azimuth );
1238 /*********************************************************************/
1242 /*********************************************************************/
1243 float RSL_get_value(Volume *v, float elev, float azimuth, float range)
1245 /* Locate 'elev' and 'azimuth' and '<range' in the Volume v
1246 * by a simple +/- epsilon on the elevation angle and azimuth angle
1251 * 1. Locate sweep using 'elev'.
1252 * 2. Call RSL_get_value_from_sweep
1254 return RSL_get_value_from_sweep ( RSL_get_sweep (v, elev), azimuth, range );
1257 /*********************************************************************/
1259 /* RSL_get_ray_above */
1261 /* Updated 5/15/95, Dennis Flanigan, Jr. */
1262 /*********************************************************************/
1263 Ray *RSL_get_ray_above(Volume *v, Ray *current_ray)
1267 if (v == NULL) return NULL;
1268 if (current_ray == NULL) return NULL;
1270 /* Find index of current Sweep */
1271 if(( i = get_closest_sweep_index(v,current_ray->h.elev)) < 0) return NULL;
1274 while( i < v->h.nsweeps)
1276 if(v->sweep[i] != NULL) break;
1280 if(i >= v->h.nsweeps) return NULL;
1282 return RSL_get_ray_from_sweep(v->sweep[i], current_ray->h.azimuth);
1286 /*********************************************************************/
1288 /* RSL_get_ray_below */
1290 /* Updated 5/15/95, Dennis Flanigan, Jr. */
1291 /*********************************************************************/
1292 Ray *RSL_get_ray_below(Volume *v, Ray *current_ray)
1296 if (v == NULL) return NULL;
1297 if (current_ray == NULL) return NULL;
1299 /* Find index of current Sweep */
1300 if(( i = get_closest_sweep_index(v,current_ray->h.elev)) < 0) return NULL;
1305 if(v->sweep[i] != NULL) break;
1309 if(i < 0) return NULL;
1311 return RSL_get_ray_from_sweep(v->sweep[i], current_ray->h.azimuth);
1314 /*********************************************************************/
1316 /* RSL_get_matching_ray */
1318 /*********************************************************************/
1319 Ray *RSL_get_matching_ray(Volume *v, Ray *ray)
1323 * Locate the closest matching ray in the Volume 'v' to 'ray'.
1324 * Typically, use this function when finding a similiar ray in another
1327 if (v == NULL) return NULL;
1328 if (ray == NULL) return NULL;
1330 return RSL_get_ray(v, ray->h.elev, ray->h.azimuth);
1333 /*********************************************************************/
1335 /* RSL_get_first_ray_of_sweep */
1336 /* RSL_get_first_ray_of_volume */
1338 /*********************************************************************/
1339 Ray *RSL_get_first_ray_of_sweep(Sweep *s)
1341 /* Because a sorting of azimuth angles may have been performed,
1342 * we need to test on the ray_num member and look for the smallest
1347 int smallest_ray_num;
1350 smallest_ray_num = 9999999;
1351 if (s == NULL) return r;
1352 for (i=0; i<s->h.nrays; i++)
1354 if (s->ray[i]->h.ray_num <= 1) return s->ray[i];
1355 if (s->ray[i]->h.ray_num < smallest_ray_num) {
1357 smallest_ray_num = r->h.ray_num;
1363 Ray *RSL_get_first_ray_of_volume(Volume *v)
1366 if (v == NULL) return NULL;
1367 for (i=0; i<v->h.nsweeps; i++)
1368 if (v->sweep[i]) return RSL_get_first_ray_of_sweep(v->sweep[i]);
1372 /*********************************************************************/
1374 /* RSL_get_first_sweep_of_volume */
1376 /*********************************************************************/
1377 Sweep *RSL_get_first_sweep_of_volume(Volume *v)
1380 if (v == NULL) return NULL;
1381 for (i=0; i<v->h.nsweeps; i++)
1382 if (RSL_get_first_ray_of_sweep(v->sweep[i])) return v->sweep[i];
1386 #define N_SPECIAL_NAMES 2
1388 * Unfortunately in C, there is no way around initializing static
1389 * arrays by specifying repetition.
1391 * There is another solution and that is to have RSL_new_radar set
1392 * a flag indicating if the application has called 'RSL_select_fields'
1393 * prior to calling the ingest routine. I choose the static = {...}; method
1397 /* Could be static and force use of 'rsl_query_field' */
1398 int rsl_qfield[MAX_RADAR_VOLUMES] = {
1411 /*********************************************************************/
1413 /* RSL_select_fields */
1415 /*********************************************************************/
1416 void RSL_select_fields(char *field_type, ...)
1420 * field_type = Case insensitive:
1421 * "all" - default, if never this routine is never called.
1422 * "none" - No fields are ingestd. Useful for getting header
1424 * "dz" - Ingest DZ volume.
1425 * "vr" - Ingest VR volume.
1426 * ... - Just list additional fields.
1428 * The last argument must be NULL. This signals this routine
1429 * when to stop parsing the field types.
1431 * Action or side-effect:
1432 * A second call to this fuction overrides any previous settings.
1433 * In other words, multiple calls are not additive. So, to get both
1434 * DZ and VR volumes, use:
1435 * RSL_select_fields("dz", "vr"); - Read both DZ and VR.
1437 * RSL_select_fields("dz"); - Read only DZ.
1438 * RSL_select_fields("vr"); - Read only VR, no DZ.
1440 * An RSL hidden array is set to flag which fields are selected.
1441 * This array is examined inside all ingest code. It is not available
1442 * to the application.
1449 for (i=0; i<MAX_RADAR_VOLUMES; i++) rsl_qfield[i] = 0;
1451 /* # arguments, should be <= MAX_RADAR_VOLUMES, but we can handle
1452 * typo's and redundancies. Each is processed in the order they
1456 c_field = field_type;
1457 va_start(ap, field_type);
1459 if (radar_verbose_flag) fprintf(stderr,"Selected fields for ingest:");
1461 /* CHECK EACH FIELD. This is a fancier case statement than C provides. */
1462 if (radar_verbose_flag) fprintf(stderr," %s", c_field);
1463 if (strcasecmp(c_field, "all") == 0) {
1464 for (i=0; i<MAX_RADAR_VOLUMES; i++) rsl_qfield[i] = 1;
1465 } else if (strcasecmp(c_field, "none") == 0) {
1466 for (i=0; i<MAX_RADAR_VOLUMES; i++) rsl_qfield[i] = 0;
1469 for (i=0; i<MAX_RADAR_VOLUMES; i++)
1470 if (strcasecmp(c_field, RSL_ftype[i]) == 0) {
1472 break; /* Break the for loop. */
1475 if (i == MAX_RADAR_VOLUMES) {
1476 if (radar_verbose_flag)
1477 fprintf(stderr, "\nRSL_select_fields: Invalid field name <<%s>> specified.\n", c_field);
1480 c_field = va_arg(ap, char *);
1483 if (radar_verbose_flag) fprintf(stderr,"\n");
1489 int rsl_query_field(char *c_field)
1493 * RSL interface, for library code developers, to rsl ingest code,
1494 * which is intended to be part of RSL ingest code, which access
1495 * the hidden array 'rsl_qfield' and reports if that field is to
1498 * Return 1 if YES, meaning yes ingest this field type.
1503 * All ingest code is meant to use this routine to decide whether
1504 * or not to allocate memory for a field type. For data formats
1505 * that are very large, this will help optimize the ingest on
1506 * small memory machines and hopefully avoid unnessary swapping.
1508 * LASSEN is a good example where there may be 10 or 12 input field
1509 * types, but the application only wants 2 or 3 of them.
1511 * The application interface is RSL_select_fields.
1515 /* Quiet the compilier when -pedantic. :-) */
1516 RSL_f_list[0] = RSL_f_list[0];
1517 RSL_invf_list[0] = RSL_invf_list[0];
1519 for (i=0; i<MAX_RADAR_VOLUMES; i++)
1520 if (strcasecmp(c_field, RSL_ftype[i]) == 0) {
1522 break; /* Break the for loop. */
1525 if (i == MAX_RADAR_VOLUMES) { /* We should never see this message for
1526 * properly written ingest code.
1528 fprintf(stderr, "rsl_query_field: Invalid field name <<%s>> specified.\n", c_field);
1531 /* 'i' is the index. Is it set? */
1532 return rsl_qfield[i];
1536 /* Could be static and force use of 'rsl_query_sweep' */
1537 int *rsl_qsweep = NULL; /* If NULL, then read all sweeps. Otherwise,
1538 * read what is on the list.
1540 #define RSL_MAX_QSWEEP 500 /* It'll be rediculious to have more. :-) */
1541 int rsl_qsweep_max = RSL_MAX_QSWEEP;
1543 /*********************************************************************/
1545 /* RSL_read_these_sweeps */
1547 /*********************************************************************/
1548 void RSL_read_these_sweeps(char *csweep, ...)
1554 /* "all", "none", "0", "1", "2", "3", ... */
1556 /* # arguments, should be <= 'max # sweeps expected', but, what is it?
1557 * We can handle typo's and redundancies. Each is processed in the
1558 * order they appear.
1562 va_start(ap, csweep);
1564 rsl_qsweep_max = -1;
1565 if (rsl_qsweep == NULL)
1566 rsl_qsweep = (int *)calloc(RSL_MAX_QSWEEP, sizeof(int));
1568 /* else Clear the array - a second call to this function over-rides
1569 * any previous settings. This holds even if the second call has
1573 for(i = 0;i< RSL_MAX_QSWEEP; i++)
1577 if (radar_verbose_flag) fprintf(stderr,"Selected sweeps for ingest:");
1578 for (;c_sweep; c_sweep = va_arg(ap, char *))
1580 /* CHECK EACH FIELD. This is a fancier case statement than C provides. */
1581 if (radar_verbose_flag) fprintf(stderr," %s", c_sweep);
1582 if (strcasecmp(c_sweep, "all") == 0) {
1583 for (i=0; i<RSL_MAX_QSWEEP; i++) rsl_qsweep[i] = 1;
1584 rsl_qsweep_max = RSL_MAX_QSWEEP;
1585 } else if (strcasecmp(c_sweep, "none") == 0) {
1586 /* Commented this out to save runtime -GJW
1587 * rsl_qsweep[] already initialized to 0 above.
1589 * for (i=0; i<RSL_MAX_QSWEEP; i++) rsl_qsweep[i] = 0;
1590 * rsl_qsweep_max = -1;
1593 i = sscanf(c_sweep,"%d", &isweep);
1594 if (i == 0) { /* No match, bad argument. */
1595 if (radar_verbose_flag) fprintf(stderr,"\nRSL_read_these_sweeps: bad parameter %s. Ignoring.\n", c_sweep);
1599 if (isweep < 0 || isweep > RSL_MAX_QSWEEP) {
1600 if (radar_verbose_flag) fprintf(stderr,"\nRSL_read_these_sweeps: parameter %s not in [0,%d). Ignoring.\n", c_sweep, RSL_MAX_QSWEEP);
1604 if (isweep > rsl_qsweep_max) rsl_qsweep_max = isweep;
1605 rsl_qsweep[isweep] = 1;
1609 if (radar_verbose_flag) fprintf(stderr,"\n");
1614 void RSL_fix_time (Ray *ray)
1618 /* Fixes possible overflow values in month, day, year, hh, mm, ss */
1619 /* Normally, ss should be the overflow. This code ensures end of
1620 * month, year and century are handled correctly by using the Unix
1623 if (ray == NULL) return;
1624 memset(&the_time, 0, sizeof(struct tm));
1625 the_time.tm_sec = ray->h.sec;
1626 fsec = ray->h.sec - the_time.tm_sec;
1627 the_time.tm_min = ray->h.minute;
1628 the_time.tm_hour = ray->h.hour;
1629 the_time.tm_mon = ray->h.month - 1;
1630 the_time.tm_year = ray->h.year - 1900;
1631 the_time.tm_mday = ray->h.day;
1632 the_time.tm_isdst = -1;
1633 (void) mktime(&the_time);
1634 /* The time is fixed. */
1635 ray->h.sec = the_time.tm_sec;
1637 ray->h.minute = the_time.tm_min;
1638 ray->h.hour = the_time.tm_hour;
1639 ray->h.month = the_time.tm_mon + 1;
1640 ray->h.year = the_time.tm_year + 1900;
1641 ray->h.day = the_time.tm_mday;
1645 /*********************************************************************/
1647 /* RSL_add_dbz_offset_to_ray */
1649 /*********************************************************************/
1651 Add the calibration factor 'dbz_offset' to each ray bin which
1652 contains a valid value.
1654 void RSL_add_dbz_offset_to_ray(Ray *r, float dbz_offset)
1659 if (r == NULL) return;
1660 for (ibin=0; ibin<r->h.nbins; ibin++)
1662 val = r->h.f(r->range[ibin]);
1663 if ( val >= (float)NOECHO ) continue; /* Invalid value */
1664 r->range[ibin] = r->h.invf(val + dbz_offset);
1668 /*********************************************************************/
1670 /* RSL_add_dbz_offset_to_sweep */
1672 /*********************************************************************/
1673 void RSL_add_dbz_offset_to_sweep(Sweep *s, float dbz_offset)
1676 if (s == NULL) return;
1677 for (iray=0; iray<s->h.nrays; iray++)
1678 RSL_add_dbz_offset_to_ray(s->ray[iray], dbz_offset);
1681 /*********************************************************************/
1683 /* RSL_add_dbz_offset_to_volume */
1685 /*********************************************************************/
1686 void RSL_add_dbz_offset_to_volume(Volume *v, float dbz_offset)
1689 if (v == NULL) return;
1690 for (isweep=0; isweep<v->h.nsweeps; isweep++)
1691 RSL_add_dbz_offset_to_sweep(v->sweep[isweep], dbz_offset);