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 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
  21     Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
  22     Boston, MA  02110-1301, USA.
  23 */
  24
  25
  26 /*
  27  * This file contains routines for processing Message Type 31, the digital
  28  * radar message type introduced in WSR-88D Level II Build 10.  For more
  29  * information, see the "Interface Control Document for the RDA/RPG" at the
  30  * WSR-88D Radar Operations Center web site.
  31  */
  32
  33 #include "rsl.h"
  34 #include "wsr88d.h"
  35 #include <string.h>
  36
  37 /* Data descriptions in the following data structures are from the "Interface
  38  * Control Document for the RDA/RPG", Build 10.0 Draft, WSR-88D Radar
  39  * Operations Center.
  40  */
  41
  42 typedef struct {
  43     short rpg[6]; /* 12 bytes inserted by RPG Communications Mgr. Ignored. */
  44     unsigned short msg_size;  /* Message size for this segment, in halfwords */
  45     unsigned char  channel;   /* RDA Redundant Channel */
  46     unsigned char  msg_type;  /* Message type. For example, 31 */
  47     unsigned short id_seq;    /* Msg seq num = 0 to 7FFF, then roll over to 0 */
  48     unsigned short msg_date;  /* Modified Julian date from 1/1/70 */
  49     unsigned int   msg_time;  /* Packet generation time in ms past midnight */
  50     unsigned short num_segs;  /* Number of segments for this message */
  51     unsigned short seg_num;   /* Number of this segment */
  52 } Wsr88d_msg_hdr;
  53
  54 typedef struct {
  55     char radar_id[4];
  56     unsigned int   ray_time; /* Data collection time in milliseconds past midnight GMT */
  57     unsigned short ray_date; /* Julian date - 2440586.5 (1/01/1970) */
  58     unsigned short azm_num ; /* Radial number within elevation scan */
  59     float azm;      /* Azimuth angle in degrees (0 to 359.956055) */
  60     unsigned char compression_code; /* 0 = uncompressed, 1 = BZIP2, 2 = zlib */
  61     unsigned char spare; /* for word alignment */
  62     unsigned short radial_len; /* radial length in bytes, including data header block */
  63     unsigned char azm_res;
  64     unsigned char radial_status;
  65     unsigned char elev_num;
  66     unsigned char sector_cut_num;
  67     float elev;  /* Elevation angle in degrees (-7.0 to 70.0) */
  68     unsigned char radial_spot_blanking;
  69     unsigned char azm_indexing_mode;
  70     unsigned short data_block_count;
  71     /* Data Block Indexes */
  72     unsigned int vol_const;
  73     unsigned int elev_const;
  74     unsigned int radial_const;
  75     unsigned int field1;
  76     unsigned int field2;
  77     unsigned int field3;
  78     unsigned int field4;
  79     unsigned int field5;
  80     unsigned int field6;
  81 } Ray_header_m31;  /* Called Data Header Block in RDA/RPG document. */
  82
  83 typedef struct {
  84     char dataname[4];
  85     unsigned int  reserved;
  86     unsigned short ngates;
  87     short range_first_gate;
  88     short range_samp_interval;
  89     short thresh_not_overlayed;
  90     short snr_thresh;
  91     unsigned char controlflag;
  92     unsigned char datasize_bits;
  93     float scale;
  94     float offset;
  95 } Data_moment_hdr;
  96
  97 #define MAX_RADIAL_LENGTH 14288
  98
  99 typedef struct {
 100     Ray_header_m31 ray_hdr;
 101     float unamb_rng;
 102     float nyq_vel;
 103     unsigned char data[MAX_RADIAL_LENGTH];
 104 } Wsr88d_ray_m31;
 105
 106
 107
 108 void wsr88d_swap_m31_hdr(Wsr88d_msg_hdr *msghdr)
 109 {
 110     swap_2_bytes(&msghdr->msg_size);
 111     swap_2_bytes(&msghdr->id_seq);
 112     swap_2_bytes(&msghdr->msg_date);
 113     swap_4_bytes(&msghdr->msg_time);
 114     swap_2_bytes(&msghdr->num_segs);
 115     swap_2_bytes(&msghdr->seg_num);
 116 }
 117
 118
 119 void wsr88d_swap_m31_ray_hdr(Ray_header_m31 *ray_hdr)
 120 {
 121     int *data_ptr;
 122
 123     swap_4_bytes(&ray_hdr->ray_time);
 124     swap_2_bytes(&ray_hdr->ray_date);
 125     swap_2_bytes(&ray_hdr->azm_num);
 126     swap_4_bytes(&ray_hdr->azm);
 127     swap_2_bytes(&ray_hdr->radial_len);
 128     swap_4_bytes(&ray_hdr->elev);
 129     swap_2_bytes(&ray_hdr->data_block_count);
 130     data_ptr = (int *) &ray_hdr->vol_const;
 131     for (; data_ptr <= (int *) &ray_hdr->field6; data_ptr++)
 132         swap_4_bytes(data_ptr);
 133 }
 134
 135
 136 void wsr88d_swap_data_hdr(Data_moment_hdr *this_field)
 137 {
 138     short *halfword;
 139     halfword = (short *) &this_field->ngates;
 140     for (; halfword < (short *) &this_field->controlflag; halfword++)
 141         swap_2_bytes(halfword);
 142     swap_4_bytes(&this_field->scale);
 143     swap_4_bytes(&this_field->offset);
 144 }
 145
 146
 147 float wsr88d_get_angle(short bitfield)
 148 {
 149     short mask = 1;
 150     int i;
 151     float angle = 0.;
 152     float value[13] = {0.043945, 0.08789, 0.17578, 0.35156, .70313, 1.40625,
 153         2.8125, 5.625, 11.25, 22.5, 45., 90., 180.};
 154
 155     /* Find which bits are set and sum corresponding values to get angle. */
 156
 157     bitfield = bitfield >> 3;  /* 3 least significant bits aren't used. */
 158     for (i = 0; i < 13; i++) {
 159         if (bitfield & mask) angle += value[i];
 160         bitfield = bitfield >> 1;
 161     }
 162     return angle;
 163 }
 164
 165
 166 float wsr88d_get_azim_rate(short bitfield)
 167 {
 168     short mask = 1;
 169     int i;
 170     float rate = 0.;
 171     float value[12] = {0.0109863, 0.021972656, 0.043945, 0.08789, 0.17578,
 172         0.35156, .70313, 1.40625, 2.8125, 5.625, 11.25, 22.5};
 173
 174     /* Find which bits are set and sum corresponding values to get rate. */
 175
 176     bitfield = bitfield >> 3;  /* 3 least significant bits aren't used. */
 177     for (i = 0; i < 12; i++) {
 178         if (bitfield & mask) rate += value[i];
 179         bitfield = bitfield >> 1;
 180     }
 181     if (bitfield >> 15) rate = -rate;
 182     return rate;
 183 }
 184
 185 #define WSR88D_MAX_SWEEPS 20
 186
 187 typedef struct {
 188     int vcp;
 189     int num_cuts;
 190     float vel_res;
 191     float fixed_angle[WSR88D_MAX_SWEEPS];
 192     float azim_rate[WSR88D_MAX_SWEEPS];
 193     int waveform[WSR88D_MAX_SWEEPS];
 194     int super_res_ctrl[WSR88D_MAX_SWEEPS];
 195     int surveil_prf_num[WSR88D_MAX_SWEEPS];
 196     int doppler_prf_num[WSR88D_MAX_SWEEPS];
 197 } VCP_data;
 198
 199 static VCP_data vcp_data;
 200
 201 void wsr88d_get_vcp_data(short *msgtype5)
 202 {
 203     short azim_rate, fixed_angle, vel_res;
 204     short sres_and_survprf; /* super res ctrl and surveil prf, one byte each */
 205     short chconf_and_waveform;
 206     int i;
 207
 208     vcp_data.vcp = (unsigned short) msgtype5[2];
 209     vcp_data.num_cuts = msgtype5[3];
 210     if (little_endian()) {
 211         swap_2_bytes(&vcp_data.vcp);
 212         swap_2_bytes(&vcp_data.num_cuts);
 213     }
 214     vel_res = msgtype5[5];
 215     if (little_endian()) swap_2_bytes(&vel_res);
 216     vel_res = vel_res >> 8;
 217     if (vel_res == 2) vcp_data.vel_res = 0.5;
 218     else if (vel_res == 4) vcp_data.vel_res = 1.0;
 219     else vcp_data.vel_res = 0.0;
 220     /* Get elevation related information for each sweep. */
 221     for (i=0; i < vcp_data.num_cuts; i++) {
 222         fixed_angle = msgtype5[11 + i*23];
 223         azim_rate = msgtype5[15 + i*23];
 224         chconf_and_waveform = msgtype5[12 + i*23];
 225         sres_and_survprf = msgtype5[13 + i*23];
 226         vcp_data.doppler_prf_num[i] = msgtype5[23 + i*23];
 227         if (little_endian()) {
 228             swap_2_bytes(&fixed_angle);
 229             swap_2_bytes(&azim_rate);
 230             swap_2_bytes(&chconf_and_waveform);
 231             swap_2_bytes(&sres_and_survprf);
 232             swap_2_bytes(&vcp_data.doppler_prf_num[i]);
 233         }
 234         vcp_data.fixed_angle[i] = wsr88d_get_angle(fixed_angle);
 235         vcp_data.azim_rate[i] = wsr88d_get_azim_rate(azim_rate);
 236         vcp_data.waveform[i] = chconf_and_waveform & 0xff;
 237         vcp_data.super_res_ctrl[i] = sres_and_survprf >> 8;
 238         vcp_data.surveil_prf_num[i] = sres_and_survprf & 0xff;
 239     }
 240 }
 241
 242
 243 void get_wsr88d_unamb_and_nyq_vel(Wsr88d_ray_m31 *wsr88d_ray, float *unamb_rng,
 244         float *nyq_vel)
 245 {
 246     int dindex, found;
 247     short nyq_vel_sh, unamb_rng_sh;
 248
 249     found = 0;
 250     dindex = wsr88d_ray->ray_hdr.radial_const;
 251     if (strncmp((char *) &wsr88d_ray->data[dindex], "RRAD", 4) == 0) found = 1;
 252     else {
 253         dindex = wsr88d_ray->ray_hdr.elev_const;
 254         if (strncmp((char *) &wsr88d_ray->data[dindex], "RRAD", 4) == 0)
 255             found = 1;
 256         else {
 257             dindex = wsr88d_ray->ray_hdr.vol_const;
 258             if (strncmp((char *) &wsr88d_ray->data[dindex], "RRAD", 4) == 0)
 259                 found = 1;
 260         }
 261     }
 262     if (found) {
 263         memcpy(&unamb_rng_sh, &wsr88d_ray->data[dindex+6], 2);
 264         memcpy(&nyq_vel_sh, &wsr88d_ray->data[dindex+16], 2);
 265         if (little_endian()) {
 266             swap_2_bytes(&unamb_rng_sh);
 267             swap_2_bytes(&nyq_vel_sh);
 268         }
 269         *unamb_rng = unamb_rng_sh / 10.;
 270         *nyq_vel = nyq_vel_sh / 100.;
 271     } else {
 272         *unamb_rng = 0.;
 273         *nyq_vel = 0.;
 274     }
 275 }
 276
 277
 278 int read_wsr88d_ray_m31(Wsr88d_file *wf, int msg_size,
 279         Wsr88d_ray_m31 *wsr88d_ray)
 280 {
 281     int n;
 282     float nyq_vel, unamb_rng;
 283
 284     /* Read wsr88d ray. */
 285
 286     n = fread(wsr88d_ray->data, msg_size, 1, wf->fptr);
 287     if (n < 1) {
 288         fprintf(stderr,"read_wsr88d_ray_m31: Read failed.\n");
 289         return 0;
 290     }
 291
 292     /* Copy data header block to ray header structure. */
 293     memcpy(&wsr88d_ray->ray_hdr, &wsr88d_ray->data, sizeof(Ray_header_m31));
 294
 295     if (little_endian()) wsr88d_swap_m31_ray_hdr(&wsr88d_ray->ray_hdr);
 296
 297     /* Retrieve unambiguous range and Nyquist velocity here so that we don't
 298      * have to do it for each data moment later.
 299      */
 300     get_wsr88d_unamb_and_nyq_vel(wsr88d_ray, &unamb_rng, &nyq_vel);
 301     wsr88d_ray->unamb_rng = unamb_rng;
 302     wsr88d_ray->nyq_vel = nyq_vel;
 303
 304     return 1;
 305 }
 306
 307
 308 void wsr88d_load_ray_hdr(Wsr88d_ray_m31 *wsr88d_ray, Ray *ray)
 309 {
 310     int month, day, year, hour, minute, sec;
 311     float fsec;
 312     Wsr88d_ray m1_ray;
 313     Ray_header_m31 ray_hdr;
 314
 315     ray_hdr = wsr88d_ray->ray_hdr;
 316     m1_ray.ray_date = ray_hdr.ray_date;
 317     m1_ray.ray_time = ray_hdr.ray_time;
 318
 319     wsr88d_get_date(&m1_ray, &month, &day, &year);
 320     wsr88d_get_time(&m1_ray, &hour, &minute, &sec, &fsec);
 321     ray->h.year = year + 1900;
 322     ray->h.month = month;
 323     ray->h.day = day;
 324     ray->h.hour = hour;
 325     ray->h.minute = minute;
 326     ray->h.sec = sec + fsec;
 327     ray->h.azimuth = ray_hdr.azm;
 328     ray->h.ray_num = ray_hdr.azm_num;
 329     ray->h.elev = ray_hdr.elev;
 330     ray->h.elev_num = ray_hdr.elev_num;
 331     ray->h.unam_rng = wsr88d_ray->unamb_rng;
 332     ray->h.nyq_vel = wsr88d_ray->nyq_vel;
 333     int elev_index;
 334     elev_index = ray_hdr.elev_num - 1;
 335     ray->h.azim_rate = vcp_data.azim_rate[elev_index];
 336     ray->h.fix_angle = vcp_data.fixed_angle[elev_index];
 337     ray->h.vel_res = vcp_data.vel_res;
 338     if (ray_hdr.azm_res != 1)
 339         ray->h.beam_width = 1.0;
 340     else ray->h.beam_width = 0.5;
 341
 342     /* For convenience, use message type 1 routines to get some values.
 343      * First load VCP and elevation numbers into a msg 1 ray.
 344      */
 345     m1_ray.vol_cpat = vcp_data.vcp;
 346     m1_ray.elev_num = ray_hdr.elev_num;
 347     m1_ray.unam_rng = (short) (wsr88d_ray->unamb_rng * 10.);
 348     /* Get values from message type 1 routines. */
 349     ray->h.frequency = wsr88d_get_frequency(&m1_ray);
 350     ray->h.pulse_width = wsr88d_get_pulse_width(&m1_ray);
 351     ray->h.pulse_count = wsr88d_get_pulse_count(&m1_ray);
 352     ray->h.prf = (int) wsr88d_get_prf(&m1_ray);
 353     ray->h.wavelength = 0.1071;
 354 }
 355
 356
 357 int wsr88d_get_vol_index(char* dataname)
 358 {
 359     if (strncmp(dataname, "DREF", 4) == 0) return DZ_INDEX;
 360     if (strncmp(dataname, "DVEL", 4) == 0) return VR_INDEX;
 361     if (strncmp(dataname, "DSW",  3) == 0) return SW_INDEX;
 362     if (strncmp(dataname, "DZDR", 4) == 0) return DR_INDEX;
 363     if (strncmp(dataname, "DPHI", 4) == 0) return PH_INDEX;
 364     if (strncmp(dataname, "DRHO", 4) == 0) return RH_INDEX;
 365
 366     return -1;
 367 }
 368
 369 #define MAXRAYS_M31 800
 370 #define MAXSWEEPS 20
 371
 372 void wsr88d_load_ray_into_radar(Wsr88d_ray_m31 *wsr88d_ray, int isweep,
 373         Radar *radar)
 374 {
 375     /* Load data into ray structure for each data field. */
 376
 377     int data_index;
 378     int *field_offset;
 379     int ifield, nfields;
 380     int iray;
 381     const nconstblocks = 3;
 382
 383     Data_moment_hdr data_hdr;
 384     int ngates, do_swap;
 385     int i, hdr_size;
 386     unsigned short item;
 387     float value, scale, offset;
 388     unsigned char *data;
 389     Range (*invf)(float x);
 390     float (*f)(Range x);
 391     Ray *ray;
 392     int vol_index, waveform;
 393     char *type_str;
 394
 395     extern int rsl_qfield[]; /* See RSL_select_fields in volume.c */
 396
 397     enum waveforms {surveillance=1, doppler_w_amb_res, doppler_no_amb_res,
 398         batch};
 399
 400     int merging_split_cuts;
 401
 402     merging_split_cuts =  wsr88d_merge_split_cuts_is_set();
 403     nfields = wsr88d_ray->ray_hdr.data_block_count - nconstblocks;
 404     field_offset = (int *) &wsr88d_ray->ray_hdr.radial_const;
 405     do_swap = little_endian();
 406     iray = wsr88d_ray->ray_hdr.azm_num - 1;
 407
 408     for (ifield=0; ifield < nfields; ifield++) {
 409         field_offset++;
 410         data_index = *field_offset;
 411         /* Get data moment header. */
 412         hdr_size = sizeof(data_hdr);
 413         memcpy(&data_hdr, &wsr88d_ray->data[data_index], hdr_size);
 414         if (do_swap) wsr88d_swap_data_hdr(&data_hdr);
 415         data_index += hdr_size;
 416
 417         vol_index = wsr88d_get_vol_index(data_hdr.dataname);
 418         if (vol_index < 0) {
 419             fprintf(stderr,"wsr88d_load_ray_into_radar: Unknown dataname %s.  "
 420                     "isweep = %d, iray = %d.\n", data_hdr.dataname, isweep,
 421                     iray);
 422             return;
 423         }
 424
 425         /* Is this field in the selected fields list? */
 426         if (!rsl_qfield[vol_index]) continue;
 427
 428         switch (vol_index) {
 429             case DZ_INDEX: f = DZ_F; invf = DZ_INVF;
 430                  type_str = "Reflectivity"; break;
 431             case VR_INDEX: f = VR_F; invf = VR_INVF;
 432                  type_str = "Velocity"; break;
 433             case SW_INDEX: f = SW_F; invf = SW_INVF;
 434                  type_str = "Spectrum width"; break;
 435             case DR_INDEX: f = DR_F; invf = DR_INVF;
 436                  type_str = "Diff. Reflectivity"; break;
 437             case PH_INDEX: f = PH_F; invf = PH_INVF;
 438                  type_str = "Diff. Phase"; break;
 439             case RH_INDEX: f = RH_F; invf = RH_INVF;
 440                  type_str = "Correlation Coef (Rho)"; break;
 441         }
 442
 443         waveform = vcp_data.waveform[isweep];
 444
 445         /* If this field is reflectivity, check to see if it's from the velocity
 446          * sweep in a split cut.  If so, we normally skip it since we already
 447          * have reflectivity from the surveillance sweep.  It is kept only when
 448          * the user has turned off merging of split cuts.  We skip over this
 449          * field if all of the following are true: surveillance PRF number is 0,
 450          * waveform is Contiguous Doppler with Ambiguity Resolution (range
 451          * unfolding), and we're merging split cuts.
 452          */
 453         if (vol_index == DZ_INDEX && (vcp_data.surveil_prf_num[isweep] == 0 &&
 454                     vcp_data.waveform[isweep] == doppler_w_amb_res &&
 455                     merging_split_cuts))
 456             continue;
 457
 458         /* Load the data for this field. */
 459         if (radar->v[vol_index] == NULL) {
 460             radar->v[vol_index] = RSL_new_volume(MAXSWEEPS);
 461             radar->v[vol_index]->h.f = f;
 462             radar->v[vol_index]->h.invf = invf;
 463             radar->v[vol_index]->h.type_str = type_str;
 464         }
 465         if (radar->v[vol_index]->sweep[isweep] == NULL) {
 466             radar->v[vol_index]->sweep[isweep] = RSL_new_sweep(MAXRAYS_M31);
 467             radar->v[vol_index]->sweep[isweep]->h.f = f;
 468             radar->v[vol_index]->sweep[isweep]->h.invf = invf;
 469         }
 470         ngates = data_hdr.ngates;
 471         ray = RSL_new_ray(ngates);
 472
 473         /* Convert data to float, then use range function to store in ray.
 474          * Note: data range is 2-255. 0 means signal is below threshold, and 1
 475          * means range folded.
 476          */
 477
 478         offset = data_hdr.offset;
 479         scale = data_hdr.scale;
 480         if (data_hdr.scale == 0) scale = 1.0;
 481         data = &wsr88d_ray->data[data_index];
 482         for (i = 0; i < ngates; i++) {
 483             if (data_hdr.datasize_bits != 16) {
 484                 item = *data;
 485                 data++;
 486             } else {
 487                 item = *(unsigned short *)data;
 488                 if (do_swap) swap_2_bytes(&item);
 489                 data += 2;
 490             }
 491             if (item > 1)
 492                 value = (item - offset) / scale;
 493             else value = (item == 0) ? BADVAL : RFVAL;
 494             ray->range[i] = invf(value);
 495             ray->h.f = f;
 496             ray->h.invf = invf;
 497         }
 498         wsr88d_load_ray_hdr(wsr88d_ray, ray);
 499         ray->h.range_bin1 = data_hdr.range_first_gate;
 500         ray->h.gate_size = data_hdr.range_samp_interval;
 501         ray->h.nbins = ngates;
 502         radar->v[vol_index]->sweep[isweep]->ray[iray] = ray;
 503         radar->v[vol_index]->sweep[isweep]->h.nrays = iray+1;
 504     } /* for each data field */
 505 }
 506
 507
 508 void wsr88d_load_sweep_header(Radar *radar, int isweep)
 509 {
 510     int ivolume, nrays;
 511     Sweep *sweep;
 512     Ray *last_ray;
 513
 514     for (ivolume=0; ivolume < MAX_RADAR_VOLUMES; ivolume++) {
 515         if (radar->v[ivolume] != NULL &&
 516                 radar->v[ivolume]->sweep[isweep] != NULL) {
 517             sweep = radar->v[ivolume]->sweep[isweep];
 518             nrays = sweep->h.nrays;
 519             if (nrays == 0) continue;
 520             last_ray = sweep->ray[nrays-1];
 521             sweep->h.sweep_num = last_ray->h.elev_num;
 522             sweep->h.elev = vcp_data.fixed_angle[isweep];
 523             sweep->h.beam_width = last_ray->h.beam_width;
 524             sweep->h.vert_half_bw = sweep->h.beam_width / 2.;
 525             sweep->h.horz_half_bw = sweep->h.beam_width / 2.;
 526         }
 527     }
 528 }
 529
 530
 531 Radar *wsr88d_load_m31_into_radar(Wsr88d_file *wf)
 532 {
 533     Wsr88d_msg_hdr msghdr;
 534     Wsr88d_ray_m31 wsr88d_ray;
 535     short non31_seg_remainder[1202]; /* Remainder after message header */
 536     int end_of_vos = 0, isweep = 0;
 537     int msg_hdr_size, msg_size, n;
 538     int prev_elev_num = 1, prev_raynum = 0, raynum = 0;
 539     Radar *radar = NULL;
 540     enum radial_status {START_OF_ELEV, INTERMED_RADIAL, END_OF_ELEV, BEGIN_VOS,
 541         END_VOS};
 542
 543
 544     /* Message type 31 is a variable length message.  All other types consist of
 545      * 1 or more segments of length 2432 bytes.  To handle all types, we read
 546      * the message header and check the type.  If not 31, then simply read
 547      * the remainder of the 2432-byte segment.  If it is 31, use the size given
 548      * in message header to determine how many bytes to read.
 549      */
 550
 551     n = fread(&msghdr, sizeof(Wsr88d_msg_hdr), 1, wf->fptr);
 552
 553     /* printf("msgtype = %d\n", msghdr.msg_type); */
 554     msg_hdr_size = sizeof(Wsr88d_msg_hdr) - sizeof(msghdr.rpg);
 555
 556     radar = RSL_new_radar(MAX_RADAR_VOLUMES);
 557     memset(&wsr88d_ray, 0, sizeof(Wsr88d_ray_m31)); /* Initialize to be safe. */
 558
 559     while (! end_of_vos) {
 560         if (msghdr.msg_type == 31) {
 561             if (little_endian()) wsr88d_swap_m31_hdr(&msghdr);
 562
 563             /* Get size of the remainder of message.  The given size is in
 564              * halfwords; convert it to bytes.
 565              */
 566             msg_size = (int) msghdr.msg_size * 2 - msg_hdr_size;
 567
 568             n = read_wsr88d_ray_m31(wf, msg_size, &wsr88d_ray);
 569             if (n <= 0) return NULL;
 570             raynum = wsr88d_ray.ray_hdr.azm_num;
 571             if (raynum > MAXRAYS_M31) {
 572                 fprintf(stderr,"Error: raynum = %d, exceeds MAXRAYS_M31"
 573                         " (%d)\n", raynum, MAXRAYS_M31);
 574                 fprintf(stderr,"isweep = %d\n", isweep);
 575                 RSL_free_radar(radar);
 576                 return NULL;
 577             }
 578
 579             /* Check for an unexpected start of new elevation, and issue a
 580              * warning if this has occurred.  This condition usually means
 581              * less rays then expected in the sweep that just ended.
 582              */
 583             if (wsr88d_ray.ray_hdr.radial_status == START_OF_ELEV &&
 584                     wsr88d_ray.ray_hdr.elev_num-1 > isweep) {
 585                 fprintf(stderr,"Warning: Radial status is Start-of-Elevation, "
 586                         "but End-of-Elevation was not\n"
 587                         "issued for elevation number %d.  Number of rays = %d"
 588                         "\n", prev_elev_num, prev_raynum);
 589                 wsr88d_load_sweep_header(radar, isweep);
 590                 isweep++;
 591                 prev_elev_num = wsr88d_ray.ray_hdr.elev_num - 1;
 592             }
 593
 594             /* Load ray into radar structure. */
 595             wsr88d_load_ray_into_radar(&wsr88d_ray, isweep, radar);
 596             prev_raynum = raynum;
 597
 598             /* Check for end of sweep */
 599             if (wsr88d_ray.ray_hdr.radial_status == END_OF_ELEV) {
 600                 wsr88d_load_sweep_header(radar, isweep);
 601                 isweep++;
 602                 prev_elev_num = wsr88d_ray.ray_hdr.elev_num;
 603             }
 604         }
 605         else { /* msg_type not 31 */
 606             n = fread(&non31_seg_remainder, sizeof(non31_seg_remainder), 1,
 607                     wf->fptr);
 608             if (n < 1) {
 609                 fprintf(stderr,"Warning: load_wsr88d_m31_into_radar: ");
 610                 if (feof(wf->fptr) != 0)
 611                     fprintf(stderr, "Unexpected end of file.\n");
 612                 else
 613                     fprintf(stderr,"Read failed.\n");
 614                 fprintf(stderr,"Current sweep index: %d\n"
 615                         "Last ray read: %d\n", isweep, prev_raynum);
 616                 wsr88d_load_sweep_header(radar, isweep);
 617                 return radar;
 618             }
 619             if (msghdr.msg_type == 5) {
 620                 wsr88d_get_vcp_data(non31_seg_remainder);
 621                 radar->h.vcp = vcp_data.vcp;
 622             }
 623         }
 624
 625         /* If not at end of volume scan, read next message header. */
 626         if (wsr88d_ray.ray_hdr.radial_status != END_VOS) {
 627             n = fread(&msghdr, sizeof(Wsr88d_msg_hdr), 1, wf->fptr);
 628             if (n < 1) {
 629                 fprintf(stderr,"Warning: load_wsr88d_m31_into_radar: ");
 630                 if (feof(wf->fptr) != 0)
 631                     fprintf(stderr,"Unexpected end of file.\n");
 632                 else fprintf(stderr,"Failed reading msghdr.\n");
 633                 fprintf(stderr,"Current sweep index: %d\n"
 634                         "Last ray read: %d\n", isweep, prev_raynum);
 635                 wsr88d_load_sweep_header(radar, isweep);
 636                 end_of_vos = 1;
 637             }
 638         }
 639         else {
 640             end_of_vos = 1;
 641             wsr88d_load_sweep_header(radar, isweep);
 642         }
 643     }  /* while not end of vos */
 644
 645     return radar;
 646 }