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