Pileus Git - ~andy/rsl/blob - 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
  70 static float (*f)(Range x);
  71 static Range (*invf)(float x);
  72
  73 FILE *file;
  74
  75 void  get_extended_header_info(NSIG_Sweep **nsig_sweep, int xh_size, int iray,
  76                                int nparams,
  77                                int *msec, float *azm, float *elev,
  78                                float *pitch, float *roll, float *heading,
  79                                float *azm_rate, float *elev_rate,
  80                                float *pitch_rate, float *roll_rate, float *heading_rate,
  81                                float *lat, float *lon, int *alt, float *rvc,
  82                                float *vel_east, float *vel_north, float *vel_up)
  83 {
  84   static NSIG_Ext_header_ver1 xh;
  85   int data_type, itype;
  86
  87   *msec = *azm = *elev = *pitch = *roll = *heading =
  88         *azm_rate = *elev_rate = *pitch_rate = *roll_rate = *heading_rate =
  89         *lat = *lon = *alt = *rvc = 0;
  90
  91   /* Determine where 'itype' for extended header is. */
  92   for (itype = 0; itype<nparams; itype++) {
  93         data_type = NSIG_I2(nsig_sweep[itype]->idh.data_type);
  94         if (data_type == NSIG_DTB_EXH) break;
  95   }
  96   /*  printf("...extended header itype=%d, nparams=%d\n", itype, nparams); */
  97   if (itype == nparams) return;  /* No extended header. */
  98
  99   /* Version 1. */
 100   if (nsig_sweep[itype]->ray[iray] == NULL) return;
 101   if (nsig_sweep[itype]->ray[iray]->range == NULL) return;
 102   memmove(&xh, nsig_sweep[itype]->ray[iray]->range, sizeof(xh));
 103   *msec = NSIG_I4(xh.msec);
 104   /*   printf("...extended header msec= %d\n", *msec); */
 105   if (xh_size <= 20) /* Stop, only have version 0. */
 106         return;
 107
 108   /* Version 1 processing. */
 109   *azm = nsig_from_bang(xh.azm);
 110   *elev = nsig_from_bang(xh.elev);
 111   *pitch = nsig_from_bang(xh.pitch);
 112   *roll   = nsig_from_bang(xh.roll);
 113   *heading = nsig_from_bang(xh.heading);
 114   *azm_rate = nsig_from_bang(xh.azm_rate);
 115   *elev_rate = nsig_from_bang(xh.elev_rate);
 116   *pitch_rate = nsig_from_bang(xh.pitch_rate);
 117   *roll_rate   = nsig_from_bang(xh.roll_rate);
 118 #ifdef NSIG_VER2
 119
 120 #else
 121   *heading_rate = nsig_from_bang(xh.heading_rate);
 122 #endif
 123   *lat = nsig_from_fourb_ang(xh.lat);
 124   *lon = nsig_from_fourb_ang(xh.lon);
 125   if(*lat > 180.0) *lat -= 360.0;
 126   if(*lon > 180.0) *lon -= 360.0;
 127   *alt       = NSIG_I2(xh.alt);
 128   *rvc       = NSIG_I2(xh.rad_vel_cor)/100.0; /* cm to m */
 129   *vel_east  = NSIG_I2(xh.vel_e)/100.0; /* cm to m */
 130   *vel_north = NSIG_I2(xh.vel_n)/100.0; /* cm to m */
 131   *vel_up    = NSIG_I2(xh.vel_u)/100.0; /* cm to m */
 132   return;
 133 }
 134
 135 /** Main code **/
 136 Radar *
 137 #ifdef NSIG_VER2
 138 RSL_nsig2_to_radar
 139 #else
 140 RSL_nsig_to_radar
 141 #endif
 142 (char *filename)
 143 {
 144    FILE *fp;
 145    /* RSL structures */
 146    Radar                    *radar;
 147    Ray                      *ray;
 148
 149    int i, j, k, n;
 150    int year, month, day;
 151    int hour, minute, sec;
 152    int numbins, numsweep;
 153    int num_rays, sea_lvl_hgt;
 154    int radar_number, num_samples;
 155    int latd, latm, lats, lond, lonm, lons;
 156    int data_type;
 157    int bin_num;
 158    int sweep_year, sweep_day, sweep_month;
 159    int sweep_hour, sweep_minute, sweep_second;
 160    int sweep_sec;
 161    int z_flag_unc, z_flag_cor, v_flag, w_flag, speckle;
 162    int ant_scan_mode;
 163    float second;
 164    float pw;
 165    float bin_space;
 166    float prf, wave, beam_width;
 167    float vert_half_bw, horz_half_bw;
 168    float rng_last_bin;
 169    float rng_first_bin, freq;
 170    float max_vel, sweep_rate, azim_rate;
 171    float ray_data;
 172    float az1, az2;
 173    double tmp;
 174    float sqi, log, csr, sig, cal_dbz;
 175    char radar_type[50], state[2], city[15];
 176    char site_name[16];
 177   NSIG_Product_file *prod_file;
 178   short id;
 179   int data_mask, nrays;
 180   int nparams, nsweeps;
 181   NSIG_Sweep **nsig_sweep;
 182   NSIG_Ray *ray_p;
 183   int itype, ifield;
 184   Sweep *sweep;
 185   int msec;
 186   float azm, elev, pitch, roll, heading, azm_rate, elev_rate,
 187         pitch_rate, roll_rate, heading_rate,
 188         lat, lon;
 189   int alt;  /* Altitude */
 190   float rvc;  /* Radial correction velocity m/s */
 191   float vel_east, vel_north, vel_up; /* Platform velocity vectors m/sec */
 192   int xh_size;
 193   extern int *rsl_qsweep; /* See RSL_read_these_sweeps in volume.c */
 194   extern int rsl_qsweep_max;
 195   extern float rsl_kdp_wavelen;
 196
 197   radar = NULL;
 198   if (radar_verbose_flag)
 199         fprintf(stderr, "open file: %s\n", filename);
 200
 201   /** Opening nsig file **/
 202   if((fp = nsig_open(filename)) == NULL) return NULL;
 203
 204 #ifdef NSIG_VER2
 205   sprintf(radar_type, "nsig2");
 206   radar_number = 22;  /** Arbitrary number given to nsig2 data **/
 207 #else
 208   sprintf(radar_type, "nsig");
 209   radar_number = 21;  /* What are these suppose to be? */
 210 #endif
 211   sprintf(state,"NA");
 212   sprintf(city,"NA");
 213
 214   /* MAINLINE CODE */
 215
 216   prod_file = (NSIG_Product_file *)calloc(1, sizeof(NSIG_Product_file));
 217
 218   n = nsig_read_record(fp, (char *)&prod_file->rec1);
 219   nsig_endianess(&prod_file->rec1);
 220   if (radar_verbose_flag)
 221         fprintf(stderr, "Read %d bytes for rec1.\n", n);
 222
 223   id = NSIG_I2(prod_file->rec1.struct_head.id);
 224   if (radar_verbose_flag)
 225         fprintf(stderr, "ID = %d\n", (int)id);
 226   if (id != 7 && id != 27) { /* testing: Use 27 for Version 2 data */
 227         fprintf(stderr, "File is not a SIGMET version 1 nor version 2 raw product file.\n");
 228         return NULL;
 229   }
 230
 231   n = nsig_read_record(fp, (char *)&prod_file->rec2);
 232   if (radar_verbose_flag)
 233         fprintf(stderr, "Read %d bytes for rec2.\n", n);
 234
 235    /** Test for scan mode -- If scan is a RHI will return NULL  **/
 236    /** because RSL can't handle RHI's.  In the future, replace  **/
 237    /** NULL will a routine to convert RHI's to RSL Format       **/
 238    ant_scan_mode =NSIG_I2(prod_file->rec2.task_config.scan_info.ant_scan_mode);
 239    if(ant_scan_mode == 2)
 240           {
 241           if (radar_verbose_flag)
 242           fprintf(stderr, "RHI scan detected. Unable to process, returning NULL.\n");
 243           /*      return NULL; */
 244           }
 245
 246   /* Count the bits set in 'data_mask' to determine the number
 247    * of parameters present.
 248    */
 249   xh_size = NSIG_I2(prod_file->rec2.ingest_head.size_ext_ray_headers);
 250   nrays = NSIG_I2(prod_file->rec2.ingest_head.num_rays);
 251   if (radar_verbose_flag)
 252         fprintf(stderr, "Expecting %d rays in each sweep.\n", nrays);
 253   memmove(&data_mask, prod_file->rec2.task_config.dsp_info.data_mask, sizeof(fourb));
 254   for (nparams=i=0; i<32; i++)
 255         nparams += (data_mask >> i) & 0x1;
 256
 257   /* Number of sweeps */
 258   nsweeps = NSIG_I2(prod_file->rec2.task_config.scan_info.num_swp);
 259
 260
 261
 262    memmove(site_name, prod_file->rec1.prod_end.site_name, sizeof(prod_file->rec1.prod_end.site_name));
 263    site_name[sizeof(site_name)-1] = '\0';
 264   if (radar_verbose_flag) {
 265         fprintf(stderr, "nparams = %d, nsweeps = %d\n", nparams, nsweeps);
 266         fprintf(stderr, "Site name = <%s>\n", site_name);
 267   }
 268
 269     /* nsig_sweep = nsig_read_sweep(fp, prod_file)
 270      *
 271      * Use: nsig_sweep[i]->ray[j]->range
 272      *
 273      * where 'range' is [0..nbins-1]
 274      */
 275
 276     /*
 277      * All the information you need is in:
 278      *    prod_file->rec1
 279      *        .struct_head, .prod_config .prod_end
 280      *    prod_file->rec2
 281      *        .struct_head, .ingest_head, .task_config .device_stat,
 282      *        .dsp1, .dsp2
 283      *    nsig_sweep[0..nparams-1]  'nparams' is the true number
 284      *                              of parameters present.  You
 285      *                              must check the 'id' (or type)
 286      *                              to determine the field type.
 287      *                              So far seen, nparams <= 6.
 288      *    nsig_sweep[i]->bhdr     <NSIG_Raw_prod_bhdr>
 289      *    nsig_sweep[i]->idh      <NSIG_Ingest_data_header>
 290      *    nsig_sweep[i]->ray[j]   <NSIG_Ray *>
 291      *
 292      * Note:
 293      *    For extended header access, you'll typically use nsig_sweep[0]
 294      *    (double check the id) and the ray data allocated (nsig_ray->range)
 295      *    is a pointer to the extended header, either v0 or v1.
 296      *    You can typecast the pointer to NSIG_Ext_header_ver0 or
 297      *    NSIG_Ext_header_ver1, as you like.  To determine which
 298      *    version of the extended headers you have use:
 299      *      xh_size <= 20 for version 0, else version 1.
 300      *    Access:
 301      *      xh_size = NSIG_I2(prod_file->rec2.ingest_head.size_ext_ray_headers)
 302      *
 303      * Functions:
 304      *    NSIG_I2(nsig_sweep[i]->idh.num_rays_act);   -- # of rays. (j)
 305      *    NSIG_I2(nsig_sweep[i]->ray[j]->h.num_bins); -- # of bins in a ray.
 306      *
 307      *    NSIG_I2(x), NSIG_I4(x)   - Convert data, x, to floating point.
 308      *
 309      *    IMPORTANT NOTE: It must be known whether or not to perform
 310      *                    byte-swapping.  To determine this, call
 311      *                    'nsig_endianess'.  It returns 0 for no-swapping
 312      *                    and 1 for swapping.  Additionally, it transparently
 313      *                    initializes the nsig library to automatically
 314      *                    swap when using NSIG_I2 or NSIG_I4.
 315      *                    The function 'nsig_read_sweep' automatically
 316      *                    calls 'nsig_endianess', too.
 317      */
 318
 319    sea_lvl_hgt = NSIG_I2(prod_file->rec1.prod_end.grnd_sea_ht);
 320
 321    if (radar_verbose_flag)
 322          fprintf(stderr, "sea: %d\n", sea_lvl_hgt);
 323    if (radar_verbose_flag)
 324          fprintf(stderr, "site_name: %s", site_name);
 325
 326    /** Determine beamwidth from input variables (not saved in nsig file) **/
 327    if(strncmp(site_name,"mit",3) == 0 || strncmp(site_name,"MIT",3) == 0)
 328          beam_width = MIT_BEAMWIDTH;
 329    else if(strncmp(site_name,"tog",3) == 0 || strncmp(site_name,"TOG",3) == 0)
 330          beam_width = TOG_BEAMWIDTH;
 331    else if(strncmp(site_name,"kwa",3) == 0 || strncmp(site_name,"KWA",3) == 0)
 332          beam_width = KWA_BEAMWIDTH;
 333    else
 334          beam_width = DEFAULT_BEAMWIDTH;
 335
 336    if (radar_verbose_flag)
 337          fprintf(stderr, "beamwidth: %f\n", beam_width);
 338
 339    vert_half_bw = beam_width/2.0;
 340    horz_half_bw = beam_width/2.0;
 341
 342    /** Reading date and time **/
 343    month = NSIG_I2(prod_file->rec2.ingest_head.start_time.month);
 344    year = NSIG_I2(prod_file->rec2.ingest_head.start_time.year);
 345    day = NSIG_I2(prod_file->rec2.ingest_head.start_time.day);
 346    sec = NSIG_I4(prod_file->rec2.ingest_head.start_time.sec);
 347
 348    /* converting seconds since mid to time of day */
 349    tmp = sec/3600.0;
 350    hour = (int)tmp;
 351    tmp = (tmp - hour) * 60.0;
 352    minute = (int)tmp;
 353    second = (tmp - minute) * 60.0;
 354
 355    /** records of the nsig file.                             **/
 356    num_rays = 0;
 357    pw = (NSIG_I4(prod_file->rec1.prod_end.pulse_wd))/100.0; /* pulse width */
 358    prf = NSIG_I4(prod_file->rec1.prod_end.prf);  /* pulse repetition frequency */
 359    wave = (NSIG_I4(prod_file->rec1.prod_end.wavelen))/100.0; /* wavelength (cm) */
 360    rsl_kdp_wavelen = wave;  /* EXTERNAL (volume.c) This sets KD_F and KD_INVF
 361                              * to operate with the proper wavelength.
 362                              */
 363    numbins = NSIG_I4(prod_file->rec1.prod_end.num_bin);   /* # bins in ray */
 364    rng_first_bin = (float)NSIG_I4(prod_file->rec1.prod_end.rng_f_bin)/100.0;
 365    rng_last_bin = (float)NSIG_I4(prod_file->rec1.prod_end.rng_l_bin)/100.0;
 366    bin_space = ((rng_last_bin-rng_first_bin)/numbins); /*rng res (m)*/
 367
 368    numsweep = NSIG_I2(prod_file->rec2.task_config.scan_info.num_swp); /* # sweeps in volume */
 369    num_samples = NSIG_I2(prod_file->rec1.prod_end.num_samp);
 370    sweep_rate = 3.0; /** Approximate value -- info not stored **/
 371    azim_rate = sweep_rate*360.0/60.0;
 372    max_vel = wave*prf/(100.0*4.0);
 373    freq = (299793000.0/wave)*1.0e-4; /** freq in MHZ **/
 374
 375    sqi = NSIG_I2(prod_file->rec2.task_config.calib_info.sqi)/256.0;
 376    log = NSIG_I2(prod_file->rec2.task_config.calib_info.noise)/16.0;
 377    csr = NSIG_I2(prod_file->rec2.task_config.calib_info.clutr_corr)/(-16.0);
 378    sig = NSIG_I2(prod_file->rec2.task_config.calib_info.power)/16.0;
 379    cal_dbz = NSIG_I2(prod_file->rec2.task_config.calib_info.cal_ref)/16.0;
 380    z_flag_unc = NSIG_I2(prod_file->rec2.task_config.calib_info.z_flag_unc);
 381    z_flag_cor = NSIG_I2(prod_file->rec2.task_config.calib_info.z_flag_cor);
 382    v_flag = NSIG_I2(prod_file->rec2.task_config.calib_info.v_flag);
 383    w_flag = NSIG_I2(prod_file->rec2.task_config.calib_info.w_flag);
 384    speckle = NSIG_I2(prod_file->rec2.task_config.calib_info.speckle);
 385
 386    /** Verbose calibration information **/
 387    if (radar_verbose_flag)
 388           {
 389           fprintf(stderr, "LOG = %5.2f\n", log);
 390           fprintf(stderr, "SQI = %5.2f\n", sqi);
 391           fprintf(stderr, "CSR = %5.2f\n", csr);
 392           fprintf(stderr, "SIG = %5.2f\n", sig);
 393           fprintf(stderr, "Calibration reflectivity: %5.2f dBZ\n", cal_dbz);
 394           fprintf(stderr, "ZT flags: %d\n", z_flag_unc);  /** can find these **/
 395           fprintf(stderr, "DZ flags: %d\n", z_flag_cor);  /** defn in the    **/
 396           fprintf(stderr, "VR flags: %d\n", v_flag);      /** SIGMET Doc     **/
 397           fprintf(stderr, "SW flags: %d\n", w_flag);
 398           fprintf(stderr, "Flags: -3856  = SQI thresholding\n");
 399           fprintf(stderr, "       -21846 = LOG thresholding\n");
 400           fprintf(stderr, "       -24416 = LOG & SQI thresholding\n");
 401           fprintf(stderr, "       -24516 = LOG & SQI & SIG thresholding\n");
 402           fprintf(stderr, "speckle remover: %d\n", speckle);
 403           }
 404
 405    if (radar_verbose_flag)
 406          fprintf(stderr, "vel: %f prf: %f\n", max_vel, prf);
 407
 408    /** Extracting Latitude and Longitude from nsig file **/
 409    lat = nsig_from_fourb_ang(prod_file->rec2.ingest_head.lat_rad);
 410    lon = nsig_from_fourb_ang(prod_file->rec2.ingest_head.lon_rad);
 411    if(lat > 180.0) lat -= 360.0;
 412    if(lon > 180.0) lon -= 360.0;
 413    if (radar_verbose_flag)
 414          fprintf(stderr, "nsig_to_radar: lat %f, lon %f\n", lat, lon);
 415    /** Latitude deg, min, sec **/
 416    latd = (int)lat;
 417    tmp = (lat - latd) * 60.0;
 418    latm = (int)tmp;
 419    lats = (int)((tmp - latm) * 60.0);
 420    /** Longitude deg, min, sec **/
 421    lond = (int)lon;
 422    tmp = (lon - lond) * 60.0;
 423    lonm = (int)tmp;
 424    lons = (int)((tmp - lonm) * 60.0);
 425
 426    /** Allocating memory for radar structure **/
 427    radar = RSL_new_radar(MAX_RADAR_VOLUMES);
 428    if (radar == NULL)
 429           {
 430           fprintf(stderr, "nsig_to_radar: radar is NULL\n");
 431           return NULL;
 432           }
 433
 434    /** Filling Radar Header **/
 435    radar->h.month = month;
 436    radar->h.day = day;
 437    radar->h.year = year; /* Year 2000 compliant. */
 438    radar->h.hour = hour;
 439    radar->h.minute = minute;
 440    radar->h.sec = second;
 441    sprintf(radar->h.radar_type, "%s", radar_type);
 442    radar->h.number = radar_number;
 443    memmove(radar->h.name, site_name, sizeof(radar->h.name));
 444    memmove(radar->h.radar_name, site_name, sizeof(radar->h.radar_name));
 445    memmove(radar->h.city, city, sizeof(radar->h.city));
 446    memmove(radar->h.state, state, sizeof(radar->h.state));
 447    radar->h.latd = latd;
 448    radar->h.latm = latm;
 449    radar->h.lats = lats;
 450    radar->h.lond = lond;
 451    radar->h.lonm = lonm;
 452    radar->h.lons = lons;
 453    radar->h.height = (int)sea_lvl_hgt;
 454    radar->h.spulse = (int)(pw*1000);
 455    radar->h.lpulse = (int)(pw*1000);
 456
 457    if (radar_verbose_flag) {
 458 #ifdef NSIG_VER2
 459          fprintf(stderr, "\nSIGMET version 2 raw product file.\n");
 460 #else
 461          fprintf(stderr, "\nSIGMET version 1 raw product file.\n");
 462 #endif
 463          fprintf(stderr, "Date: %2.2d/%2.2d/%4.4d %2.2d:%2.2d:%f\n",
 464                          radar->h.month, radar->h.day, radar->h.year,
 465                          radar->h.hour, radar->h.minute, radar->h.sec);
 466          fprintf(stderr, "Name: ");
 467          for (i=0; i<sizeof(radar->h.name); i++)
 468            fprintf(stderr, "%c", radar->h.name[i]);
 469          fprintf(stderr, "\n");
 470          fprintf(stderr, "Lat/lon (%d %d' %d'', %d %d' %d'')\n",
 471                          radar->h.latd, radar->h.latm, radar->h.lats,
 472                          radar->h.lond, radar->h.lonm, radar->h.lons);
 473    }
 474
 475    /** Converting data **/
 476    if (radar_verbose_flag) fprintf(stderr, "Expecting %d sweeps.\n", numsweep);
 477    for(i = 0; i < numsweep; i++)
 478           {
 479                 nsig_sweep = nsig_read_sweep(fp, prod_file);
 480                 if (nsig_sweep == NULL) { /* EOF possibility */
 481                   if (feof(fp)) break;
 482                   else continue;
 483                 }
 484                 if (rsl_qsweep != NULL) {
 485                   if (i > rsl_qsweep_max) break;
 486                   if (rsl_qsweep[i] == 0) continue;
 487                 }
 488                 if (radar_verbose_flag)
 489                   fprintf(stderr, "Read sweep # %d\n", i);
 490         /* The whole sweep is 'nsig_sweep' ... pretty slick.
 491          *
 492          * nsig_sweep[itype]  -- [0..nparams], if non-null.
 493          */
 494         for (itype=0; itype<nparams; itype++) {
 495           if (nsig_sweep[itype] == NULL) continue;
 496
 497           /** Reading date and time **/
 498           sweep_month = NSIG_I2(nsig_sweep[itype]->idh.time.month);
 499           sweep_year = NSIG_I2(nsig_sweep[itype]->idh.time.year);
 500           sweep_day = NSIG_I2(nsig_sweep[itype]->idh.time.day);
 501           sweep_sec = NSIG_I4(nsig_sweep[itype]->idh.time.sec);
 502 #ifdef NSIG_VER2
 503           msec      = NSIG_I2(nsig_sweep[itype]->idh.time.msec);
 504           /*      printf("....... msec == %d\n", msec); */
 505 #endif
 506           /* converting seconds since mid to time of day */
 507           tmp = sweep_sec/3600.0;
 508           sweep_hour = (int)tmp;
 509           tmp = (tmp - sweep_hour) * 60.0;
 510           sweep_minute = (int)tmp;
 511           sweep_second = sweep_sec - (sweep_hour*3600 + sweep_minute*60);
 512
 513           num_rays = NSIG_I2(nsig_sweep[itype]->idh.num_rays_exp);
 514
 515           data_type = NSIG_I2(nsig_sweep[itype]->idh.data_type);
 516
 517           ifield = 0;
 518           switch (data_type) {
 519           case NSIG_DTB_EXH:
 520                   ifield = -1;
 521                   break;
 522           case NSIG_DTB_UCR:
 523                 ifield = ZT_INDEX;
 524                 f      = ZT_F;
 525                 invf   = ZT_INVF;
 526                 break;
 527           case NSIG_DTB_CR:
 528                 ifield = DZ_INDEX;
 529                 f      = DZ_F;
 530                 invf   = DZ_INVF;
 531                 break;
 532           case NSIG_DTB_VEL:
 533                 ifield = VR_INDEX;
 534                 f      = VR_F;
 535                 invf   = VR_INVF;
 536                 break;
 537           case NSIG_DTB_WID:
 538                 ifield = SW_INDEX;
 539                 f      = SW_F;
 540                 invf   = SW_INVF;
 541                 break;
 542           case NSIG_DTB_ZDR:
 543                 ifield = DR_INDEX;
 544                 f      = DR_F;
 545                 invf   = DR_INVF;
 546                 break;
 547           case NSIG_DTB_KDP:
 548                 ifield = KD_INDEX;
 549                 f      = KD_F;
 550                 invf   = KD_INVF;
 551                 break;
 552           case NSIG_DTB_PHIDP:     /* SRB 990127 */
 553                 ifield = PH_INDEX;
 554                 f      = PH_F;
 555                 invf   = PH_INVF;
 556                 break;
 557           case NSIG_DTB_RHOHV:     /* SRB 000414 */
 558                 ifield = RH_INDEX;
 559                 f      = RH_F;
 560                 invf   = RH_INVF;
 561                 break;
 562           case NSIG_DTB_SQI:
 563                 ifield = SQ_INDEX;
 564                 f      = SQ_F;
 565                 invf   = SQ_INVF;
 566                 break;
 567           default:
 568                 fprintf(stderr,"Unknown field type: %d  Skipping it.\n", data_type);
 569                 continue;
 570           }
 571
 572           if (radar_verbose_flag)
 573             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));
 574
 575           if (data_type != NSIG_DTB_EXH) {
 576                 if ((radar->v[ifield] == NULL)) {
 577                   if (rsl_qfield[ifield]) {
 578                          radar->v[ifield] = RSL_new_volume(numsweep);
 579                          radar->v[ifield]->h.f = f;
 580                          radar->v[ifield]->h.invf = invf;
 581                    } else {
 582                          /* Skip this field, because, the user does not want it. */
 583                          continue;
 584                    }
 585                 }
 586                  if (radar->v[ifield]->sweep[i] == NULL)
 587                    radar->v[ifield]->sweep[i] = RSL_new_sweep(num_rays);
 588                  }
 589           else
 590           continue;    /* Skip the actual extended header processing.
 591                         * This is different than getting it, so that
 592                         * the information is available for the other
 593                         * fields when filling the RSL ray headers.
 594                         */
 595
 596           /** DATA conversion time **/
 597           sweep = radar->v[ifield]->sweep[i];
 598           sweep->h.f = f;
 599           sweep->h.invf = invf;
 600           sweep->h.sweep_num = i;
 601           sweep->h.beam_width = beam_width;
 602           sweep->h.vert_half_bw = vert_half_bw;
 603           sweep->h.horz_half_bw = horz_half_bw;
 604           elev = nsig_from_bang(nsig_sweep[itype]->idh.fix_ang);
 605           sweep->h.elev = elev;
 606
 607           for(j = 0; j < num_rays; j++)
 608                 {
 609                   ray_p = nsig_sweep[itype]->ray[j];
 610                   if (ray_p == NULL) continue;
 611                   bin_num = NSIG_I2(ray_p->h.num_bins);
 612
 613                   /* Load extended header information, if available.
 614                    * We need to pass the entire nsig_sweep and search for
 615                    * the extended header field (it may not be data_type==0).
 616                    */
 617                   get_extended_header_info(nsig_sweep, xh_size, j, nparams,
 618                                            &msec, &azm, &elev,
 619                                            &pitch, &roll, &heading,
 620                                            &azm_rate, &elev_rate,
 621                                            &pitch_rate, &roll_rate, &heading_rate,
 622                                            &lat, &lon, &alt, &rvc,
 623                                            &vel_east, &vel_north, &vel_up);
 624
 625
 626                   if (radar->v[ifield]->sweep[i]->ray[j] == NULL)
 627                         radar->v[ifield]->sweep[i]->ray[j] = RSL_new_ray(bin_num);
 628                   ray = radar->v[ifield]->sweep[i]->ray[j];
 629                   ray->h.f = f;
 630                   ray->h.invf = invf;
 631                   /** Ray is at nsig_sweep[itype].ray->... **/
 632                   /** Loading nsig data into data structure **/
 633
 634                   ray->h.month  = sweep_month;
 635                   ray->h.day    = sweep_day;
 636                   ray->h.year   = sweep_year; /* Year 2000 compliant. */
 637                   ray->h.hour   = sweep_hour;
 638                   ray->h.minute = sweep_minute;
 639                   if (msec == 0) { /* No extended header */
 640                         ray->h.sec  = NSIG_I2(ray_p->h.sec) + sweep_second;
 641                         elev = sweep->h.elev;
 642                   } else
 643                         ray->h.sec  = sweep_second + msec/1000.0;
 644
 645                   /* add time ... handles end of min,hour,month,year and century. */
 646                   if (ray->h.sec >= 60) /* Should I fix the time no matter what? */
 647                         RSL_fix_time(ray);  /* Repair second overflow. */
 648
 649                   ray->h.ray_num    = j;
 650                   ray->h.elev_num   = i;
 651                   ray->h.range_bin1 = (int)rng_first_bin;
 652                   ray->h.gate_size  = (int)(bin_space+.5); /* Nearest int */
 653                   ray->h.vel_res    = bin_space;
 654                   ray->h.sweep_rate = sweep_rate;
 655                   ray->h.prf        = (int)prf;
 656                         if (prf != 0)
 657                           ray->h.unam_rng = 299793000.0 / (2.0 * prf * 1000.0);  /* km */
 658                         else
 659                           ray->h.unam_rng = 0.0;
 660                         ray->h.fix_angle = (float)sweep->h.elev;
 661                   ray->h.azim_rate  = azim_rate;
 662                   ray->h.pulse_count = (float)num_samples;
 663                   ray->h.pulse_width = pw;
 664                   ray->h.beam_width  = beam_width;
 665                   ray->h.frequency   = freq / 1000.0;  /* GHz */
 666                   ray->h.wavelength  = wave/100.0;     /* meters */
 667                   ray->h.nyq_vel     = max_vel;        /* m/s */
 668                   if (elev == 0.) elev = sweep->h.elev;
 669                   ray->h.elev        = elev;
 670                   /* Compute mean azimuth angle for ray. */
 671                   az1 = nsig_from_bang(ray_p->h.beg_azm);
 672                   az2 = nsig_from_bang(ray_p->h.end_azm);
 673                   /*              printf("az1, %f, az2 %f\n", az1, az2); */
 674                   if(az1 > az2)
 675                         if((az1 - az2) > 180.0) az2 += 360.0;
 676                         else
 677                           ;
 678                   else
 679                         if((az2 - az1) > 180.0) az1 += 360.0;
 680
 681                   az1 = (az1 + az2) / 2.0;
 682                   if (az1 > 360) az1 -= 360;
 683                   ray->h.azimuth     = az1;
 684
 685                   /* From the extended header information, we learn the following. */
 686                   ray->h.pitch        = pitch;
 687                   ray->h.roll         = roll;
 688                   ray->h.heading      = heading;
 689                   ray->h.pitch_rate   = pitch_rate;
 690                   ray->h.roll_rate    = roll_rate;
 691                   ray->h.heading_rate = heading_rate;
 692                   ray->h.lat          = lat;
 693                   ray->h.lon          = lon;
 694                   ray->h.alt          = alt;
 695                   ray->h.rvc          = rvc;
 696                   ray->h.vel_east     = vel_east;
 697                   ray->h.vel_north    = vel_north;
 698                   ray->h.vel_up       = vel_up;
 699
 700                   /*              printf("Processing sweep[%d]->ray[%d]: %d %f %f %f %f %f %f %f %f %d nbins=%d, bin1=%d gate=%d\n",
 701                                  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);
 702                                  */
 703                   if (data_type == NSIG_DTB_EXH) continue;
 704                   ray_data = 0;
 705                   for(k = 0; k < bin_num; k++) {
 706                         switch(data_type) {
 707                         case NSIG_DTB_UCR:
 708                         case NSIG_DTB_CR:
 709                           if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO;
 710                           else ray_data = (float)((ray_p->range[k]-64.0)/2.0);
 711                           break;
 712
 713                         case NSIG_DTB_VEL:
 714                           if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO;
 715                           else ray_data = (float)((ray_p->range[k]*max_vel/127.0)+
 716                                                    max_vel*(1.0-255.0/127.0));
 717                           break;
 718
 719                         case NSIG_DTB_WID:
 720                           if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO;
 721                           else ray_data =(float)((ray_p->range[k])/256.0)*max_vel;
 722                           break;
 723
 724                         case NSIG_DTB_ZDR:
 725                           if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO;
 726                           else ray_data = (float)((ray_p->range[k]-128.0)/16.0);
 727                           break;
 728
 729                           /*
 730                            * Special optimization note:
 731                            * For KDP, PHIDP, RHOHV we skip the float conversion,
 732                            * and carry the native sigmet data values into RSL storage.
 733                            */
 734                         case NSIG_DTB_KDP:
 735                                 ray_data = ray_p->range[k];
 736                                 /* F_OFFSET *must* match the definition in volume.c */
 737 #define F_OFFSET 4
 738                                 if (ray_data == 0 || ray_data == 255) ray_data = NSIG_NO_ECHO;
 739                                 else ray_data += F_OFFSET;
 740                                 break;
 741
 742                         case NSIG_DTB_RHOHV:
 743                         case NSIG_DTB_PHIDP:
 744                           if (ray_p->range[k] <= 0 || ray_p->range[k] >= 255)
 745                             ray_data = NSIG_NO_ECHO;
 746                           else
 747                             ray_data = ray_p->range[k];
 748                           break;
 749                         case NSIG_DTB_SQI:
 750                           if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO;
 751                           else ray_data =
 752                               (float)sqrt((ray_p->range[k]-1.0)/253.0);
 753                         }
 754
 755                         if (ray_data == NSIG_NO_ECHO)
 756                           ray->range[k] = ray->h.invf(BADVAL);
 757                         else
 758                           if ( (data_type == NSIG_DTB_KDP) ||
 759                                    (data_type == NSIG_DTB_RHOHV) ||
 760                                    (data_type == NSIG_DTB_PHIDP) )
 761                                 ray->range[k] = (Range)ray_data;
 762                           else
 763                                 ray->range[k] = ray->h.invf(ray_data);
 764
 765                         /*
 766                         if (data_type == NSIG_DTB_KDP)
 767                         printf("v[%d]->sweep[%d]->ray[%d]->range[%d] = %f, %d, %f\n",
 768                                    ifield, i, j, k, ray->h.f(ray->range[k]),
 769                                    (int)ray_p->range[k], ray_data);
 770                         */
 771                   }
 772                 }
 773             }
 774         nsig_free_sweep(nsig_sweep);
 775           }
 776
 777    /* Do not reset radar->h.nvolumes. It is already set properly. */
 778    if (radar_verbose_flag)
 779          fprintf(stderr, "Max index of radar->v[0..%d]\n", radar->h.nvolumes);
 780
 781
 782    /** close nsig file **/
 783    nsig_close(fp);
 784
 785    radar = RSL_prune_radar(radar);
 786    /** return radar pointer **/
 787    return radar;
 788 }