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

   1 #include <string.h>
   2 #include <math.h>
   3 #include <stdlib.h>
   4 #ifndef _rsl_h
   5 #include "rsl.h"
   6 #endif
   7
   8 void replace_ray_gaps_with_null(Sweep *sweep)
   9 {
  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.
  13      */
  14
  15     int i, j, maxrays, nrays;
  16     Sweep *newsweep;
  17
  18     if (sweep == NULL) {
  19         fprintf(stderr,"replace_ray_gaps_with_null: Null sweep.\n");
  20         return;
  21     }
  22
  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.
  26      */
  27     if (sweep->ray[nrays-1] && sweep->ray[nrays-1]->h.ray_num == nrays) return;
  28
  29     /* Determine expected maximum rays from beamwidth. */
  30     maxrays = 720;
  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);
  34
  35     /* Copy rays to new sweep based on ray number.  New sweep will contain
  36      * NULL where rays are missing.
  37      */
  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++];
  41
  42     /* Copy rays back to original sweep. */
  43     for (i=0; i<maxrays; i++) sweep->ray[i] = newsweep->ray[i];
  44
  45     sweep->h.nrays = maxrays;
  46     free(newsweep);
  47 }
  48
  49
  50 int get_ray_index(Sweep *sweep, Ray *ray)
  51 {
  52     /* Find the array index for the given ray in the sweep structure. */
  53
  54     float target_azimuth, this_azimuth;
  55     int iray, incr, maxrays, nchk;
  56
  57     if (ray == NULL) return -1;
  58
  59     target_azimuth = ray->h.azimuth;
  60
  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)
  64         return iray;
  65
  66     /* Search for ray with matching azimuth by incrementing ray index forward
  67      * or backward, depending on azimuth.
  68      */
  69     this_azimuth = sweep->ray[iray]->h.azimuth;
  70     if (this_azimuth < target_azimuth) incr = 1;
  71     else incr = -1;
  72     iray+=incr;
  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;
  78         iray+=incr;
  79     }
  80     /* Control shouldn't end up here.  Ray came from this sweep, so there
  81      * has to be a match.
  82      */
  83     fprintf(stderr,"get_ray_index: Didn't find matching azimuth, but there "
  84             "should be one.\n");
  85     fprintf(stderr,"Current values:\n");
  86     fprintf(stderr, "target_azimuth: %.2f\n"
  87             "this_azimuth:   %.2f\n"
  88             "iray:           %d\n"
  89             "incr:           %d\n"
  90             "maxrays:        %d\n"
  91             "nchk:           %d\n",
  92             target_azimuth, this_azimuth, iray, incr, maxrays, nchk);
  93     return -1;
  94 }
  95
  96
  97 int get_azimuth_match(Sweep *sweep1, Sweep *sweep2, int *iray1, int *iray2)
  98 {
  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
 101      * argument iray2.
 102      *
 103      * Return value: function returns 1 for successful match, 0 otherwise.
 104      */
 105
 106     Ray *ray2;
 107     float matchlimit;
 108     const notfound = 0, found = 1;
 109
 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;
 113
 114     ray2 = RSL_get_closest_ray_from_sweep(sweep2,
 115             sweep1->ray[*iray1]->h.azimuth, matchlimit);
 116
 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);
 120
 121     if (ray2 == NULL) {
 122         *iray2 = -1;
 123         return notfound;
 124     }
 125
 126     *iray2 = get_ray_index(sweep2, ray2);
 127     return found;
 128 }
 129
 130
 131 int get_first_azimuth_match(Sweep *dzsweep, Sweep *vrsweep, int *dz_ray_index,
 132         int *vr_ray_index)
 133 {
 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.
 137      */
 138
 139     Ray *dz_ray, *vr_ray;
 140     int iray_dz = 0, iray_vr = 0, match_found = 0;
 141     const notfound = 0, found = 1;
 142
 143     /* Get first non-null VR ray. */
 144     vr_ray = vrsweep->ray[iray_vr];
 145     while (vr_ray  == NULL && iray_vr < vrsweep->h.nrays)
 146         iray_vr++;
 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);
 151         return notfound;
 152     }
 153
 154     /* Get DZ ray with azimuth closest to first VR ray. */
 155
 156     while (!match_found) {
 157         if (get_azimuth_match(vrsweep, dzsweep, &iray_vr, &iray_dz))
 158             match_found = 1;
 159         else {
 160             /* If still no match, get next non-null VR ray. */
 161             iray_vr++;
 162             while ((vr_ray = vrsweep->ray[iray_vr]) == NULL) {
 163                 iray_vr++;
 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);
 171                     return notfound;
 172                 }
 173             }
 174         }
 175     }
 176
 177     *vr_ray_index = iray_vr;
 178     *dz_ray_index = iray_dz;
 179     return found;
 180 }
 181
 182
 183 void reorder_rays(Sweep *dzsweep, Sweep *vrsweep, Sweep *swsweep)
 184 {
 185     /* Reorder rays of the VR and SW sweeps so that they match by azimuth
 186      * the rays in the DZ sweep.
 187      */
 188     int i, iray_dz, iray_vr, iray_dz_start, ray_num;
 189     int maxrays;
 190     int more_rays;
 191     Sweep *new_vrsweep, *new_swsweep;
 192
 193     /* Allocate new VR and SW sweeps.  Allocate maximum number of rays based
 194      * on beam width.
 195      */
 196
 197     if (vrsweep != NULL) {
 198         maxrays = (vrsweep->h.beam_width < 0.8) ? 720 : 360;
 199         if (vrsweep->h.nrays > maxrays) maxrays = vrsweep->h.nrays;
 200     }
 201     else {
 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);
 205         return;
 206     }
 207
 208     replace_ray_gaps_with_null(dzsweep);
 209     replace_ray_gaps_with_null(vrsweep);
 210     if (swsweep) replace_ray_gaps_with_null(swsweep);
 211
 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);
 217         return;
 218     }
 219
 220     /* If DZ and VR azimuths already match, no need to do ray alignment. */
 221     if (iray_dz == iray_vr) return;
 222
 223     new_vrsweep = RSL_new_sweep(maxrays);
 224     memcpy(&new_vrsweep->h, &vrsweep->h, sizeof(Sweep_header));
 225     if (swsweep) {
 226         new_swsweep = RSL_new_sweep(maxrays);
 227         memcpy(&new_swsweep->h, &swsweep->h, sizeof(Sweep_header));
 228     }
 229
 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.
 232      */
 233
 234     iray_dz_start = iray_dz;
 235     ray_num = dzsweep->ray[iray_dz]->h.ray_num;
 236     more_rays = 1;
 237     while (more_rays) {
 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;
 241         }
 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;
 245         }
 246         iray_dz++;
 247         iray_vr++;
 248         ray_num++;
 249
 250         /* Handle sweep wraparound. */
 251         if (iray_dz == dzsweep->h.nrays || iray_dz == maxrays) {
 252             iray_dz = 0;
 253             ray_num = 1;
 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);
 260             }
 261         }
 262         if (iray_vr == vrsweep->h.nrays) iray_vr = 0;
 263
 264         /* If we're back at first match, we're done. */
 265         if (iray_dz == iray_dz_start) more_rays = 0;
 266     }
 267
 268     /* Copy ray pointers in new order back to original sweeps. */
 269     /* TODO: Can we simply replace the sweeps themselves with the new sweeps? */
 270
 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];
 274     }
 275
 276     /* TODO:
 277      * Update nrays in sweep header with ray number of last ray
 278      * (ray count includes NULL rays).
 279      */
 280     vrsweep->h.nrays = maxrays;
 281     if (swsweep) swsweep->h.nrays = maxrays;
 282 }
 283
 284
 285 void align_this_split_cut(Radar *radar, int iswp)
 286 {
 287     int vrindex[] = {VR_INDEX, V2_INDEX, V3_INDEX};
 288     int swindex[] = {SW_INDEX, S2_INDEX, S3_INDEX};
 289     int ivr, nvrfields;
 290     Sweep *dzsweep, *vrsweep, *swsweep;
 291
 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);
 296         return;
 297     }
 298
 299     if (radar->h.vcp != 121) nvrfields = 1;
 300     else {
 301         if (iswp < 4) nvrfields = 3;
 302         else nvrfields = 2;
 303     }
 304
 305     for (ivr = 0; ivr < nvrfields; ivr++) {
 306         if (radar->v[vrindex[ivr]]) {
 307             vrsweep = radar->v[vrindex[ivr]]->sweep[iswp];
 308         }
 309         else {
 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);
 313             return;
 314         }
 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 "
 318                     "VR don't match:\n"
 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);
 323             return;
 324         }
 325         if (radar->v[swindex[ivr]])
 326             swsweep = radar->v[swindex[ivr]]->sweep[iswp];
 327         else
 328             swsweep = NULL;
 329         reorder_rays(dzsweep, vrsweep, swsweep);
 330     }
 331 }
 332
 333
 334 static Radar *copy_radar_for_split_cut_adjustments(Radar *r)
 335 {
 336   /* Return a copy of the radar structure for modification at split cuts. */
 337
 338   Radar *r_new;
 339   int i,j,k;
 340
 341   r_new = RSL_new_radar(MAX_RADAR_VOLUMES);
 342   r_new->h = r->h;
 343   /* Copy the contents of VR and SW volumes, which will be modified. For other
 344    * fields, just copy pointers.
 345    */
 346   for (k=0; k < MAX_RADAR_VOLUMES; k++) {
 347     if (r->v[k] == NULL) continue;
 348     switch(k) {
 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];
 359             }
 360         }
 361         break;
 362       default:
 363         /* For fields other than VR or SW types, simply copy volume pointer.
 364          */
 365         r_new->v[k] = r->v[k];
 366         break;
 367     }
 368   }
 369   return r_new;
 370 }
 371
 372
 373 Radar *wsr88d_align_split_cut_rays(Radar *radar_in)
 374 {
 375 /* This is the driver function for aligning rays in split cuts by azimuth. */
 376
 377     int iswp;
 378     int nsplit, vcp;
 379     Radar *radar;
 380
 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");
 388         return radar_in;
 389     }
 390
 391     vcp = radar_in->h.vcp;
 392
 393     /* How many split cuts are in this VCP? */
 394     nsplit = 0;
 395     switch (vcp) {
 396         case 11: case 21: case 32: case 211: case 221:
 397             nsplit = 2;
 398             break;
 399         case 12: case 31: case 212:
 400             nsplit = 3;
 401             break;
 402         case 121:
 403             nsplit = 5;
 404             break;
 405     }
 406
 407     if (nsplit == 0) {
 408         fprintf(stderr, "wsr88d_align_split_cut_rays: Unknown VCP: %d\n", vcp);
 409         return radar_in;
 410     }
 411
 412     /* Make a copy of radar.  It is the copy that we'll modify. */
 413     radar = copy_radar_for_split_cut_adjustments(radar_in);
 414
 415     /* Align rays in each split cut. */
 416
 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);
 422             return radar;
 423         }
 424         align_this_split_cut(radar, iswp);
 425     }
 426
 427     return radar;
 428 }