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

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