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

   1 /*
   2     NASA/TRMM, Code 910.1.
   3     This is the TRMM Office Radar Software Library.
   4     Copyright (C) 1996, 1997
   5             John H. Merritt
   6             Space Applications Corporation
   7             Vienna, Virginia
   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 Free
  21     Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  22 */
  23 #include <stdio.h>
  24 #include <string.h>
  25 #include <stdlib.h>
  26 #include <unistd.h>
  27
  28 #include "rsl.h"
  29
  30 extern int radar_verbose_flag;
  31 typedef short UF_buffer[16384]; /* Some UF files are bigger than 4096
  32                                  * that the UF doc's specify.
  33                                  */
  34
  35 #define UF_MORE 0
  36 #define UF_DONE 1
  37
  38 static float (*f)(Range x);
  39 static Range (*invf)(float x);
  40
  41 Volume *reset_nsweeps_in_volume(Volume *volume)
  42 {
  43   int i;
  44   if (volume == NULL) return NULL;
  45   for (i=volume->h.nsweeps; i>0; i--)
  46         if (volume->sweep[i-1] != NULL) {
  47           volume->h.nsweeps = i;
  48           break;
  49         }
  50   return volume;
  51 }
  52 Radar *reset_nsweeps_in_all_volumes(Radar *radar)
  53 {
  54   int i;
  55   if (radar == NULL) return NULL;
  56   for (i=0; i<radar->h.nvolumes; i++)
  57         radar->v[i] = reset_nsweeps_in_volume(radar->v[i]);
  58   return radar;
  59 }
  60
  61 Volume *copy_sweeps_into_volume(Volume *new_volume, Volume *old_volume)
  62 {
  63   int i;
  64   int nsweeps;
  65   if (old_volume == NULL) return new_volume;
  66   if (new_volume == NULL) return new_volume;
  67   nsweeps = new_volume->h.nsweeps; /* Save before copying old header. */
  68   new_volume->h = old_volume->h;
  69   new_volume->h.nsweeps = nsweeps;
  70   for (i=0; i<old_volume->h.nsweeps; i++)
  71         new_volume->sweep[i] = old_volume->sweep[i]; /* Just copy pointers. */
  72   /* Free the old sweep array, not the pointers to sweeps. */
  73   free(old_volume->sweep);
  74   return new_volume;
  75 }
  76
  77 void swap2(short *buf, int n)
  78 {
  79   short *end_addr;
  80   end_addr = buf + n;
  81   while (buf < end_addr) {
  82         swap_2_bytes(buf);
  83         buf++;
  84   }
  85 }
  86
  87
  88 /********************************************************************/
  89 /*********************************************************************/
  90 /*                                                                   */
  91 /*                  uf_into_radar                                    */
  92 /*                                                                   */
  93 /*  By: John Merritt                                                 */
  94 /*      Space Applications Corporation                               */
  95 /*      August 26, 1994                                              */
  96 /*********************************************************************/
  97 int uf_into_radar(UF_buffer uf, Radar **the_radar)
  98 {
  99
 100 /* Missing data flag : -32768 when a signed short. */
 101 #define UF_NO_DATA 0X8000
 102
 103   /* Any convensions may be observed, however, Radial Velocity must be VE. */
 104   /* Typically:
 105    *    DM = Reflectivity (dB(mW)).
 106    *    DZ = Reflectivity (dBZ).
 107    *    VR = Radial Velocity.
 108    *    CZ = Corrected Reflectivity. (Quality controlled: AP removed, etc.)
 109    *    SW = Spectrum Width.
 110    *    DR = Differential Reflectivity.
 111    *    XZ = X-band Reflectivity.
 112    *
 113    * These fields may appear in any order in the UF file.
 114    *
 115    * RETURN:
 116    *   UF_DONE if we're done with the UF ingest.
 117    *   UF_MORE if we need more UF data.
 118    */
 119
 120   /* These are pointers to various locations within the UF buffer 'uf'.
 121    * They are used to index the different components of the UF structure in
 122    * a manor consistant with the UF documentation.  For instance, uf_ma[1]
 123    * will be equivalenced to the second word (2 bytes/each) of the UF
 124    * buffer.
 125    */
 126   short *uf_ma;  /* Mandatory header block. */
 127   short *uf_lu;  /* Local Use header block.  */
 128   short *uf_dh;  /* Data header.  */
 129   short *uf_fh;  /* Field header. */
 130   short *uf_data; /* Data. */
 131
 132   /* The length of each header. */
 133   int len_data, len_lu;
 134
 135   int current_fh_index;
 136   float scale_factor;
 137
 138   int nfields, isweep, ifield, iray, i, m;
 139   static int pulled_time_from_first_ray = 0;
 140   char *field_type; /* For printing the field type upon error. */
 141         short proj_name[4];
 142   Ray *ray;
 143   Sweep *sweep;
 144   Radar *radar;
 145   float x;
 146   short missing_data;
 147   Volume *new_volume;
 148   int nbins;
 149   extern int rsl_qfield[];
 150   extern int *rsl_qsweep; /* See RSL_read_these_sweeps in volume.c */
 151   extern int rsl_qsweep_max;
 152
 153   radar = *the_radar;
 154
 155 /*
 156  * The organization of the Radar structure is by volumes, then sweeps, then
 157  * rays, then gates.  This is different from the UF file organization.
 158  * The UF format is sweeps, rays, then gates for all field types (volumes).
 159  */
 160
 161
 162 /* Set up all the UF pointers. */
 163   uf_ma = uf;
 164   uf_lu = uf + uf_ma[3] - 1;
 165   uf_dh = uf + uf_ma[4] - 1;
 166
 167   nfields =  uf_dh[0];
 168   isweep = uf_ma[9] - 1;
 169
 170   if (rsl_qsweep != NULL) {
 171         if (isweep > rsl_qsweep_max) return UF_DONE;
 172         if (rsl_qsweep[isweep] == 0) return UF_MORE;
 173   }
 174
 175
 176 /* Here is a sticky part.  We must make sure that if we encounter any
 177  * additional fields that were not previously present, that we are able
 178  * to load them.  This will require us to copy the entire radar structure
 179  * and whack off the old one.  But, we must be sure that it really is a
 180  * new field.  This is not so trivial as a couple of lines of code; I will
 181  * have to think about this a little bit more.  See STICKYSOLVED below.
 182  */
 183 #ifdef STICKYSOLVED
 184   if (radar == NULL) radar = RSL_new_radar(nfields);
 185   /* Sticky solution here. */
 186 #else
 187   if (radar == NULL) {
 188         radar = RSL_new_radar(MAX_RADAR_VOLUMES);
 189         *the_radar = radar;
 190         pulled_time_from_first_ray = 0;
 191         for (i=0; i<MAX_RADAR_VOLUMES; i++)
 192           if (rsl_qfield[i]) /* See RSL_select_fields in volume.c */
 193                 radar->v[i] = RSL_new_volume(20);
 194   }
 195
 196 #endif
 197 /* For LITTLE ENDIAN:
 198  * WE "UNSWAP" character strings.  Because there are so few of them,
 199  * it is easier to catch them here.  The entire UF buffer is swapped prior
 200  * to entry to here, therefore, undo-ing these swaps; sets the
 201  * character strings right.
 202  */
 203
 204   for (i=0; i<nfields; i++) {
 205         if (little_endian()) swap_2_bytes(&uf_dh[3+2*i]); /* Unswap. */
 206         if (strncmp((char *)&uf_dh[3+2*i], "DZ", 2) == 0) ifield = DZ_INDEX;
 207         else if (strncmp((char *)&uf_dh[3+2*i], "VR", 2) == 0) ifield = VR_INDEX;
 208         else if (strncmp((char *)&uf_dh[3+2*i], "SW", 2) == 0) ifield = SW_INDEX;
 209         else if (strncmp((char *)&uf_dh[3+2*i], "CZ", 2) == 0) ifield = CZ_INDEX;
 210         else if (strncmp((char *)&uf_dh[3+2*i], "UZ", 2) == 0) ifield = ZT_INDEX;
 211         else if (strncmp((char *)&uf_dh[3+2*i], "ZT", 2) == 0) ifield = ZT_INDEX;
 212         else if (strncmp((char *)&uf_dh[3+2*i], "DR", 2) == 0) ifield = DR_INDEX;
 213         else if (strncmp((char *)&uf_dh[3+2*i], "ZD", 2) == 0) ifield = ZD_INDEX;
 214         else if (strncmp((char *)&uf_dh[3+2*i], "DM", 2) == 0) ifield = DM_INDEX;
 215         else if (strncmp((char *)&uf_dh[3+2*i], "RH", 2) == 0) ifield = RH_INDEX;
 216         else if (strncmp((char *)&uf_dh[3+2*i], "PH", 2) == 0) ifield = PH_INDEX;
 217         else if (strncmp((char *)&uf_dh[3+2*i], "XZ", 2) == 0) ifield = XZ_INDEX;
 218         else if (strncmp((char *)&uf_dh[3+2*i], "CD", 2) == 0) ifield = CD_INDEX;
 219         else if (strncmp((char *)&uf_dh[3+2*i], "MZ", 2) == 0) ifield = MZ_INDEX;
 220         else if (strncmp((char *)&uf_dh[3+2*i], "MD", 2) == 0) ifield = MD_INDEX;
 221         else if (strncmp((char *)&uf_dh[3+2*i], "ZE", 2) == 0) ifield = ZE_INDEX;
 222         else if (strncmp((char *)&uf_dh[3+2*i], "VE", 2) == 0) ifield = VE_INDEX;
 223         else if (strncmp((char *)&uf_dh[3+2*i], "KD", 2) == 0) ifield = KD_INDEX;
 224         else if (strncmp((char *)&uf_dh[3+2*i], "TI", 2) == 0) ifield = TI_INDEX;
 225         else if (strncmp((char *)&uf_dh[3+2*i], "DX", 2) == 0) ifield = DX_INDEX;
 226         else if (strncmp((char *)&uf_dh[3+2*i], "CH", 2) == 0) ifield = CH_INDEX;
 227         else if (strncmp((char *)&uf_dh[3+2*i], "AH", 2) == 0) ifield = AH_INDEX;
 228         else if (strncmp((char *)&uf_dh[3+2*i], "CV", 2) == 0) ifield = CV_INDEX;
 229         else if (strncmp((char *)&uf_dh[3+2*i], "AV", 2) == 0) ifield = AV_INDEX;
 230         else if (strncmp((char *)&uf_dh[3+2*i], "SQ", 2) == 0) ifield = SQ_INDEX;
 231         else { /* etc.  DON'T know how to handle this yet. */
 232           field_type = (char *)&uf_dh[3+2*i];
 233           fprintf(stderr, "Unknown field type %c%c\n", (char)field_type[0], (char)field_type[1]);
 234           continue;
 235         }
 236         switch (ifield) {
 237         case DZ_INDEX: f = DZ_F; invf = DZ_INVF; break;
 238         case VR_INDEX: f = VR_F; invf = VR_INVF; break;
 239         case SW_INDEX: f = SW_F; invf = SW_INVF; break;
 240         case CZ_INDEX: f = CZ_F; invf = CZ_INVF; break;
 241         case ZT_INDEX: f = ZT_F; invf = ZT_INVF; break;
 242         case DR_INDEX: f = DR_F; invf = DR_INVF; break;
 243         case ZD_INDEX: f = ZD_F; invf = ZD_INVF; break;
 244         case DM_INDEX: f = DM_F; invf = DM_INVF; break;
 245         case RH_INDEX: f = RH_F; invf = RH_INVF; break;
 246         case PH_INDEX: f = PH_F; invf = PH_INVF; break;
 247         case XZ_INDEX: f = XZ_F; invf = XZ_INVF; break;
 248         case CD_INDEX: f = CD_F; invf = CD_INVF; break;
 249         case MZ_INDEX: f = MZ_F; invf = MZ_INVF; break;
 250         case MD_INDEX: f = MD_F; invf = MD_INVF; break;
 251         case ZE_INDEX: f = ZE_F; invf = ZE_INVF; break;
 252         case VE_INDEX: f = VE_F; invf = VE_INVF; break;
 253         case KD_INDEX: f = KD_F; invf = KD_INVF; break;
 254         case TI_INDEX: f = TI_F; invf = TI_INVF; break;
 255         case DX_INDEX: f = DX_F; invf = DX_INVF; break;
 256         case CH_INDEX: f = CH_F; invf = CH_INVF; break;
 257         case AH_INDEX: f = AH_F; invf = AH_INVF; break;
 258         case CV_INDEX: f = CV_F; invf = CV_INVF; break;
 259         case AV_INDEX: f = AV_F; invf = AV_INVF; break;
 260         case SQ_INDEX: f = SQ_F; invf = SQ_INVF; break;
 261         default: f = DZ_F; invf = DZ_INVF; break;
 262         }
 263
 264         /* Do we place the data into this volume? */
 265         if (radar->v[ifield] == NULL) continue; /* Nope. */
 266
 267         if (isweep >= radar->v[ifield]->h.nsweeps) { /* Exceeded sweep limit.
 268                                                       * Allocate more sweeps.
 269                                                       * Copy all previous sweeps.
 270                                                       */
 271           if (radar_verbose_flag)
 272                 fprintf(stderr,"Exceeded sweep allocation of %d. Adding 20 more.\n", isweep);
 273           new_volume = RSL_new_volume(radar->v[ifield]->h.nsweeps+20);
 274           new_volume = copy_sweeps_into_volume(new_volume, radar->v[ifield]);
 275           radar->v[ifield] = new_volume;
 276         }
 277
 278         if (radar->v[ifield]->sweep[isweep] == NULL) {
 279           if (radar_verbose_flag)
 280                 fprintf(stderr,"Allocating new sweep for field %d, isweep %d\n", ifield, isweep);
 281           radar->v[ifield]->sweep[isweep] = RSL_new_sweep(1000);
 282           radar->v[ifield]->sweep[isweep]->h.nrays = 0; /* Increment this for each
 283                                                          * ray encountered.
 284                                                          */
 285           radar->v[ifield]->h.f = f;
 286           radar->v[ifield]->h.invf = invf;
 287           radar->v[ifield]->sweep[isweep]->h.f = f;
 288           radar->v[ifield]->sweep[isweep]->h.invf = invf;
 289                 radar->v[ifield]->sweep[isweep]->h.sweep_num = uf_ma[9];
 290                 radar->v[ifield]->sweep[isweep]->h.elev = uf_ma[35] / 64.0;
 291         }
 292
 293
 294
 295         current_fh_index = uf_dh[4+2*i];
 296         uf_fh = uf + current_fh_index - 1;
 297         sweep = radar->v[ifield]->sweep[isweep];
 298         iray =  sweep->h.nrays;
 299         nbins = uf_fh[5];
 300         radar->v[ifield]->sweep[isweep]->ray[iray] = RSL_new_ray(nbins);
 301         ray   = radar->v[ifield]->sweep[isweep]->ray[iray];
 302         sweep->h.nrays += 1;
 303
 304
 305         if (ray) {
 306                 /*
 307                  * ---- Beginning of MANDATORY HEADER BLOCK.
 308                  */
 309           ray->h.ray_num = uf_ma[7];
 310           if (little_endian()) swap2(&uf_ma[10], 8);
 311           memcpy(radar->h.radar_name, &uf_ma[10], 8);
 312           if (little_endian()) swap2(&uf_ma[10], 8/2);
 313           memcpy(radar->h.name, &uf_ma[14], 8);
 314           if (little_endian()) swap2(&uf_ma[14], 8/2);
 315
 316       /* All components of lat/lon are the same sign.  If not, then
 317        * what ever wrote the UF was in error.  A simple RSL program
 318        * can repair the damage, however, not here.
 319        */
 320           ray->h.lat = uf_ma[18] + uf_ma[19]/60.0 + uf_ma[20]/64.0/3600;
 321           ray->h.lon = uf_ma[21] + uf_ma[22]/60.0 + uf_ma[23]/64.0/3600;
 322           ray->h.alt      = uf_ma[24];
 323           ray->h.year     = uf_ma[25];
 324           if (ray->h.year < 1900) {
 325             ray->h.year += 1900;
 326             if (ray->h.year < 1980) ray->h.year += 100; /* Year >= 2000. */
 327           }
 328           ray->h.month    = uf_ma[26];
 329           ray->h.day      = uf_ma[27];
 330           ray->h.hour     = uf_ma[28];
 331           ray->h.minute   = uf_ma[29];
 332           ray->h.sec      = uf_ma[30];
 333           ray->h.azimuth  = uf_ma[32] / 64.0;
 334
 335           /* If Local Use Header is present and contains azimuth, use that
 336            * azimuth for VR and SW. This is for WSR-88D, which runs separate
 337            * scans for DZ and VR/SW at the lower elevations, which means DZ
 338            * VR/SW and have different azimuths in the "same" ray.
 339            */
 340           len_lu = uf_ma[4] - uf_ma[3];
 341           if (len_lu == 2 && (ifield == VR_INDEX || ifield == SW_INDEX)) {
 342               if (strncmp((char *)uf_lu,"ZA",2) == 0 ||
 343                       strncmp((char *)uf_lu,"AZ",2) == 0)
 344                   ray->h.azimuth = uf_lu[1] / 64.0;
 345           }
 346           if (ray->h.azimuth < 0.) ray->h.azimuth += 360.; /* make it 0 to 360. */
 347           ray->h.elev     = uf_ma[33] / 64.0;
 348           ray->h.elev_num = sweep->h.sweep_num;
 349           ray->h.fix_angle  = sweep->h.elev = uf_ma[35] / 64.0;
 350           ray->h.azim_rate  = uf_ma[36] / 64.0;
 351           ray->h.sweep_rate = ray->h.azim_rate * (60.0/360.0);
 352           missing_data      = uf_ma[44];
 353
 354           if (pulled_time_from_first_ray == 0) {
 355             radar->h.height = uf_ma[24];
 356                 radar->h.latd = uf_ma[18];
 357                 radar->h.latm = uf_ma[19];
 358                 radar->h.lats = uf_ma[20] / 64.0;
 359                 radar->h.lond = uf_ma[21];
 360                 radar->h.lonm = uf_ma[22];
 361                 radar->h.lons = uf_ma[23] / 64.0;
 362                 radar->h.year  = ray->h.year;
 363                 radar->h.month = ray->h.month;
 364                 radar->h.day   = ray->h.day;
 365                 radar->h.hour  = ray->h.hour;
 366                 radar->h.minute = ray->h.minute;
 367                 radar->h.sec    = ray->h.sec;
 368                 strcpy(radar->h.radar_type, "uf");
 369                 pulled_time_from_first_ray = 1;
 370           }
 371           /*
 372            * ---- End of MANDATORY HEADER BLOCK.
 373            */
 374
 375           /* ---- Optional header used for MCTEX files. */
 376                 /* If this is a MCTEX file, the first 4 words following the
 377                          mandatory header contain the string 'MCTEX'. */
 378           memcpy(proj_name, (short *)(uf + uf_ma[2] - 1), 8);
 379           if (little_endian()) swap2(proj_name, 4);
 380
 381
 382           /* ---- Local Use Header (if present) was checked during Mandatory
 383            *      Header processing above.
 384            */
 385
 386           /* ---- Begining of FIELD HEADER. */
 387           uf_fh = uf+current_fh_index - 1;
 388           scale_factor      = uf_fh[1];
 389           ray->h.range_bin1 = uf_fh[2] * 1000.0 + uf_fh[3];
 390           ray->h.gate_size  = uf_fh[4];
 391
 392           ray->h.nbins      = uf_fh[5];
 393           ray->h.pulse_width  = uf_fh[6]/(RSL_SPEED_OF_LIGHT/1.0e6);
 394
 395                 if (strncmp((char *)proj_name, "MCTEX", 5) == 0)  /* MCTEX? */
 396                 {
 397                         /* The beamwidth values are not correct in Mctex UF files. */
 398                         ray->h.beam_width = 1.0;
 399                         sweep->h.beam_width = ray->h.beam_width;
 400                         sweep->h.horz_half_bw = ray->h.beam_width/2.0;
 401                         sweep->h.vert_half_bw = ray->h.beam_width/2.0;
 402                 }
 403                 else  /* Not MCTEX */
 404                 {
 405                         ray->h.beam_width = uf_fh[7] / 64.0;
 406                         sweep->h.beam_width = uf_fh[7]  / 64.0;
 407                         sweep->h.horz_half_bw = uf_fh[7] / 128.0; /* DFF 4/4/95 */
 408                         sweep->h.vert_half_bw = uf_fh[8] / 128.0; /* DFF 4/4/95 */
 409                 }
 410                 /*              fprintf (stderr, "uf_fh[7] = %d, [8] = %d\n", (int)uf_fh[7], (int)uf_fh[8]); */
 411                 if((int)uf_fh[7] == -32768) {
 412                         ray->h.beam_width     = 1;
 413                         sweep->h.beam_width   = 1;
 414                         sweep->h.horz_half_bw = .5;
 415                         sweep->h.vert_half_bw = .5;
 416                 }
 417
 418           ray->h.frequency    = uf_fh[9]  / 64.0;
 419           ray->h.wavelength   = uf_fh[11] / 64.0 / 100.0;  /* cm to m. */
 420           ray->h.pulse_count  = uf_fh[12];
 421           if (ifield == DZ_INDEX || ifield == ZT_INDEX) {
 422             radar->v[ifield]->h.calibr_const  = uf_fh[16] / 100.0;
 423                                                     /* uf value scaled by 100 */
 424           }
 425           else {
 426             radar->v[ifield]->h.calibr_const  = 0.0;
 427           }
 428           if (uf_fh[17] == (short)UF_NO_DATA) x = 0;
 429           else x = uf_fh[17] / 1000000.0;  /* PRT in seconds. */
 430           if (x != 0) {
 431             ray->h.prf = 1/x;
 432             ray->h.unam_rng = RSL_SPEED_OF_LIGHT / (2.0 * ray->h.prf * 1000.0);
 433           }
 434           else {
 435             ray->h.prf = 0.0;
 436             ray->h.unam_rng = 0.0;
 437           }
 438
 439           if (VR_INDEX == ifield || VE_INDEX == ifield) {
 440                 ray->h.nyq_vel = uf_fh[19] / scale_factor;
 441           }
 442
 443           /* ---- End of FIELD HEADER. */
 444
 445           ray->h.f = f;
 446           ray->h.invf = invf;
 447
 448           /* ---- Begining of FIELD DATA. */
 449           uf_data = uf+uf_fh[0] - 1;
 450
 451           len_data = ray->h.nbins;  /* Known because of RSL_new_ray. */
 452           for (m=0; m<len_data; m++) {
 453                 if (uf_data[m] == (short)UF_NO_DATA)
 454                   ray->range[m] = invf(BADVAL); /* BADVAL */
 455                 else {
 456                   if(uf_data[m] == missing_data)
 457                         ray->range[m] = invf(NOECHO); /* NOECHO */
 458                   else
 459                         ray->range[m] = invf((float)uf_data[m]/scale_factor);
 460                 }
 461           }
 462         }
 463   }
 464   return UF_MORE;
 465 }
 466
 467
 468 /*********************************************************************/
 469 /*                                                                   */
 470 /*                  swap_uf_buffer                                   */
 471 /*                                                                   */
 472 /*  By: John Merritt                                                 */
 473 /*      Space Applications Corporation                               */
 474 /*      October 4, 1994                                              */
 475 /*********************************************************************/
 476 void swap_uf_buffer(UF_buffer uf)
 477 {
 478   short *addr_end;
 479
 480   addr_end = uf + sizeof(UF_buffer)/sizeof(short);
 481   while (uf < addr_end)
 482         swap_2_bytes(uf++);
 483 }
 484
 485 enum UF_type {NOT_UF, TRUE_UF, TWO_BYTE_UF, FOUR_BYTE_UF};
 486
 487
 488 /*********************************************************************/
 489 /*                                                                   */
 490 /*                  RSL_uf_to_radar_fp                               */
 491 /*                                                                   */
 492 /*  By: John Merritt                                                 */
 493 /*      Space Applications Corporation                               */
 494 /*      September 22, 1995                                           */
 495 /*********************************************************************/
 496 Radar *RSL_uf_to_radar_fp(FILE *fp)
 497 {
 498   union {
 499         char buf[6];
 500         short sword;
 501         int word;
 502   } magic;
 503   Radar *radar;
 504   int nbytes;
 505   short sbytes;
 506   UF_buffer uf;
 507   enum UF_type uf_type;
 508 #define NEW_BUFSIZ 16384
 509
 510
 511   radar = NULL;
 512   setvbuf(fp,NULL,_IOFBF,(size_t)NEW_BUFSIZ); /* Faster i/o? */
 513   if (fread(magic.buf, sizeof(char), 6, fp) <= 0) return NULL;
 514 /*
 515  * Check for fortran record length delimeters, NCAR kludge.
 516  */
 517   if (strncmp("UF", magic.buf, 2) == 0) uf_type = TRUE_UF;
 518   else if (strncmp("UF", &magic.buf[2], 2) == 0) uf_type = TWO_BYTE_UF;
 519   else if (strncmp("UF", &magic.buf[4], 2) == 0) uf_type = FOUR_BYTE_UF;
 520   else uf_type = NOT_UF;
 521
 522   switch (uf_type) {
 523   case FOUR_BYTE_UF:
 524         if (radar_verbose_flag) fprintf(stderr,"UF file with 4 byte FORTRAN record delimeters.\n");
 525         /* Handle first record specially, since we needed magic information. */
 526         nbytes = magic.word;
 527         if (little_endian()) swap_4_bytes(&nbytes);
 528         memcpy(uf, &magic.buf[4], 2);
 529         (void)fread(&uf[1], sizeof(char), nbytes-2, fp);
 530         if (little_endian()) swap_uf_buffer(uf);
 531         (void)fread(&nbytes, sizeof(int), 1, fp);
 532         if (uf_into_radar(uf, &radar) == UF_DONE) break;
 533         /* Now the rest of the file. */
 534         while(fread(&nbytes, sizeof(int), 1, fp) > 0) {
 535           if (little_endian()) swap_4_bytes(&nbytes);
 536
 537           (void)fread(uf, sizeof(char), nbytes, fp);
 538           if (little_endian()) swap_uf_buffer(uf);
 539
 540           (void)fread(&nbytes, sizeof(int), 1, fp);
 541
 542           if (uf_into_radar(uf, &radar) == UF_DONE) break;
 543         }
 544         break;
 545
 546   case TWO_BYTE_UF:
 547         if (radar_verbose_flag) fprintf(stderr,"UF file with 2 byte FORTRAN record delimeters.\n");
 548         /* Handle first record specially, since we needed magic information. */
 549         sbytes = magic.sword;
 550         if (little_endian()) swap_2_bytes(&sbytes);
 551         memcpy(uf, &magic.buf[2], 4);
 552         (void)fread(&uf[2], sizeof(char), sbytes-4, fp);
 553         if (little_endian()) swap_uf_buffer(uf);
 554         (void)fread(&sbytes, sizeof(short), 1, fp);
 555         uf_into_radar(uf, &radar);
 556         /* Now the rest of the file. */
 557         while(fread(&sbytes, sizeof(short), 1, fp) > 0) {
 558           if (little_endian()) swap_2_bytes(&sbytes);
 559
 560           (void)fread(uf, sizeof(char), sbytes, fp);
 561           if (little_endian()) swap_uf_buffer(uf);
 562
 563           (void)fread(&sbytes, sizeof(short), 1, fp);
 564
 565           if (uf_into_radar(uf, &radar) == UF_DONE) break;
 566         }
 567         break;
 568
 569   case TRUE_UF:
 570         if (radar_verbose_flag) fprintf(stderr,"UF file with no FORTRAN record delimeters.  Good.\n");
 571         /* Handle first record specially, since we needed magic information. */
 572         memcpy(&sbytes, &magic.buf[2], 2); /* Record length is in word #2. */
 573         if (little_endian()) swap_2_bytes(&sbytes); /* # of 2 byte words. */
 574
 575         memcpy(uf, &magic.buf[0], 6);
 576         (void)fread(&uf[3], sizeof(short), sbytes-3, fp);
 577         if (little_endian()) swap_uf_buffer(uf);
 578         uf_into_radar(uf, &radar);
 579         /* Now the rest of the file. */
 580         while(fread(uf, sizeof(short), 2, fp) > 0) {
 581           memcpy(&sbytes, &uf[1], 2);  /* Record length is in word #2. */
 582           if (little_endian()) swap_2_bytes(&sbytes);
 583
 584           (void)fread(&uf[2], sizeof(short), sbytes-2, fp);  /* Have words 1,2. */
 585           if (little_endian()) swap_uf_buffer(uf);
 586
 587           if (uf_into_radar(uf, &radar) == UF_DONE) break;
 588         }
 589         break;
 590
 591   case NOT_UF: return NULL; break;
 592   }
 593   radar = reset_nsweeps_in_all_volumes(radar);
 594   radar = RSL_prune_radar(radar);
 595
 596   return radar;
 597 }
 598
 599
 600 /*********************************************************************/
 601 /*                                                                   */
 602 /*                  RSL_uf_to_radar                                  */
 603 /*                                                                   */
 604 /*  By: John Merritt                                                 */
 605 /*      Space Applications Corporation                               */
 606 /*      September 22, 1995                                           */
 607 /*********************************************************************/
 608 Radar *RSL_uf_to_radar(char *infile)
 609 {
 610 /*
 611  * This routine ingests a UF file and fills the Radar structure.
 612  * This routine allocates space via the system routine malloc.
 613  *
 614  * If *infile is NULL, read from stdin.
 615  */
 616   FILE *fp;
 617   Radar *radar;
 618
 619   radar = NULL;
 620   if (infile == NULL) {
 621         int save_fd;
 622         save_fd = dup(0);
 623         fp = fdopen(save_fd, "r");
 624   }  else if ((fp = fopen(infile, "r")) == NULL) {
 625         perror(infile);
 626         return radar;
 627   }
 628   fp = uncompress_pipe(fp); /* Transparently gunzip. */
 629   radar = RSL_uf_to_radar_fp(fp);
 630   rsl_pclose(fp);
 631
 632   return radar;
 633 }