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