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