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