3 This is the TRMM Office Radar Software Library.
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.
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.
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/>.
25 * This file contains routines for processing Message Type 31, the digital
26 * radar message type introduced in WSR-88D Level II Build 10.
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
38 short rpg[6]; /* Unused. Not really part of message header, but is
39 * inserted by RPG Communications Manager for each message.
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;
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. */
82 unsigned int reserved;
83 unsigned short ngates;
84 short range_first_gate;
85 short range_samp_interval;
86 short thresh_not_overlayed;
88 unsigned char controlflag;
89 unsigned char datasize_bits;
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.
103 Ray_header_m31 ray_hdr;
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.
116 #define MAXRAYS_M31 800
119 enum {START_OF_ELEV, INTERMED_RADIAL, END_OF_ELEV, BEGIN_VOS, END_VOS};
122 void wsr88d_swap_m31_hdr(Wsr88d_msg_hdr *msghdr)
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);
133 void wsr88d_swap_m31_ray(Ray_header_m31 *wsr88d_ray)
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);
150 void wsr88d_swap_data_moment(Data_moment_hdr *this_field)
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);
161 void testprt(Ray_header_m31 ray_hdr)
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);
180 float wsr88d_get_angle(short bitfield)
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.};
188 /* find which bits are set and sum corresponding values to get angle. */
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;
199 float wsr88d_get_azim_rate(short bitfield)
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};
207 /* Find which bits are set and sum corresponding values to get rate. */
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;
214 if (bitfield >> 15) rate = -rate;
218 #define WSR88D_MAX_SWEEPS 20
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];
232 static VCP_data vcp_data;
234 void wsr88d_get_vcp_data(short *msgtype5)
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;
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);
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]);
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;
276 Data_moment *read_data_moment(unsigned char *radial, int begin_data_block)
278 Data_moment *this_field;
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));
288 this_field = malloc(sizeof(Data_moment));
289 hdr_size = sizeof(Data_moment_hdr);
290 /* printf("hdr_size = %d\n", hdr_size); */
292 dptr = radial + begin_data_block;
293 memcpy(&this_field->data_hdr, dptr, hdr_size);
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;
304 int read_wsr88d_ray_m31(Wsr88d_file *wf, int msg_size,
305 Wsr88d_ray_m31 *wsr88d_ray)
307 Ray_header_m31 ray_hdr;
308 int n, begin_data_block, bytes_to_read, ray_hdr_len;
309 unsigned char *radial;
311 /* Read ray data header block */
312 n = fread(&wsr88d_ray->ray_hdr, sizeof(Ray_header_m31), 1, wf->fptr);
314 fprintf(stderr,"read_wsr88d_ray_m31: Read failed.\n");
317 if (little_endian()) wsr88d_swap_m31_ray(&wsr88d_ray->ray_hdr);
319 ray_hdr = wsr88d_ray->ray_hdr;
322 if (ray_hdr.azm_num == 1) testprt(wsr88d_ray->ray_hdr);
325 ray_hdr_len = sizeof(ray_hdr);
328 * Use header_size offset with data pointers
329 * Copy data into structures
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);
336 fprintf(stderr,"read_wsr88d_ray_m31: Read failed.\n");
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);
351 wsr88d_ray->unamb_rng = 0;
352 wsr88d_ray->nyq_vel = 0;
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);
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);
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);
369 /* For testing: print reflectivity data. */
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]);
386 void wsr88d_load_ray_data(Data_moment *data_block, Ray *ray)
388 Data_moment_hdr data_hdr;
391 float value, scale, offset;
393 Range (*invf)(float x);
396 data_hdr = data_block->data_hdr;
397 data = data_block->data;
400 printf("wsr88d_load_ray_data: ");
401 printf(" DataName: %s\n", data_hdr.dataname);
404 ngates = data_hdr.ngates;
406 printf(" ngates = %d\n", ngates);
407 printf(" scale = %f\n", data_hdr.scale);
408 printf(" offset = %f\n", data_hdr.offset);
410 offset = data_hdr.offset;
411 scale = data_hdr.scale;
413 /* Note: data range is 2-255. 0 means signal is below threshold, and 1
414 * means range folded.
418 for (i=0; i < 10; i++) {
419 value = (data[i] - offset) / scale;
420 printf(" bytevalue = %d, value = %f\n", data[i], value);
425 if (strncmp(data_hdr.dataname, "DREF", 4) == 0) {
426 /* convert data to float
427 * convert float to range and put in ray.range
431 for (i = 0; i < ngates; i++) {
433 value = (data[i] - offset) / scale;
434 else value = (data[i] == 0) ? BADVAL : RFVAL;
435 ray->range[i] = invf(value);
442 if (strncmp(data_hdr.dataname, "DVEL", 4) == 0) {
443 /* convert data to float
444 * convert float to range and put in ray.range
448 for (i = 0; i < ngates; i++) {
450 value = (data[i] - offset) / scale;
451 else value = (data[i] == 0) ? BADVAL : RFVAL;
452 ray->range[i] = invf(value);
459 if (strncmp(data_hdr.dataname, "DSW", 3) == 0) {
460 /* convert data to float
461 * convert float to range and put in ray.range
465 for (i = 0; i < ngates; i++) {
467 value = (data[i] - offset) / scale;
468 else value = (data[i] == 0) ? BADVAL : RFVAL;
469 ray->range[i] = invf(value);
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;
481 void wsr88d_load_ray_hdr(Wsr88d_ray_m31 wsr88d_ray, Ray *ray)
483 int month, day, year, hour, minute, sec;
486 Ray_header_m31 ray_hdr;
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;
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;
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.;
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;
512 /* Get some values using message type 1 routines.
513 * First load VCP and elevation numbers into msg 1 ray.
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;
529 int wsr88d_get_vol_index(char* dataname)
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. */
542 void wsr88d_load_ray_into_radar(Wsr88d_ray_m31 wsr88d_ray, int isweep, int iray,
546 Range (*invf)(float x);
548 int vol_index, waveform;
550 enum waveforms {surveillance=1, doppler_ambres, doppler_no_ambres, batch};
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
557 put this data moment in ray
558 endfor each data moment
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.
572 if (wsr88d_ray.ray_hdr.dbptr_ref > 0) {
573 vol_index = wsr88d_get_vol_index(wsr88d_ray.ref->data_hdr.dataname);
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;
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.
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))
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;
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;
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;
608 if (wsr88d_ray.ray_hdr.dbptr_vel > 0) {
609 vol_index = wsr88d_get_vol_index(wsr88d_ray.vel->data_hdr.dataname);
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;
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;
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;
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;
633 if (wsr88d_ray.ray_hdr.dbptr_sw > 0) {
634 vol_index = wsr88d_get_vol_index(wsr88d_ray.sw->data_hdr.dataname);
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;
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;
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;
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;
662 void wsr88d_load_sweep_header(Radar *radar, int isweep,
663 Wsr88d_ray_m31 wsr88d_ray)
666 Ray_header_m31 ray_hdr;
670 ray_hdr = wsr88d_ray.ray_hdr;
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.;
688 Radar *load_wsr88d_m31_into_radar(Wsr88d_file *wf)
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;
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.
709 n = fread(&msghdr, sizeof(Wsr88d_msg_hdr), 1, wf->fptr);
711 /* printf("msgtype = %d\n", msghdr.msg_type); */
712 msg_hdr_size = sizeof(Wsr88d_msg_hdr) - sizeof(msghdr.rpg);
715 radar = RSL_new_radar(MAX_RADAR_VOLUMES);
717 while (! end_of_vos) {
718 if (msghdr.msg_type == 31) {
719 if (little_endian()) wsr88d_swap_m31_hdr(&msghdr);
721 /* Get size in bytes of message following message header.
722 * The size given in message header is in halfwords, so double it.
724 msg_size = (int) msghdr.msg_size * 2 - msg_hdr_size;
726 n = read_wsr88d_ray_m31(wf, msg_size, &wsr88d_ray);
727 if (n <= 0) return NULL;
729 /* Load this ray into radar structure ray. */
730 wsr88d_load_ray_into_radar(wsr88d_ray, isweep, iray, radar);
733 else { /* msg_type not 31 */
734 n = fread(&non31_seg_remainder, sizeof(non31_seg_remainder), 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);
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); */
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);
760 if (wsr88d_ray.ray_hdr.radial_status != END_VOS) {
761 n = fread(&msghdr, sizeof(Wsr88d_msg_hdr), 1, wf->fptr);
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);
775 wsr88d_load_sweep_header(radar, isweep, wsr88d_ray);
777 if (feof(wf->fptr) != 0) end_of_vos = 1;