Pileus Git - ~andy/rsl/blob - src/nsig_to_radar.c

   1 /*
   2     NASA/TRMM, Code 910.1.
   3     This is the TRMM Office Radar Software Library.
   4     Copyright (C) 1996  Paul A. Kucera of Applied Research Corporation,
   5                         Landover, Maryland, a NASA/GSFC on-site contractor.
   6
   7     This library is free software; you can redistribute it and/or
   8     modify it under the terms of the GNU Library General Public
   9     License as published by the Free Software Foundation; either
  10     version 2 of the License, or (at your option) any later version.
  11
  12     This library is distributed in the hope that it will be useful,
  13     but WITHOUT ANY WARRANTY; without even the implied warranty of
  14     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15     Library General Public License for more details.
  16
  17     You should have received a copy of the GNU Library General Public
  18     License along with this library; if not, write to the Free
  19     Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  20         */
  21 /*************************************************************/
  22 /*                                                           */
  23 /*    Function: nsig_to_radar.c                              */
  24 /*                                                           */
  25 /*    Paul A. Kucera                                         */
  26 /*    Applied Research Corporation                           */
  27 /*    NASA/GSFC                                              */
  28 /*    TRMM/Code 910.1                                        */
  29 /*                                                           */
  30 /*    Modifications by:                                      */
  31 /*    John H. Merritt                                        */
  32 /*    Space Applications Corporation                         */
  33 /*    NASA/GSFC                                              */
  34 /*    TRMM/Code 910.1                                        */
  35 /*                                                           */
  36 /*    Started: 08 AUG 96                                     */
  37 /*                                                           */
  38 /*    Derived from Paul Kucera's nsig_to_radar.c             */
  39 /*    Copyright 1996                                         */
  40 /*************************************************************/
  41
  42 #include<stdio.h>
  43 #include<unistd.h>
  44 #include<sys/types.h>
  45 #include<sys/stat.h>
  46 #include<fcntl.h>
  47 #include<stdlib.h>
  48 #include<math.h>
  49 #include<string.h>
  50
  51 #include"nsig.h"
  52 #include"rsl.h"
  53
  54 extern int radar_verbose_flag;
  55 extern int rsl_qfield[]; /* See RSL_select_fields */
  56
  57    /*  We need this entry for various things esp in Ray_header  */
  58 #define MISSING_HEADER_DATA -9999
  59    /*  The following is speed of light _in_ _air_ !  */
  60 #define SPEED_OF_LIGHT 299702547
  61 #define MIT_BEAMWIDTH      1.65
  62 #define TOG_BEAMWIDTH      1.65
  63 #define KWA_BEAMWIDTH      1.0
  64 #define DEFAULT_BEAMWIDTH  1.0
  65 #define NSIG_NO_DATA       -1
  66 #define MAX_NSIG_SWEEPS    30
  67 #define MAX_NSIG_RAYS      400
  68 #define NSIG_NO_ECHO       -32.0
  69 #define NSIG_NO_ECHO2     -999.0
  70
  71 static float (*f)(Range x);
  72 static Range (*invf)(float x);
  73
  74 extern FILE *file;
  75
  76 void  get_extended_header_info(NSIG_Sweep **nsig_sweep, int xh_size, int iray,
  77                                int nparams,
  78                                int *msec, float *azm, float *elev,
  79                                float *pitch, float *roll, float *heading,
  80                                float *azm_rate, float *elev_rate,
  81                                float *pitch_rate, float *roll_rate,
  82                                float *heading_rate,
  83                                float *lat, float *lon, int *alt, float *rvc,
  84                                float *vel_east, float *vel_north, float *vel_up)
  85 {
  86   static NSIG_Ext_header_ver1 xh;
  87   int data_type, itype;
  88
  89   *msec = *azm = *elev = *pitch = *roll = *heading =
  90         *azm_rate = *elev_rate = *pitch_rate = *roll_rate = *heading_rate =
  91         *lat = *lon = *alt = *rvc = 0;
  92
  93   /* Determine where 'itype' for extended header is. */
  94   for (itype = 0; itype<nparams; itype++) {
  95         data_type = NSIG_I2(nsig_sweep[itype]->idh.data_type);
  96     if (data_type == NSIG_DTB_EXH) break;
  97   }
  98   /*  printf("...extended header itype=%d, nparams=%d\n", itype, nparams); */
  99   if (itype == nparams) return;  /* No extended header. */
 100
 101   /* Version 1. */
 102   if (nsig_sweep[itype]->ray[iray] == NULL) return;
 103   if (nsig_sweep[itype]->ray[iray]->range == NULL) return;
 104   memmove(&xh, nsig_sweep[itype]->ray[iray]->range, sizeof(xh));
 105   *msec = NSIG_I4(xh.msec);
 106   /*   printf("...extended header msec= %d\n", *msec); */
 107   if (xh_size <= 20) /* Stop, only have version 0. */
 108     return;
 109
 110   /* Version 1 processing. */
 111   *azm = nsig_from_bang(xh.azm);
 112   *elev = nsig_from_bang(xh.elev);
 113   *pitch = nsig_from_bang(xh.pitch);
 114   *roll   = nsig_from_bang(xh.roll);
 115   *heading = nsig_from_bang(xh.heading);
 116   *azm_rate = nsig_from_bang(xh.azm_rate);
 117   *elev_rate = nsig_from_bang(xh.elev_rate);
 118   *pitch_rate = nsig_from_bang(xh.pitch_rate);
 119   *roll_rate   = nsig_from_bang(xh.roll_rate);
 120 #ifdef NSIG_VER2
 121
 122 #else
 123   *heading_rate = nsig_from_bang(xh.heading_rate);
 124 #endif
 125   *lat = nsig_from_fourb_ang(xh.lat);
 126   *lon = nsig_from_fourb_ang(xh.lon);
 127   if(*lat > 180.0) *lat -= 360.0;
 128   if(*lon > 180.0) *lon -= 360.0;
 129   *alt       = NSIG_I2(xh.alt);
 130   *rvc       = NSIG_I2(xh.rad_vel_cor)/100.0; /* cm to m */
 131   *vel_east  = NSIG_I2(xh.vel_e)/100.0; /* cm to m */
 132   *vel_north = NSIG_I2(xh.vel_n)/100.0; /* cm to m */
 133   *vel_up    = NSIG_I2(xh.vel_u)/100.0; /* cm to m */
 134   return;
 135 }
 136
 137 /** Main code **/
 138 Radar *
 139 #ifdef NSIG_VER2
 140 RSL_nsig2_to_radar
 141 #else
 142 RSL_nsig_to_radar
 143 #endif
 144 (char *filename)
 145 {
 146   FILE *fp;
 147   /* RSL structures */
 148   Radar                    *radar;
 149   Ray                      *ray;
 150
 151   int i, j, k, n;
 152   int year, month, day;
 153   int hour, minute, sec;
 154   int numbins, numsweep;
 155   int num_rays, sea_lvl_hgt;
 156   int radar_number, num_samples;
 157   int latd, latm, lats, lond, lonm, lons;
 158   int data_type;
 159   int bin_num;
 160   int sweep_year, sweep_day, sweep_month;
 161   int sweep_hour, sweep_minute, sweep_second;
 162   int sweep_sec;
 163   int z_flag_unc, z_flag_cor, v_flag, w_flag, speckle;
 164   int ant_scan_mode;
 165   float second;
 166   float pw;
 167   float bin_space;
 168   float prf, prf2, wave, beam_width;
 169   int prf_mode;
 170   float prf_modes[] = {1.0, 2.0/3.0, 3.0/4.0, 4.0/5.0};
 171   float vert_half_bw, horz_half_bw;
 172   float rng_last_bin;
 173   float rng_first_bin, freq;
 174   float max_vel, sweep_rate, azim_rate;
 175   float ray_data;
 176   float az1, az2;
 177   float elev1, elev2;
 178   double tmp;
 179   float sqi, log, csr, sig, cal_dbz;
 180   char radar_type[50], state[4], city[15];
 181   char site_name[16];
 182   NSIG_Product_file *prod_file;
 183   short id;
 184   int data_mask, nrays;
 185 #ifdef NSIG_VER2
 186   int masks[5];
 187 #endif
 188   int nparams, nsweeps;
 189   NSIG_Sweep **nsig_sweep;
 190   NSIG_Ray *ray_p;
 191   int itype, ifield;
 192   unsigned short nsig_2byte;   /* New for 2-byte data types, Aug 2009 */
 193   twob nsig_twob;
 194   Sweep *sweep;
 195   int msec;
 196   const unsigned short low10bits = 0x3ff;
 197   float azm, elev, pitch, roll, heading, azm_rate, elev_rate,
 198     pitch_rate, roll_rate, heading_rate,
 199     lat, lon;
 200   float fix_angle;
 201   int alt;  /* Altitude */
 202   float rvc;  /* Radial correction velocity m/s */
 203   float vel_east, vel_north, vel_up; /* Platform velocity vectors m/sec */
 204   int xh_size;
 205   float incr;
 206   extern int *rsl_qsweep; /* See RSL_read_these_sweeps in volume.c */
 207   extern int rsl_qsweep_max;
 208   extern float rsl_kdp_wavelen;
 209
 210   radar = NULL;
 211   if (radar_verbose_flag)
 212     fprintf(stderr, "open file: %s\n", filename);
 213
 214   /** Opening nsig file **/
 215   if((fp = nsig_open(filename)) == NULL) return NULL;
 216
 217 #ifdef NSIG_VER2
 218   sprintf(radar_type, "nsig2");
 219   radar_number = 22;  /** Arbitrary number given to nsig2 data **/
 220 #else
 221   sprintf(radar_type, "nsig");
 222   radar_number = 21;  /* What are these suppose to be? */
 223 #endif
 224   sprintf(state,"NA");
 225   sprintf(city,"NA");
 226
 227   /* MAINLINE CODE */
 228
 229   prod_file = (NSIG_Product_file *)calloc(1, sizeof(NSIG_Product_file));
 230
 231   n = nsig_read_record(fp, (char *)&prod_file->rec1);
 232   nsig_endianess(&prod_file->rec1);
 233   if (radar_verbose_flag)
 234     fprintf(stderr, "Read %d bytes for rec1.\n", n);
 235
 236   id = NSIG_I2(prod_file->rec1.struct_head.id);
 237   if (radar_verbose_flag)
 238     fprintf(stderr, "ID = %d\n", (int)id);
 239   if (id != 7 && id != 27) { /* testing: Use 27 for Version 2 data */
 240     fprintf(stderr, "File is not a SIGMET version 1 nor version 2 raw product file.\n");
 241     return NULL;
 242   }
 243
 244   n = nsig_read_record(fp, (char *)&prod_file->rec2);
 245   if (radar_verbose_flag)
 246     fprintf(stderr, "Read %d bytes for rec2.\n", n);
 247
 248   /* Count the bits set in 'data_mask' to determine the number
 249    * of parameters present.
 250    */
 251   xh_size = NSIG_I2(prod_file->rec2.ingest_head.size_ext_ray_headers);
 252   nrays = NSIG_I2(prod_file->rec2.ingest_head.num_rays);
 253   if (radar_verbose_flag)
 254     fprintf(stderr, "Expecting %d rays in each sweep.\n", nrays);
 255 #ifdef NSIG_VER2
 256   memmove(&masks[0], prod_file->rec2.task_config.dsp_info.data_mask_cur.mask_word_0,
 257     sizeof(fourb));
 258   memmove(&masks[1], &prod_file->rec2.task_config.dsp_info.data_mask_cur.mask_word_1,
 259     4*sizeof(fourb));
 260   nparams = 0;
 261   for (j=0; j < 5; j++) {
 262     data_mask = masks[j];
 263     for (i=0; i<32; i++)
 264       nparams += (data_mask >> i) & 0x1;
 265   }
 266 #else
 267   memmove(&data_mask, prod_file->rec2.task_config.dsp_info.data_mask, sizeof(fourb));
 268   for (nparams=i=0; i<32; i++)
 269     nparams += (data_mask >> i) & 0x1;
 270 #endif
 271
 272   /* Number of sweeps */
 273   nsweeps = NSIG_I2(prod_file->rec2.task_config.scan_info.num_swp);
 274
 275
 276
 277    memmove(site_name, prod_file->rec1.prod_end.site_name, sizeof(prod_file->rec1.prod_end.site_name));
 278    site_name[sizeof(site_name)-1] = '\0';
 279   if (radar_verbose_flag) {
 280     fprintf(stderr, "nparams = %d, nsweeps = %d\n", nparams, nsweeps);
 281     fprintf(stderr, "Site name = <%s>\n", site_name);
 282   }
 283
 284     /* nsig_sweep = nsig_read_sweep(fp, prod_file)
 285      *
 286      * Use: nsig_sweep[i]->ray[j]->range
 287      *
 288      * where 'range' is [0..nbins-1]
 289      */
 290
 291     /*
 292      * All the information you need is in:
 293      *    prod_file->rec1
 294      *        .struct_head, .prod_config .prod_end
 295      *    prod_file->rec2
 296      *        .struct_head, .ingest_head, .task_config .device_stat,
 297      *        .dsp1, .dsp2
 298      *    nsig_sweep[0..nparams-1]  'nparams' is the true number
 299      *                              of parameters present.  You
 300      *                              must check the 'id' (or type)
 301      *                              to determine the field type.
 302      *                              So far seen, nparams <= 6.
 303      *    nsig_sweep[i]->bhdr     <NSIG_Raw_prod_bhdr>
 304      *    nsig_sweep[i]->idh      <NSIG_Ingest_data_header>
 305      *    nsig_sweep[i]->ray[j]   <NSIG_Ray *>
 306      *
 307      * Note:
 308      *    For extended header access, you'll typically use nsig_sweep[0]
 309      *    (double check the id) and the ray data allocated (nsig_ray->range)
 310      *    is a pointer to the extended header, either v0 or v1.
 311      *    You can typecast the pointer to NSIG_Ext_header_ver0 or
 312      *    NSIG_Ext_header_ver1, as you like.  To determine which
 313      *    version of the extended headers you have use:
 314      *      xh_size <= 20 for version 0, else version 1.
 315      *    Access:
 316      *      xh_size = NSIG_I2(prod_file->rec2.ingest_head.size_ext_ray_headers)
 317      *
 318      * Functions:
 319      *    NSIG_I2(nsig_sweep[i]->idh.num_rays_act);   -- # of rays. (j)
 320      *    NSIG_I2(nsig_sweep[i]->ray[j]->h.num_bins); -- # of bins in a ray.
 321      *
 322      *    NSIG_I2(x), NSIG_I4(x)   - Convert data, x, to floating point.
 323      *
 324      *    IMPORTANT NOTE: It must be known whether or not to perform
 325      *                    byte-swapping.  To determine this, call
 326      *                    'nsig_endianess'.  It returns 0 for no-swapping
 327      *                    and 1 for swapping.  Additionally, it transparently
 328      *                    initializes the nsig library to automatically
 329      *                    swap when using NSIG_I2 or NSIG_I4.
 330      *                    The function 'nsig_read_sweep' automatically
 331      *                    calls 'nsig_endianess', too.
 332      */
 333
 334    sea_lvl_hgt = NSIG_I2(prod_file->rec1.prod_end.grnd_sea_ht);
 335
 336    if (radar_verbose_flag)
 337      fprintf(stderr, "sea: %d\n", sea_lvl_hgt);
 338    if (radar_verbose_flag)
 339      fprintf(stderr, "site_name: %s", site_name);
 340
 341    /** Determine beamwidth from input variables (not saved in nsig file) **/
 342    if(strncmp(site_name,"mit",3) == 0 || strncmp(site_name,"MIT",3) == 0)
 343      beam_width = MIT_BEAMWIDTH;
 344    else if(strncmp(site_name,"tog",3) == 0 || strncmp(site_name,"TOG",3) == 0)
 345      beam_width = TOG_BEAMWIDTH;
 346    else if(strncmp(site_name,"kwa",3) == 0 || strncmp(site_name,"KWA",3) == 0)
 347      beam_width = KWA_BEAMWIDTH;
 348    else
 349      beam_width = DEFAULT_BEAMWIDTH;
 350
 351    if (radar_verbose_flag)
 352      fprintf(stderr, "beamwidth: %f\n", beam_width);
 353
 354    vert_half_bw = beam_width/2.0;
 355    horz_half_bw = beam_width/2.0;
 356
 357    /** Reading date and time **/
 358    month = NSIG_I2(prod_file->rec2.ingest_head.start_time.month);
 359    year = NSIG_I2(prod_file->rec2.ingest_head.start_time.year);
 360    day = NSIG_I2(prod_file->rec2.ingest_head.start_time.day);
 361    sec = NSIG_I4(prod_file->rec2.ingest_head.start_time.sec);
 362
 363    /* converting seconds since mid to time of day */
 364    tmp = sec/3600.0;
 365    hour = (int)tmp;
 366    tmp = (tmp - hour) * 60.0;
 367    minute = (int)tmp;
 368    second = (tmp - minute) * 60.0;
 369
 370    /** records of the nsig file.                             **/
 371    num_rays = 0;
 372    pw = (NSIG_I4(prod_file->rec1.prod_end.pulse_wd))/100.0; /* pulse width */
 373    prf = NSIG_I4(prod_file->rec1.prod_end.prf);  /* pulse repetition frequency */
 374    prf_mode = NSIG_I2(prod_file->rec2.task_config.dsp_info.prf_mode);
 375    prf2 = prf * prf_modes[prf_mode];
 376    wave = (NSIG_I4(prod_file->rec1.prod_end.wavelen))/100.0; /* wavelength (cm) */
 377    rsl_kdp_wavelen = wave;  /* EXTERNAL (volume.c) This sets KD_F and KD_INVF
 378                              * to operate with the proper wavelength.
 379                              */
 380    numbins = NSIG_I4(prod_file->rec1.prod_end.num_bin);   /* # bins in ray */
 381    rng_first_bin = (float)NSIG_I4(prod_file->rec1.prod_end.rng_f_bin)/100.0;
 382    rng_last_bin = (float)NSIG_I4(prod_file->rec1.prod_end.rng_l_bin)/100.0;
 383    bin_space = ((rng_last_bin-rng_first_bin)/numbins); /*rng res (m)*/
 384
 385    numsweep = NSIG_I2(prod_file->rec2.task_config.scan_info.num_swp); /* # sweeps in volume */
 386    num_samples = NSIG_I2(prod_file->rec1.prod_end.num_samp);
 387    sweep_rate = 3.0; /** Approximate value -- info not stored **/
 388    azim_rate = sweep_rate*360.0/60.0;
 389    if (prf_mode != 0)
 390    {
 391         float max_vel1 = wave*prf/(100.0*4.0);
 392         float max_vel2 = wave*prf2/(100.0*4.0);
 393
 394         max_vel = (max_vel1 * max_vel2)/(max_vel1-max_vel2);
 395    }
 396    else
 397    {
 398         max_vel = wave*prf/(100.0*4.0);
 399    }
 400
 401    freq = (299793000.0/wave)*1.0e-4; /** freq in MHZ **/
 402
 403    sqi = NSIG_I2(prod_file->rec2.task_config.calib_info.sqi)/256.0;
 404    log = NSIG_I2(prod_file->rec2.task_config.calib_info.noise)/16.0;
 405    csr = NSIG_I2(prod_file->rec2.task_config.calib_info.clutr_corr)/(-16.0);
 406    sig = NSIG_I2(prod_file->rec2.task_config.calib_info.power)/16.0;
 407    cal_dbz = NSIG_I2(prod_file->rec2.task_config.calib_info.cal_ref)/16.0;
 408    z_flag_unc = NSIG_I2(prod_file->rec2.task_config.calib_info.z_flag_unc);
 409    z_flag_cor = NSIG_I2(prod_file->rec2.task_config.calib_info.z_flag_cor);
 410    v_flag = NSIG_I2(prod_file->rec2.task_config.calib_info.v_flag);
 411    w_flag = NSIG_I2(prod_file->rec2.task_config.calib_info.w_flag);
 412    speckle = NSIG_I2(prod_file->rec2.task_config.calib_info.speckle);
 413
 414    /** Verbose calibration information **/
 415    if (radar_verbose_flag)
 416       {
 417       fprintf(stderr, "LOG = %5.2f\n", log);
 418       fprintf(stderr, "SQI = %5.2f\n", sqi);
 419       fprintf(stderr, "CSR = %5.2f\n", csr);
 420       fprintf(stderr, "SIG = %5.2f\n", sig);
 421       fprintf(stderr, "Calibration reflectivity: %5.2f dBZ\n", cal_dbz);
 422       fprintf(stderr, "ZT flags: %d\n", z_flag_unc);  /** can find these **/
 423       fprintf(stderr, "DZ flags: %d\n", z_flag_cor);  /** defn in the    **/
 424       fprintf(stderr, "VR flags: %d\n", v_flag);      /** SIGMET Doc     **/
 425       fprintf(stderr, "SW flags: %d\n", w_flag);
 426       fprintf(stderr, "Flags: -3856  = SQI thresholding\n");
 427       fprintf(stderr, "       -21846 = LOG thresholding\n");
 428       fprintf(stderr, "       -24416 = LOG & SQI thresholding\n");
 429       fprintf(stderr, "       -24516 = LOG & SQI & SIG thresholding\n");
 430       fprintf(stderr, "speckle remover: %d\n", speckle);
 431       }
 432
 433    if (radar_verbose_flag)
 434      fprintf(stderr, "vel: %f prf: %f prf2: %f\n", max_vel, prf, prf2);
 435
 436    /** Extracting Latitude and Longitude from nsig file **/
 437    lat = nsig_from_fourb_ang(prod_file->rec2.ingest_head.lat_rad);
 438    lon = nsig_from_fourb_ang(prod_file->rec2.ingest_head.lon_rad);
 439    if(lat > 180.0) lat -= 360.0;
 440    if(lon > 180.0) lon -= 360.0;
 441    if (radar_verbose_flag)
 442      fprintf(stderr, "nsig_to_radar: lat %f, lon %f\n", lat, lon);
 443    /** Latitude deg, min, sec **/
 444    latd = (int)lat;
 445    tmp = (lat - latd) * 60.0;
 446    latm = (int)tmp;
 447    lats = (int)((tmp - latm) * 60.0);
 448    /** Longitude deg, min, sec **/
 449    lond = (int)lon;
 450    tmp = (lon - lond) * 60.0;
 451    lonm = (int)tmp;
 452    lons = (int)((tmp - lonm) * 60.0);
 453
 454    /** Allocating memory for radar structure **/
 455    radar = RSL_new_radar(MAX_RADAR_VOLUMES);
 456    if (radar == NULL)
 457       {
 458       fprintf(stderr, "nsig_to_radar: radar is NULL\n");
 459       return NULL;
 460       }
 461
 462    /** Filling Radar Header **/
 463    radar->h.month = month;
 464    radar->h.day = day;
 465    radar->h.year = year; /* Year 2000 compliant. */
 466    radar->h.hour = hour;
 467    radar->h.minute = minute;
 468    radar->h.sec = second;
 469    sprintf(radar->h.radar_type, "%s", radar_type);
 470    radar->h.number = radar_number;
 471    memmove(radar->h.name, site_name, sizeof(radar->h.name));
 472    memmove(radar->h.radar_name, site_name, sizeof(radar->h.radar_name));
 473    memmove(radar->h.city, city, sizeof(radar->h.city));
 474    memmove(radar->h.state, state, sizeof(radar->h.state));
 475    radar->h.latd = latd;
 476    radar->h.latm = latm;
 477    radar->h.lats = lats;
 478    radar->h.lond = lond;
 479    radar->h.lonm = lonm;
 480    radar->h.lons = lons;
 481    radar->h.height = (int)sea_lvl_hgt;
 482    radar->h.spulse = (int)(pw*1000);
 483    radar->h.lpulse = (int)(pw*1000);
 484    ant_scan_mode = NSIG_I2(prod_file->rec2.task_config.scan_info.ant_scan_mode);
 485    if(ant_scan_mode == 2 || ant_scan_mode == 7) radar->h.scan_mode = RHI;
 486    else radar->h.scan_mode = PPI;
 487
 488    if (radar_verbose_flag) {
 489 #ifdef NSIG_VER2
 490      fprintf(stderr, "\nSIGMET version 2 raw product file.\n");
 491 #else
 492      fprintf(stderr, "\nSIGMET version 1 raw product file.\n");
 493 #endif
 494      fprintf(stderr, "Date: %2.2d/%2.2d/%4.4d %2.2d:%2.2d:%f\n",
 495              radar->h.month, radar->h.day, radar->h.year,
 496              radar->h.hour, radar->h.minute, radar->h.sec);
 497      fprintf(stderr, "Name: ");
 498      for (i=0; i<sizeof(radar->h.name); i++)
 499        fprintf(stderr, "%c", radar->h.name[i]);
 500      fprintf(stderr, "\n");
 501      fprintf(stderr, "Lat/lon (%d %d' %d'', %d %d' %d'')\n",
 502              radar->h.latd, radar->h.latm, radar->h.lats,
 503              radar->h.lond, radar->h.lonm, radar->h.lons);
 504    }
 505
 506    /** Converting data **/
 507    if (radar_verbose_flag) fprintf(stderr, "Expecting %d sweeps.\n", numsweep);
 508    for(i = 0; i < numsweep; i++)
 509       {
 510         nsig_sweep = nsig_read_sweep(fp, prod_file);
 511         if (nsig_sweep == NULL) { /* EOF possibility */
 512           if (feof(fp)) break;
 513           else continue;
 514         }
 515         if (rsl_qsweep != NULL) {
 516           if (i > rsl_qsweep_max) break;
 517           if (rsl_qsweep[i] == 0) continue;
 518         }
 519         if (radar_verbose_flag)
 520           fprintf(stderr, "Read sweep # %d\n", i);
 521     /* The whole sweep is 'nsig_sweep' ... pretty slick.
 522          *
 523          * nsig_sweep[itype]  -- [0..nparams], if non-null.
 524          */
 525     for (itype=0; itype<nparams; itype++) {
 526       if (nsig_sweep[itype] == NULL) continue;
 527
 528       /** Reading date and time **/
 529       sweep_month = NSIG_I2(nsig_sweep[itype]->idh.time.month);
 530       sweep_year = NSIG_I2(nsig_sweep[itype]->idh.time.year);
 531       sweep_day = NSIG_I2(nsig_sweep[itype]->idh.time.day);
 532       sweep_sec = NSIG_I4(nsig_sweep[itype]->idh.time.sec);
 533 #ifdef NSIG_VER2
 534       msec = NSIG_I2(nsig_sweep[itype]->idh.time.msec) & low10bits;
 535       /*      printf("....... msec == %d\n", msec); */
 536 #endif
 537       /* converting seconds since mid to time of day */
 538       tmp = sweep_sec/3600.0;
 539       sweep_hour = (int)tmp;
 540       tmp = (tmp - sweep_hour) * 60.0;
 541       sweep_minute = (int)tmp;
 542       sweep_second = sweep_sec - (sweep_hour*3600 + sweep_minute*60);
 543
 544       num_rays = NSIG_I2(nsig_sweep[itype]->idh.num_rays_exp);
 545
 546       data_type = NSIG_I2(nsig_sweep[itype]->idh.data_type);
 547
 548       ifield = 0;
 549       switch (data_type) {
 550       case NSIG_DTB_EXH:
 551           ifield = -1;
 552           break;
 553       case NSIG_DTB_UCR:
 554       case NSIG_DTB_UCR2:
 555         ifield = ZT_INDEX;
 556         f      = ZT_F;
 557         invf   = ZT_INVF;
 558         break;
 559       case NSIG_DTB_CR:
 560       case NSIG_DTB_CR2:
 561         ifield = DZ_INDEX;
 562         f      = DZ_F;
 563         invf   = DZ_INVF;
 564         break;
 565       case NSIG_DTB_VEL:
 566       case NSIG_DTB_VEL2:
 567         ifield = VR_INDEX;
 568         f      = VR_F;
 569         invf   = VR_INVF;
 570         break;
 571       case NSIG_DTB_WID:
 572       case NSIG_DTB_WID2:
 573         ifield = SW_INDEX;
 574         f      = SW_F;
 575         invf   = SW_INVF;
 576         break;
 577       case NSIG_DTB_ZDR:
 578       case NSIG_DTB_ZDR2:
 579         ifield = DR_INDEX;
 580         f      = DR_F;
 581         invf   = DR_INVF;
 582         break;
 583       case NSIG_DTB_KDP:
 584         ifield = KD_INDEX;
 585         f      = KD_F;
 586         invf   = KD_INVF;
 587         break;
 588       case NSIG_DTB_PHIDP:     /* SRB 990127 */
 589         ifield = PH_INDEX;
 590         f      = PH_F;
 591         invf   = PH_INVF;
 592         break;
 593       case NSIG_DTB_RHOHV:     /* SRB 000414 */
 594         ifield = RH_INDEX;
 595         f      = RH_F;
 596         invf   = RH_INVF;
 597         break;
 598       case NSIG_DTB_VELC:
 599       case NSIG_DTB_VELC2:
 600         ifield = VC_INDEX;
 601         f      = VC_F;
 602         invf   = VC_INVF;
 603         break;
 604       case NSIG_DTB_KDP2:
 605         ifield = KD_INDEX;
 606         f      = KD_F;
 607         invf   = KD_INVF;
 608         break;
 609       case NSIG_DTB_PHIDP2:
 610         ifield = PH_INDEX;
 611         f      = PH_F;
 612         invf   = PH_INVF;
 613         break;
 614       case NSIG_DTB_RHOHV2:
 615         ifield = RH_INDEX;
 616         f      = RH_F;
 617         invf   = RH_INVF;
 618         break;
 619       case NSIG_DTB_SQI:
 620       case NSIG_DTB_SQI2:
 621         ifield = SQ_INDEX;
 622         f      = SQ_F;
 623         invf   = SQ_INVF;
 624         break;
 625       case NSIG_DTB_HCLASS:
 626       case NSIG_DTB_HCLASS2:
 627         ifield = HC_INDEX;
 628         f      = HC_F;
 629         invf   = HC_INVF;
 630         break;
 631       case NSIG_DTB_DBZ2:
 632         ifield = CZ_INDEX;
 633         f      = CZ_F;
 634         invf   = CZ_INVF;
 635         break;
 636       case NSIG_DTB_ZDRC2:
 637         ifield = ZD_INDEX;
 638         f      = ZD_F;
 639         invf   = ZD_INVF;
 640         break;
 641       case NSIG_DTB_DBTE8:
 642         ifield = ET_INDEX;
 643         f      = ZT_F;
 644         invf   = ZT_INVF;
 645         break;
 646       case NSIG_DTB_DBZE8:
 647         ifield = EZ_INDEX;
 648         f      = DZ_F;
 649         invf   = DZ_INVF;
 650         break;
 651       default:
 652         fprintf(stderr,"Unknown field type: %d  Skipping it.\n", data_type);
 653         continue;
 654       }
 655
 656       if (radar_verbose_flag)
 657         fprintf(stderr, "     nsig_sweep[%d], data_type = %d, rays(expected) = %d, nrays(actual) = %d\n", itype, data_type, num_rays, NSIG_I2(nsig_sweep[itype]->idh.num_rays_act));
 658
 659       if (data_type != NSIG_DTB_EXH) {
 660         if ((radar->v[ifield] == NULL)) {
 661           if (rsl_qfield[ifield]) {
 662              radar->v[ifield] = RSL_new_volume(numsweep);
 663              radar->v[ifield]->h.f = f;
 664              radar->v[ifield]->h.invf = invf;
 665            } else {
 666              /* Skip this field, because, the user does not want it. */
 667              continue;
 668            }
 669         }
 670          if (radar->v[ifield]->sweep[i] == NULL)
 671            radar->v[ifield]->sweep[i] = RSL_new_sweep(num_rays);
 672          }
 673       else
 674       continue;    /* Skip the actual extended header processing.
 675                     * This is different than getting it, so that
 676                     * the information is available for the other
 677                     * fields when filling the RSL ray headers.
 678                     */
 679
 680       /** DATA conversion time **/
 681       sweep = radar->v[ifield]->sweep[i];
 682       sweep->h.f = f;
 683       sweep->h.invf = invf;
 684       sweep->h.sweep_num = i;
 685       sweep->h.beam_width = beam_width;
 686       sweep->h.vert_half_bw = vert_half_bw;
 687       sweep->h.horz_half_bw = horz_half_bw;
 688       fix_angle = nsig_from_bang(nsig_sweep[itype]->idh.fix_ang);
 689       if (radar->h.scan_mode == PPI) sweep->h.elev = fix_angle;
 690       else sweep->h.azimuth = fix_angle;
 691
 692       for(j = 0; j < num_rays; j++)
 693         {
 694           ray_p = nsig_sweep[itype]->ray[j];
 695           if (ray_p == NULL) continue;
 696           bin_num = NSIG_I2(ray_p->h.num_bins);
 697
 698           /* Load extended header information, if available.
 699            * We need to pass the entire nsig_sweep and search for
 700            * the extended header field (it may not be data_type==0).
 701            */
 702           get_extended_header_info(nsig_sweep, xh_size, j, nparams,
 703                        &msec, &azm, &elev,
 704                        &pitch, &roll, &heading,
 705                        &azm_rate, &elev_rate,
 706                        &pitch_rate, &roll_rate, &heading_rate,
 707                        &lat, &lon, &alt, &rvc,
 708                        &vel_east, &vel_north, &vel_up);
 709
 710
 711           if (radar->v[ifield]->sweep[i]->ray[j] == NULL)
 712             radar->v[ifield]->sweep[i]->ray[j] = RSL_new_ray(bin_num);
 713           ray = radar->v[ifield]->sweep[i]->ray[j];
 714           ray->h.f = f;
 715           ray->h.invf = invf;
 716           /** Ray is at nsig_sweep[itype].ray->... **/
 717           /** Loading nsig data into data structure **/
 718
 719           ray->h.month  = sweep_month;
 720           ray->h.day    = sweep_day;
 721           ray->h.year   = sweep_year; /* Year 2000 compliant. */
 722           ray->h.hour   = sweep_hour;
 723           ray->h.minute = sweep_minute;
 724           if (msec == 0) { /* No extended header */
 725             ray->h.sec  = NSIG_I2(ray_p->h.sec) + sweep_second;
 726             elev = sweep->h.elev;
 727           } else
 728             ray->h.sec  = sweep_second + msec/1000.0;
 729
 730           /* add time ... handles end of min,hour,month,year and century. */
 731           if (ray->h.sec >= 60) /* Should I fix the time no matter what? */
 732             RSL_fix_time(ray);  /* Repair second overflow. */
 733
 734           ray->h.ray_num    = j;
 735           ray->h.elev_num   = i;
 736           ray->h.range_bin1 = (int)rng_first_bin;
 737           ray->h.gate_size  = (int)(bin_space+.5); /* Nearest int */
 738           ray->h.vel_res    = bin_space;
 739           ray->h.sweep_rate = sweep_rate;
 740           ray->h.prf        = (int)prf;
 741             if (prf != 0)
 742               ray->h.unam_rng = 299793000.0 / (2.0 * prf * 1000.0);  /* km */
 743             else
 744               ray->h.unam_rng = 0.0;
 745           ray->h.prf2 = (int) prf2;
 746           ray->h.fix_angle = fix_angle;
 747           ray->h.azim_rate  = azim_rate;
 748           ray->h.pulse_count = num_samples;
 749           ray->h.pulse_width = pw;
 750           ray->h.beam_width  = beam_width;
 751           ray->h.frequency   = freq / 1000.0;  /* GHz */
 752           ray->h.wavelength  = wave/100.0;     /* meters */
 753           ray->h.nyq_vel     = max_vel;        /* m/s */
 754           if (elev == 0.) elev = sweep->h.elev;
 755
 756           /* Compute mean azimuth angle for ray. */
 757           az1 = nsig_from_bang(ray_p->h.beg_azm);
 758           az2 = nsig_from_bang(ray_p->h.end_azm);
 759           /*          printf("az1, %f, az2 %f\n", az1, az2); */
 760           if(az1 > az2)
 761             if((az1 - az2) > 180.0) az2 += 360.0;
 762             else
 763               ;
 764           else
 765             if((az2 - az1) > 180.0) az1 += 360.0;
 766
 767           az1 = (az1 + az2) / 2.0;
 768           if (az1 > 360) az1 -= 360;
 769           ray->h.azimuth     = az1;
 770
 771           /* Compute mean elevation angle for ray. */
 772           elev1 = nsig_from_bang(ray_p->h.beg_elev);
 773           elev2 = nsig_from_bang(ray_p->h.end_elev);
 774           /*          printf("elev1, %f, elev2 %f\n", elev1, elev2); */
 775           if(elev1 > elev2)
 776             if((elev1 - elev2) > 180.0) elev2 += 360.0;
 777             else
 778               ;
 779           else
 780             if((elev2 - elev1) > 180.0) elev2 -= 360.0;
 781
 782           elev1 = (elev1 + elev2) / 2.0;
 783           if (elev1 > 360) elev1 -= 360;
 784           ray->h.elev = elev1;
 785
 786           /* From the extended header information, we learn the following. */
 787           ray->h.pitch        = pitch;
 788           ray->h.roll         = roll;
 789           ray->h.heading      = heading;
 790           ray->h.pitch_rate   = pitch_rate;
 791           ray->h.roll_rate    = roll_rate;
 792           ray->h.heading_rate = heading_rate;
 793           ray->h.lat          = lat;
 794           ray->h.lon          = lon;
 795           ray->h.alt          = alt;
 796           ray->h.rvc          = rvc;
 797           ray->h.vel_east     = vel_east;
 798           ray->h.vel_north    = vel_north;
 799           ray->h.vel_up       = vel_up;
 800
 801           /*          printf("Processing sweep[%d]->ray[%d]: %d %f %f %f %f %f %f %f %f %d nbins=%d, bin1=%d gate=%d\n",
 802                  i, j, msec, ray->h.sec, ray->h.azimuth, ray->h.elev, ray->h.pitch, ray->h.roll, ray->h.heading, ray->h.lat, ray->h.lon, ray->h.alt, ray->h.nbins, ray->h.range_bin1, ray->h.gate_size);
 803                  */
 804           /* TODO: ingest data header contains a value for bits-per-bin.
 805            * This might be of use to allocate an array for ray->range with
 806            * either 1-byte or 2-byte elements.  Then there's no need for
 807            * memmove() whenever we need 2 bytes.
 808            */
 809
 810           if (data_type == NSIG_DTB_EXH) continue;
 811           ray_data = 0;
 812           for(k = 0; k < bin_num; k++) {
 813             switch(data_type) {
 814             case NSIG_DTB_UCR:
 815             case NSIG_DTB_CR:
 816             case NSIG_DTB_DBTE8:
 817             case NSIG_DTB_DBZE8:
 818               if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO;
 819               else ray_data = (float)((ray_p->range[k]-64.0)/2.0);
 820               break;
 821             /* Simplified the velocity conversion for NSIG_DTB_VEL, using
 822              * formula from IRIS Programmer's Manual. BLK, Oct 9 2009.
 823              */
 824             case NSIG_DTB_VEL:
 825               if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO;
 826               else ray_data = (float)((ray_p->range[k]-128.0)/127.0)*max_vel;
 827               break;
 828
 829             case NSIG_DTB_WID:
 830               if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO;
 831               else ray_data =(float)((ray_p->range[k])/256.0)*max_vel;
 832               break;
 833
 834             case NSIG_DTB_ZDR:
 835               if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO;
 836               else ray_data = (float)((ray_p->range[k]-128.0)/16.0);
 837               break;
 838
 839             case NSIG_DTB_KDP:
 840                 if (ray_p->range[k] == 0 || ray_p->range[k] == 255 ||
 841                     rsl_kdp_wavelen == 0.0) {
 842                   ray_data = NSIG_NO_ECHO;
 843                   break;
 844                 }
 845                 if (ray_p->range[k] < 128)
 846                   ray_data = (-0.25 *
 847                     pow((double)600.0,(double)((127-ray_p->range[k])/126.0))) /
 848                       rsl_kdp_wavelen;
 849                 else if (ray_p->range[k] > 128)
 850                   ray_data = (0.25 *
 851                     pow((double)600.0,(double)((ray_p->range[k]-129)/126.0))) /
 852                       rsl_kdp_wavelen;
 853                 else
 854                   ray_data = 0.0;
 855                 break;
 856
 857             case NSIG_DTB_PHIDP:
 858               if (ray_p->range[k] == 0 || ray_p->range[k] == 255)
 859                 ray_data = NSIG_NO_ECHO;
 860               else
 861                 ray_data = 180.0*((ray_p->range[k]-1.0)/254.0);
 862               break;
 863
 864             case NSIG_DTB_RHOHV:
 865               if (ray_p->range[k] == 0 || ray_p->range[k] == 255)
 866                 ray_data = NSIG_NO_ECHO;
 867               else
 868                 ray_data = sqrt((double)((ray_p->range[k]-1.0)/253.0));
 869               break;
 870
 871             case NSIG_DTB_HCLASS:
 872               if (ray_p->range[k] == 0 || ray_p->range[k] == 255)
 873                 ray_data = NSIG_NO_ECHO;
 874               else
 875                 ray_data = ray_p->range[k];
 876               break;
 877
 878             case NSIG_DTB_SQI:
 879               if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO;
 880               else ray_data = (float)sqrt((ray_p->range[k]-1.0)/253.0);
 881               break;
 882
 883             case NSIG_DTB_VELC:
 884               if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO;
 885               else {
 886                 incr=75./127.;  /*  (+|- 75m/s) / 254 values */
 887                 ray_data = (float)(ray_p->range[k]-128)*incr;
 888               }
 889               break;
 890
 891             case NSIG_DTB_UCR2:
 892             case NSIG_DTB_CR2:
 893             case NSIG_DTB_VEL2:
 894             case NSIG_DTB_VELC2:
 895             case NSIG_DTB_ZDR2:
 896             case NSIG_DTB_KDP2:
 897               memmove(nsig_twob, &ray_p->range[2*k], 2);
 898               nsig_2byte = NSIG_I2(nsig_twob);
 899               if (nsig_2byte == 0 || nsig_2byte == 65535)
 900                 ray_data = NSIG_NO_ECHO2;
 901               else ray_data = (float)(nsig_2byte-32768)/100.;
 902               break;
 903
 904             case NSIG_DTB_WID2:
 905               memmove(nsig_twob, &ray_p->range[2*k], 2);
 906               nsig_2byte = NSIG_I2(nsig_twob);
 907               if (nsig_2byte == 0 || nsig_2byte == 65535)
 908                 ray_data = NSIG_NO_ECHO2;
 909               else ray_data = (float)nsig_2byte/100.;
 910               break;
 911
 912             case NSIG_DTB_PHIDP2:
 913               memmove(nsig_twob, &ray_p->range[2*k], 2);
 914               nsig_2byte = NSIG_I2(nsig_twob);
 915               if (nsig_2byte == 0 || nsig_2byte == 65535)
 916                 ray_data = NSIG_NO_ECHO;
 917               else
 918                 ray_data = 360.*(nsig_2byte-1)/65534.;
 919               break;
 920
 921             case NSIG_DTB_SQI2:
 922             case NSIG_DTB_RHOHV2:
 923               memmove(nsig_twob, &ray_p->range[2*k], 2);
 924               nsig_2byte = NSIG_I2(nsig_twob);
 925               if (nsig_2byte == 0 || nsig_2byte == 65535)
 926                 ray_data = NSIG_NO_ECHO2;
 927               else ray_data = (float)(nsig_2byte-1)/65533.;
 928               break;
 929
 930             case NSIG_DTB_HCLASS2:
 931               memmove(nsig_twob, &ray_p->range[2*k], 2);
 932               nsig_2byte = NSIG_I2(nsig_twob);
 933               if (nsig_2byte == 0 || nsig_2byte == 65535)
 934                 ray_data = NSIG_NO_ECHO2;
 935               else
 936                 ray_data = nsig_2byte;
 937             }
 938
 939             if (ray_data == NSIG_NO_ECHO || ray_data == NSIG_NO_ECHO2)
 940               ray->range[k] = ray->h.invf(BADVAL);
 941             else
 942               ray->range[k] = ray->h.invf(ray_data);
 943
 944             /*
 945             if (data_type == NSIG_DTB_KDP)
 946             printf("v[%d]->sweep[%d]->ray[%d]->range[%d] = %f, %d, %f\n",
 947                    ifield, i, j, k, ray->h.f(ray->range[k]),
 948                    (int)ray_p->range[k], ray_data);
 949             */
 950           }
 951         }
 952         }
 953         nsig_free_sweep(nsig_sweep);
 954       }
 955
 956    /* Do not reset radar->h.nvolumes. It is already set properly. */
 957    if (radar_verbose_flag)
 958      fprintf(stderr, "Max index of radar->v[0..%d]\n", radar->h.nvolumes);
 959
 960
 961    /** close nsig file **/
 962    nsig_close(fp);
 963
 964    radar = RSL_prune_radar(radar);
 965    /** return radar pointer **/
 966    return radar;
 967 }