Pileus Git - ~andy/rsl/blob - wsr88d.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  *------------------------------------------------------------
  25  * v1.14 5/12/95
  26  *------------------------------------------------------------
  27  *  Procedures:
  28  *   wsr88d_open
  29  *   wsr88d_close
  30  *   wsr88d_read_file_header
  31  *   wsr88d_read_tape_header
  32  *   wsr88d_read_sweep
  33  *   wsr88d_read_ray
  34  *   wsr88d_perror
  35  *   wsr88d_ray_to_float
  36  *
  37  *  Functions:
  38  *   wsr88d_get_nyquist
  39  *   wsr88d_get_atmos_atten_factor
  40  *   wsr88d_get_velocity_resolution
  41  *   wsr88d_get_volume_coverage
  42  *   wsr88d_get_elevation_angle
  43  *   wsr88d_get_azimuth
  44  *   wsr88d_get_range
  45  *   wsr88d_get_data
  46  *   wsr88d_get_time
  47  *   wsr88d_get_vcp_info(int vcp_num,int el_num)
  48  *   wsr88d_get_fix_angle(Wsr88d_ray *ray)
  49  *   wsr88d_get_pulse_count(Wsr88d_ray *ray)
  50  *   wsr88d_get_azimuth_rate(Wsr88d_ray *ray)
  51  *   wsr88d_get_pulse_width(Wsr88d_ray *ray)
  52  *   wsr88d_get_prt(Wsr88d_ray *ray)
  53  *   wsr88d_get_prf(Wsr88d_ray *ray)
  54  *   wsr88d_get_wavelength(Wsr88d_ray *ray)
  55  *   wsr88d_get_frequency(Wsr88d_ray *ray)
  56  *
  57  *  Misc. routines: (canidates for possible inclusion into the library)
  58  *   print_head
  59  *   print_packet_info
  60  *   free_and_clear_sweep
  61  *   clear_sweep
  62  *   print_sweep_info
  63  *   Cvt_date                    <- From Dan Austin
  64  *   Cvt_time                    <- From Dan Austin
  65  *
  66  */
  67 #include <stdio.h>
  68 #include <stdlib.h>
  69 #include <string.h>
  70 #include <unistd.h>
  71
  72 #include "wsr88d.h"
  73
  74 static int little_endian(void)
  75 {
  76   union {
  77     unsigned char byte[4];
  78     int val;
  79   } word;
  80   word.val = 0;
  81   word.byte[3] = 0x1;
  82   return word.val != 1;
  83 }
  84
  85
  86 static void swap_4_bytes(void *word)
  87 {
  88   unsigned char *byte;
  89   unsigned char temp;
  90   byte    = word;
  91   temp    = byte[0];
  92   byte[0] = byte[3];
  93   byte[3] = temp;
  94   temp    = byte[1];
  95   byte[1] = byte[2];
  96   byte[2] = temp;
  97 }
  98
  99 static void swap_2_bytes(void *word)
 100 {
 101   unsigned char *byte;
 102   unsigned char temp;
 103   byte    = word;
 104   temp    = byte[0];
 105   byte[0] = byte[1];
 106   byte[1] = temp;
 107 }
 108
 109 /**********************************************************************/
 110 /*   D E B U G G I N G     R O U T I N E S    F O L L O W             */
 111 /**********************************************************************/
 112 /************************************************
 113  * Cvt_date-  convert the date in days since 1/1/70 (Julian) to mm/dd/yy
 114  * parameters:
 115  * long int date_in - input julian date
 116  * returns: output date
 117  * calls from: Cvt_pckt_hdr
 118  * calls to: none
 119  ************************************************/
 120
 121 #include <time.h>
 122 int Cvt_date(long int date_in)
 123 {
 124   int mm, dd, yy;
 125   time_t itime;
 126   struct tm *tm_time;
 127   itime = date_in - 1;
 128   itime *= 24*60*60; /* Seconds/day * days. */
 129
 130   tm_time = gmtime(&itime);
 131   mm = tm_time->tm_mon+1;
 132   dd = tm_time->tm_mday;
 133   yy = tm_time->tm_year;
 134
 135   return 10000.0*yy+100.0*mm+dd;
 136 }
 137
 138 /************************************************
 139  * Cvt_time- converts 24 hr time in msecs after midnight to hhmmss
 140  * parameters:
 141  * long int time_in - input time in msecs after midnight
 142  * returns: double *time_out - output time
 143  * calls from: Cvt_pckt_hdr
 144  * calls to: none
 145  ************************************************/
 146 float Cvt_time(long int time_in)
 147 {
 148     double t;
 149     int hh,mm;
 150
 151     t = time_in;
 152     t /= 1000.0;
 153     hh = t/3600;
 154     t -= hh*3600;
 155     mm = t/60;
 156     t -= mm*60;
 157
 158     return hh*10000 + mm*100 + (float) t;
 159 }
 160 /**********************************************************************/
 161 /*                                                                    */
 162 /*  done 2/28        print_head                                       */
 163 /*                                                                    */
 164 /**********************************************************************/
 165 void print_head(Wsr88d_file_header h)
 166 {
 167   int i;
 168   fprintf(stderr,"Filename : ");
 169   for (i=0;i<9;i++) fprintf(stderr,"%c", h.title.filename[i]);   printf("\n");
 170
 171   fprintf(stderr,"Extension: ");
 172   for (i=0;i<3;i++) fprintf(stderr,"%c", h.title.ext[i]);   printf("\n");
 173
 174   fprintf(stderr,"Julian date: %d\n", Cvt_date(h.title.file_date));
 175   fprintf(stderr,"       time: %f\n", Cvt_time(h.title.file_time));
 176
 177
 178 }
 179
 180 void print_packet_info(Wsr88d_packet *p)
 181 {
 182   fprintf(stderr,"%5hd %5hd %5hd %5hd %5hd %5hd %5hd %10.3f %6d\n",
 183          p->msg_type, p->id_seq, p->azm, p->ray_num, p->ray_status, p->elev, p->elev_num,
 184          Cvt_time((int)p->ray_time), Cvt_date((int)p->ray_date));
 185 }
 186
 187
 188
 189
 190 /**********************************************************************/
 191 /* End of debug routines.                                             */
 192 /**********************************************************************/
 193
 194 void free_and_clear_sweep(Wsr88d_sweep *s, int low, int high)
 195 {
 196 /* Frees and sets the ray pointers to NULL.
 197  * Assumes that rays pointers have been allocated.
 198  */
 199   int i;
 200   for (i=low; i<high; i++)
 201     if (s->ray[i] != NULL) {
 202       free(s->ray[i]);
 203       s->ray[i] = NULL;
 204     }
 205 }
 206
 207 void clear_sweep(Wsr88d_sweep *s, int low, int high)
 208 {
 209 /*
 210  * Simply set all sweep pointers to NULL.
 211  */
 212   int i;
 213   for (i=low; i<high; i++)
 214     s->ray[i] = NULL;
 215 }
 216
 217 void wsr88d_print_sweep_info(Wsr88d_sweep *s)
 218 {
 219   int i;
 220
 221   fprintf(stderr,"Mtype    ID  azim  ray# rstat  elev elev#       time   date\n");
 222   fprintf(stderr,"----- ----- ----- ----- ----- ----- ----- ---------- ------\n");
 223
 224   for (i=0; i<MAX_RAYS_IN_SWEEP; i++) {
 225     if (s->ray[i] != NULL)
 226       print_packet_info((Wsr88d_packet *) s->ray[i]);
 227   }
 228 }
 229
 230 /**********************************************************************/
 231 /*                                                                    */
 232 /*  done 2/28             wsr88d_open                                 */
 233 /*                                                                    */
 234 /**********************************************************************/
 235
 236 Wsr88d_file *wsr88d_open(char *filename)
 237 {
 238   Wsr88d_file *wf = (Wsr88d_file *)malloc(sizeof(Wsr88d_file));
 239   int save_fd;
 240
 241   if ( strcmp(filename, "stdin") == 0 ) {
 242     save_fd = dup(0);
 243     wf->fptr = fdopen(save_fd,"r");
 244   } else {
 245     wf->fptr = fopen(filename, "r");
 246   }
 247
 248   if (wf->fptr == NULL) return NULL;
 249   wf->fptr = uncompress_pipe(wf->fptr);
 250 #define NEW_BUFSIZ 16384
 251   setvbuf(wf->fptr,NULL,_IOFBF,(size_t)NEW_BUFSIZ); /* Faster i/o? */
 252   return wf;
 253 }
 254
 255
 256 /**********************************************************************/
 257 /*                                                                    */
 258 /*  done 2/28             wsr88d_perror                               */
 259 /*                                                                    */
 260 /**********************************************************************/
 261 int wsr88d_perror(char *message)
 262 {
 263 /*
 264  * I want to use a global 'wsr88d_errno' and
 265  * have this routine print an appropriate message.
 266  */
 267
 268   /* This is a simple model now. */
 269   fprintf(stderr, "wsr88d_error: ");
 270   perror(message);
 271   return 0;
 272 }
 273
 274 /**********************************************************************/
 275 /*                                                                    */
 276 /*  done 2/28             wsr88d_close                                */
 277 /*                                                                    */
 278 /**********************************************************************/
 279 int wsr88d_close(Wsr88d_file *wf)
 280 {
 281   int rc;
 282   rc = rsl_pclose(wf->fptr);
 283   free(wf);
 284   return rc;
 285 }
 286
 287
 288 /**********************************************************************/
 289 /*                                                                    */
 290 /*                     wsr88d_swap_file_header                        */
 291 /*                                                                    */
 292 /**********************************************************************/
 293 void wsr88d_swap_file_header(Wsr88d_file_header *header)
 294 {
 295   swap_4_bytes(&header->title.file_date);
 296   swap_4_bytes(&header->title.file_time);
 297 }
 298
 299 /**********************************************************************/
 300 /*                                                                    */
 301 /*  done 2/28          wsr88d_read_file_header                        */
 302 /*                                                                    */
 303 /**********************************************************************/
 304 int wsr88d_read_file_header(Wsr88d_file *wf,
 305                             Wsr88d_file_header *wsr88d_file_header)
 306 {
 307   int n;
 308   n = fread(&wsr88d_file_header->title,
 309             sizeof(wsr88d_file_header->title), 1, wf->fptr);
 310   if (little_endian())
 311     wsr88d_swap_file_header(wsr88d_file_header);
 312   return n;
 313 }
 314
 315 /**********************************************************************/
 316 /*                                                                    */
 317 /*  done 8/18          wsr88d_read_tape_header                        */
 318 /*                                                                    */
 319 /**********************************************************************/
 320 int wsr88d_read_tape_header(char *first_file,
 321                             Wsr88d_tape_header *wsr88d_tape_header)
 322 {
 323   FILE *fp;
 324   int n;
 325   char c;
 326
 327   if ((fp = fopen(first_file, "r")) == NULL) {
 328     perror(first_file);
 329     return 0;
 330   }
 331
 332   n = fread(wsr88d_tape_header, sizeof(Wsr88d_tape_header), 1, fp);
 333   if (n == 0) {
 334     fprintf(stderr, "WARNING: %s is smaller than 31616 bytes.  It is not a tape header file.\n", first_file);
 335   } else {
 336     /* Try to read one more character.  If we can, then this is not a
 337      * tape header file.  I suppose that we could look for '.' as the
 338      * 9-th character and if it were there, then too this is not a tape
 339      * header file.
 340      */
 341     if (fread(&c, sizeof(char), 1, fp) > 0) {
 342       fprintf(stderr, "WARNING: %s is larger than 31616 bytes.  It is not a tape header file.\n", first_file);
 343       memset(wsr88d_tape_header, 0, sizeof(Wsr88d_tape_header));
 344       n = 0;
 345     } else { /* Ok so far. Now check the first 8 bytes for "ARCHIVE2" */
 346       if (strncmp(wsr88d_tape_header->archive2, "ARCHIVE2", 8) != 0) {
 347         fprintf(stderr, "WARNING: %s is 31616 bytes.  However, the first 8 bytes are not 'ARCHIVE2'.\n", first_file);
 348         memset(wsr88d_tape_header, 0, sizeof(Wsr88d_tape_header));
 349         n = 0;
 350       }
 351     }
 352
 353   }
 354   fclose(fp);
 355   return n;
 356 }
 357
 358
 359
 360 /**********************************************************************/
 361 /*                                                                    */
 362 /*  not done N/A       wsr88d_read_header                             */
 363 /*                                                                    */
 364 /**********************************************************************/
 365 int wsr88d_read_header(Wsr88d_file *wf, Wsr88d_header *wsr88d_header)
 366 {
 367   fprintf(stderr,"Routine: wsr88d_read_header\n");
 368   return 0;
 369 }
 370
 371
 372 /**********************************************************************/
 373 /*                                                                    */
 374 /*  done 3/2           wsr88d_read_sweep                              */
 375 /*                                                                    */
 376 /**********************************************************************/
 377 int wsr88d_read_sweep(Wsr88d_file *wf, Wsr88d_sweep *wsr88d_sweep)
 378 {
 379   int n;
 380   Wsr88d_ray wsr88d_ray;
 381   int nrec;
 382   int end_of_volume;
 383   int ray_num;
 384
 385 /* One sweep is defined as staying at the same RDA elevation number. */
 386 /* We can read the file and check for that, however, we will need to
 387  * buffer our input.  The solution is to read the file and check the
 388  * radial status.  If it is '2' then we have reached the END OF ELEVATION.
 389  * Here is a complete list of radial status codes:
 390  *    0 = Start of new elevation.
 391  *    1 = Intermediate radial.
 392  *    2 = End of elevation.
 393  *    3 = Beginning of volume scan.
 394  *    4 = End of volume scan.
 395  */
 396
 397 /* Algorithm steps:
 398  *  1. Skip packets until.  Start of new elevation or
 399  *     Beginning of Volume scan.  STAT=0 or 3.
 400  *  2. Read until End of elevation.  STAT=2 or 4.  Skip message type != 1.
 401  */
 402
 403   nrec = 0;
 404   ray_num = 0;
 405   n = wsr88d_read_ray(wf, &wsr88d_ray);
 406
 407 /* Step 1. */
 408   while ((wsr88d_ray.msg_type & 15) != 1 && n > 0) {
 409     /*
 410     fprintf(stderr,"SKIPPING packet: type %d, radial status %d\n",
 411            wsr88d_ray.msg_type, wsr88d_ray.ray_status);
 412     */
 413     n = wsr88d_read_ray(wf, &wsr88d_ray);
 414   }
 415
 416   if (n <= 0) return n; /* Read failure. */
 417   end_of_volume = 0;
 418 /* Step 2. */
 419   while ( ! end_of_volume ) {
 420     if ((wsr88d_ray.msg_type & 15) != 1) {
 421       /*
 422       fprintf(stderr,"SKIPPING (amid a sweep) packet: type %d, "
 423             "radial status %d\n",
 424              wsr88d_ray.msg_type, wsr88d_ray.ray_status);
 425        */
 426
 427     } else {
 428       /* Load this ray into the sweep. */
 429       ray_num = wsr88d_ray.ray_num;
 430       /* Double check against #  records we've seen. */
 431       /* It is possible that a reset occurs and we begin to overwrite
 432        * previously loaded rays.  I've seen this occur, rarely, in the
 433        * WSR88D data.  I must trust 'ray_num'.
 434        */
 435       /*
 436       if (nrec+1 != ray_num) {
 437         fprintf(stderr, "Data says %d is ray_num, but, I've seen %d "
 438                 "records.\n", ray_num, nrec+1);
 439       }
 440       */
 441       if (wsr88d_sweep->ray[ray_num] == NULL) {
 442         wsr88d_sweep->ray[ray_num] = (Wsr88d_ray *) malloc (sizeof(Wsr88d_ray));
 443       }
 444       memcpy(wsr88d_sweep->ray[ray_num], &wsr88d_ray, sizeof(Wsr88d_ray));
 445     }
 446     n = wsr88d_read_ray(wf, &wsr88d_ray);
 447     if (n > 0) nrec++;
 448     end_of_volume = wsr88d_ray.ray_status == 2 ||
 449                     wsr88d_ray.ray_status == 4 ||
 450                     n <= 0;
 451   }
 452
 453   /* Process the last packet of the input data. */
 454   if ((wsr88d_ray.ray_status == 2 || wsr88d_ray.ray_status == 4) &&
 455           (wsr88d_ray.msg_type & 15) == 1) {
 456     /* Load this ray into the sweep. */
 457     ray_num = wsr88d_ray.ray_num;
 458     if (wsr88d_sweep->ray[ray_num] == NULL) {
 459       wsr88d_sweep->ray[ray_num] = (Wsr88d_ray *) malloc (sizeof(Wsr88d_ray));
 460     }
 461     memcpy(wsr88d_sweep->ray[ray_num], &wsr88d_ray, sizeof(Wsr88d_ray));
 462   }
 463
 464   /* Just to be safe, clear all ray pointers left in this sweep to
 465    * the maximum MAX_RAYS_IN_SWEEP.  This is required when the
 466    * wsr88d_sweep is reused and not cleared.
 467    */
 468   free_and_clear_sweep(wsr88d_sweep, ray_num+1, MAX_RAYS_IN_SWEEP);
 469
 470 /*
 471   fprintf(stderr,"Processed %d records for elevation number %d\n",
 472          nrec+1, wsr88d_ray.elev_num);
 473   wsr88d_print_sweep_info(wsr88d_sweep);
 474 */
 475   return nrec;
 476 }
 477
 478 /**********************************************************************/
 479 /*                                                                    */
 480 /*                      wsr88d_swap_ray                               */
 481 /*                                                                    */
 482 /**********************************************************************/
 483 void wsr88d_swap_ray(Wsr88d_ray *wsr88d_ray)
 484 {
 485   short *half_word;
 486   half_word = (short *)wsr88d_ray;
 487   for (; half_word<(short *)&wsr88d_ray->msg_time; half_word++)
 488     swap_2_bytes(half_word);
 489
 490   swap_4_bytes(&wsr88d_ray->msg_time);
 491   swap_2_bytes(&wsr88d_ray->num_seg);
 492   swap_2_bytes(&wsr88d_ray->seg_num);
 493   swap_4_bytes(&wsr88d_ray->ray_time);
 494
 495   half_word = (short *) &wsr88d_ray->ray_date;
 496   for (; half_word<(short *)&wsr88d_ray->sys_cal; half_word++)
 497     swap_2_bytes(half_word);
 498
 499   swap_4_bytes(&wsr88d_ray->sys_cal);
 500
 501   half_word = (short *) &wsr88d_ray->refl_ptr;
 502   for (; half_word<(short *)&wsr88d_ray->data[0]; half_word++)
 503     swap_2_bytes(half_word);
 504
 505 }
 506
 507 /**********************************************************************/
 508 /*                                                                    */
 509 /*  done 2/28           wsr88d_read_ray                               */
 510 /*                                                                    */
 511 /**********************************************************************/
 512 int wsr88d_read_ray(Wsr88d_file *wf, Wsr88d_ray *wsr88d_ray)
 513 {
 514   int n;
 515   n = fread(wsr88d_ray, sizeof(Wsr88d_ray), 1, wf->fptr);
 516 /*  if (n > 0) print_packet_info(wsr88d_ray); */
 517   if (little_endian())
 518     wsr88d_swap_ray(wsr88d_ray);
 519
 520   return n;
 521 }
 522
 523 /**********************************************************************/
 524 /*                                                                    */
 525 /*  not done N/A     wsr88d_read_ray_header                           */
 526 /*                                                                    */
 527 /**********************************************************************/
 528 int wsr88d_read_ray_header(Wsr88d_file *wf,
 529                            Wsr88d_ray_header *wsr88d_ray_header)
 530 {
 531   fprintf(stderr,"Stub routine: wsr88d_read_ray_header.\n");
 532   return 0;
 533 }
 534
 535 /**********************************************************************/
 536 /*                                                                    */
 537 /*  done 3/3         wsr88d_ray_to_float                              */
 538 /*                                                                    */
 539 /**********************************************************************/
 540 int wsr88d_ray_to_float(Wsr88d_ray *ray,
 541                         int THE_DATA_WANTED, float v[], int *n)
 542 {
 543 /*
 544  *  Input: *ray             -  WSR-88D packet
 545  * Output: THE_DATA_WANTED  -  Indicates which field to convert.  Fields:
 546  *                             WSR88D_DZ, WSR88D_VR, WSR88D_SW
 547  *         v[]              -  The output vector of float values.
 548  *         n                -  Length of the output vector.  0 = no data.
 549  *
 550  * Returns n.
 551  *
 552  * No allocation of space for the output vector performed here.
 553  */
 554
 555 /* Code from Dan Austin (cvt_pckt_data.c) was the template for this. */
 556
 557   /* declarations   */
 558   int num_ref_gates,num_vel_gates,num_spec_gates;
 559   int refl_ptr, vel_ptr,spec_ptr,res_flag;
 560   int ival;
 561   int i;
 562
 563   *n = 0;
 564   num_ref_gates  = ray->num_refl;
 565   num_vel_gates  = ray->num_dop;
 566   num_spec_gates = ray->num_dop;  /* 'num_dop', this is not a typo. */
 567
 568 /* The data pointers are specified from the begining of the
 569  * 'Digital Radar Data (Message) Header'.  Since we have a structure
 570  * that defines all the header variables and a member called 'data'.
 571  * we must subtract the length of the 'message header' from the data
 572  * pointer.  Hopefully, the reflecivity pointer will be 0 meaning the
 573  * first element of the 'data' member; ray->data[0];
 574  */
 575 #define LENGTH_OF_MESSAGE 100
 576   if (num_ref_gates > 0) refl_ptr = ray->refl_ptr - LENGTH_OF_MESSAGE;
 577   else refl_ptr = 0;
 578
 579   vel_ptr  = ray->vel_ptr - LENGTH_OF_MESSAGE;
 580   spec_ptr = ray->spc_ptr - LENGTH_OF_MESSAGE;
 581
 582   res_flag = ray->vel_res;
 583
 584
 585 /*
 586   fprintf(stderr,"refl_ptr = %d  #g = %d, ", refl_ptr, num_ref_gates);
 587   fprintf(stderr," vel_ptr = %d  #g = %d, ", vel_ptr, num_vel_gates);
 588   fprintf(stderr,"spec_ptr = %d  #g = %d, ", spec_ptr, num_spec_gates);
 589   fprintf(stderr,"res_flag = %d\n", res_flag);
 590 */
 591
 592   if (THE_DATA_WANTED == WSR88D_DZ) {
 593     /* do the reflectivity data  (dbZ)*/
 594     if (refl_ptr+num_ref_gates > 2300)
 595       fprintf(stderr, "WARNING: # refl index (%d) exceeds maximum (2300)\n",
 596               refl_ptr+num_ref_gates);
 597     else {
 598     for(i=0; i<num_ref_gates; i++) {
 599       ival = ray->data[refl_ptr+i];
 600       if(ival > 1)
 601           v[i] = (((ival-2.0)/2.0)-32.0);
 602       else if (ival == 1)
 603         v[i] = WSR88D_RFVAL;
 604       else /* ival = 0 */
 605         v[i] = WSR88D_BADVAL;
 606     }
 607     *n = num_ref_gates;
 608     }
 609
 610   } else if (THE_DATA_WANTED == WSR88D_VR) {
 611     /* do the velocity data  (M/S) */
 612     if (vel_ptr+num_vel_gates > 2300)
 613       fprintf(stderr, "WARNING: # vel index (%d) exceeds maximum (2300)\n",
 614               vel_ptr+num_vel_gates);
 615     else {
 616     for(i=0; i<num_vel_gates;i++)   {
 617       ival = ray->data[vel_ptr+i];
 618       if(ival > 1)
 619         if (res_flag == 2) /* High resolution: 0.5 m/s */
 620           v[i] = (((ival-2.0)/2.0)-63.5);
 621         else
 622           v[i] = ((ival-2.0)-127.0);
 623       else if (ival == 1)
 624         v[i] = WSR88D_RFVAL;
 625       else /* ival = 0 */
 626         v[i] = WSR88D_BADVAL;
 627     }
 628     *n = num_vel_gates;
 629     }
 630
 631   } else if (THE_DATA_WANTED == WSR88D_SW) {
 632     /* now do the spectrum width data (M/S)*/
 633     if (spec_ptr+num_spec_gates > 2300)
 634       fprintf(stderr, "WARNING: # spec index (%d) exceeds maximum (2300)\n",
 635               spec_ptr+num_spec_gates);
 636     else {
 637     for(i=0;i<num_spec_gates;i++) {
 638       ival = ray->data[spec_ptr+i];
 639         if(ival > 1)
 640           v[i] = (((ival-2)/2.0)-63.5);
 641         else if (ival == 1)
 642           v[i] = WSR88D_RFVAL;
 643         else /* ival = 0 */
 644           v[i] = WSR88D_BADVAL;
 645     }
 646     *n = num_spec_gates;
 647     }
 648   }
 649
 650   return *n;
 651 }
 652
 653
 654
 655 /**********************************************************************/
 656 /*        Functions that convert some message header values.          */
 657 /**********************************************************************/
 658 /**********************************************************************/
 659 /*                                                                    */
 660 /*  done 3/3   float wsr88d_get_nyquist                               */
 661 /*  done 3/3   float wsr88d_get_atmos_atten_factor                    */
 662 /*  done 3/3   float wsr88d_get_velocity_resolution                   */
 663 /*  done 3/3   int   wsr88d_get_volume_coverage                       */
 664 /*  done 3/3   float wsr88d_get_elevation_angle                       */
 665 /*  done 3/3   float wsr88d_get_azimuth                               */
 666 /*  done 3/3   float wsr88d_get_range                                 */
 667 /*  done 3/3   void  wsr88d_get_date                                  */
 668 /*  done 3/3   void  wsr88d_get_time                                  */
 669 /*  done 5/20  int  *wsr88d_get_vcp_info                              */
 670 /*  done 5/20  float wsr88d_get_fix_angle                             */
 671 /*  done 5/20  int   wsr88d_get_pulse_count                           */
 672 /*  done 5/20  float wsr88d_get_azimuth_rate                          */
 673 /*  done 5/20  float wsr88d_get_pulse_width                           */
 674 /*  done 5/20  float wsr88d_get_prf                                   */
 675 /*  done 5/20  float wsr88d_get_prt                                   */
 676 /*  done 5/20  float wsr88d_get_wavelength                            */
 677 /*  done 5/20  float wsr88d_get_frequency                             */
 678 /*                                                                    */
 679 /**********************************************************************/
 680 float wsr88d_get_nyquist(Wsr88d_ray *ray)
 681 {
 682   return ray->nyq_vel/100.0;
 683 }
 684
 685 float wsr88d_get_atmos_atten_factor(Wsr88d_ray *ray)
 686 {
 687   return ray->atm_att/1000.0;
 688 }
 689
 690 float wsr88d_get_velocity_resolution(Wsr88d_ray *ray)
 691 {
 692   if (ray->vel_res == 2) return 0.5;
 693   return 0.0;
 694 }
 695
 696 int   wsr88d_get_volume_coverage(Wsr88d_ray *ray)
 697 {
 698   if (ray->vol_cpat == 11) return 11;
 699   if (ray->vol_cpat == 12) return 12;
 700   if (ray->vol_cpat == 21) return 21;
 701   if (ray->vol_cpat == 31) return 31;
 702   if (ray->vol_cpat == 32) return 32;
 703   if (ray->vol_cpat == 121) return 121;
 704   if (ray->vol_cpat == 211) return 211;
 705   if (ray->vol_cpat == 212) return 212;
 706   if (ray->vol_cpat == 221) return 221;
 707   return 0;
 708 }
 709
 710 float wsr88d_get_elevation_angle(Wsr88d_ray *ray)
 711 {
 712   return ray->elev/8.0*(180.0/4096.0);
 713 }
 714
 715 float wsr88d_get_azimuth(Wsr88d_ray *ray)
 716 {
 717   return ray->azm/8.0*(180.0/4096.0);
 718 }
 719
 720 float wsr88d_get_range(Wsr88d_ray *ray)
 721 {
 722   return ray->unam_rng/10.0;
 723 }
 724
 725 #include <time.h>
 726 void wsr88d_get_date(Wsr88d_ray *ray, int *mm, int *dd, int *yy)
 727 {
 728 /*
 729  * mm (1-12)
 730  * dd (1-31)
 731  * yy (ex. 93)
 732  */
 733   time_t itime;
 734   struct tm *tm_time;
 735   if (ray == NULL) {
 736     *mm = *dd = *yy = 0;
 737     return;
 738   }
 739
 740   itime = ray->ray_date - 1;
 741   itime *= 24*60*60; /* Seconds/day * days. */
 742
 743   tm_time = gmtime(&itime);
 744   *mm = tm_time->tm_mon+1;
 745   *dd = tm_time->tm_mday;
 746   *yy = tm_time->tm_year;
 747 }
 748
 749 void wsr88d_get_time(Wsr88d_ray *ray, int *hh, int *mm, int *ss, float *fsec)
 750 {
 751   /*
 752    * hh (0-23)
 753    * mm (0-59)
 754    * ss (0-59)
 755    * fsec (fraction of second)
 756    */
 757   double t;
 758
 759   if (ray == NULL) {
 760     *hh = *mm = *ss = *fsec = 0;
 761     return;
 762   }
 763   t = ray->ray_time;
 764   t /= 1000.0;
 765   *hh = t/3600;
 766   t -= *hh*3600;
 767   *mm = t/60;
 768   t -= *mm*60;
 769   *ss = (int)t;
 770   *fsec = t - *ss;
 771 }
 772
 773
 774
 775 /*
 776  * Get_vcp_info - gets info about the volume coverage pattern for this scan
 777  * parameters:
 778  * int vcp_num - volume coverage pattern number
 779  * int el_num - elevation number w/in vcp
 780  * returns: int *vcp_info - ptr to array w/vcp info
 781  * calls from: Nexrad2uf
 782  * calls to: none
 783  */
 784
 785 /* this database contains volume coverage patterns & associated info:   */
 786 /* (0)= vcp # (1)=pulse width for vcp  "Id$"                            */
 787 /* line[1-n]: (n,0)= elev. # for vcp (n,1)= (fixed angle)*8*(4096/180)  */
 788 /* (n,2)= pulse count (n,3)= (azimuthal sweep rate)*8*(4096/45)         */
 789
 790 static int vcp11[68] ={11,514,88,17,13600,88,0,14000,264,16,12664,264,0,14000,440,6,11736,608,6,24760,784,6,24760,952,10,12712,1128,10,12720,1368,0,18328,1584,0,18496,1824,0,18512,2184,0,18544,2552,0,18576,3040,0,18640,3552,0,18712};
 791
 792 static int vcp12[53]={12,514,91,15,15401,91,0,18204,164,15,15401,164,0,18204,237,15,15401,237,0,18204,328,3,19297,437,3,20393,564,3,20393,728,3,20393,928,3,20393,1165,0,20680,1456,0,20680,1820,0,21033,2276,0,20929,2840,0,20929,3550,0,20929};
 793
 794 static int vcp21[48]={21,514,88,28,8256,88,0,8272,264,28,8256,264,0,8272,440,8,7888,608,8,7888,784,8,8160,1096,12,8160,1800,0,10640,2656,0,10432,3552,0,10496};
 795
 796 static int vcp31[36]={31,516,88,63,3672,88,0,3688,272,63,3672,272,0,3688,456,63,3672,456,0,3688,640,0,3688,816,0,3688};
 797
 798 static int vcp32[32]={32,514,88,64,3616,88,0,3312,272,64,3616,272,0,3312,456,11,2960,640,11,2960,816,11,2960};
 799
 800 static int vcp121[62]={121,514,91,11,21336,91,0,21696,91,0,19952,91,0,15584,264,11,21336,264,0,21696,264,0,19952,264,0,15584,437,6,13985,437,0,19952,437,0,15584,610,6,15729,610,0,19952,610,0,15584,783,6,11872,783,0,21481,1092,6,14712,1802,0,21481,2658,0,21696,3550,0,21696};
 801
 802 static int vcp300[20]={300,514,88,28,8256,88,0,8272,440,8,8160,1800,0,10384};
 803
 804 int *wsr88d_get_vcp_info(int vcp_num,int el_num)
 805 {
 806 /*
 807  * This routine from Dan Austin.  Program component of nex2uf.
 808  */
 809     static int vcp_info[4];
 810     int fix_angle;
 811     int pulse_cnt;
 812     int az_rate;
 813     int pulse_width;
 814
 815     /* case statement to get vcp info */
 816     switch(vcp_num) {
 817     case 11:
 818       fix_angle =   vcp11[(3*el_num)-1];
 819       pulse_cnt =   vcp11[(3*el_num)];
 820       az_rate =     vcp11[(3*el_num)+1];
 821       pulse_width = vcp11[1];
 822       break;
 823     case 12:
 824       fix_angle =   vcp12[(3*el_num)-1];
 825       pulse_cnt =   vcp12[(3*el_num)];
 826       az_rate =     vcp12[(3*el_num)+1];
 827       pulse_width = vcp12[1];
 828       break;
 829     case 21:
 830       fix_angle =   vcp21[(3*el_num)-1];
 831       pulse_cnt =   vcp21[(3*el_num)];
 832       az_rate =     vcp21[(3*el_num)+1];
 833       pulse_width = vcp21[1];
 834       break;
 835     case 31:
 836       fix_angle =   vcp31[(3*el_num)-1];
 837       pulse_cnt =   vcp31[(3*el_num)];
 838       az_rate =     vcp31[(3*el_num)+1];
 839       pulse_width = vcp31[1];
 840       break;
 841     case 32:
 842       fix_angle =   vcp32[(3*el_num)-1];
 843       pulse_cnt =   vcp32[(3*el_num)];
 844       az_rate =     vcp32[(3*el_num)+1];
 845       pulse_width = vcp32[1];
 846       break;
 847     case 300:
 848       fix_angle =   vcp300[(3*el_num)-1];
 849       pulse_cnt =   vcp300[(3*el_num)];
 850       az_rate =     vcp300[(3*el_num)+1];
 851       pulse_width = vcp300[1];
 852       break;
 853     case 121:
 854       fix_angle =   vcp121[(3*el_num)-1];
 855       pulse_cnt =   vcp121[(3*el_num)];
 856       az_rate =     vcp121[(3*el_num)+1];
 857       pulse_width = vcp121[1];
 858       break;
 859     case 211:
 860       fix_angle =   vcp11[(3*el_num)-1];
 861       pulse_cnt =   vcp11[(3*el_num)];
 862       az_rate =     vcp11[(3*el_num)+1];
 863       pulse_width = vcp11[1];
 864       break;
 865     case 212:
 866       fix_angle =   vcp12[(3*el_num)-1];
 867       pulse_cnt =   vcp12[(3*el_num)];
 868       az_rate =     vcp12[(3*el_num)+1];
 869       pulse_width = vcp12[1];
 870       break;
 871     case 221:
 872       fix_angle =   vcp21[(3*el_num)-1];
 873       pulse_cnt =   vcp21[(3*el_num)];
 874       az_rate =     vcp21[(3*el_num)+1];
 875       pulse_width = vcp21[1];
 876       break;
 877     default:
 878       fix_angle  = 0;
 879       pulse_cnt  = 0;
 880       az_rate    = 0;
 881       pulse_width= 0;
 882       break;
 883     }
 884
 885     /* get array for output */
 886     vcp_info[0]=fix_angle;
 887     vcp_info[1]=pulse_cnt;
 888     vcp_info[2]=az_rate;
 889     vcp_info[3]=pulse_width;
 890
 891
 892     /* return the value array   */
 893     return(vcp_info);
 894 }
 895
 896
 897 float wsr88d_get_fix_angle(Wsr88d_ray *ray)
 898 {
 899   int *vcp_info;
 900   vcp_info = wsr88d_get_vcp_info(ray->vol_cpat, ray->elev_num);
 901   return vcp_info[0]/8.0*180./4096.0;
 902 }
 903 int wsr88d_get_pulse_count(Wsr88d_ray *ray)
 904 {
 905   int *vcp_info;
 906   vcp_info = wsr88d_get_vcp_info(ray->vol_cpat, ray->elev_num);
 907   return vcp_info[1];
 908 }
 909 float wsr88d_get_azimuth_rate(Wsr88d_ray *ray)
 910 {
 911   int *vcp_info;
 912   vcp_info = wsr88d_get_vcp_info(ray->vol_cpat, ray->elev_num);
 913   return vcp_info[2]/8.0*45./4096.0;
 914 }
 915 float wsr88d_get_pulse_width(Wsr88d_ray *ray)
 916 {
 917   int *vcp_info;
 918   vcp_info = wsr88d_get_vcp_info(ray->vol_cpat, ray->elev_num);
 919   return vcp_info[3]/299.792458;
 920 }
 921
 922 float wsr88d_get_prf(Wsr88d_ray *ray)
 923 {
 924   float prf;
 925   float c = 299792458.0;
 926   float range;
 927
 928   range = wsr88d_get_range(ray)*1000.0;
 929   if (range != 0) prf = c/(2*range);
 930   else prf = 0.0;
 931
 932   return prf;
 933 }
 934
 935 float wsr88d_get_prt(Wsr88d_ray *ray)
 936 {
 937   float prf;
 938   float prt;
 939
 940   prf = wsr88d_get_prf(ray);
 941   if (prf != 0) prt = 1.0/prf;
 942   else prt = 0;
 943   return prt;
 944 }
 945
 946 /* Note: wsr88d_get_wavelength() below is no longer used because of differences
 947  * in wavelength for velocity and reflectivity.  The function computes
 948  * wavelength when Nyquist is present, but returns a constant wavelength
 949  * otherwise.  Nyquist is present for velocity, but not for reflectivity.  The
 950  * fact is that WSR-88D radars use a constant wavelength, 10.7 cm., which is
 951  * the value now used where this function was formerly called in
 952  * wsr88d_load_sweep_into_volume().
 953  */
 954
 955 float wsr88d_get_wavelength(Wsr88d_ray *ray)
 956 {
 957   float wavelength;
 958   float prf;
 959   float nyquist;
 960
 961   prf = wsr88d_get_prf(ray);
 962   nyquist = wsr88d_get_nyquist(ray);
 963     /* If required info to determine wavelength does not exist,
 964          just use 10 cm. All wsr88d radars are 10cm. MJK */
 965   if ((prf == 0) || (nyquist == 0.0)) wavelength = 0.10;
 966   else wavelength = 4*nyquist/prf;
 967   return wavelength;
 968 }
 969
 970 float wsr88d_get_frequency(Wsr88d_ray *ray)
 971 {
 972   float freq;
 973   float c = 299792458.0;
 974
 975     /* Carrier freq (GHz). Revised 12 Jun 97. MJK */
 976     freq = (c / wsr88d_get_wavelength(ray)) * 1.0e-9;
 977   return freq;
 978 }
 979