8 void replace_ray_gaps_with_null(Sweep *sweep)
10 /* Check for ray gaps and replace any with null. A ray gap means that
11 * the difference in ray numbers in a pair of consecutive rays is greater
12 * than 1. When such a case is found, nulls are inserted for missing rays.
15 int i, j, maxrays, nrays;
19 fprintf(stderr,"replace_ray_gaps_with_null: Null sweep.\n");
23 nrays = sweep->h.nrays;
24 /* If last ray of sweep is non-null, and its ray number is same as nrays,
25 * then there are no gaps.
27 if (sweep->ray[nrays-1] && sweep->ray[nrays-1]->h.ray_num == nrays) return;
29 /* Determine expected maximum rays from beamwidth. */
31 if (sweep->h.beam_width > 0.9) maxrays = 360;
32 if (sweep->h.nrays > maxrays) maxrays = sweep->h.nrays;
33 newsweep = RSL_new_sweep(maxrays);
35 /* Copy rays to new sweep based on ray number. New sweep will contain
36 * NULL where rays are missing.
38 for (i=0, j=0; i<maxrays; i++)
39 if (sweep->ray[j] && sweep->ray[j]->h.ray_num == i+1)
40 newsweep->ray[i] = sweep->ray[j++];
42 /* Copy rays back to original sweep. */
43 for (i=0; i<maxrays; i++) sweep->ray[i] = newsweep->ray[i];
45 sweep->h.nrays = maxrays;
50 int get_ray_index(Sweep *sweep, Ray *ray)
52 /* Find the array index for the given ray in the sweep structure. */
54 float target_azimuth, this_azimuth;
55 int iray, incr, maxrays, nchk;
57 if (ray == NULL) return -1;
59 target_azimuth = ray->h.azimuth;
61 /* Start with ray_num - 1. That will normally be a match. */
62 iray = ray->h.ray_num - 1;
63 if (sweep->ray[iray]->h.azimuth == target_azimuth)
66 /* Search for ray with matching azimuth by incrementing ray index forward
67 * or backward, depending on azimuth.
69 this_azimuth = sweep->ray[iray]->h.azimuth;
70 if (this_azimuth < target_azimuth) incr = 1;
73 maxrays = sweep->h.nrays;
74 for (nchk = 1; nchk <= maxrays; nchk++) {
75 if (iray < 0) iray = maxrays - 1;
76 if (iray >= maxrays) iray = 0;
77 if (sweep->ray[iray]->h.azimuth == target_azimuth) return iray;
80 /* Control shouldn't end up here. Ray came from this sweep, so there
83 fprintf(stderr,"get_ray_index: Didn't find matching azimuth, but there "
85 fprintf(stderr,"Current values:\n");
86 fprintf(stderr, "target_azimuth: %.2f\n"
87 "this_azimuth: %.2f\n"
92 target_azimuth, this_azimuth, iray, incr, maxrays, nchk);
97 int get_azimuth_match(Sweep *sweep1, Sweep *sweep2, int *iray1, int *iray2)
99 /* Using the azimuth of sweep1->ray[iray1], find the ray in sweep2 with
100 * matching azimuth. The index of matching ray is returned through
103 * Return value: function returns 1 for successful match, 0 otherwise.
108 const notfound = 0, found = 1;
110 matchlimit = sweep2->h.horz_half_bw;
111 /* Assure that matchlimit is wide enough--a problem with pre-Build10. */
112 if (sweep2->h.beam_width > .9) matchlimit = 0.5;
114 ray2 = RSL_get_closest_ray_from_sweep(sweep2,
115 sweep1->ray[*iray1]->h.azimuth, matchlimit);
117 /* If no match, try again with a slightly wider limit. */
118 if (ray2 == NULL) ray2 = RSL_get_closest_ray_from_sweep(sweep2,
119 sweep1->ray[*iray1]->h.azimuth, matchlimit + .04 * matchlimit);
126 *iray2 = get_ray_index(sweep2, ray2);
131 int get_first_azimuth_match(Sweep *dzsweep, Sweep *vrsweep, int *dz_ray_index,
134 /* Find the DZ ray whose azimuth matches the first ray of VR sweep.
135 * Indexes for the matching rays are returned through arguments.
136 * Return value: 1 for successful match, 0 otherwise.
139 Ray *dz_ray, *vr_ray;
140 int iray_dz = 0, iray_vr = 0, match_found = 0;
141 const notfound = 0, found = 1;
143 /* Get first non-null VR ray. */
144 vr_ray = vrsweep->ray[iray_vr];
145 while (vr_ray == NULL && iray_vr < vrsweep->h.nrays)
147 if (vr_ray == NULL) {
148 fprintf(stderr,"get_first_azimuth_match: All rays in VR sweep number %d"
149 " are NULL.\nMerge split cut not done for this sweep.\n",
150 vrsweep->h.sweep_num);
154 /* Get DZ ray with azimuth closest to first VR ray. */
156 while (!match_found) {
157 if (get_azimuth_match(vrsweep, dzsweep, &iray_vr, &iray_dz))
160 /* If still no match, get next non-null VR ray. */
162 while ((vr_ray = vrsweep->ray[iray_vr]) == NULL) {
164 if (iray_vr >= vrsweep->h.nrays) {
165 fprintf(stderr,"get_first_azimuth_match: Couldn't find a "
166 "matching azimuth.\nSweeps compared:\n");
167 fprintf(stderr,"DZ: sweep number %d, elev %f\n"
168 "VR: sweep number %d, elev %f\n",
169 dzsweep->h.sweep_num, dzsweep->h.elev,
170 vrsweep->h.sweep_num, vrsweep->h.elev);
177 *vr_ray_index = iray_vr;
178 *dz_ray_index = iray_dz;
183 void reorder_rays(Sweep *dzsweep, Sweep *vrsweep, Sweep *swsweep)
185 /* Reorder rays of the VR and SW sweeps so that they match by azimuth
186 * the rays in the DZ sweep.
188 int i, iray_dz, iray_vr, iray_dz_start, ray_num;
191 Sweep *new_vrsweep, *new_swsweep;
193 /* Allocate new VR and SW sweeps. Allocate maximum number of rays based
197 if (vrsweep != NULL) {
198 maxrays = (vrsweep->h.beam_width < 0.8) ? 720 : 360;
199 if (vrsweep->h.nrays > maxrays) maxrays = vrsweep->h.nrays;
202 fprintf(stderr,"reorder_rays: Velocity sweep number %d is NULL.\n"
203 "Merge split cut not done for this sweep.\n",
204 vrsweep->h.sweep_num);
208 replace_ray_gaps_with_null(dzsweep);
209 replace_ray_gaps_with_null(vrsweep);
210 if (swsweep) replace_ray_gaps_with_null(swsweep);
212 /* Get ray indices for first match of DZ and VR azimuths. */
213 if (get_first_azimuth_match(dzsweep, vrsweep, &iray_dz, &iray_vr) == 0) {
214 fprintf(stderr,"reorder_rays: Merge split cut not done for DZ sweep "
215 "number %d and VR sweep %d.\n",
216 dzsweep->h.sweep_num, vrsweep->h.sweep_num);
220 /* If DZ and VR azimuths already match, no need to do ray alignment. */
221 if (iray_dz == iray_vr) return;
223 new_vrsweep = RSL_new_sweep(maxrays);
224 memcpy(&new_vrsweep->h, &vrsweep->h, sizeof(Sweep_header));
226 new_swsweep = RSL_new_sweep(maxrays);
227 memcpy(&new_swsweep->h, &swsweep->h, sizeof(Sweep_header));
230 /* Copy VR rays to new VR sweep such that their azimuths match those of
231 * DZ rays at the same indices. Do the same for SW.
234 iray_dz_start = iray_dz;
235 ray_num = dzsweep->ray[iray_dz]->h.ray_num;
238 if (vrsweep->ray[iray_vr]) {
239 new_vrsweep->ray[iray_dz] = vrsweep->ray[iray_vr];
240 new_vrsweep->ray[iray_dz]->h.ray_num = ray_num;
242 if (swsweep && swsweep->ray[iray_vr]) {
243 new_swsweep->ray[iray_dz] = swsweep->ray[iray_vr];
244 new_swsweep->ray[iray_dz]->h.ray_num = ray_num;
250 /* Handle sweep wraparound. */
251 if (iray_dz == dzsweep->h.nrays || iray_dz == maxrays) {
254 /* Check if azimuth rematch is needed. */
255 if (dzsweep->ray[iray_dz] && vrsweep->ray[iray_vr] &&
256 fabsf(dzsweep->ray[iray_dz]->h.azimuth -
257 vrsweep->ray[iray_vr]->h.azimuth)
258 > dzsweep->h.horz_half_bw) {
259 get_azimuth_match(dzsweep, vrsweep, &iray_dz, &iray_vr);
262 if (iray_vr == vrsweep->h.nrays) iray_vr = 0;
264 /* If we're back at first match, we're done. */
265 if (iray_dz == iray_dz_start) more_rays = 0;
268 /* Copy ray pointers in new order back to original sweeps. */
269 /* TODO: Can we simply replace the sweeps themselves with the new sweeps? */
271 for (i=0; i < maxrays; i++) {
272 vrsweep->ray[i] = new_vrsweep->ray[i];
273 if (swsweep) swsweep->ray[i] = new_swsweep->ray[i];
277 * Update nrays in sweep header with ray number of last ray
278 * (ray count includes NULL rays).
280 vrsweep->h.nrays = maxrays;
281 if (swsweep) swsweep->h.nrays = maxrays;
285 void align_this_split_cut(Radar *radar, int iswp)
287 int vrindex[] = {VR_INDEX, V2_INDEX, V3_INDEX};
288 int swindex[] = {SW_INDEX, S2_INDEX, S3_INDEX};
290 Sweep *dzsweep, *vrsweep, *swsweep;
292 dzsweep = radar->v[DZ_INDEX]->sweep[iswp];
293 if (dzsweep == NULL) {
294 fprintf(stderr,"align_this_split_cut: DZ sweep at index %d is NULL.\n"
295 "Merge split cut not done for this sweep.\n", iswp);
299 if (radar->h.vcp != 121) nvrfields = 1;
301 if (iswp < 4) nvrfields = 3;
305 for (ivr = 0; ivr < nvrfields; ivr++) {
306 if (radar->v[vrindex[ivr]]) {
307 vrsweep = radar->v[vrindex[ivr]]->sweep[iswp];
310 fprintf(stderr,"align_this_split_cut: No VR at INDEX %d.\n"
311 "No ray alignment done for this VR index at DZ sweep "
312 "index %d.\n", vrindex[ivr], iswp);
315 /* Make sure vrsweep and dzsweep elevs are approximately the same. */
316 if (fabsf(dzsweep->h.elev - vrsweep->h.elev) > 0.1) {
317 fprintf(stderr,"align_this_split_cut: elevation angles for DZ and "
319 "radar->v[DZ_INDEX=%d]->sweep[%d]->h.elev = %f\n"
320 "radar->v[VR_INDEX=%d]->sweep[%d]->h.elev = %f\n",
321 DZ_INDEX, iswp, dzsweep->h.elev,
322 VR_INDEX, iswp, vrsweep->h.elev);
325 if (radar->v[swindex[ivr]])
326 swsweep = radar->v[swindex[ivr]]->sweep[iswp];
329 reorder_rays(dzsweep, vrsweep, swsweep);
334 static Radar *copy_radar_for_split_cut_adjustments(Radar *r)
336 /* Return a copy of the radar structure for modification at split cuts. */
341 r_new = RSL_new_radar(MAX_RADAR_VOLUMES);
343 /* Copy the contents of VR and SW volumes, which will be modified. For other
344 * fields, just copy pointers.
346 for (k=0; k < MAX_RADAR_VOLUMES; k++) {
347 if (r->v[k] == NULL) continue;
349 case VR_INDEX: case SW_INDEX: case V2_INDEX: case S2_INDEX:
350 case V3_INDEX: case S3_INDEX:
351 r_new->v[k] = RSL_new_volume(r->v[k]->h.nsweeps);
352 r_new->v[k]->h = r->v[k]->h;
353 for (i=0; i < r->v[k]->h.nsweeps; i++) {
354 /* Add new sweeps to volume and copy ray pointers. */
355 r_new->v[k]->sweep[i] = RSL_new_sweep(r->v[k]->sweep[i]->h.nrays);
356 r_new->v[k]->sweep[i]->h = r->v[k]->sweep[i]->h;
357 for (j=0; j < r->v[k]->sweep[i]->h.nrays; j++) {
358 r_new->v[k]->sweep[i]->ray[j] = r->v[k]->sweep[i]->ray[j];
363 /* For fields other than VR or SW types, simply copy volume pointer.
365 r_new->v[k] = r->v[k];
373 Radar *wsr88d_align_split_cut_rays(Radar *radar_in)
375 /* This is the driver function for aligning rays in split cuts by azimuth. */
381 /* If VR is not present, nothing more to do. */
382 if (radar_in->v[VR_INDEX] == NULL) {
383 /* If SW without VR, notify user that VR must be included. */
384 if (radar_in->v[SW_INDEX])
385 fprintf(stderr,"wsr88d_align_split_cut_rays: VR must be selected in"
386 " order to do ray alignment\n"
387 " on SW. No ray alignment for SW in split cuts.\n");
391 vcp = radar_in->h.vcp;
393 /* How many split cuts are in this VCP? */
396 case 11: case 21: case 32: case 211: case 221:
399 case 12: case 31: case 212:
408 fprintf(stderr, "wsr88d_align_split_cut_rays: Unknown VCP: %d\n", vcp);
412 /* Make a copy of radar. It is the copy that we'll modify. */
413 radar = copy_radar_for_split_cut_adjustments(radar_in);
415 /* Align rays in each split cut. */
417 for (iswp = 0; iswp < nsplit; iswp++) {
418 if (radar->v[DZ_INDEX]->sweep[iswp] == NULL) {
419 fprintf(stderr, "wsr88d_align_split_cut_rays: Encountered "
420 "unexpected NULL DZ sweep\n");
421 fprintf(stderr, "iswp (sweep index) = %d\n", iswp);
424 align_this_split_cut(radar, iswp);