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

   1 /*
   2     NASA/TRMM, Code 613
   3     This is the TRMM Office Radar Software Library.
   4     Copyright (C) 2008
   5             Bart Kelley
   6             SSAI
   7             Lanham, Maryland
   8
   9     This program is free software: you can redistribute it and/or modify
  10     it under the terms of the GNU General Public License as published by
  11     the Free Software Foundation, either version 3 of the License, or
  12     (at your option) any later version.
  13
  14     This program 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
  17     GNU General Public License for more details.
  18
  19     You should have received a copy of the GNU General Public License
  20     along with this program.  If not, see <http://www.gnu.org/licenses/>.
  21 */
  22
  23
  24 /*
  25  * This file contains routines for processing Message Type 31, the digital
  26  * radar message type introduced in WSR-88D Level II Build 10.
  27  */
  28
  29 #include "rsl.h"
  30 #include "wsr88d.h"
  31
  32 /* Data descriptions in the following data structures are from the "Interface
  33  * Control Document for the RDA/RPG", Build 10.0 Draft, WSR-88D Radar
  34  * Operations Center.
  35  */
  36
  37 typedef struct {
  38     short rpg[6]; /* Unused.  Not really part of message header, but is
  39                    * inserted by RPG Communications Manager for each message.
  40                    */
  41     unsigned short msg_size;  /* for this segment, in halfwords */
  42     unsigned char  channel;   /* RDA Redundant Channel */
  43     unsigned char  msg_type;  /* Message Type */
  44     unsigned short id_seq;    /* I.d. Seq = 0 to 7FFF, then roll over to 0 */
  45     unsigned short msg_date;  /* modified Julian date from 1/1/70 */
  46     unsigned int   msg_time;  /* packet generation time in ms past midnite */
  47     unsigned short num_segs;
  48     unsigned short seg_num;
  49 } Wsr88d_msg_hdr;
  50
  51 typedef struct {
  52     char radar_id[4];
  53     unsigned int   ray_time; /* Data collection time in milliseconds past midnight GMT */
  54     unsigned short ray_date; /* Julian date - 2440586.5 (1/01/1970) */
  55     unsigned short azm_num ; /* Radial number within elevation scan */
  56     float azm;      /* Azimuth angle in degrees (0 to 359.956055) */
  57     unsigned char compression_code; /* 0 = uncompressed, 1 = BZIP2, 2 = zlib */
  58     unsigned char spare; /* for word alignment */
  59     unsigned short radial_len; /* radial length in bytes, including data header block */
  60     unsigned char azm_res;
  61     unsigned char radial_status;
  62     unsigned char elev_num;
  63     unsigned char sector_cut_num;
  64     float elev;  /* Elevation angle in degrees (-7.0 to 70.0) */
  65     unsigned char radial_spot_blanking;
  66     unsigned char azm_indexing_mode;
  67     unsigned short data_block_count;
  68     /* Data Block Pointers */
  69     unsigned int dbptr_vol_const;
  70     unsigned int dbptr_elev_const;
  71     unsigned int dbptr_radial_const;
  72     unsigned int dbptr_ref;
  73     unsigned int dbptr_vel;
  74     unsigned int dbptr_sw;
  75     unsigned int dbptr_zdr;
  76     unsigned int dbptr_phi;
  77     unsigned int dbptr_rho;
  78 } Ray_header_m31;  /* Called Data Header Block in RDA/RPG document. */
  79
  80 typedef struct {
  81     char dataname[4];
  82     unsigned int  reserved;
  83     unsigned short ngates;
  84     short range_first_gate;
  85     short range_samp_interval;
  86     short thresh_not_overlayed;
  87     short snr_thresh;
  88     unsigned char controlflag;
  89     unsigned char datasize_bits;
  90     float scale;
  91     float offset;
  92 } Data_moment_hdr;
  93
  94 typedef struct {
  95     Data_moment_hdr data_hdr;
  96     /* TODO: data type will need to be changed to "void" to handle Dual Pole:
  97      * some data are in 2-byte words.
  98      */
  99     unsigned char *data;
 100 } Data_moment;
 101
 102 typedef struct {
 103     Ray_header_m31 ray_hdr;
 104     short unamb_rng;
 105     short nyq_vel;
 106     /* Pointers to data */
 107     /* TODO: Replace data moment pointers with single pointer to radial.
 108      *       Can simply make unsigned char array as in MSG1 version.
 109      *       Maximum radial length is 14288 bytes.
 110      */
 111     Data_moment *ref;
 112     Data_moment *vel;
 113     Data_moment *sw;
 114 } Wsr88d_ray_m31;
 115
 116 #define MAXRAYS_M31 800
 117 #define MAXSWEEPS 20
 118
 119 enum {START_OF_ELEV, INTERMED_RADIAL, END_OF_ELEV, BEGIN_VOS, END_VOS};
 120
 121
 122 void wsr88d_swap_m31_hdr(Wsr88d_msg_hdr *msghdr)
 123 {
 124     swap_2_bytes(&msghdr->msg_size);
 125     swap_2_bytes(&msghdr->id_seq);
 126     swap_2_bytes(&msghdr->msg_date);
 127     swap_4_bytes(&msghdr->msg_time);
 128     swap_2_bytes(&msghdr->num_segs);
 129     swap_2_bytes(&msghdr->seg_num);
 130 }
 131
 132
 133 void wsr88d_swap_m31_ray(Ray_header_m31 *wsr88d_ray)
 134 {
 135     int *fullword;
 136
 137     swap_4_bytes(&wsr88d_ray->ray_time);
 138     swap_2_bytes(&wsr88d_ray->ray_date);
 139     swap_2_bytes(&wsr88d_ray->azm_num);
 140     swap_4_bytes(&wsr88d_ray->azm);
 141     swap_2_bytes(&wsr88d_ray->radial_len);
 142     swap_4_bytes(&wsr88d_ray->elev);
 143     swap_2_bytes(&wsr88d_ray->data_block_count);
 144     fullword = (int *) &wsr88d_ray->dbptr_vol_const;
 145     for (; fullword <= (int *) &wsr88d_ray->dbptr_rho; fullword++)
 146         swap_4_bytes(fullword);
 147 }
 148
 149
 150 void wsr88d_swap_data_moment(Data_moment_hdr *this_field)
 151 {
 152     short *halfword;
 153     halfword = (short *) &this_field->ngates;
 154     for (; halfword < (short *) &this_field->controlflag; halfword++)
 155         swap_2_bytes(halfword);
 156     swap_4_bytes(&this_field->scale);
 157     swap_4_bytes(&this_field->offset);
 158 }
 159
 160
 161 void testprt(Ray_header_m31 ray_hdr)
 162 {
 163     /* For testing: Print some values from data block header. */
 164     printf("\nray_time: %d\n",ray_hdr.ray_time);
 165     printf("ray_date: %d\n",ray_hdr.ray_date);
 166     printf("azm_num: %d\n",ray_hdr.azm_num);
 167     printf("azm: %f\n",ray_hdr.azm);
 168     printf("radial_len: %d\n",ray_hdr.radial_len);
 169     printf("elev: %f\n",ray_hdr.elev);
 170     printf("data block count: %d\n",ray_hdr.data_block_count);
 171     printf("dbptr_vol_const: %d\n",ray_hdr.dbptr_vol_const);
 172     printf("dbptr_elev_const: %d\n",ray_hdr.dbptr_elev_const);
 173     printf("dbptr_radial_const: %d\n",ray_hdr.dbptr_radial_const);
 174     printf("dbptr_ref: %d\n",ray_hdr.dbptr_ref);
 175     printf("dbptr_vel: %d\n",ray_hdr.dbptr_vel);
 176     printf("dbptr_sw: %d\n",ray_hdr.dbptr_sw);
 177 }
 178
 179
 180 float wsr88d_get_angle(short bitfield)
 181 {
 182     short mask = 1;
 183     int i;
 184     float angle = 0.;
 185     float value[13] = {0.043945, 0.08789, 0.17578, 0.35156, .70313, 1.40625,
 186         2.8125, 5.625, 11.25, 22.5, 45., 90., 180.};
 187
 188     /* find which bits are set and sum corresponding values to get angle. */
 189
 190     bitfield = bitfield >> 3;  /* 3 least significant bits aren't used. */
 191     for (i = 0; i < 13; i++) {
 192         if (bitfield & mask) angle += value[i];
 193         bitfield = bitfield >> 1;
 194     }
 195     return angle;
 196 }
 197
 198
 199 float wsr88d_get_azim_rate(short bitfield)
 200 {
 201     short mask = 1;
 202     int i;
 203     float rate = 0.;
 204     float value[12] = {0.0109863, 0.021972656, 0.043945, 0.08789, 0.17578,
 205         0.35156, .70313, 1.40625, 2.8125, 5.625, 11.25, 22.5};
 206
 207     /* Find which bits are set and sum corresponding values to get rate. */
 208
 209     bitfield = bitfield >> 3;  /* 3 least significant bits aren't used. */
 210     for (i = 0; i < 12; i++) {
 211         if (bitfield & mask) rate += value[i];
 212         bitfield = bitfield >> 1;
 213     }
 214     if (bitfield >> 15) rate = -rate;
 215     return rate;
 216 }
 217
 218 #define WSR88D_MAX_SWEEPS 20
 219
 220 typedef struct {
 221     int vcp;
 222     int num_cuts;
 223     float vel_res;
 224     float fixed_angle[WSR88D_MAX_SWEEPS];
 225     float azim_rate[WSR88D_MAX_SWEEPS];
 226     int waveform[WSR88D_MAX_SWEEPS];
 227     int super_res_ctrl[WSR88D_MAX_SWEEPS];
 228     int surveil_prf_num[WSR88D_MAX_SWEEPS];
 229     int doppler_prf_num[WSR88D_MAX_SWEEPS];
 230 } VCP_data;
 231
 232 static VCP_data vcp_data;
 233
 234 void wsr88d_get_vcp_data(short *msgtype5)
 235 {
 236     short azim_rate, fixed_angle, vel_res;
 237     short sres_and_survprf; /* super res ctrl and surveil prf, one byte each */
 238     short chconf_and_waveform;
 239     int i;
 240
 241     vcp_data.vcp = (unsigned short) msgtype5[2];
 242     vcp_data.num_cuts = msgtype5[3];
 243     if (little_endian()) {
 244         swap_2_bytes(&vcp_data.vcp);
 245         swap_2_bytes(&vcp_data.num_cuts);
 246     }
 247     vel_res = msgtype5[5];
 248     if (little_endian()) swap_2_bytes(&vel_res);
 249     vel_res = vel_res >> 8;
 250     if (vel_res == 2) vcp_data.vel_res = 0.5;
 251     else if (vel_res == 4) vcp_data.vel_res = 1.0;
 252     else vcp_data.vel_res = 0.0;
 253     /* Get elevation related information for each sweep. */
 254     for (i=0; i < vcp_data.num_cuts; i++) {
 255         fixed_angle = msgtype5[11 + i*23];
 256         azim_rate = msgtype5[15 + i*23];
 257         chconf_and_waveform = msgtype5[12 + i*23];
 258         sres_and_survprf = msgtype5[13 + i*23];
 259         vcp_data.doppler_prf_num[i] = msgtype5[23 + i*23];
 260         if (little_endian()) {
 261             swap_2_bytes(&fixed_angle);
 262             swap_2_bytes(&azim_rate);
 263             swap_2_bytes(&chconf_and_waveform);
 264             swap_2_bytes(&sres_and_survprf);
 265             swap_2_bytes(&vcp_data.doppler_prf_num[i]);
 266         }
 267         vcp_data.fixed_angle[i] = wsr88d_get_angle(fixed_angle);
 268         vcp_data.azim_rate[i] = wsr88d_get_azim_rate(azim_rate);
 269         vcp_data.waveform[i] = chconf_and_waveform & 0xff;
 270         vcp_data.super_res_ctrl[i] = sres_and_survprf >> 8;
 271         vcp_data.surveil_prf_num[i] = sres_and_survprf & 0xff;
 272     }
 273 }
 274
 275
 276 Data_moment *read_data_moment(unsigned char *radial, int begin_data_block)
 277 {
 278     Data_moment *this_field;
 279     unsigned char* data;
 280     unsigned char* dptr;
 281     int hdr_size, n;
 282
 283     /*
 284     printf("read_data_moment: ");
 285     printf("sizeof(Data_moment) = %d\n", sizeof(Data_moment));
 286     printf("sizeof(Data_moment_hdr) = %d\n", sizeof(Data_moment_hdr));
 287     */
 288     this_field = malloc(sizeof(Data_moment));
 289     hdr_size = sizeof(Data_moment_hdr);
 290     /* printf("hdr_size = %d\n", hdr_size); */
 291
 292     dptr = radial + begin_data_block;
 293     memcpy(&this_field->data_hdr, dptr, hdr_size);
 294     dptr += hdr_size;
 295
 296     if (little_endian()) wsr88d_swap_data_moment(&this_field->data_hdr);
 297     data = (unsigned char *) malloc(this_field->data_hdr.ngates);
 298     memcpy(data, dptr, this_field->data_hdr.ngates);
 299     this_field->data = data;
 300     return this_field;
 301 }
 302
 303
 304 int read_wsr88d_ray_m31(Wsr88d_file *wf, int msg_size,
 305         Wsr88d_ray_m31 *wsr88d_ray)
 306 {
 307     Ray_header_m31 ray_hdr;
 308     int n, begin_data_block, bytes_to_read, ray_hdr_len;
 309     unsigned char *radial;
 310
 311     /* Read ray data header block */
 312     n = fread(&wsr88d_ray->ray_hdr, sizeof(Ray_header_m31), 1, wf->fptr);
 313     if (n < 1) {
 314         fprintf(stderr,"read_wsr88d_ray_m31: Read failed.\n");
 315         return 0;
 316     }
 317     if (little_endian()) wsr88d_swap_m31_ray(&wsr88d_ray->ray_hdr);
 318
 319     ray_hdr = wsr88d_ray->ray_hdr;
 320
 321     /*
 322     if (ray_hdr.azm_num == 1) testprt(wsr88d_ray->ray_hdr);
 323     */
 324
 325     ray_hdr_len = sizeof(ray_hdr);
 326
 327     /* Read in radial
 328      * Use header_size offset with data pointers
 329      * Copy data into structures
 330      */
 331
 332     bytes_to_read = msg_size - ray_hdr_len;
 333     radial = (unsigned char *) malloc(bytes_to_read);
 334     n = fread(radial, bytes_to_read, 1, wf->fptr);
 335     if (n < 1) {
 336         fprintf(stderr,"read_wsr88d_ray_m31: Read failed.\n");
 337         return 0;
 338     }
 339
 340
 341     if (ray_hdr.dbptr_radial_const != 0) {
 342         begin_data_block = ray_hdr.dbptr_radial_const - ray_hdr_len;
 343         memcpy(&wsr88d_ray->unamb_rng, &radial[begin_data_block+6], 2);
 344         memcpy(&wsr88d_ray->nyq_vel, &radial[begin_data_block+16], 2);
 345         if (little_endian()) {
 346             swap_2_bytes(&wsr88d_ray->unamb_rng);
 347             swap_2_bytes(&wsr88d_ray->nyq_vel);
 348         }
 349     }
 350     else {
 351         wsr88d_ray->unamb_rng = 0;
 352         wsr88d_ray->nyq_vel = 0;
 353     }
 354
 355     if (ray_hdr.dbptr_ref != 0) {
 356         begin_data_block = ray_hdr.dbptr_ref - ray_hdr_len;
 357         wsr88d_ray->ref = read_data_moment(radial, begin_data_block);
 358     }
 359     if (ray_hdr.dbptr_vel != 0) {
 360         begin_data_block = ray_hdr.dbptr_vel - ray_hdr_len;
 361         wsr88d_ray->vel = read_data_moment(radial, begin_data_block);
 362     }
 363     if (ray_hdr.dbptr_sw != 0) {
 364         begin_data_block = ray_hdr.dbptr_sw - ray_hdr_len;
 365         wsr88d_ray->sw = read_data_moment(radial, begin_data_block);
 366     }
 367
 368
 369     /* For testing: print reflectivity data. */
 370     int prtdata = 0;
 371     { int i;
 372         if (prtdata) {
 373             for (i=0; i < wsr88d_ray->ref->data_hdr.ngates; i++) {
 374                 if (i % 10 == 0) printf("\n");
 375                 printf(" %d", wsr88d_ray->ref->data[i]);
 376             }
 377             printf("\n");
 378         }
 379     }
 380
 381     free(radial);
 382     return 1;
 383 }
 384
 385
 386 void wsr88d_load_ray_data(Data_moment *data_block, Ray *ray)
 387 {
 388     Data_moment_hdr data_hdr;
 389     int ngates;
 390     int i;
 391     float value, scale, offset;
 392     unsigned char *data;
 393     Range (*invf)(float x);
 394     float (*f)(Range x);
 395
 396     data_hdr = data_block->data_hdr;
 397     data = data_block->data;
 398
 399     /*
 400     printf("wsr88d_load_ray_data: ");
 401     printf("  DataName: %s\n", data_hdr.dataname);
 402     */
 403
 404     ngates = data_hdr.ngates;
 405     /*
 406     printf("  ngates = %d\n", ngates);
 407     printf("  scale = %f\n", data_hdr.scale);
 408     printf("  offset = %f\n", data_hdr.offset);
 409     */
 410     offset = data_hdr.offset;
 411     scale = data_hdr.scale;
 412
 413     /* Note: data range is 2-255. 0 means signal is below threshold, and 1
 414      * means range folded.
 415      */
 416
 417     /* Test print
 418     for (i=0; i < 10; i++) {
 419     value = (data[i] - offset) / scale;
 420     printf("  bytevalue = %d, value = %f\n", data[i], value);
 421     }
 422     */
 423
 424     /* Reflectivity */
 425     if (strncmp(data_hdr.dataname, "DREF", 4) == 0) {
 426         /* convert data to float
 427          * convert float to range and put in ray.range
 428          */
 429         f = DZ_F;
 430         invf = DZ_INVF;
 431         for (i = 0; i < ngates; i++) {
 432             if (data[i] > 1)
 433                 value = (data[i] - offset) / scale;
 434             else value = (data[i] == 0) ? BADVAL : RFVAL;
 435             ray->range[i] = invf(value);
 436             ray->h.f = f;
 437             ray->h.invf = invf;
 438         }
 439     }
 440
 441     /* Velocity */
 442     if (strncmp(data_hdr.dataname, "DVEL", 4) == 0) {
 443         /* convert data to float
 444          * convert float to range and put in ray.range
 445          */
 446         f = VR_F;
 447         invf = VR_INVF;
 448         for (i = 0; i < ngates; i++) {
 449             if (data[i] > 1)
 450                 value = (data[i] - offset) / scale;
 451             else value = (data[i] == 0) ? BADVAL : RFVAL;
 452             ray->range[i] = invf(value);
 453             ray->h.f = f;
 454             ray->h.invf = invf;
 455         }
 456     }
 457
 458     /* Spectrum Width */
 459     if (strncmp(data_hdr.dataname, "DSW", 3) == 0) {
 460         /* convert data to float
 461          * convert float to range and put in ray.range
 462          */
 463         f = SW_F;
 464         invf = SW_INVF;
 465         for (i = 0; i < ngates; i++) {
 466             if (data[i] > 1)
 467                 value = (data[i] - offset) / scale;
 468             else value = (data[i] == 0) ? BADVAL : RFVAL;
 469             ray->range[i] = invf(value);
 470             ray->h.f = f;
 471             ray->h.invf = invf;
 472         }
 473     }
 474     ray->h.range_bin1 = data_hdr.range_first_gate;
 475     ray->h.gate_size = data_hdr.range_samp_interval;
 476     ray->h.nbins = ngates;
 477 }
 478
 479
 480
 481 void wsr88d_load_ray_hdr(Wsr88d_ray_m31 wsr88d_ray, Ray *ray)
 482 {
 483     int month, day, year, hour, minute, sec;
 484     float fsec;
 485     Wsr88d_ray m1_ray;
 486     Ray_header_m31 ray_hdr;
 487
 488     ray_hdr = wsr88d_ray.ray_hdr;
 489     m1_ray.ray_date = ray_hdr.ray_date;
 490     m1_ray.ray_time = ray_hdr.ray_time;
 491
 492     wsr88d_get_date(&m1_ray, &month, &day, &year);
 493     wsr88d_get_time(&m1_ray, &hour, &minute, &sec, &fsec);
 494     ray->h.year = year + 1900;
 495     ray->h.month = month;
 496     ray->h.day = day;
 497     ray->h.hour = hour;
 498     ray->h.minute = minute;
 499     ray->h.sec = sec + fsec;
 500     ray->h.azimuth = ray_hdr.azm;
 501     ray->h.ray_num = ray_hdr.azm_num;
 502     ray->h.elev = ray_hdr.elev;
 503     ray->h.elev_num = ray_hdr.elev_num;
 504     ray->h.unam_rng = wsr88d_ray.unamb_rng / 10.;
 505     ray->h.nyq_vel = wsr88d_ray.nyq_vel / 100.;
 506     int elev_index;
 507     elev_index = ray_hdr.elev_num - 1;
 508     ray->h.azim_rate = vcp_data.azim_rate[elev_index];
 509     ray->h.fix_angle = vcp_data.fixed_angle[elev_index];
 510     ray->h.vel_res = vcp_data.vel_res;
 511
 512     /* Get some values using message type 1 routines.
 513      * First load VCP and elevation numbers into msg 1 ray.
 514      */
 515     m1_ray.vol_cpat = vcp_data.vcp;
 516     m1_ray.elev_num = ray_hdr.elev_num;
 517     m1_ray.unam_rng = wsr88d_ray.unamb_rng;
 518     if (ray_hdr.azm_res != 1)
 519         ray->h.beam_width = 1.0;
 520     else ray->h.beam_width = 0.5;
 521     ray->h.frequency = wsr88d_get_frequency(&m1_ray);
 522     ray->h.pulse_width = wsr88d_get_pulse_width(&m1_ray);
 523     ray->h.pulse_count = wsr88d_get_pulse_count(&m1_ray);
 524     ray->h.prf = wsr88d_get_prf(&m1_ray);
 525     ray->h.wavelength = 0.1071;
 526 }
 527
 528
 529 int wsr88d_get_vol_index(char* dataname)
 530 {
 531     int vol_index;
 532
 533     if (strncmp(dataname, "DREF", 4) == 0) vol_index = DZ_INDEX;
 534     if (strncmp(dataname, "DVEL", 4) == 0) vol_index = VR_INDEX;
 535     if (strncmp(dataname, "DSW", 3) == 0) vol_index = SW_INDEX;
 536     /* TODO: Add the other data moments. */
 537
 538     return vol_index;
 539 }
 540
 541
 542 void wsr88d_load_ray_into_radar(Wsr88d_ray_m31 wsr88d_ray, int isweep, int iray,
 543         Radar *radar)
 544 {
 545     Ray *ray;
 546     Range (*invf)(float x);
 547     float (*f)(Range x);
 548     int vol_index, waveform;
 549
 550     enum waveforms {surveillance=1, doppler_ambres, doppler_no_ambres, batch};
 551
 552 /*  for each data moment do
 553         get the name of the data moment
 554         get new volume if needed
 555         get new sweep if needed
 556         get new ray
 557         put this data moment in ray
 558     endfor each data moment
 559 */
 560
 561     /* Note: The data moment type can only be determined by the name within the
 562      * individual data moment block.  The name of the pointer is not reliable.
 563      * For example, in legacy resolution mode, dbptr_ref points to the velocity
 564      * data moment in velocity split cuts.  Apparently the location of the first
 565      * data moment is always stored in the reflectivity pointer (dbptr_ref), and
 566      * sometimes this is velocity.  When this occurs, the velocity data pointer
 567      * (dbptr_vel) then points to spectrum width.
 568      * With super resolution, there actually is reflectivity in the velocity
 569      * split cut.  It's there for use with the recombination algorithm.
 570      */
 571
 572     if (wsr88d_ray.ray_hdr.dbptr_ref > 0) {
 573         vol_index = wsr88d_get_vol_index(wsr88d_ray.ref->data_hdr.dataname);
 574         switch (vol_index) {
 575             case DZ_INDEX: f = DZ_F; invf = DZ_INVF; break;
 576             case VR_INDEX: f = VR_F; invf = VR_INVF; break;
 577             case SW_INDEX: f = SW_F; invf = SW_INVF; break;
 578             default: f = DZ_F; invf = DZ_INVF; break;
 579         }
 580         /* If this is reflectivity, check the waveform type to make sure
 581          * it isn't from a Doppler split cut.
 582          * We only keep reflectivity if the waveform type is surveillance or
 583          * batch, or the elevation is above the split cut elevations.
 584          */
 585         waveform = vcp_data.waveform[isweep];
 586         if (vol_index != DZ_INDEX ||
 587                 (waveform == surveillance || waveform == batch ||
 588                  vcp_data.fixed_angle[isweep] >= 6.0))
 589         {
 590             if (radar->v[vol_index] == NULL) {
 591                 radar->v[vol_index] = RSL_new_volume(MAXSWEEPS);
 592                 radar->v[vol_index]->h.f = f;
 593                 radar->v[vol_index]->h.invf = invf;
 594             }
 595             if (radar->v[vol_index]->sweep[isweep] == NULL) {
 596                 radar->v[vol_index]->sweep[isweep] = RSL_new_sweep(MAXRAYS_M31);
 597                 radar->v[vol_index]->sweep[isweep]->h.f = f;
 598                 radar->v[vol_index]->sweep[isweep]->h.invf = invf;
 599             }
 600             ray = RSL_new_ray(wsr88d_ray.ref->data_hdr.ngates);
 601             wsr88d_load_ray_data(wsr88d_ray.ref, ray);
 602             wsr88d_load_ray_hdr(wsr88d_ray, ray);
 603             radar->v[vol_index]->sweep[isweep]->ray[iray] = ray;
 604             radar->v[vol_index]->sweep[isweep]->h.nrays = iray+1;
 605         }
 606     }
 607
 608     if (wsr88d_ray.ray_hdr.dbptr_vel > 0) {
 609         vol_index = wsr88d_get_vol_index(wsr88d_ray.vel->data_hdr.dataname);
 610         switch (vol_index) {
 611             case DZ_INDEX: f = DZ_F; invf = DZ_INVF; break;
 612             case VR_INDEX: f = VR_F; invf = VR_INVF; break;
 613             case SW_INDEX: f = SW_F; invf = SW_INVF; break;
 614             default: f = DZ_F; invf = DZ_INVF; break;
 615         }
 616         if (radar->v[vol_index] == NULL) {
 617             radar->v[vol_index] = RSL_new_volume(MAXSWEEPS);
 618             radar->v[vol_index]->h.f = f;
 619             radar->v[vol_index]->h.invf = invf;
 620         }
 621         if (radar->v[vol_index]->sweep[isweep] == NULL) {
 622             radar->v[vol_index]->sweep[isweep] = RSL_new_sweep(MAXRAYS_M31);
 623             radar->v[vol_index]->sweep[isweep]->h.f = f;
 624             radar->v[vol_index]->sweep[isweep]->h.invf = invf;
 625         }
 626         ray = RSL_new_ray(wsr88d_ray.vel->data_hdr.ngates);
 627         wsr88d_load_ray_data(wsr88d_ray.vel, ray);
 628         wsr88d_load_ray_hdr(wsr88d_ray, ray);
 629         radar->v[vol_index]->sweep[isweep]->ray[iray] = ray;
 630         radar->v[vol_index]->sweep[isweep]->h.nrays = iray+1;
 631     }
 632
 633     if (wsr88d_ray.ray_hdr.dbptr_sw > 0) {
 634         vol_index = wsr88d_get_vol_index(wsr88d_ray.sw->data_hdr.dataname);
 635         switch (vol_index) {
 636             case DZ_INDEX: f = DZ_F; invf = DZ_INVF; break;
 637             case VR_INDEX: f = VR_F; invf = VR_INVF; break;
 638             case SW_INDEX: f = SW_F; invf = SW_INVF; break;
 639             default: f = DZ_F; invf = DZ_INVF; break;
 640         }
 641         if (radar->v[vol_index] == NULL) {
 642             radar->v[vol_index] = RSL_new_volume(MAXSWEEPS);
 643             radar->v[vol_index]->h.f = f;
 644             radar->v[vol_index]->h.invf = invf;
 645         }
 646         if (radar->v[vol_index]->sweep[isweep] == NULL) {
 647             radar->v[vol_index]->sweep[isweep] = RSL_new_sweep(MAXRAYS_M31);
 648             radar->v[vol_index]->sweep[isweep]->h.f = f;
 649             radar->v[vol_index]->sweep[isweep]->h.invf = invf;
 650         }
 651         ray = RSL_new_ray(wsr88d_ray.sw->data_hdr.ngates);
 652         wsr88d_load_ray_data(wsr88d_ray.sw, ray);
 653         wsr88d_load_ray_hdr(wsr88d_ray, ray);
 654         radar->v[vol_index]->sweep[isweep]->ray[iray] = ray;
 655         radar->v[vol_index]->sweep[isweep]->h.nrays = iray+1;
 656     }
 657
 658 }
 659
 660
 661
 662 void wsr88d_load_sweep_header(Radar *radar, int isweep,
 663         Wsr88d_ray_m31 wsr88d_ray)
 664 {
 665     int ivolume;
 666     Ray_header_m31 ray_hdr;
 667     Sweep *sweep;
 668     int vcp;
 669
 670     ray_hdr = wsr88d_ray.ray_hdr;
 671
 672     for (ivolume=0; ivolume < MAX_RADAR_VOLUMES; ivolume++) {
 673         if (radar->v[ivolume] != NULL && radar->v[ivolume]->sweep[isweep] != NULL) {
 674             sweep = radar->v[ivolume]->sweep[isweep];
 675             radar->v[ivolume]->sweep[isweep]->h.sweep_num = ray_hdr.elev_num;
 676             radar->v[ivolume]->sweep[isweep]->h.elev =
 677                 vcp_data.fixed_angle[isweep];
 678             if (ray_hdr.azm_res != 1)
 679                 sweep->h.beam_width = 1.0;
 680             else sweep->h.beam_width = 0.5;
 681             sweep->h.vert_half_bw = sweep->h.beam_width / 2.;
 682             sweep->h.horz_half_bw = sweep->h.beam_width / 2.;
 683         }
 684     }
 685 }
 686
 687
 688 Radar *load_wsr88d_m31_into_radar(Wsr88d_file *wf)
 689 {
 690     Wsr88d_msg_hdr msghdr;
 691     Wsr88d_ray_m31 wsr88d_ray;
 692     short non31_seg_remainder[1202]; /* Remainder after message header */
 693     int end_of_vos = 0, isweep = 0, iray = 0;
 694     int msg_hdr_size, msg_size, n;
 695     Radar *radar = NULL;
 696
 697 /* Message type 31 is a variable length message, whereas all other message
 698  * types are made up of 2432 byte segments.  To handle this, we read the
 699  * message header and check the message type.  If it is not 31, then we read
 700  * the remainder of the constant size segment.  If message type is 31, we read
 701  * the remainder of the message by the size given in message header.
 702  * When reading the message header, we must include 12 bytes inserted
 703  * by RPG, which we ignore, followed by the 8 halfwords (short integers) which
 704  * make up the actual message header.
 705  * For more information, see the "Interface Control Document for the RDA/RPG"
 706  * at the WSR-88D Radar Operations Center web site.
 707  */
 708
 709     n = fread(&msghdr, sizeof(Wsr88d_msg_hdr), 1, wf->fptr);
 710
 711     /* printf("msgtype = %d\n", msghdr.msg_type); */
 712     msg_hdr_size = sizeof(Wsr88d_msg_hdr) - sizeof(msghdr.rpg);
 713
 714
 715     radar = RSL_new_radar(MAX_RADAR_VOLUMES);
 716
 717     while (! end_of_vos) {
 718         if (msghdr.msg_type == 31) {
 719             if (little_endian()) wsr88d_swap_m31_hdr(&msghdr);
 720
 721             /* Get size in bytes of message following message header.
 722              * The size given in message header is in halfwords, so double it.
 723              */
 724             msg_size = (int) msghdr.msg_size * 2 - msg_hdr_size;
 725
 726             n = read_wsr88d_ray_m31(wf, msg_size, &wsr88d_ray);
 727             if (n <= 0) return NULL;
 728
 729             /* Load this ray into radar structure ray. */
 730             wsr88d_load_ray_into_radar(wsr88d_ray, isweep, iray, radar);
 731             iray++;
 732         }
 733         else { /* msg_type not 31 */
 734             n = fread(&non31_seg_remainder, sizeof(non31_seg_remainder), 1,
 735                     wf->fptr);
 736             if (n < 1) {
 737                 fprintf(stderr,"Warning: load_wsr88d_m31_into_radar: ");
 738                 if (feof(wf->fptr) != 0)
 739                     fprintf(stderr, "Unexpected end of file.\n");
 740                 else fprintf(stderr,"Read failed.\n");
 741                 fprintf(stderr,"Current sweep number: %d\n"
 742                         "Last ray read: %d\n", isweep+1, iray);
 743                 wsr88d_load_sweep_header(radar, isweep, wsr88d_ray);
 744                 return radar;
 745             }
 746             if (msghdr.msg_type == 5) {
 747                 wsr88d_get_vcp_data(non31_seg_remainder);
 748                 radar->h.vcp = vcp_data.vcp;
 749                 /* printf("VCP = %d\n", vcp_data.vcp); */
 750             }
 751         }
 752
 753         /* Check for end of sweep */
 754         if (wsr88d_ray.ray_hdr.radial_status == END_OF_ELEV) {
 755             wsr88d_load_sweep_header(radar, isweep, wsr88d_ray);
 756             isweep++;
 757             iray = 0;
 758         }
 759
 760         if (wsr88d_ray.ray_hdr.radial_status != END_VOS) {
 761             n = fread(&msghdr, sizeof(Wsr88d_msg_hdr), 1, wf->fptr);
 762             if (n < 1) {
 763                 fprintf(stderr,"Warning: load_wsr88d_m31_into_radar: ");
 764                 if (feof(wf->fptr) != 0) fprintf(stderr,
 765                         "Unexpected end of file.\n");
 766                 else fprintf(stderr,"Failed reading msghdr.\n");
 767                 fprintf(stderr,"Current sweep number: %d\n"
 768                         "Last ray read: %d\n", isweep+1, iray);
 769                 wsr88d_load_sweep_header(radar, isweep, wsr88d_ray);
 770                 return radar;
 771             }
 772         }
 773         else {
 774             end_of_vos = 1;
 775             wsr88d_load_sweep_header(radar, isweep, wsr88d_ray);
 776         }
 777         if (feof(wf->fptr) != 0) end_of_vos = 1;
 778     }
 779
 780     return radar;
 781 }