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