Pileus Git - ~andy/rsl/blob - volume.c

   1 /*
   2     NASA/TRMM, Code 910.1.
   3     This is the TRMM Office Radar Software Library.
   4     Copyright (C) 1996, 1997
   5             John H. Merritt
   6             Space Applications Corporation
   7             Vienna, Virginia
   8
   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.
  13
  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.
  18
  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.
  22 */
  23 /*
  24  * Volume functions coded in this file:
  25  *
  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
  50  *
  51  * See image_gen.c for the Volume image generation functions.
  52  *
  53  * See doc/rsl_index.html for HTML formatted documentation.
  54  *
  55  * All routines herein coded, unless otherwise stated, by:
  56  *     John Merritt
  57  *     Space Applications Corporation
  58  *
  59  *  Contributions by
  60  *  Dennis Flanigan, Jr.
  61  *  flanigan@lance.colostate.edu
  62  */
  63 #include <stdlib.h>
  64 #include <stdio.h>
  65 #include <math.h>
  66 #include <string.h>
  67 #include <stdarg.h>
  68 int strcasecmp(const char *s1, const char *s2);
  69
  70 #define USE_RSL_VARS
  71 #include "rsl.h"
  72
  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)
  76
  77 extern int radar_verbose_flag;
  78
  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,
  87  * APFLAG, NOECHO.
  88  *
  89  * The conversion functions may NOT be macros.
  90  */
  91
  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
  96 #else
  97 #define F_FACTOR 2.0
  98 #define F_DR_FACTOR 10.0
  99 #define F_DZ_RANGE_OFFSET 32
 100 #endif
 101
 102 /* #define F_VR_OFFSET 63.5 */
 103 #define F_VR_OFFSET 127.0
 104 #define F_DR_OFFSET 12.0
 105
 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.
 109  *
 110  * You must change nsig_to_radar.c where F_OFFSET is used for optimization.
 111  */
 112 #define F_OFFSET 4
 113
 114
 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. */
 123 }
 124
 125 float VR_F(Range x) {
 126   float val;
 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); */
 130         return val;
 131   }
 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. */
 137 }
 138
 139 float DR_F(Range x) {    /* Differential reflectivity */
 140   float val;
 141   if (x >= F_OFFSET) { /* This test works when Range is unsigned. */
 142         val = (((float)x-F_OFFSET)/F_DR_FACTOR - F_DR_OFFSET);
 143         return val;
 144   }
 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. */
 150 }
 151
 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;
 159   return BADVAL;
 160  }
 161
 162 float HC_F(Range x) {  /* HydroClass (Sigmet) */
 163   if (x == 0) return BADVAL;
 164   return (float)x;
 165 }
 166
 167 /****************************
 168  Sigmet RhoHV : one_byte values
 169 > RohHV = sqrt((N-1)/253)
 170 > 0     bad/no value
 171 > 1     0.0000
 172 > 2     0.0629
 173 > 128   0.7085
 174 > 253   0.9980
 175 > 254   1.0000
 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.;
 181 }
 182
 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.
 187 > 0     bad/no value
 188 > 1     0.00 deg
 189 > 2     0.71 deg
 190 > 101   70.87 deg
 191 > 254   179.29 deg
 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.;
 197  }
 198
 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.
 202                               */
 203 /* KD_F for 1 or 2 byte. */
 204 float KD_F(Range x)
 205 {
 206 /****** Commented-out code for 1-byte Sigmet native data format.
 207   if (x >= F_OFFSET) {
 208     x -= F_OFFSET;
 209     if (rsl_kdp_wavelen == 0.0) return BADVAL;
 210     if (x < 128)
 211       return (float)(
 212                      -0.25 * pow((double)600.0,(double)((127-x)/126.0))
 213                      )/rsl_kdp_wavelen;
 214     else if (x > 128)
 215       return (float)(
 216                       0.25 * pow((double)600.0,(double)((x-129)/126.0))
 217                      )/rsl_kdp_wavelen;
 218     else
 219       return 0.0;
 220   }
 221   if (x == 1) return RFVAL;
 222   if (x == 2) return APFLAG;
 223   if (x == 3) return NOECHO;
 224   return BADVAL;
 225 ******/
 226
 227   if (x == 0) return BADVAL;
 228   return (x-32768.)/100.;
 229 }
 230
 231 float NP_F(Range x) {  /* Normalized Coherent Power (DORADE) */
 232   if (x == 0) return BADVAL;
 233   return (float)x / 100.;
 234 }
 235
 236 /* Signal Quality Index */
 237 float SQ_F(Range x) {
 238   if (x == 0) return BADVAL;
 239   return (float)(x-1) / 65533.;
 240 }
 241
 242 float TI_F(Range x)
 243 {
 244   if (x >= F_OFFSET) return (float)(x);
 245   if (x == 0) return BADVAL;
 246   if (x == 1) return RFVAL;
 247   if (x == 2) return APFLAG;
 248   if (x == 3) return NOECHO;
 249   return BADVAL;
 250 }
 251
 252 float SW_F(Range x) {  return VR_F(x); }
 253 float CZ_F(Range x) {  return DZ_F(x); }
 254 float ZT_F(Range x) {  return DZ_F(x); }
 255 float ZD_F(Range x) {  return DR_F(x); } /* Differential reflectivity */
 256 float CD_F(Range x) {  return DR_F(x); } /* Differential reflectivity */
 257 float XZ_F(Range x) {  return DZ_F(x); }
 258 float MZ_F(Range x) {  return (float)x; } /* DZ Mask */
 259 float MD_F(Range x) {  return MZ_F(x); }  /* ZD Mask */
 260 float ZE_F(Range x) {  return DZ_F(x); }
 261 float VE_F(Range x) {  return VR_F(x); }
 262 float DM_F(Range x) {  return DZ_F(x); }
 263 float DX_F(Range x) {  return DZ_F(x); }
 264 float CH_F(Range x) {  return DZ_F(x); }
 265 float AH_F(Range x) {  return DZ_F(x); }
 266 float CV_F(Range x) {  return DZ_F(x); }
 267 float AV_F(Range x) {  return DZ_F(x); }
 268 float VS_F(Range x) {  return VR_F(x); }
 269 float VL_F(Range x) {  return VR_F(x); }
 270 float VG_F(Range x) {  return VR_F(x); }
 271 float VT_F(Range x) {  return VR_F(x); }
 272 float VC_F(Range x) {  return VR_F(x); }
 273
 274
 275
 276
 277 /* Unfortunately, floats are stored differently than ints/shorts.  So,
 278  * we cannot simply set up a switch statement and we must test for
 279  * all the special cases first.  We must test for exactness.
 280  */
 281 Range DZ_INVF(float x)
 282 {
 283   if (x == BADVAL) return (Range)0;
 284   if (x == RFVAL)  return (Range)1;
 285   if (x == APFLAG)  return (Range)2;
 286   if (x == NOECHO)  return (Range)3;
 287   if (x < -F_DZ_RANGE_OFFSET) return (Range)0;
 288   return (Range)(F_FACTOR*(x+F_DZ_RANGE_OFFSET)+.5 + F_OFFSET); /* Default wsr88d. */
 289 }
 290
 291 Range VR_INVF(float x)
 292 {
 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_VR_OFFSET)    return (Range)0;
 298   return (Range)(F_FACTOR*(x+F_VR_OFFSET)+.5 + F_OFFSET); /* Default wsr88d coding. */
 299 }
 300
 301 Range DR_INVF(float x)     /* Differential reflectivity */
 302 {
 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_DR_OFFSET)    return (Range)0;
 308   return (Range)(F_DR_FACTOR*(x + F_DR_OFFSET) + F_OFFSET + 0.5);
 309 }
 310
 311 Range HC_INVF(float x)  /* HydroClass (Sigmet) */
 312 {
 313   if (x == BADVAL) return (Range)0;
 314   return (Range)(x + 0.5); /* Round */
 315 }
 316
 317 Range LR_INVF(float x) /* MCTEX */
 318 {
 319   if (x == BADVAL) return (Range)0;
 320   if (x == RFVAL)  return (Range)1;
 321   if (x == APFLAG)  return (Range)2;
 322   if (x == NOECHO)  return (Range)3;
 323   return (Range)((6.*x + 250) + 0.5); /* Round */
 324 }
 325
 326 /**************************
 327  Sigmet RhoHV : one_byte values
 328 > RohHV = sqrt((N-1)/253)
 329 > 0     bad/no value
 330 > 1     0.0000
 331 > 2     0.0629
 332 > 128   0.7085
 333 > 253   0.9980
 334 > 254   1.0000
 335 ****************************/
 336 /* RH_INVF for 1 or 2 byte data. */
 337 Range RH_INVF(float x) {
 338   if (x == BADVAL) return (Range)0;
 339   /* return (Range)(x * x * 253.0 + 1.0 + 0.5); */
 340   return (Range)(x * 65533. + 1. +.5);
 341 }
 342
 343 /******************************
 344  Sigmet PhiDP : one_byte values
 345 > PhiDP (mod 180) = 180 * ((N-1)/254)
 346 > The range is from 0 to 180 degrees in steps of 0.71 as follows.
 347 > 0     bad/no value
 348 > 1     0.00 deg
 349 > 2     0.71 deg
 350 > 101   70.87 deg
 351 > 254   179.29 deg
 352 *******************************/
 353 Range PH_INVF(float x) {
 354   if (x == BADVAL) return (Range)0;
 355   /* return (Range)((x / 180.0) * 254.0 + 1.0 + 0.5); */
 356   return (Range)(x*65534./360. + 1.0 + 0.5);
 357 }
 358
 359
 360 /* KD_INVF for 1 or 2 byte data. */
 361 Range KD_INVF(float x) {
 362   if (x == BADVAL) return (Range)0;
 363 /****** Commented-out code for 1-byte Sigmet native data format.
 364   if (x == RFVAL)  return (Range)1;
 365   if (x == APFLAG)  return (Range)2;
 366   if (x == NOECHO)  return (Range)3;
 367   if (rsl_kdp_wavelen == 0.0) return (Range)0;
 368   if (x < 0) {
 369         x = 127 -
 370           126 * (log((double)-x) - log((double)(0.25/rsl_kdp_wavelen))) /
 371           log((double)600.0) +
 372           0.5;
 373   } else if (x > 0) {
 374         x = 129 +
 375           126 * (log((double)x) - log((double)0.25/rsl_kdp_wavelen)) /
 376           log((double)600.0) +
 377           0.5;
 378   } else {
 379         x = 128;
 380   }
 381   x += F_OFFSET;
 382 ******/
 383   return (Range)(x * 100. + 32768. + 0.5);
 384
 385 }
 386
 387 /* Signal Quality Index */
 388 Range SQ_INVF(float x)
 389 {
 390   if (x == BADVAL) return (Range)0;
 391   return (Range)(x * 65533. + 1. +.5);
 392 }
 393
 394
 395 Range NP_INVF(float x) /* Normalized Coherent Power (DORADE) */
 396 {
 397   if (x == BADVAL) return (0);
 398   return (Range)(x * 100.);
 399 }
 400
 401 Range TI_INVF(float x) /* MCTEX */
 402 {
 403   if (x == BADVAL) return (Range)0;
 404   if (x == RFVAL)  return (Range)1;
 405   if (x == APFLAG)  return (Range)2;
 406   if (x == NOECHO)  return (Range)3;
 407   return (Range)(x);
 408 }
 409
 410
 411 Range SW_INVF(float x) {  return VR_INVF(x); }
 412 Range CZ_INVF(float x) {  return DZ_INVF(x); }
 413 Range ZT_INVF(float x) {  return DZ_INVF(x); }
 414 Range ZD_INVF(float x) {  return DR_INVF(x); } /* Differential reflectivity */
 415 Range CD_INVF(float x) {  return DR_INVF(x); } /* Differential reflectivity */
 416 Range XZ_INVF(float x) {  return DZ_INVF(x); }
 417 Range MZ_INVF(float x) {  return (Range)x;   } /* DZ Mask */
 418 Range MD_INVF(float x) {  return MZ_INVF(x); } /* ZD Mask */
 419 Range ZE_INVF(float x) {  return DZ_INVF(x); }
 420 Range VE_INVF(float x) {  return VR_INVF(x); }
 421 Range DM_INVF(float x) {  return DZ_INVF(x); }
 422 Range DX_INVF(float x) {  return DZ_INVF(x); }
 423 Range CH_INVF(float x) {  return DZ_INVF(x); }
 424 Range AH_INVF(float x) {  return DZ_INVF(x); }
 425 Range CV_INVF(float x) {  return DZ_INVF(x); }
 426 Range AV_INVF(float x) {  return DZ_INVF(x); }
 427 Range VS_INVF(float x) {  return VR_INVF(x); }
 428 Range VL_INVF(float x) {  return VR_INVF(x); }
 429 Range VG_INVF(float x) {  return VR_INVF(x); }
 430 Range VT_INVF(float x) {  return VR_INVF(x); }
 431 Range VC_INVF(float x) {  return VR_INVF(x); }
 432
 433
 434
 435 /**********************************************************************/
 436 /*        M E M O R Y    M A N A G E M E N T   R O U T I N E S        */
 437 /**********************************************************************/
 438 /**********************************************************************/
 439 /*                                                                    */
 440 /*                        new_volume                                  */
 441 /*                        new_sweep                                   */
 442 /*                        new_ray                                     */
 443 /*                                                                    */
 444 /**********************************************************************/
 445 Volume *RSL_new_volume(int max_sweeps)
 446 {
 447   /*
 448    * A volume consists of a header section and an array of sweeps.
 449    */
 450   Volume *v;
 451   v = (Volume *)calloc(1, sizeof(Volume));
 452   if (v == NULL) perror("RSL_new_volume");
 453   v->sweep = (Sweep **) calloc(max_sweeps, sizeof(Sweep*));
 454   if (v->sweep == NULL) perror("RSL_new_volume, Sweep*");
 455   v->h.nsweeps = max_sweeps; /* A default setting. */
 456   return v;
 457 }
 458
 459 /*
 460  * The 'Sweep_list' structure is internal to RSL.  It maintains a list
 461  * of sweeps allocated and it contains pointers to a hash table of Rays
 462  * separately for each sweep.  There is no reason to access this internal
 463  * structure except when optimizing new RSL routines that access Rays.
 464  * Otherwise, the RSL interfaces should suffice.
 465  *
 466  * The hash table is a means of finding rays, by azimuth, quickly.
 467  * To find a ray is simple:  use the hash function to get close
 468  * to the ray, if not right on it the first time.  Collisions of rays in
 469  * the hash table are handled by a link list of rays from a hash entry.
 470  * Typically, the first ray of the sweep is not the ray with the smallest
 471  * azimuth angle.  We are confident that the order of Rays in the Sweep
 472  * is by azimuth angle, but that cannot be guarenteed.  Therefore, this
 473  * hash scheme is required.
 474  *
 475  * The 'Sweep_list' contains the address of the memory allocated to
 476  * sweep.  The list is sorted by addresses.  There is no
 477  * memory limit to the number of sweeps.  If the number of sweeps exceeds
 478  * the current allocation for the Sweep_list, then a new Sweep_list is
 479  * allocated, which is bigger, and the old list copied to it.
 480  *
 481  * Sweep_list is at least as long as the number of sweeps allocated.
 482  */
 483
 484 typedef struct {
 485   Sweep *s_addr;
 486   Hash_table *hash;
 487 } Sweep_list;
 488
 489 /*
 490  * By design of RSL, this should be "#define STATIC static"
 491  *
 492  * It is OK to "#define STATIC static", but, if you do, then
 493  * the examples (run by run_tests in examples/) will fail for
 494  * those programs that test these data structures.  I normally,
 495  * don't set this #define, for that reason.
 496  */
 497
 498 #define STATIC
 499 STATIC int RSL_max_sweeps = 0; /* Initial allocation for sweep_list.
 500                                                         * RSL_new_sweep will allocate the space first
 501                                                         * time around.
 502                                                         */
 503 STATIC int RSL_nsweep_addr = 0; /* A count of sweeps in the table. */
 504 STATIC Sweep_list *RSL_sweep_list = NULL;
 505 STATIC int RSL_nextents = 0;
 506
 507 void FREE_HASH_NODE(Azimuth_hash *node)
 508 {
 509   if (node == NULL) return;
 510   FREE_HASH_NODE(node->next); /* Tail recursive link list removal. */
 511   free(node);
 512 }
 513
 514 void FREE_HASH_TABLE(Hash_table *table)
 515 {
 516   int i;
 517   if (table == NULL) return;
 518   for (i=0; i<table->nindexes; i++)
 519         FREE_HASH_NODE(table->indexes[i]); /* A possible linked list of Rays. */
 520   free(table->indexes);
 521   free(table);
 522 }
 523
 524 void REMOVE_SWEEP(Sweep *s)
 525 {
 526   int i;
 527   int j;
 528   /* Find where it goes, split the list and slide the tail down one. */
 529   for (i=0; i<RSL_nsweep_addr; i++)
 530         if (s == RSL_sweep_list[i].s_addr) break;
 531
 532   if (i == RSL_nsweep_addr) return; /* Not found. */
 533   /* This sweep is at 'i'. */
 534   /* Deallocate the memory for the hash table. */
 535   FREE_HASH_TABLE(RSL_sweep_list[i].hash);
 536
 537   RSL_nsweep_addr--;
 538   for (j=i; j<RSL_nsweep_addr; j++)
 539         RSL_sweep_list[j] = RSL_sweep_list[j+1];
 540
 541   RSL_sweep_list[RSL_nsweep_addr].s_addr = NULL;
 542   RSL_sweep_list[RSL_nsweep_addr].hash = NULL;
 543 }
 544
 545
 546 int INSERT_SWEEP(Sweep *s)
 547 {
 548   Sweep_list *new_list;
 549   int i,j;
 550
 551   if (RSL_nsweep_addr >= RSL_max_sweeps) { /* Current list is too small. */
 552         RSL_nextents++;
 553         new_list = (Sweep_list *) calloc(100*RSL_nextents, sizeof(Sweep_list));
 554         if (new_list == NULL) {
 555           perror("INSERT_SWEEP");
 556           exit(2);
 557         }
 558     /* Copy the old list to the new one. */
 559         for (i=0; i<RSL_max_sweeps; i++) new_list[i] = RSL_sweep_list[i];
 560         RSL_max_sweeps = 100*RSL_nextents;
 561         free(RSL_sweep_list);
 562         RSL_sweep_list = new_list;
 563   }
 564   /* Find where it goes, split the list and slide the tail down one. */
 565   for (i=0; i<RSL_nsweep_addr; i++)
 566         if (s < RSL_sweep_list[i].s_addr) break;
 567
 568   /* This sweep goes at 'i'. But first we must split the list. */
 569   for (j=RSL_nsweep_addr; j>i; j--)
 570         RSL_sweep_list[j] = RSL_sweep_list[j-1];
 571
 572   RSL_sweep_list[i].s_addr = s;
 573   RSL_sweep_list[i].hash = NULL;
 574   RSL_nsweep_addr++;
 575   return i;
 576 }
 577
 578 int SWEEP_INDEX(Sweep *s)
 579 {
 580   /* Locate the sweep in the RSL_sweep_list.  Return the index. */
 581   /* Simple linear search; but this will be a binary search. */
 582   int i;
 583   for (i=0; i<RSL_nsweep_addr; i++)
 584         if (s == RSL_sweep_list[i].s_addr) return i;
 585   return -1;
 586 }
 587
 588 Sweep *RSL_new_sweep(int max_rays)
 589 {
 590   /*
 591    * A sweep consists of a header section and an array of rays.
 592    */
 593   Sweep *s;
 594   s = (Sweep  *)calloc(1, sizeof(Sweep));
 595   if (s == NULL) perror("RSL_new_sweep");
 596   INSERT_SWEEP(s);
 597   s->ray = (Ray **) calloc(max_rays, sizeof(Ray*));
 598   if (s->ray == NULL) perror("RSL_new_sweep, Ray*");
 599   s->h.nrays = max_rays; /* A default setting. */
 600   return s;
 601 }
 602
 603 Ray *RSL_new_ray(int max_bins)
 604 {
 605   /*
 606    * A ray consists of a header section and an array of Range types (floats).
 607    */
 608   Ray *r;
 609   r = (Ray *)calloc(1, sizeof(Ray));
 610   if (r == NULL) perror("RSL_new_ray");
 611   r->range = (Range *) calloc(max_bins, sizeof(Range));
 612   if (r->range == NULL) perror("RSL_new_ray, Range");
 613   r->h.nbins = max_bins; /* A default setting. */
 614 /*  fprintf(stderr,"range[0] = %x, range[%d] = %x\n", &r->range[0], max_bins-1, &r->range[max_bins-1]);*/
 615   return r;
 616 }
 617
 618 /**********************************************************************/
 619 /*                                                                    */
 620 /*                      clear_ray                                     */
 621 /*                      clear_sweep                                   */
 622 /*                      clear_volume                                  */
 623 /*                                                                    */
 624 /**********************************************************************/
 625 Ray *RSL_clear_ray(Ray *r)
 626 {
 627   if (r == NULL) return r;
 628   memset(r->range, 0, sizeof(Range)*r->h.nbins);
 629   return r;
 630 }
 631 Sweep *RSL_clear_sweep(Sweep *s)
 632 {
 633   int i;
 634   if (s == NULL) return s;
 635   for (i=0; i<s->h.nrays; i++) {
 636         RSL_clear_ray(s->ray[i]);
 637   }
 638   return s;
 639 }
 640 Volume *RSL_clear_volume(Volume *v)
 641 {
 642   int i;
 643   if (v == NULL) return v;
 644   for (i=0; i<v->h.nsweeps; i++) {
 645         RSL_clear_sweep(v->sweep[i]);
 646   }
 647   return v;
 648 }
 649 /**********************************************************************/
 650 /*                                                                    */
 651 /*                       free_ray                                     */
 652 /*                       free_sweep                                   */
 653 /*                       free_volume                                  */
 654 /*                                                                    */
 655 /**********************************************************************/
 656 void RSL_free_ray(Ray *r)
 657 {
 658   if (r == NULL) return;
 659   if (r->range) free(r->range);
 660   free(r);
 661 }
 662 void RSL_free_sweep(Sweep *s)
 663 {
 664   int i;
 665   if (s == NULL) return;
 666   for (i=0; i<s->h.nrays; i++) {
 667         RSL_free_ray(s->ray[i]);
 668   }
 669   if (s->ray) free(s->ray);
 670   REMOVE_SWEEP(s); /* Remove from internal Sweep list. */
 671   free(s);
 672 }
 673 void RSL_free_volume(Volume *v)
 674 {
 675   int i;
 676   if (v == NULL) return;
 677
 678   for (i=0; i<v->h.nsweeps; i++)
 679          {
 680          RSL_free_sweep(v->sweep[i]);
 681          }
 682   if (v->sweep) free(v->sweep);
 683   free(v);
 684 }
 685
 686 /**********************************************************************/
 687 /*                                                                    */
 688 /*                       copy_ray                                     */
 689 /*                       copy_sweep                                   */
 690 /*                       copy_volume                                  */
 691 /*                                                                    */
 692 /**********************************************************************/
 693 Ray *RSL_copy_ray(Ray *r)
 694 {
 695   Ray *new_ray;
 696
 697   if (r == NULL) return NULL;
 698   new_ray = RSL_new_ray(r->h.nbins);
 699   new_ray->h = r->h;
 700   memcpy(new_ray->range, r->range, r->h.nbins*sizeof(Range));
 701   return new_ray;
 702 }
 703 Sweep *RSL_copy_sweep(Sweep *s)
 704 {
 705   int i;
 706   Sweep *n_sweep;
 707
 708   if (s == NULL) return NULL;
 709   n_sweep = RSL_new_sweep(s->h.nrays);
 710   if (n_sweep == NULL) return NULL;
 711   n_sweep->h = s->h;
 712
 713   for (i=0; i<s->h.nrays; i++) {
 714         n_sweep->ray[i] = RSL_copy_ray(s->ray[i]);
 715   }
 716   return n_sweep;
 717 }
 718
 719
 720
 721 Volume *RSL_copy_volume(Volume *v)
 722 {
 723   int i;
 724   Volume *new_vol;
 725
 726   if (v == NULL) return NULL;
 727   new_vol = RSL_new_volume(v->h.nsweeps);
 728   new_vol->h = v->h;
 729
 730   for (i=0; i<v->h.nsweeps; i++) {
 731         new_vol->sweep[i] = RSL_copy_sweep(v->sweep[i]);
 732   }
 733   return new_vol;
 734 }
 735
 736
 737 /**********************************************************************/
 738 /**********************************************************************/
 739 /*              G E N E R A L    F U N C T I O N S                    */
 740 /**********************************************************************/
 741 /**********************************************************************/
 742
 743 double angle_diff(float x, float y)
 744 {
 745   double d;
 746   d = fabs((double)(x - y));
 747   if (d > 180) d = 360 - d;
 748   return d;
 749 }
 750
 751 /**********************************************************************/
 752 /*                                                                    */
 753 /*                 RSL_get_next_cwise_ray                             */
 754 /* Dennis Flanigan                                                    */
 755 /* Mods by John Merritt 10/20/95                                      */
 756 /**********************************************************************/
 757 Ray *RSL_get_next_cwise_ray(Sweep *s, Ray *ray)
 758 {
 759   /* The fastest way to do this is to gain access to the hash table
 760    * which maintains a linked list of sorted rays.
 761    */
 762   Hash_table *hash_table;
 763   Azimuth_hash *closest;
 764   int hindex;
 765   float ray_angle;
 766
 767   if (s == NULL) return NULL;
 768   if (ray == NULL) return NULL;
 769   /* Find a non-NULL index close to hindex that we want. */
 770   hash_table = hash_table_for_sweep(s);
 771   if (hash_table == NULL) return NULL; /* Nada. */
 772   ray_angle = ray->h.azimuth;
 773   hindex = hash_bin(hash_table,ray_angle);
 774
 775   /* Find hash entry with closest Ray */
 776   closest = the_closest_hash(hash_table->indexes[hindex],ray_angle);
 777
 778   return closest->ray_high->ray;
 779 }
 780
 781 /**********************************************************************/
 782 /*                                                                    */
 783 /*                 RSL_get_next_ccwise_ray                            */
 784 /* JHM 10/20/95                                                       */
 785 /**********************************************************************/
 786 Ray *RSL_get_next_ccwise_ray(Sweep *s, Ray *ray)
 787 {
 788   /* The fastest way to do this is to gain access to the hash table
 789    * which maintains a linked list of sorted rays.
 790    */
 791   Hash_table *hash_table;
 792   Azimuth_hash *closest;
 793   int hindex;
 794   float ray_angle;
 795
 796   if (s == NULL) return NULL;
 797   if (ray == NULL) return NULL;
 798   /* Find a non-NULL index close to hindex that we want. */
 799   hash_table = hash_table_for_sweep(s);
 800   if (hash_table == NULL) return NULL; /* Nada. */
 801   ray_angle = ray->h.azimuth;
 802   hindex = hash_bin(hash_table,ray_angle);
 803
 804   /* Find hash entry with closest Ray */
 805   closest = the_closest_hash(hash_table->indexes[hindex],ray_angle);
 806
 807   return closest->ray_low->ray;
 808 }
 809
 810
 811 /******************************************
 812  *                                        *
 813  * cwise_angle_diff                       *
 814  *                                        *
 815  * Dennis Flanigan,Jr. 5/17/95            *
 816  ******************************************/
 817 double cwise_angle_diff(float x,float y)
 818    {
 819    /* Returns the clockwise angle difference of x to y.
 820         * If x = 345 and y = 355 return 10.
 821         * If x = 345 and y = 335 return 350
 822         */
 823    double d;
 824
 825    d = (double)(y - x);
 826    if (d < 0) d += 360;
 827    return d;
 828    }
 829
 830 /******************************************
 831  *                                        *
 832  * ccwise_angle_diff                      *
 833  *                                        *
 834  * Dennis Flanigan,Jr. 5/17/95            *
 835  ******************************************/
 836 double ccwise_angle_diff(float x,float y)
 837    {
 838    /* Returns the counterclockwise angle differnce of x to y.
 839         * If x = 345 and y = 355 return 350.
 840         * If x = 345 and y = 335 return 10
 841         */
 842    double d;
 843
 844    d = (double)(x - y);
 845    if (d < 0) d += 360;
 846    return d;
 847    }
 848
 849 /*****************************************
 850  *                                       *
 851  * the_closest_hash                      *
 852  *                                       *
 853  * Dennis Flanigan,Jr. 4/29/95           *
 854  *****************************************/
 855 Azimuth_hash *the_closest_hash(Azimuth_hash *hash, float ray_angle)
 856    {
 857    /* Return the hash pointer with the minimum ray angle difference. */
 858
 859    double clow,chigh,cclow;
 860    Azimuth_hash *high,*low;
 861
 862    if (hash == NULL) return NULL;
 863
 864    /* Set low pointer to hash index with ray angle just below
 865         * requested angle and high pointer to just above requested
 866         * angle.
 867         */
 868
 869    /* set low and high pointers to initial search locations*/
 870    low = hash;
 871    high = hash->ray_high;
 872
 873    /* Search until clockwise angle to high is less then clockwise
 874         * angle to low.
 875         */
 876
 877    clow   = cwise_angle_diff(ray_angle,low->ray->h.azimuth);
 878    chigh  = cwise_angle_diff(ray_angle,high->ray->h.azimuth);
 879    cclow  = ccwise_angle_diff(ray_angle,low->ray->h.azimuth);
 880
 881    while((chigh > clow) && (clow != 0))
 882           {
 883           if (clow < cclow)
 884                  {
 885                  low  = low->ray_low;
 886                  high = low->ray_high;  /* Not the same low as line before ! */
 887                  }
 888           else
 889                  {
 890                  low  = low->ray_high;
 891                  high = low->ray_high;  /* Not the same low as line before ! */
 892                  }
 893
 894           clow   = cwise_angle_diff(ray_angle,low->ray->h.azimuth);
 895           chigh  = cwise_angle_diff(ray_angle,high->ray->h.azimuth);
 896           cclow  = ccwise_angle_diff(ray_angle,low->ray->h.azimuth);
 897           }
 898
 899    if(chigh <= cclow)
 900           {
 901           return high;
 902           }
 903    else
 904           {
 905           return low;
 906           }
 907    }
 908
 909
 910 /*******************************************************************/
 911 /*                                                                 */
 912 /*       get_closest_sweep_index                                   */
 913 /*                                                                 */
 914 /* Dennis Flanigan, Jr.  5/15/95                                   */
 915 /*******************************************************************/
 916 int get_closest_sweep_index(Volume *v,float sweep_angle)
 917    {
 918    Sweep *s;
 919    int i,ci;
 920    float delta_angle = 91;
 921    float check_angle;
 922
 923    if(v == NULL) return -1;
 924
 925    ci = 0;
 926
 927    for (i=0; i<v->h.nsweeps; i++)
 928       {
 929       s = v->sweep[i];
 930       if (s == NULL) continue;
 931       check_angle = fabs((double)(s->h.elev - sweep_angle));
 932
 933       if(check_angle <= delta_angle)
 934          {
 935          delta_angle = check_angle;
 936          ci = i;
 937          }
 938       }
 939
 940    return ci;
 941    }
 942
 943
 944
 945
 946
 947 /********************************************************************/
 948 /*                                                                  */
 949 /*     RSL_get_closest_sweep                                        */
 950 /*                                                                  */
 951 /*  Dennis Flanigan, Jr. 5/15/95                                    */
 952 /********************************************************************/
 953 Sweep *RSL_get_closest_sweep(Volume *v,float sweep_angle,float limit)
 954    {
 955    /* Find closest sweep to requested angle.  Assume PPI sweep for
 956     * now. Meaning: sweep_angle represents elevation angle from
 957     * 0->90 degrees
 958     */
 959    Sweep *s;
 960    float delta_angle;
 961    int ci;
 962
 963    if (v == NULL) return NULL;
 964
 965    if((ci = get_closest_sweep_index(v,sweep_angle)) < 0)
 966       {
 967       return NULL;
 968       }
 969
 970    s = v->sweep[ci];
 971
 972    delta_angle = fabs((double)(s->h.elev - sweep_angle));
 973
 974    if( delta_angle <= limit)
 975       {
 976       return s;
 977       }
 978    else
 979       {
 980       return NULL;
 981       }
 982    }
 983
 984 /**********************************************************************/
 985 /*     These are more specific routines to make coding hierarchical.  */
 986 /*                                                                    */
 987 /* done 4/7/95    Ray   *RSL_get_ray_from_sweep                       */
 988 /* done 3/31      float  RSL_get_value_from_sweep                     */
 989 /* done 3/31      float  RSL_get_value_from_ray                       */
 990 /* done 4/1       float  RSL_get_value_at_h                           */
 991 /*                                                                    */
 992 /**********************************************************************/
 993 Ray *RSL_get_ray_from_sweep(Sweep *s, float ray_angle)
 994 {
 995   /* Locate the Ray * for ray_angle in the sweep. */
 996
 997   /* Sanity checks. */
 998   if (s == NULL) return NULL;
 999   if (ray_angle < 0) ray_angle += 360.0; /* Only positive angles. */
1000   if (ray_angle >= 360) ray_angle -= 360;
1001
1002   return RSL_get_closest_ray_from_sweep(s,ray_angle,s->h.horz_half_bw);
1003 }
1004
1005 /**********************************************
1006  *                                            *
1007  *           hash_bin                         *
1008  *                                            *
1009  * Dennis Flanigan, Jr. 4/27/95               *
1010  **********************************************/
1011 int hash_bin(Hash_table *table,float angle)
1012 {
1013   /* Internal Routine to calculate the hashing bin index
1014    * of a given angle.
1015    */
1016   int hash;
1017   float res;
1018
1019   res = 360.0/table->nindexes;
1020   hash = (int)(angle/res + res/2.0);/*Centered about bin.*/
1021
1022   if(hash >= table->nindexes) hash = hash - table->nindexes;
1023
1024   /* Could test see which direction is closer, but
1025    * why bother?
1026    */
1027   while(table->indexes[hash] == NULL) {
1028         hash++;
1029         if(hash >= table->nindexes) hash = 0;
1030   }
1031
1032   return hash;
1033 }
1034
1035 Hash_table *hash_table_for_sweep(Sweep *s)
1036 {
1037   int i;
1038
1039   i = SWEEP_INDEX(s);
1040   if (i==-1) { /* Obviously, an unregistered sweep.  Most likely the
1041                                 * result of pointer assignments.
1042                                 */
1043         i = INSERT_SWEEP(s);
1044   }
1045
1046   if (RSL_sweep_list[i].hash == NULL) { /* First time.  Construct the table. */
1047         RSL_sweep_list[i].hash = construct_sweep_hash_table(s);
1048   }
1049
1050   return RSL_sweep_list[i].hash;
1051 }
1052
1053 /*********************************************************************/
1054 /*                                                                   */
1055 /*                    RSL_get_closest_ray_from_sweep                 */
1056 /*                                                                   */
1057 /* Dennis Flanigan  4/30/95                                          */
1058 /*********************************************************************/
1059 Ray *RSL_get_closest_ray_from_sweep(Sweep *s,float ray_angle, float limit)
1060    {
1061    /*
1062     * Return closest Ray in Sweep within limit (angle) specified
1063     * in parameter list.  Assume PPI mode.
1064         */
1065    int hindex;
1066    Hash_table *hash_table;
1067    Azimuth_hash *closest;
1068    double close_diff;
1069
1070    if (s == NULL) return NULL;
1071    /* Find a non-NULL index close to hindex that we want. */
1072    hash_table = hash_table_for_sweep(s);
1073    if (hash_table == NULL) return NULL; /* Nada. */
1074
1075    hindex = hash_bin(hash_table,ray_angle);
1076
1077    /* Find hash entry with closest Ray */
1078    closest = the_closest_hash(hash_table->indexes[hindex],ray_angle);
1079
1080    /* Is closest ray within limit parameter ? If
1081         * so return ray, else return NULL.
1082         */
1083
1084    close_diff = angle_diff(ray_angle,closest->ray->h.azimuth);
1085
1086    if(close_diff <= limit) return closest->ray;
1087
1088    return NULL;
1089  }
1090
1091
1092 /*********************************************************************/
1093 /*                                                                   */
1094 /*                    Rsl_get_value_from_sweep                       */
1095 /*                                                                   */
1096 /*********************************************************************/
1097 float RSL_get_value_from_sweep(Sweep *s, float azim, float r)
1098 {
1099  /* Locate the polar point (r,azim) in the sweep. */
1100   Ray *ray;
1101   if (s == NULL) return BADVAL;
1102   ray = RSL_get_ray_from_sweep(s, azim);
1103   if (ray == NULL) return BADVAL;
1104   return RSL_get_value_from_ray(ray, r);
1105 }
1106
1107
1108 /*********************************************************************/
1109 /*                                                                   */
1110 /*                 RSL_get_range_of_range_index                      */
1111 /* D. Flanigan 8/18/95                                               */
1112 /*********************************************************************/
1113 float RSL_get_range_of_range_index(Ray *ray, int index)
1114    {
1115    if (ray == NULL) return 0.0;
1116    if (index >= ray->h.nbins) return 0.0;
1117    return ray->h.range_bin1/1000.0 + index*ray->h.gate_size/1000.0;
1118    }
1119
1120
1121 /************************************/
1122 /*  RSL_get_value_from_ray          */
1123 /*                                  */
1124 /*  Updated 4/4/95  D. Flanigan     */
1125 /*                                  */
1126 /************************************/
1127 float RSL_get_value_from_ray(Ray *ray, float r)
1128    {
1129    int bin_index;
1130    float rm;
1131
1132    rm = r * 1000;
1133
1134    if (ray == NULL) return BADVAL;
1135
1136    if(ray->h.gate_size == 0)
1137           {
1138           if(radar_verbose_flag)
1139                  {
1140                  fprintf(stderr,"RSL_get_value_from_ray: ray->h.gate_size == 0\n");
1141                  }
1142           return BADVAL;
1143           }
1144
1145    /* range_bin1 is range to center of first bin */
1146    bin_index = (int)(((rm - ray->h.range_bin1)/ray->h.gate_size) + 0.5);
1147
1148    /* Bin indexes go from 0 to nbins - 1 */
1149    if (bin_index >= ray->h.nbins || bin_index < 0) return BADVAL;
1150
1151    return ray->h.f(ray->range[bin_index]);
1152    }
1153
1154
1155 /*********************************************************************/
1156 /*                                                                   */
1157 /*                    RSL_get_value_at_h                             */
1158 /*                                                                   */
1159 /*********************************************************************/
1160 float RSL_get_value_at_h(Volume *v, float azim, float grnd_r, float h)
1161 {
1162   float elev, r;
1163
1164   RSL_get_slantr_and_elev(grnd_r, h, &r, &elev);
1165   return RSL_get_value(v, elev, azim, r);
1166 }
1167
1168
1169 /**********************************************************************/
1170 /*     These take a Volume and return the appropriate structure.      */
1171 /*                                                                    */
1172 /* done 4/21/95  Sweep  *RSL_get_sweep                                */
1173 /* done 4/1      Ray    *RSL_get_ray                                  */
1174 /* done 4/1      float  *RSL_get_value                                */
1175 /* done 5/3      Ray    *RSL_get_ray_above                            */
1176 /* done 5/3      Ray    *RSL_get_ray_below                            */
1177 /* done 5/12     Ray    *RSL_get_ray_from_other_volume                */
1178 /*                                                                    */
1179 /**********************************************************************/
1180
1181
1182
1183 /*********************************************************************/
1184 /*                                                                   */
1185 /*                    RSL_get_sweep                                  */
1186 /*                                                                   */
1187 /* Updated 5/15/95 Dennis Flanigan, Jr.                              */
1188 /*********************************************************************/
1189 Sweep *RSL_get_sweep(Volume *v, float sweep_angle)
1190    {
1191    /* Return a sweep with +/- 1/2 beam_width of 'elev', if found. */
1192    int   i = 0;
1193
1194    if (v == NULL) return NULL;
1195    while(v->sweep[i] == NULL) i++;
1196
1197    return RSL_get_closest_sweep(v,sweep_angle,v->sweep[i]->h.vert_half_bw);
1198    }
1199
1200
1201 /*********************************************************************/
1202 /*                                                                   */
1203 /*                    RSL_get_ray                                    */
1204 /*                                                                   */
1205 /*********************************************************************/
1206 Ray *RSL_get_ray(Volume *v, float elev, float azimuth)
1207 {
1208   /* Locate 'elev' and 'azimuth' in the Volume v by a simple +/- epsilon on
1209    * the elevation angle and azimuth angle.
1210    */
1211
1212   /*
1213    * 1. Locate sweep using azimuth; call RSL_get_sweep.
1214    * 2. Call RSL_get_ray_from_sweep
1215    */
1216
1217   return  RSL_get_ray_from_sweep( RSL_get_sweep( v, elev ), azimuth );
1218 }
1219
1220 /*********************************************************************/
1221 /*                                                                   */
1222 /*                    RSL_get_value                                  */
1223 /*                                                                   */
1224 /*********************************************************************/
1225 float RSL_get_value(Volume *v, float elev, float azimuth, float range)
1226 {
1227   /* Locate 'elev' and 'azimuth' and '<range' in the Volume v
1228    *  by a simple +/- epsilon on the elevation angle and azimuth angle
1229    *  and range.
1230    */
1231
1232   /*
1233    * 1. Locate sweep using 'elev'.
1234    * 2. Call RSL_get_value_from_sweep
1235    */
1236   return RSL_get_value_from_sweep ( RSL_get_sweep (v, elev), azimuth, range );
1237 }
1238
1239 /*********************************************************************/
1240 /*                                                                   */
1241 /*                    RSL_get_ray_above                              */
1242 /*                                                                   */
1243 /* Updated 5/15/95, Dennis Flanigan, Jr.                             */
1244 /*********************************************************************/
1245 Ray *RSL_get_ray_above(Volume *v, Ray *current_ray)
1246    {
1247    int i;
1248
1249    if (v == NULL) return NULL;
1250    if (current_ray == NULL) return NULL;
1251
1252    /* Find index of current Sweep */
1253    if(( i = get_closest_sweep_index(v,current_ray->h.elev)) < 0) return NULL;
1254
1255    i++;
1256    while( i < v->h.nsweeps)
1257           {
1258           if(v->sweep[i] != NULL) break;
1259           i++;
1260           }
1261
1262    if(i >= v->h.nsweeps) return NULL;
1263
1264    return RSL_get_ray_from_sweep(v->sweep[i], current_ray->h.azimuth);
1265    }
1266
1267
1268 /*********************************************************************/
1269 /*                                                                   */
1270 /*                    RSL_get_ray_below                              */
1271 /*                                                                   */
1272 /* Updated 5/15/95, Dennis Flanigan, Jr.                             */
1273 /*********************************************************************/
1274 Ray *RSL_get_ray_below(Volume *v, Ray *current_ray)
1275    {
1276    int i;
1277
1278    if (v == NULL) return NULL;
1279    if (current_ray == NULL) return NULL;
1280
1281    /* Find index of current Sweep */
1282    if(( i = get_closest_sweep_index(v,current_ray->h.elev)) < 0) return NULL;
1283
1284    i--;
1285    while( i >= 0)
1286           {
1287           if(v->sweep[i] != NULL) break;
1288           i--;
1289           }
1290
1291    if(i < 0) return NULL;
1292
1293    return RSL_get_ray_from_sweep(v->sweep[i], current_ray->h.azimuth);
1294    }
1295
1296 /*********************************************************************/
1297 /*                                                                   */
1298 /*                    RSL_get_matching_ray                           */
1299 /*                                                                   */
1300 /*********************************************************************/
1301 Ray *RSL_get_matching_ray(Volume *v, Ray *ray)
1302 {
1303
1304   /*
1305    * Locate the closest matching ray in the Volume 'v' to 'ray'.
1306    * Typically, use this function when finding a similiar ray in another
1307    * volume.
1308    */
1309   if (v == NULL) return NULL;
1310   if (ray == NULL) return NULL;
1311
1312   return RSL_get_ray(v, ray->h.elev, ray->h.azimuth);
1313 }
1314
1315 /*********************************************************************/
1316 /*                                                                   */
1317 /*                 RSL_get_first_ray_of_sweep                        */
1318 /*                 RSL_get_first_ray_of_volume                       */
1319 /*                                                                   */
1320 /*********************************************************************/
1321 Ray *RSL_get_first_ray_of_sweep(Sweep *s)
1322 {
1323 /* Because a sorting of azimuth angles may have been performed,
1324  * we need to test on the ray_num member and look for the smallest
1325  * one.
1326  */
1327   int i;
1328   Ray *r;
1329   int smallest_ray_num;
1330
1331   r = NULL;
1332   smallest_ray_num = 9999999;
1333   if (s == NULL) return r;
1334   for (i=0; i<s->h.nrays; i++)
1335         if (s->ray[i]) {
1336           if (s->ray[i]->h.ray_num <= 1) return s->ray[i];
1337           if (s->ray[i]->h.ray_num < smallest_ray_num) {
1338                 r = s->ray[i];
1339                 smallest_ray_num = r->h.ray_num;
1340           }
1341         }
1342   return r;
1343 }
1344
1345 Ray *RSL_get_first_ray_of_volume(Volume *v)
1346 {
1347   int i;
1348   if (v == NULL) return NULL;
1349   for (i=0; i<v->h.nsweeps; i++)
1350         if (v->sweep[i]) return RSL_get_first_ray_of_sweep(v->sweep[i]);
1351   return NULL;
1352 }
1353
1354 /*********************************************************************/
1355 /*                                                                   */
1356 /*                 RSL_get_first_sweep_of_volume                     */
1357 /*                                                                   */
1358 /*********************************************************************/
1359 Sweep *RSL_get_first_sweep_of_volume(Volume *v)
1360 {
1361   int i;
1362   if (v == NULL) return NULL;
1363   for (i=0; i<v->h.nsweeps; i++)
1364         if (RSL_get_first_ray_of_sweep(v->sweep[i])) return v->sweep[i];
1365   return NULL;
1366 }
1367
1368 #define N_SPECIAL_NAMES 2
1369 /*
1370  * Unfortunately in C, there is no way around initializing static
1371  * arrays by specifying repetition.
1372  *
1373  * There is another solution and that is to have RSL_new_radar set
1374  * a flag indicating if the application has called 'RSL_select_fields'
1375  * prior to calling the ingest routine.  I choose the static = {...}; method
1376  * for now.
1377  */
1378
1379 /* Could be static and force use of 'rsl_query_field' */
1380 int rsl_qfield[MAX_RADAR_VOLUMES] = {
1381   1, 1, 1, 1, 1,
1382   1, 1, 1, 1, 1,
1383   1, 1, 1, 1, 1,
1384   1, 1, 1, 1, 1,
1385   1, 1, 1, 1, 1,
1386   1, 1, 1, 1, 1,
1387   1, 1
1388  };
1389
1390
1391 /*********************************************************************/
1392 /*                                                                   */
1393 /*                 RSL_select_fields                                 */
1394 /*                                                                   */
1395 /*********************************************************************/
1396 void RSL_select_fields(char *field_type, ...)
1397 {
1398   /*
1399    * 10/15/96
1400    * field_type = Case insensitive:
1401    *           "all"   - default, if never this routine is never called.
1402    *           "none"  - No fields are ingestd. Useful for getting header
1403    *                     information.
1404    *           "dz"    - Ingest DZ volume.
1405    *           "vr"    - Ingest VR volume.
1406    *            ...    - Just list additional fields.
1407    *
1408    * The last argument must be NULL.  This signals this routine
1409    * when to stop parsing the field types.
1410    *
1411    * Action or side-effect:
1412    *   A second call to this fuction overrides any previous settings.
1413    *   In other words, multiple calls are not additive.  So, to get both
1414    *   DZ and VR volumes, use:
1415    *      RSL_select_fields("dz", "vr");  - Read both DZ and VR.
1416    *   and not:
1417    *      RSL_select_fields("dz");  -  Read only DZ.
1418    *      RSL_select_fields("vr");  -  Read only VR, no DZ.
1419    *
1420    * An RSL hidden array is set to flag which fields are selected.
1421    * This array is examined inside all ingest code.  It is not available
1422    * to the application.
1423    */
1424
1425   va_list ap;
1426   char *c_field;
1427   int i;
1428
1429   for (i=0; i<MAX_RADAR_VOLUMES; i++) rsl_qfield[i] = 0;
1430
1431   /* # arguments, should be <= MAX_RADAR_VOLUMES, but we can handle
1432    * typo's and redundancies.  Each is processed in the order they
1433    * appear.
1434    */
1435
1436   c_field = field_type;
1437   va_start(ap, field_type);
1438
1439   if (radar_verbose_flag) fprintf(stderr,"Selected fields for ingest:");
1440   while (c_field) {
1441         /* CHECK EACH FIELD. This is a fancier case statement than C provides. */
1442         if (radar_verbose_flag) fprintf(stderr," %s", c_field);
1443         if (strcasecmp(c_field, "all") == 0) {
1444           for (i=0; i<MAX_RADAR_VOLUMES; i++) rsl_qfield[i] = 1;
1445         } else if (strcasecmp(c_field, "none") == 0) {
1446           for (i=0; i<MAX_RADAR_VOLUMES; i++) rsl_qfield[i] = 0;
1447         } else {
1448
1449           for (i=0; i<MAX_RADAR_VOLUMES; i++)
1450                 if (strcasecmp(c_field, RSL_ftype[i]) == 0) {
1451                   rsl_qfield[i] = 1;
1452                   break; /* Break the for loop. */
1453                 }
1454
1455           if (i == MAX_RADAR_VOLUMES) {
1456                 if (radar_verbose_flag)
1457                   fprintf(stderr, "\nRSL_select_fields: Invalid field name <<%s>> specified.\n", c_field);
1458           }
1459         }
1460         c_field = va_arg(ap, char *);
1461   }
1462
1463   if (radar_verbose_flag) fprintf(stderr,"\n");
1464   va_end(ap);
1465
1466 }
1467
1468
1469 int rsl_query_field(char *c_field)
1470 {
1471   /*
1472    * 10/15/96
1473    * RSL interface, for library code developers, to rsl ingest code,
1474    * which is intended to be part of RSL ingest code, which access
1475    * the hidden array 'rsl_qfield' and reports if that field is to
1476    * be ingested.
1477    *
1478    * Return 1 if YES, meaning yes ingest this field type.
1479    * Return 0 if NO.
1480    */
1481
1482   /*
1483    * All ingest code is meant to use this routine to decide whether
1484    * or not to allocate memory for a field type.  For data formats
1485    * that are very large, this will help optimize the ingest on
1486    * small memory machines and hopefully avoid unnessary swapping.
1487    *
1488    * LASSEN is a good example where there may be 10 or 12 input field
1489    * types, but the application only wants 2 or 3 of them.
1490    *
1491    * The application interface is RSL_select_fields.
1492    */
1493   int i;
1494
1495   /* Quiet the compilier when -pedantic. :-) */
1496   RSL_f_list[0] = RSL_f_list[0];
1497   RSL_invf_list[0] = RSL_invf_list[0];
1498
1499   for (i=0; i<MAX_RADAR_VOLUMES; i++)
1500         if (strcasecmp(c_field, RSL_ftype[i]) == 0) {
1501           rsl_qfield[i] = 1;
1502           break; /* Break the for loop. */
1503         }
1504
1505   if (i == MAX_RADAR_VOLUMES) { /* We should never see this message for
1506                                                                  * properly written ingest code.
1507                                                                  */
1508         fprintf(stderr, "rsl_query_field: Invalid field name <<%s>> specified.\n", c_field);
1509   }
1510
1511   /* 'i' is the index. Is it set? */
1512   return rsl_qfield[i];
1513 }
1514
1515
1516 /* Could be static and force use of 'rsl_query_sweep' */
1517 int *rsl_qsweep = NULL;  /* If NULL, then read all sweeps. Otherwise,
1518                                                   * read what is on the list.
1519                                                   */
1520 #define RSL_MAX_QSWEEP 500 /* It'll be rediculious to have more. :-) */
1521 int rsl_qsweep_max = RSL_MAX_QSWEEP;
1522
1523 /*********************************************************************/
1524 /*                                                                   */
1525 /*                 RSL_read_these_sweeps                             */
1526 /*                                                                   */
1527 /*********************************************************************/
1528 void RSL_read_these_sweeps(char *csweep, ...)
1529 {
1530   va_list ap;
1531   char *c_sweep;
1532   int i, isweep;
1533
1534   /* "all", "none", "0", "1", "2", "3", ... */
1535
1536   /* # arguments, should be <= 'max # sweeps expected', but, what is it?
1537    * We can handle typo's and redundancies.  Each is processed in the
1538    * order they appear.
1539    */
1540
1541   c_sweep = csweep;
1542   va_start(ap, csweep);
1543
1544   rsl_qsweep_max = -1;
1545   if (rsl_qsweep == NULL)
1546     rsl_qsweep = (int *)calloc(RSL_MAX_QSWEEP, sizeof(int));
1547
1548   /* else Clear the array - a second call to this function over-rides
1549    * any previous settings.  This holds even if the second call has
1550    * bad arguments.
1551    */
1552   else
1553     for(i = 0;i< RSL_MAX_QSWEEP; i++)
1554       rsl_qsweep[i] = 0;
1555
1556
1557   if (radar_verbose_flag) fprintf(stderr,"Selected sweeps for ingest:");
1558   for (;c_sweep;        c_sweep = va_arg(ap, char *))
1559         {
1560         /* CHECK EACH FIELD. This is a fancier case statement than C provides. */
1561         if (radar_verbose_flag) fprintf(stderr," %s", c_sweep);
1562         if (strcasecmp(c_sweep, "all") == 0) {
1563           for (i=0; i<RSL_MAX_QSWEEP; i++) rsl_qsweep[i] = 1;
1564           rsl_qsweep_max = RSL_MAX_QSWEEP;
1565         } else if (strcasecmp(c_sweep, "none") == 0) {
1566          /* Commented this out to save runtime -GJW
1567           * rsl_qsweep[] already initialized to 0 above.
1568           *
1569           * for (i=0; i<RSL_MAX_QSWEEP; i++) rsl_qsweep[i] = 0;
1570           * rsl_qsweep_max = -1;
1571           */
1572         } else {
1573           i = sscanf(c_sweep,"%d", &isweep);
1574           if (i == 0) { /* No match, bad argument. */
1575                 if (radar_verbose_flag) fprintf(stderr,"\nRSL_read_these_sweeps: bad parameter %s.  Ignoring.\n", c_sweep);
1576                 continue;
1577           }
1578
1579           if (isweep < 0 || isweep > RSL_MAX_QSWEEP) {
1580                 if (radar_verbose_flag) fprintf(stderr,"\nRSL_read_these_sweeps: parameter %s not in [0,%d).  Ignoring.\n", c_sweep, RSL_MAX_QSWEEP);
1581                 continue;
1582           }
1583
1584           if (isweep > rsl_qsweep_max) rsl_qsweep_max = isweep;
1585           rsl_qsweep[isweep] = 1;
1586         }
1587   }
1588
1589   if (radar_verbose_flag) fprintf(stderr,"\n");
1590   va_end(ap);
1591 }
1592
1593 #include <time.h>
1594 void RSL_fix_time (Ray *ray)
1595 {
1596   struct tm the_time;
1597   float fsec;
1598   /* Fixes possible overflow values in month, day, year, hh, mm, ss */
1599   /* Normally, ss should be the overflow.  This code ensures end of
1600    * month, year and century are handled correctly by using the Unix
1601    * time functions.
1602    */
1603   if (ray == NULL) return;
1604   memset(&the_time, 0, sizeof(struct tm));
1605   the_time.tm_sec = ray->h.sec;
1606   fsec = ray->h.sec - the_time.tm_sec;
1607   the_time.tm_min  = ray->h.minute;
1608   the_time.tm_hour = ray->h.hour;
1609   the_time.tm_mon  = ray->h.month - 1;
1610   the_time.tm_year = ray->h.year - 1900;
1611   the_time.tm_mday = ray->h.day;
1612   the_time.tm_isdst = -1;
1613   (void) mktime(&the_time);
1614   /* The time is fixed. */
1615   ray->h.sec    = the_time.tm_sec;
1616   ray->h.sec   += fsec;
1617   ray->h.minute = the_time.tm_min;
1618   ray->h.hour   = the_time.tm_hour;
1619   ray->h.month  = the_time.tm_mon + 1;
1620   ray->h.year   = the_time.tm_year + 1900;
1621   ray->h.day    = the_time.tm_mday;
1622   return;
1623 }
1624
1625 /*********************************************************************/
1626 /*                                                                   */
1627 /*                 RSL_add_dbz_offset_to_ray                         */
1628 /*                                                                   */
1629 /*********************************************************************/
1630 /*
1631   Add the calibration factor 'dbz_offset' to each ray bin which
1632   contains a valid value.
1633 */
1634 void RSL_add_dbz_offset_to_ray(Ray *r, float dbz_offset)
1635 {
1636   int ibin;
1637   float val;
1638
1639   if (r == NULL) return;
1640   for (ibin=0; ibin<r->h.nbins; ibin++)
1641   {
1642         val = r->h.f(r->range[ibin]);
1643         if ( val >= (float)NOECHO ) continue;  /* Invalid value */
1644         r->range[ibin] = r->h.invf(val + dbz_offset);
1645   }
1646 }
1647
1648 /*********************************************************************/
1649 /*                                                                   */
1650 /*                 RSL_add_dbz_offset_to_sweep                       */
1651 /*                                                                   */
1652 /*********************************************************************/
1653 void RSL_add_dbz_offset_to_sweep(Sweep *s, float dbz_offset)
1654 {
1655   int iray;
1656   if (s == NULL) return;
1657   for (iray=0; iray<s->h.nrays; iray++)
1658         RSL_add_dbz_offset_to_ray(s->ray[iray], dbz_offset);
1659 }
1660
1661 /*********************************************************************/
1662 /*                                                                   */
1663 /*                 RSL_add_dbz_offset_to_volume                      */
1664 /*                                                                   */
1665 /*********************************************************************/
1666 void RSL_add_dbz_offset_to_volume(Volume *v, float dbz_offset)
1667 {
1668   int isweep;
1669   if (v == NULL) return;
1670   for (isweep=0; isweep<v->h.nsweeps; isweep++)
1671         RSL_add_dbz_offset_to_sweep(v->sweep[isweep], dbz_offset);
1672 }