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         /* Load this ray into the sweep. */
 456         ray_num = wsr88d_ray.ray_num;
 457         if (wsr88d_sweep->ray[ray_num] == NULL) {
 458           wsr88d_sweep->ray[ray_num] = (Wsr88d_ray *) malloc (sizeof(Wsr88d_ray));
 459         }
 460         memcpy(wsr88d_sweep->ray[ray_num], &wsr88d_ray, sizeof(Wsr88d_ray));
 461   }
 462
 463   /* Just to be safe, clear all ray pointers left in this sweep to
 464    * the maximum MAX_RAYS_IN_SWEEP.  This is required when the
 465    * wsr88d_sweep is reused and not cleared.
 466    */
 467   free_and_clear_sweep(wsr88d_sweep, ray_num+1, MAX_RAYS_IN_SWEEP);
 468
 469 /*
 470   fprintf(stderr,"Processed %d records for elevation number %d\n",
 471                  nrec+1, wsr88d_ray.elev_num);
 472   wsr88d_print_sweep_info(wsr88d_sweep);
 473 */
 474   return nrec;
 475 }
 476
 477 /**********************************************************************/
 478 /*                                                                    */
 479 /*                      wsr88d_swap_ray                               */
 480 /*                                                                    */
 481 /**********************************************************************/
 482 void wsr88d_swap_ray(Wsr88d_ray *wsr88d_ray)
 483 {
 484   short *half_word;
 485   half_word = (short *)wsr88d_ray;
 486   for (; half_word<(short *)&wsr88d_ray->msg_time; half_word++)
 487         swap_2_bytes(half_word);
 488
 489   swap_4_bytes(&wsr88d_ray->msg_time);
 490   swap_2_bytes(&wsr88d_ray->num_seg);
 491   swap_2_bytes(&wsr88d_ray->seg_num);
 492   swap_4_bytes(&wsr88d_ray->ray_time);
 493
 494   half_word = (short *) &wsr88d_ray->ray_date;
 495   for (; half_word<(short *)&wsr88d_ray->sys_cal; half_word++)
 496         swap_2_bytes(half_word);
 497
 498   swap_4_bytes(&wsr88d_ray->sys_cal);
 499
 500   half_word = (short *) &wsr88d_ray->refl_ptr;
 501   for (; half_word<(short *)&wsr88d_ray->data[0]; half_word++)
 502         swap_2_bytes(half_word);
 503
 504 }
 505
 506 /**********************************************************************/
 507 /*                                                                    */
 508 /*  done 2/28           wsr88d_read_ray                               */
 509 /*                                                                    */
 510 /**********************************************************************/
 511 int wsr88d_read_ray(Wsr88d_file *wf, Wsr88d_ray *wsr88d_ray)
 512 {
 513   int n;
 514   n = fread(wsr88d_ray, sizeof(Wsr88d_ray), 1, wf->fptr);
 515 /*  if (n > 0) print_packet_info(wsr88d_ray); */
 516   if (little_endian())
 517         wsr88d_swap_ray(wsr88d_ray);
 518
 519   return n;
 520 }
 521
 522 /**********************************************************************/
 523 /*                                                                    */
 524 /*  not done N/A     wsr88d_read_ray_header                           */
 525 /*                                                                    */
 526 /**********************************************************************/
 527 int wsr88d_read_ray_header(Wsr88d_file *wf,
 528                                                    Wsr88d_ray_header *wsr88d_ray_header)
 529 {
 530   fprintf(stderr,"Stub routine: wsr88d_read_ray_header.\n");
 531   return 0;
 532 }
 533
 534 /**********************************************************************/
 535 /*                                                                    */
 536 /*  done 3/3         wsr88d_ray_to_float                              */
 537 /*                                                                    */
 538 /**********************************************************************/
 539 int wsr88d_ray_to_float(Wsr88d_ray *ray,
 540                                                 int THE_DATA_WANTED, float v[], int *n)
 541 {
 542 /*
 543  *  Input: *ray             -  WSR-88D packet
 544  * Output: THE_DATA_WANTED  -  Indicates which field to convert.  Fields:
 545  *                             WSR88D_DZ, WSR88D_VR, WSR88D_SW
 546  *         v[]              -  The output vector of float values.
 547  *         n                -  Length of the output vector.  0 = no data.
 548  *
 549  * Returns n.
 550  *
 551  * No allocation of space for the output vector performed here.
 552  */
 553
 554 /* Code from Dan Austin (cvt_pckt_data.c) was the template for this. */
 555
 556   /* declarations       */
 557   int num_ref_gates,num_vel_gates,num_spec_gates;
 558   int refl_ptr, vel_ptr,spec_ptr,res_flag;
 559   int ival;
 560   int i;
 561
 562   *n = 0;
 563   num_ref_gates  = ray->num_refl;
 564   num_vel_gates  = ray->num_dop;
 565   num_spec_gates = ray->num_dop;  /* 'num_dop', this is not a typo. */
 566
 567 /* The data pointers are specified from the begining of the
 568  * 'Digital Radar Data (Message) Header'.  Since we have a structure
 569  * that defines all the header variables and a member called 'data'.
 570  * we must subtract the length of the 'message header' from the data
 571  * pointer.  Hopefully, the reflecivity pointer will be 0 meaning the
 572  * first element of the 'data' member; ray->data[0];
 573  */
 574 #define LENGTH_OF_MESSAGE 100
 575   if (num_ref_gates > 0) refl_ptr = ray->refl_ptr - LENGTH_OF_MESSAGE;
 576   else refl_ptr = 0;
 577
 578   vel_ptr  = ray->vel_ptr - LENGTH_OF_MESSAGE;
 579   spec_ptr = ray->spc_ptr - LENGTH_OF_MESSAGE;
 580
 581   res_flag = ray->vel_res;
 582
 583
 584 /*
 585   fprintf(stderr,"refl_ptr = %d  #g = %d, ", refl_ptr, num_ref_gates);
 586   fprintf(stderr," vel_ptr = %d  #g = %d, ", vel_ptr, num_vel_gates);
 587   fprintf(stderr,"spec_ptr = %d  #g = %d, ", spec_ptr, num_spec_gates);
 588   fprintf(stderr,"res_flag = %d\n", res_flag);
 589 */
 590
 591   if (THE_DATA_WANTED == WSR88D_DZ) {
 592         /* do the reflectivity data  (dbZ)*/
 593         if (refl_ptr+num_ref_gates > 2300)
 594           fprintf(stderr, "WARNING: # refl index (%d) exceeds maximum (2300)\n",
 595                           refl_ptr+num_ref_gates);
 596         else {
 597         for(i=0; i<num_ref_gates; i++) {
 598           ival = ray->data[refl_ptr+i];
 599           if(ival > 1)
 600                   v[i] = (((ival-2.0)/2.0)-32.0);
 601           else if (ival == 1)
 602                 v[i] = WSR88D_RFVAL;
 603           else /* ival = 0 */
 604                 v[i] = WSR88D_BADVAL;
 605         }
 606         *n = num_ref_gates;
 607         }
 608
 609   } else if (THE_DATA_WANTED == WSR88D_VR) {
 610         /* do the velocity data  (M/S) */
 611         if (vel_ptr+num_vel_gates > 2300)
 612           fprintf(stderr, "WARNING: # vel index (%d) exceeds maximum (2300)\n",
 613                           vel_ptr+num_vel_gates);
 614         else {
 615         for(i=0; i<num_vel_gates;i++)   {
 616           ival = ray->data[vel_ptr+i];
 617           if(ival > 1)
 618                 if (res_flag == 2) /* High resolution: 0.5 m/s */
 619                   v[i] = (((ival-2.0)/2.0)-63.5);
 620                 else
 621                   v[i] = ((ival-2.0)-127.0);
 622           else if (ival == 1)
 623                 v[i] = WSR88D_RFVAL;
 624           else /* ival = 0 */
 625                 v[i] = WSR88D_BADVAL;
 626         }
 627         *n = num_vel_gates;
 628         }
 629
 630   } else if (THE_DATA_WANTED == WSR88D_SW) {
 631         /* now do the spectrum width data (M/S)*/
 632         if (spec_ptr+num_spec_gates > 2300)
 633           fprintf(stderr, "WARNING: # spec index (%d) exceeds maximum (2300)\n",
 634                           spec_ptr+num_spec_gates);
 635         else {
 636         for(i=0;i<num_spec_gates;i++) {
 637           ival = ray->data[spec_ptr+i];
 638                 if(ival > 1)
 639                   v[i] = (((ival-2)/2.0)-63.5);
 640                 else if (ival == 1)
 641                   v[i] = WSR88D_RFVAL;
 642                 else /* ival = 0 */
 643                   v[i] = WSR88D_BADVAL;
 644         }
 645         *n = num_spec_gates;
 646         }
 647   }
 648
 649   return *n;
 650 }
 651
 652
 653
 654 /**********************************************************************/
 655 /*        Functions that convert some message header values.          */
 656 /**********************************************************************/
 657 /**********************************************************************/
 658 /*                                                                    */
 659 /*  done 3/3   float wsr88d_get_nyquist                               */
 660 /*  done 3/3   float wsr88d_get_atmos_atten_factor                    */
 661 /*  done 3/3   float wsr88d_get_velocity_resolution                   */
 662 /*  done 3/3   int   wsr88d_get_volume_coverage                       */
 663 /*  done 3/3   float wsr88d_get_elevation_angle                       */
 664 /*  done 3/3   float wsr88d_get_azimuth                               */
 665 /*  done 3/3   float wsr88d_get_range                                 */
 666 /*  done 3/3   void  wsr88d_get_date                                  */
 667 /*  done 3/3   void  wsr88d_get_time                                  */
 668 /*  done 5/20  int  *wsr88d_get_vcp_info                              */
 669 /*  done 5/20  float wsr88d_get_fix_angle                             */
 670 /*  done 5/20  int   wsr88d_get_pulse_count                           */
 671 /*  done 5/20  float wsr88d_get_azimuth_rate                          */
 672 /*  done 5/20  float wsr88d_get_pulse_width                           */
 673 /*  done 5/20  float wsr88d_get_prf                                   */
 674 /*  done 5/20  float wsr88d_get_prt                                   */
 675 /*  done 5/20  float wsr88d_get_wavelength                            */
 676 /*  done 5/20  float wsr88d_get_frequency                             */
 677 /*                                                                    */
 678 /**********************************************************************/
 679 float wsr88d_get_nyquist(Wsr88d_ray *ray)
 680 {
 681   return ray->nyq_vel/100.0;
 682 }
 683
 684 float wsr88d_get_atmos_atten_factor(Wsr88d_ray *ray)
 685 {
 686   return ray->atm_att/1000.0;
 687 }
 688
 689 float wsr88d_get_velocity_resolution(Wsr88d_ray *ray)
 690 {
 691   if (ray->vel_res == 2) return 0.5;
 692   return 0.0;
 693 }
 694
 695 int   wsr88d_get_volume_coverage(Wsr88d_ray *ray)
 696 {
 697   if (ray->vol_cpat == 11) return 11;
 698   if (ray->vol_cpat == 12) return 12;
 699   if (ray->vol_cpat == 21) return 21;
 700   if (ray->vol_cpat == 31) return 31;
 701   if (ray->vol_cpat == 32) return 32;
 702   if (ray->vol_cpat == 121) return 121;
 703   return 0;
 704 }
 705
 706 float wsr88d_get_elevation_angle(Wsr88d_ray *ray)
 707 {
 708   return ray->elev/8.0*(180.0/4096.0);
 709 }
 710
 711 float wsr88d_get_azimuth(Wsr88d_ray *ray)
 712 {
 713   return ray->azm/8.0*(180.0/4096.0);
 714 }
 715
 716 float wsr88d_get_range(Wsr88d_ray *ray)
 717 {
 718   return ray->unam_rng/10.0;
 719 }
 720
 721 #include <time.h>
 722 void wsr88d_get_date(Wsr88d_ray *ray, int *mm, int *dd, int *yy)
 723 {
 724 /*
 725  * mm (1-12)
 726  * dd (1-31)
 727  * yy (ex. 93)
 728  */
 729   time_t itime;
 730   struct tm *tm_time;
 731   if (ray == NULL) {
 732     *mm = *dd = *yy = 0;
 733     return;
 734   }
 735
 736   itime = ray->ray_date - 1;
 737   itime *= 24*60*60; /* Seconds/day * days. */
 738
 739   tm_time = gmtime(&itime);
 740   *mm = tm_time->tm_mon+1;
 741   *dd = tm_time->tm_mday;
 742   *yy = tm_time->tm_year;
 743 }
 744
 745 void wsr88d_get_time(Wsr88d_ray *ray, int *hh, int *mm, int *ss, float *fsec)
 746 {
 747   /*
 748    * hh (0-23)
 749    * mm (0-59)
 750    * ss (0-59)
 751    * fsec (fraction of second)
 752    */
 753   double t;
 754
 755   if (ray == NULL) {
 756     *hh = *mm = *ss = *fsec = 0;
 757     return;
 758   }
 759   t = ray->ray_time;
 760   t /= 1000.0;
 761   *hh = t/3600;
 762   t -= *hh*3600;
 763   *mm = t/60;
 764   t -= *mm*60;
 765   *ss = (int)t;
 766   *fsec = t - *ss;
 767 }
 768
 769
 770
 771 /*
 772  * Get_vcp_info - gets info about the volume coverage pattern for this scan
 773  * parameters:
 774  * int vcp_num - volume coverage pattern number
 775  * int el_num - elevation number w/in vcp
 776  * returns: int *vcp_info - ptr to array w/vcp info
 777  * calls from: Nexrad2uf
 778  * calls to: none
 779  */
 780
 781 /* this database contains volume coverage patterns & associated info:   */
 782 /* (0)= vcp # (1)=pulse width for vcp  "Id$"                            */
 783 /* line[1-n]: (n,0)= elev. # for vcp (n,1)= (fixed angle)*8*(4096/180)  */
 784 /* (n,2)= pulse count (n,3)= (azimuthal sweep rate)*8*(4096/45)         */
 785
 786 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};
 787
 788 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};
 789
 790 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};
 791
 792 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};
 793
 794 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};
 795
 796 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};
 797
 798 static int vcp300[20]={300,514,88,28,8256,88,0,8272,440,8,8160,1800,0,10384};
 799
 800 int *wsr88d_get_vcp_info(int vcp_num,int el_num)
 801 {
 802 /*
 803  * This routine from Dan Austin.  Program component of nex2uf.
 804  */
 805         static int vcp_info[4];
 806         int fix_angle;
 807         int pulse_cnt;
 808         int az_rate;
 809         int pulse_width;
 810
 811         /* case statement to get vcp info */
 812         switch(vcp_num) {
 813         case 11:
 814           fix_angle =   vcp11[(3*el_num)-1];
 815           pulse_cnt =   vcp11[(3*el_num)];
 816           az_rate =     vcp11[(3*el_num)+1];
 817           pulse_width = vcp11[1];
 818           break;
 819         case 12:
 820           fix_angle =   vcp12[(3*el_num)-1];
 821           pulse_cnt =   vcp12[(3*el_num)];
 822           az_rate =     vcp12[(3*el_num)+1];
 823           pulse_width = vcp12[1];
 824           break;
 825         case 21:
 826           fix_angle =   vcp21[(3*el_num)-1];
 827           pulse_cnt =   vcp21[(3*el_num)];
 828           az_rate =     vcp21[(3*el_num)+1];
 829           pulse_width = vcp21[1];
 830           break;
 831         case 31:
 832           fix_angle =   vcp31[(3*el_num)-1];
 833           pulse_cnt =   vcp31[(3*el_num)];
 834           az_rate =     vcp31[(3*el_num)+1];
 835           pulse_width = vcp31[1];
 836           break;
 837         case 32:
 838           fix_angle =   vcp32[(3*el_num)-1];
 839           pulse_cnt =   vcp32[(3*el_num)];
 840           az_rate =     vcp32[(3*el_num)+1];
 841           pulse_width = vcp32[1];
 842           break;
 843         case 300:
 844           fix_angle =   vcp300[(3*el_num)-1];
 845           pulse_cnt =   vcp300[(3*el_num)];
 846           az_rate =     vcp300[(3*el_num)+1];
 847           pulse_width = vcp300[1];
 848           break;
 849         case 121:
 850           fix_angle =   vcp121[(3*el_num)-1];
 851           pulse_cnt =   vcp121[(3*el_num)];
 852           az_rate =     vcp121[(3*el_num)+1];
 853           pulse_width = vcp121[1];
 854           break;
 855         case 211:
 856           fix_angle =   vcp11[(3*el_num)-1];
 857           pulse_cnt =   vcp11[(3*el_num)];
 858           az_rate =     vcp11[(3*el_num)+1];
 859           pulse_width = vcp11[1];
 860           break;
 861         case 212:
 862           fix_angle =   vcp12[(3*el_num)-1];
 863           pulse_cnt =   vcp12[(3*el_num)];
 864           az_rate =     vcp12[(3*el_num)+1];
 865           pulse_width = vcp12[1];
 866           break;
 867         case 221:
 868           fix_angle =   vcp21[(3*el_num)-1];
 869           pulse_cnt =   vcp21[(3*el_num)];
 870           az_rate =     vcp21[(3*el_num)+1];
 871           pulse_width = vcp21[1];
 872           break;
 873         default:
 874           fix_angle  = 0;
 875           pulse_cnt  = 0;
 876           az_rate    = 0;
 877           pulse_width= 0;
 878           break;
 879         }
 880
 881         /* get array for output */
 882         vcp_info[0]=fix_angle;
 883         vcp_info[1]=pulse_cnt;
 884         vcp_info[2]=az_rate;
 885         vcp_info[3]=pulse_width;
 886
 887
 888         /* return the value array       */
 889         return(vcp_info);
 890 }
 891
 892
 893 float wsr88d_get_fix_angle(Wsr88d_ray *ray)
 894 {
 895   int *vcp_info;
 896   vcp_info = wsr88d_get_vcp_info(ray->vol_cpat, ray->elev_num);
 897   return vcp_info[0]/8.0*180./4096.0;
 898 }
 899 int wsr88d_get_pulse_count(Wsr88d_ray *ray)
 900 {
 901   int *vcp_info;
 902   vcp_info = wsr88d_get_vcp_info(ray->vol_cpat, ray->elev_num);
 903   return vcp_info[1];
 904 }
 905 float wsr88d_get_azimuth_rate(Wsr88d_ray *ray)
 906 {
 907   int *vcp_info;
 908   vcp_info = wsr88d_get_vcp_info(ray->vol_cpat, ray->elev_num);
 909   return vcp_info[2]/8.0*45./4096.0;
 910 }
 911 float wsr88d_get_pulse_width(Wsr88d_ray *ray)
 912 {
 913   int *vcp_info;
 914   vcp_info = wsr88d_get_vcp_info(ray->vol_cpat, ray->elev_num);
 915   return vcp_info[3]/299.792458;
 916 }
 917
 918 float wsr88d_get_prf(Wsr88d_ray *ray)
 919 {
 920   float prf;
 921   float c = 299792458.0;
 922   float range;
 923
 924   range = wsr88d_get_range(ray)*1000.0;
 925   if (range != 0) prf = c/(2*range);
 926   else prf = 0.0;
 927
 928   return prf;
 929 }
 930
 931 float wsr88d_get_prt(Wsr88d_ray *ray)
 932 {
 933   float prf;
 934   float prt;
 935
 936   prf = wsr88d_get_prf(ray);
 937   if (prf != 0) prt = 1.0/prf;
 938   else prt = 0;
 939   return prt;
 940 }
 941
 942 /* Note: wsr88d_get_wavelength() below is no longer used because of differences
 943  * in wavelength for velocity and reflectivity.  The function computes
 944  * wavelength when Nyquist is present, but returns a constant wavelength
 945  * otherwise.  Nyquist is present for velocity, but not for reflectivity.  The
 946  * fact is that WSR-88D radars use a constant wavelength, 10.7 cm., which is
 947  * the value now used where this function was formerly called in
 948  * wsr88d_load_sweep_into_volume().
 949  */
 950
 951 float wsr88d_get_wavelength(Wsr88d_ray *ray)
 952 {
 953   float wavelength;
 954   float prf;
 955   float nyquist;
 956
 957   prf = wsr88d_get_prf(ray);
 958   nyquist = wsr88d_get_nyquist(ray);
 959         /* If required info to determine wavelength does not exist,
 960                  just use 10 cm. All wsr88d radars are 10cm. MJK */
 961   if ((prf == 0) || (nyquist == 0.0)) wavelength = 0.10;
 962   else wavelength = 4*nyquist/prf;
 963   return wavelength;
 964 }
 965
 966 float wsr88d_get_frequency(Wsr88d_ray *ray)
 967 {
 968   float freq;
 969   float c = 299792458.0;
 970
 971         /* Carrier freq (GHz). Revised 12 Jun 97. MJK */
 972         freq = (c / wsr88d_get_wavelength(ray)) * 1.0e-9;
 973   return freq;
 974 }
 975