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.
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.
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.
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.
21 /*************************************************************/
23 /* Function: nsig_to_radar.c */
26 /* Applied Research Corporation */
30 /* Modifications by: */
32 /* Space Applications Corporation */
36 /* Started: 08 AUG 96 */
38 /* Derived from Paul Kucera's nsig_to_radar.c */
40 /*************************************************************/
54 extern int radar_verbose_flag;
55 extern int rsl_qfield[]; /* See RSL_select_fields */
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
71 static float (*f)(Range x);
72 static Range (*invf)(float x);
76 void get_extended_header_info(NSIG_Sweep **nsig_sweep, int xh_size, int iray,
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,
83 float *lat, float *lon, int *alt, float *rvc,
84 float *vel_east, float *vel_north, float *vel_up)
86 static NSIG_Ext_header_ver1 xh;
89 *msec = *azm = *elev = *pitch = *roll = *heading =
90 *azm_rate = *elev_rate = *pitch_rate = *roll_rate = *heading_rate =
91 *lat = *lon = *alt = *rvc = 0;
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;
98 /* printf("...extended header itype=%d, nparams=%d\n", itype, nparams); */
99 if (itype == nparams) return; /* No extended header. */
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. */
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);
123 *heading_rate = nsig_from_bang(xh.heading_rate);
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 */
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;
160 int sweep_year, sweep_day, sweep_month;
161 int sweep_hour, sweep_minute, sweep_second;
163 int z_flag_unc, z_flag_cor, v_flag, w_flag, speckle;
168 float prf, prf2, wave, beam_width;
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;
173 float rng_first_bin, freq;
174 float max_vel, sweep_rate, azim_rate;
179 float sqi, log, csr, sig, cal_dbz;
180 char radar_type[50], state[4], city[15];
182 NSIG_Product_file *prod_file;
184 int data_mask, nrays;
186 int nparams, nsweeps;
187 NSIG_Sweep **nsig_sweep;
190 unsigned short nsig_2byte; /* New for 2-byte data types, Aug 2009 */
194 const unsigned short low10bits = 0x3ff;
195 float azm, elev, pitch, roll, heading, azm_rate, elev_rate,
196 pitch_rate, roll_rate, heading_rate,
199 int alt; /* Altitude */
200 float rvc; /* Radial correction velocity m/s */
201 float vel_east, vel_north, vel_up; /* Platform velocity vectors m/sec */
204 extern int *rsl_qsweep; /* See RSL_read_these_sweeps in volume.c */
205 extern int rsl_qsweep_max;
206 extern float rsl_kdp_wavelen;
209 if (radar_verbose_flag)
210 fprintf(stderr, "open file: %s\n", filename);
212 /** Opening nsig file **/
213 if((fp = nsig_open(filename)) == NULL) return NULL;
216 sprintf(radar_type, "nsig2");
217 radar_number = 22; /** Arbitrary number given to nsig2 data **/
219 sprintf(radar_type, "nsig");
220 radar_number = 21; /* What are these suppose to be? */
227 prod_file = (NSIG_Product_file *)calloc(1, sizeof(NSIG_Product_file));
229 n = nsig_read_record(fp, (char *)&prod_file->rec1);
230 nsig_endianess(&prod_file->rec1);
231 if (radar_verbose_flag)
232 fprintf(stderr, "Read %d bytes for rec1.\n", n);
234 id = NSIG_I2(prod_file->rec1.struct_head.id);
235 if (radar_verbose_flag)
236 fprintf(stderr, "ID = %d\n", (int)id);
237 if (id != 7 && id != 27) { /* testing: Use 27 for Version 2 data */
238 fprintf(stderr, "File is not a SIGMET version 1 nor version 2 raw product file.\n");
242 n = nsig_read_record(fp, (char *)&prod_file->rec2);
243 if (radar_verbose_flag)
244 fprintf(stderr, "Read %d bytes for rec2.\n", n);
246 /* Count the bits set in 'data_mask' to determine the number
247 * of parameters present.
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);
254 memmove(&masks[0], prod_file->rec2.task_config.dsp_info.data_mask_cur.mask_word_0,
256 memmove(&masks[1], &prod_file->rec2.task_config.dsp_info.data_mask_cur.mask_word_1,
259 for (j=0; j < 5; j++) {
260 data_mask = masks[j];
262 nparams += (data_mask >> i) & 0x1;
265 memmove(&data_mask, prod_file->rec2.task_config.dsp_info.data_mask, sizeof(fourb));
266 for (nparams=i=0; i<32; i++)
267 nparams += (data_mask >> i) & 0x1;
270 /* Number of sweeps */
271 nsweeps = NSIG_I2(prod_file->rec2.task_config.scan_info.num_swp);
275 memmove(site_name, prod_file->rec1.prod_end.site_name, sizeof(prod_file->rec1.prod_end.site_name));
276 site_name[sizeof(site_name)-1] = '\0';
277 if (radar_verbose_flag) {
278 fprintf(stderr, "nparams = %d, nsweeps = %d\n", nparams, nsweeps);
279 fprintf(stderr, "Site name = <%s>\n", site_name);
282 /* nsig_sweep = nsig_read_sweep(fp, prod_file)
284 * Use: nsig_sweep[i]->ray[j]->range
286 * where 'range' is [0..nbins-1]
290 * All the information you need is in:
292 * .struct_head, .prod_config .prod_end
294 * .struct_head, .ingest_head, .task_config .device_stat,
296 * nsig_sweep[0..nparams-1] 'nparams' is the true number
297 * of parameters present. You
298 * must check the 'id' (or type)
299 * to determine the field type.
300 * So far seen, nparams <= 6.
301 * nsig_sweep[i]->bhdr <NSIG_Raw_prod_bhdr>
302 * nsig_sweep[i]->idh <NSIG_Ingest_data_header>
303 * nsig_sweep[i]->ray[j] <NSIG_Ray *>
306 * For extended header access, you'll typically use nsig_sweep[0]
307 * (double check the id) and the ray data allocated (nsig_ray->range)
308 * is a pointer to the extended header, either v0 or v1.
309 * You can typecast the pointer to NSIG_Ext_header_ver0 or
310 * NSIG_Ext_header_ver1, as you like. To determine which
311 * version of the extended headers you have use:
312 * xh_size <= 20 for version 0, else version 1.
314 * xh_size = NSIG_I2(prod_file->rec2.ingest_head.size_ext_ray_headers)
317 * NSIG_I2(nsig_sweep[i]->idh.num_rays_act); -- # of rays. (j)
318 * NSIG_I2(nsig_sweep[i]->ray[j]->h.num_bins); -- # of bins in a ray.
320 * NSIG_I2(x), NSIG_I4(x) - Convert data, x, to floating point.
322 * IMPORTANT NOTE: It must be known whether or not to perform
323 * byte-swapping. To determine this, call
324 * 'nsig_endianess'. It returns 0 for no-swapping
325 * and 1 for swapping. Additionally, it transparently
326 * initializes the nsig library to automatically
327 * swap when using NSIG_I2 or NSIG_I4.
328 * The function 'nsig_read_sweep' automatically
329 * calls 'nsig_endianess', too.
332 sea_lvl_hgt = NSIG_I2(prod_file->rec1.prod_end.grnd_sea_ht);
334 if (radar_verbose_flag)
335 fprintf(stderr, "sea: %d\n", sea_lvl_hgt);
336 if (radar_verbose_flag)
337 fprintf(stderr, "site_name: %s", site_name);
339 /** Determine beamwidth from input variables (not saved in nsig file) **/
340 if(strncmp(site_name,"mit",3) == 0 || strncmp(site_name,"MIT",3) == 0)
341 beam_width = MIT_BEAMWIDTH;
342 else if(strncmp(site_name,"tog",3) == 0 || strncmp(site_name,"TOG",3) == 0)
343 beam_width = TOG_BEAMWIDTH;
344 else if(strncmp(site_name,"kwa",3) == 0 || strncmp(site_name,"KWA",3) == 0)
345 beam_width = KWA_BEAMWIDTH;
347 beam_width = DEFAULT_BEAMWIDTH;
349 if (radar_verbose_flag)
350 fprintf(stderr, "beamwidth: %f\n", beam_width);
352 vert_half_bw = beam_width/2.0;
353 horz_half_bw = beam_width/2.0;
355 /** Reading date and time **/
356 month = NSIG_I2(prod_file->rec2.ingest_head.start_time.month);
357 year = NSIG_I2(prod_file->rec2.ingest_head.start_time.year);
358 day = NSIG_I2(prod_file->rec2.ingest_head.start_time.day);
359 sec = NSIG_I4(prod_file->rec2.ingest_head.start_time.sec);
361 /* converting seconds since mid to time of day */
364 tmp = (tmp - hour) * 60.0;
366 second = (tmp - minute) * 60.0;
368 /** records of the nsig file. **/
370 pw = (NSIG_I4(prod_file->rec1.prod_end.pulse_wd))/100.0; /* pulse width */
371 prf = NSIG_I4(prod_file->rec1.prod_end.prf); /* pulse repetition frequency */
372 prf_mode = NSIG_I2(prod_file->rec2.task_config.dsp_info.prf_mode);
373 prf2 = prf * prf_modes[prf_mode];
374 wave = (NSIG_I4(prod_file->rec1.prod_end.wavelen))/100.0; /* wavelength (cm) */
375 rsl_kdp_wavelen = wave; /* EXTERNAL (volume.c) This sets KD_F and KD_INVF
376 * to operate with the proper wavelength.
378 numbins = NSIG_I4(prod_file->rec1.prod_end.num_bin); /* # bins in ray */
379 rng_first_bin = (float)NSIG_I4(prod_file->rec1.prod_end.rng_f_bin)/100.0;
380 rng_last_bin = (float)NSIG_I4(prod_file->rec1.prod_end.rng_l_bin)/100.0;
381 bin_space = ((rng_last_bin-rng_first_bin)/numbins); /*rng res (m)*/
383 numsweep = NSIG_I2(prod_file->rec2.task_config.scan_info.num_swp); /* # sweeps in volume */
384 num_samples = NSIG_I2(prod_file->rec1.prod_end.num_samp);
385 sweep_rate = 3.0; /** Approximate value -- info not stored **/
386 azim_rate = sweep_rate*360.0/60.0;
389 float max_vel1 = wave*prf/(100.0*4.0);
390 float max_vel2 = wave*prf2/(100.0*4.0);
392 max_vel = (max_vel1 * max_vel2)/(max_vel1-max_vel2);
396 max_vel = wave*prf/(100.0*4.0);
399 freq = (299793000.0/wave)*1.0e-4; /** freq in MHZ **/
401 sqi = NSIG_I2(prod_file->rec2.task_config.calib_info.sqi)/256.0;
402 log = NSIG_I2(prod_file->rec2.task_config.calib_info.noise)/16.0;
403 csr = NSIG_I2(prod_file->rec2.task_config.calib_info.clutr_corr)/(-16.0);
404 sig = NSIG_I2(prod_file->rec2.task_config.calib_info.power)/16.0;
405 cal_dbz = NSIG_I2(prod_file->rec2.task_config.calib_info.cal_ref)/16.0;
406 z_flag_unc = NSIG_I2(prod_file->rec2.task_config.calib_info.z_flag_unc);
407 z_flag_cor = NSIG_I2(prod_file->rec2.task_config.calib_info.z_flag_cor);
408 v_flag = NSIG_I2(prod_file->rec2.task_config.calib_info.v_flag);
409 w_flag = NSIG_I2(prod_file->rec2.task_config.calib_info.w_flag);
410 speckle = NSIG_I2(prod_file->rec2.task_config.calib_info.speckle);
412 /** Verbose calibration information **/
413 if (radar_verbose_flag)
415 fprintf(stderr, "LOG = %5.2f\n", log);
416 fprintf(stderr, "SQI = %5.2f\n", sqi);
417 fprintf(stderr, "CSR = %5.2f\n", csr);
418 fprintf(stderr, "SIG = %5.2f\n", sig);
419 fprintf(stderr, "Calibration reflectivity: %5.2f dBZ\n", cal_dbz);
420 fprintf(stderr, "ZT flags: %d\n", z_flag_unc); /** can find these **/
421 fprintf(stderr, "DZ flags: %d\n", z_flag_cor); /** defn in the **/
422 fprintf(stderr, "VR flags: %d\n", v_flag); /** SIGMET Doc **/
423 fprintf(stderr, "SW flags: %d\n", w_flag);
424 fprintf(stderr, "Flags: -3856 = SQI thresholding\n");
425 fprintf(stderr, " -21846 = LOG thresholding\n");
426 fprintf(stderr, " -24416 = LOG & SQI thresholding\n");
427 fprintf(stderr, " -24516 = LOG & SQI & SIG thresholding\n");
428 fprintf(stderr, "speckle remover: %d\n", speckle);
431 if (radar_verbose_flag)
432 fprintf(stderr, "vel: %f prf: %f prf2: %f\n", max_vel, prf, prf2);
434 /** Extracting Latitude and Longitude from nsig file **/
435 lat = nsig_from_fourb_ang(prod_file->rec2.ingest_head.lat_rad);
436 lon = nsig_from_fourb_ang(prod_file->rec2.ingest_head.lon_rad);
437 if(lat > 180.0) lat -= 360.0;
438 if(lon > 180.0) lon -= 360.0;
439 if (radar_verbose_flag)
440 fprintf(stderr, "nsig_to_radar: lat %f, lon %f\n", lat, lon);
441 /** Latitude deg, min, sec **/
443 tmp = (lat - latd) * 60.0;
445 lats = (int)((tmp - latm) * 60.0);
446 /** Longitude deg, min, sec **/
448 tmp = (lon - lond) * 60.0;
450 lons = (int)((tmp - lonm) * 60.0);
452 /** Allocating memory for radar structure **/
453 radar = RSL_new_radar(MAX_RADAR_VOLUMES);
456 fprintf(stderr, "nsig_to_radar: radar is NULL\n");
460 /** Filling Radar Header **/
461 radar->h.month = month;
463 radar->h.year = year; /* Year 2000 compliant. */
464 radar->h.hour = hour;
465 radar->h.minute = minute;
466 radar->h.sec = second;
467 sprintf(radar->h.radar_type, "%s", radar_type);
468 radar->h.number = radar_number;
469 memmove(radar->h.name, site_name, sizeof(radar->h.name));
470 memmove(radar->h.radar_name, site_name, sizeof(radar->h.radar_name));
471 memmove(radar->h.city, city, sizeof(radar->h.city));
472 memmove(radar->h.state, state, sizeof(radar->h.state));
473 radar->h.latd = latd;
474 radar->h.latm = latm;
475 radar->h.lats = lats;
476 radar->h.lond = lond;
477 radar->h.lonm = lonm;
478 radar->h.lons = lons;
479 radar->h.height = (int)sea_lvl_hgt;
480 radar->h.spulse = (int)(pw*1000);
481 radar->h.lpulse = (int)(pw*1000);
482 ant_scan_mode = NSIG_I2(prod_file->rec2.task_config.scan_info.ant_scan_mode);
483 if(ant_scan_mode == 2 || ant_scan_mode == 7) radar->h.scan_mode = RHI;
484 else radar->h.scan_mode = PPI;
486 if (radar_verbose_flag) {
488 fprintf(stderr, "\nSIGMET version 2 raw product file.\n");
490 fprintf(stderr, "\nSIGMET version 1 raw product file.\n");
492 fprintf(stderr, "Date: %2.2d/%2.2d/%4.4d %2.2d:%2.2d:%f\n",
493 radar->h.month, radar->h.day, radar->h.year,
494 radar->h.hour, radar->h.minute, radar->h.sec);
495 fprintf(stderr, "Name: ");
496 for (i=0; i<sizeof(radar->h.name); i++)
497 fprintf(stderr, "%c", radar->h.name[i]);
498 fprintf(stderr, "\n");
499 fprintf(stderr, "Lat/lon (%d %d' %d'', %d %d' %d'')\n",
500 radar->h.latd, radar->h.latm, radar->h.lats,
501 radar->h.lond, radar->h.lonm, radar->h.lons);
504 /** Converting data **/
505 if (radar_verbose_flag) fprintf(stderr, "Expecting %d sweeps.\n", numsweep);
506 for(i = 0; i < numsweep; i++)
508 nsig_sweep = nsig_read_sweep(fp, prod_file);
509 if (nsig_sweep == NULL) { /* EOF possibility */
513 if (rsl_qsweep != NULL) {
514 if (i > rsl_qsweep_max) break;
515 if (rsl_qsweep[i] == 0) continue;
517 if (radar_verbose_flag)
518 fprintf(stderr, "Read sweep # %d\n", i);
519 /* The whole sweep is 'nsig_sweep' ... pretty slick.
521 * nsig_sweep[itype] -- [0..nparams], if non-null.
523 for (itype=0; itype<nparams; itype++) {
524 if (nsig_sweep[itype] == NULL) continue;
526 /** Reading date and time **/
527 sweep_month = NSIG_I2(nsig_sweep[itype]->idh.time.month);
528 sweep_year = NSIG_I2(nsig_sweep[itype]->idh.time.year);
529 sweep_day = NSIG_I2(nsig_sweep[itype]->idh.time.day);
530 sweep_sec = NSIG_I4(nsig_sweep[itype]->idh.time.sec);
532 msec = NSIG_I2(nsig_sweep[itype]->idh.time.msec) & low10bits;
533 /* printf("....... msec == %d\n", msec); */
535 /* converting seconds since mid to time of day */
536 tmp = sweep_sec/3600.0;
537 sweep_hour = (int)tmp;
538 tmp = (tmp - sweep_hour) * 60.0;
539 sweep_minute = (int)tmp;
540 sweep_second = sweep_sec - (sweep_hour*3600 + sweep_minute*60);
542 num_rays = NSIG_I2(nsig_sweep[itype]->idh.num_rays_exp);
544 data_type = NSIG_I2(nsig_sweep[itype]->idh.data_type);
586 case NSIG_DTB_PHIDP: /* SRB 990127 */
591 case NSIG_DTB_RHOHV: /* SRB 000414 */
607 case NSIG_DTB_PHIDP2:
612 case NSIG_DTB_RHOHV2:
623 case NSIG_DTB_HCLASS:
624 case NSIG_DTB_HCLASS2:
650 fprintf(stderr,"Unknown field type: %d Skipping it.\n", data_type);
654 if (radar_verbose_flag)
655 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));
657 if (data_type != NSIG_DTB_EXH) {
658 if ((radar->v[ifield] == NULL)) {
659 if (rsl_qfield[ifield]) {
660 radar->v[ifield] = RSL_new_volume(numsweep);
661 radar->v[ifield]->h.f = f;
662 radar->v[ifield]->h.invf = invf;
664 /* Skip this field, because, the user does not want it. */
668 if (radar->v[ifield]->sweep[i] == NULL)
669 radar->v[ifield]->sweep[i] = RSL_new_sweep(num_rays);
672 continue; /* Skip the actual extended header processing.
673 * This is different than getting it, so that
674 * the information is available for the other
675 * fields when filling the RSL ray headers.
678 /** DATA conversion time **/
679 sweep = radar->v[ifield]->sweep[i];
681 sweep->h.invf = invf;
682 sweep->h.sweep_num = i;
683 sweep->h.beam_width = beam_width;
684 sweep->h.vert_half_bw = vert_half_bw;
685 sweep->h.horz_half_bw = horz_half_bw;
686 fix_angle = nsig_from_bang(nsig_sweep[itype]->idh.fix_ang);
687 if (radar->h.scan_mode == PPI) sweep->h.elev = fix_angle;
688 else sweep->h.azimuth = fix_angle;
690 for(j = 0; j < num_rays; j++)
692 ray_p = nsig_sweep[itype]->ray[j];
693 if (ray_p == NULL) continue;
694 bin_num = NSIG_I2(ray_p->h.num_bins);
696 /* Load extended header information, if available.
697 * We need to pass the entire nsig_sweep and search for
698 * the extended header field (it may not be data_type==0).
700 get_extended_header_info(nsig_sweep, xh_size, j, nparams,
702 &pitch, &roll, &heading,
703 &azm_rate, &elev_rate,
704 &pitch_rate, &roll_rate, &heading_rate,
705 &lat, &lon, &alt, &rvc,
706 &vel_east, &vel_north, &vel_up);
709 if (radar->v[ifield]->sweep[i]->ray[j] == NULL)
710 radar->v[ifield]->sweep[i]->ray[j] = RSL_new_ray(bin_num);
711 ray = radar->v[ifield]->sweep[i]->ray[j];
714 /** Ray is at nsig_sweep[itype].ray->... **/
715 /** Loading nsig data into data structure **/
717 ray->h.month = sweep_month;
718 ray->h.day = sweep_day;
719 ray->h.year = sweep_year; /* Year 2000 compliant. */
720 ray->h.hour = sweep_hour;
721 ray->h.minute = sweep_minute;
722 if (msec == 0) { /* No extended header */
723 ray->h.sec = NSIG_I2(ray_p->h.sec) + sweep_second;
724 elev = sweep->h.elev;
726 ray->h.sec = sweep_second + msec/1000.0;
728 /* add time ... handles end of min,hour,month,year and century. */
729 if (ray->h.sec >= 60) /* Should I fix the time no matter what? */
730 RSL_fix_time(ray); /* Repair second overflow. */
734 ray->h.range_bin1 = (int)rng_first_bin;
735 ray->h.gate_size = (int)(bin_space+.5); /* Nearest int */
736 ray->h.vel_res = bin_space;
737 ray->h.sweep_rate = sweep_rate;
738 ray->h.prf = (int)prf;
740 ray->h.unam_rng = 299793000.0 / (2.0 * prf * 1000.0); /* km */
742 ray->h.unam_rng = 0.0;
743 ray->h.prf2 = (int) prf2;
744 ray->h.fix_angle = fix_angle;
745 ray->h.azim_rate = azim_rate;
746 ray->h.pulse_count = num_samples;
747 ray->h.pulse_width = pw;
748 ray->h.beam_width = beam_width;
749 ray->h.frequency = freq / 1000.0; /* GHz */
750 ray->h.wavelength = wave/100.0; /* meters */
751 ray->h.nyq_vel = max_vel; /* m/s */
752 if (elev == 0.) elev = sweep->h.elev;
754 /* Compute mean azimuth angle for ray. */
755 az1 = nsig_from_bang(ray_p->h.beg_azm);
756 az2 = nsig_from_bang(ray_p->h.end_azm);
757 /* printf("az1, %f, az2 %f\n", az1, az2); */
759 if((az1 - az2) > 180.0) az2 += 360.0;
763 if((az2 - az1) > 180.0) az1 += 360.0;
765 az1 = (az1 + az2) / 2.0;
766 if (az1 > 360) az1 -= 360;
767 ray->h.azimuth = az1;
769 /* Compute mean elevation angle for ray. */
770 elev1 = nsig_from_bang(ray_p->h.beg_elev);
771 elev2 = nsig_from_bang(ray_p->h.end_elev);
772 /* printf("elev1, %f, elev2 %f\n", elev1, elev2); */
774 if((elev1 - elev2) > 180.0) elev2 += 360.0;
778 if((elev2 - elev1) > 180.0) elev2 -= 360.0;
780 elev1 = (elev1 + elev2) / 2.0;
781 if (elev1 > 360) elev1 -= 360;
784 /* From the extended header information, we learn the following. */
785 ray->h.pitch = pitch;
787 ray->h.heading = heading;
788 ray->h.pitch_rate = pitch_rate;
789 ray->h.roll_rate = roll_rate;
790 ray->h.heading_rate = heading_rate;
795 ray->h.vel_east = vel_east;
796 ray->h.vel_north = vel_north;
797 ray->h.vel_up = vel_up;
799 /* printf("Processing sweep[%d]->ray[%d]: %d %f %f %f %f %f %f %f %f %d nbins=%d, bin1=%d gate=%d\n",
800 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);
802 /* TODO: ingest data header contains a value for bits-per-bin.
803 * This might be of use to allocate an array for ray->range with
804 * either 1-byte or 2-byte elements. Then there's no need for
805 * memmove() whenever we need 2 bytes.
808 if (data_type == NSIG_DTB_EXH) continue;
810 for(k = 0; k < bin_num; k++) {
816 if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO;
817 else ray_data = (float)((ray_p->range[k]-64.0)/2.0);
819 /* Simplified the velocity conversion for NSIG_DTB_VEL, using
820 * formula from IRIS Programmer's Manual. BLK, Oct 9 2009.
823 if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO;
824 else ray_data = (float)((ray_p->range[k]-128.0)/127.0)*max_vel;
828 if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO;
829 else ray_data =(float)((ray_p->range[k])/256.0)*max_vel;
833 if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO;
834 else ray_data = (float)((ray_p->range[k]-128.0)/16.0);
838 if (ray_p->range[k] == 0 || ray_p->range[k] == 255 ||
839 rsl_kdp_wavelen == 0.0) {
840 ray_data = NSIG_NO_ECHO;
843 if (ray_p->range[k] < 128)
845 pow((double)600.0,(double)((127-ray_p->range[k])/126.0))) /
847 else if (ray_p->range[k] > 128)
849 pow((double)600.0,(double)((ray_p->range[k]-129)/126.0))) /
856 if (ray_p->range[k] == 0 || ray_p->range[k] == 255)
857 ray_data = NSIG_NO_ECHO;
859 ray_data = 180.0*((ray_p->range[k]-1.0)/254.0);
863 if (ray_p->range[k] == 0 || ray_p->range[k] == 255)
864 ray_data = NSIG_NO_ECHO;
866 ray_data = sqrt((double)((ray_p->range[k]-1.0)/253.0));
869 case NSIG_DTB_HCLASS:
870 if (ray_p->range[k] == 0 || ray_p->range[k] == 255)
871 ray_data = NSIG_NO_ECHO;
873 ray_data = ray_p->range[k];
877 if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO;
878 else ray_data = (float)sqrt((ray_p->range[k]-1.0)/253.0);
882 if (ray_p->range[k] == 0) ray_data = NSIG_NO_ECHO;
884 incr=75./127.; /* (+|- 75m/s) / 254 values */
885 ray_data = (float)(ray_p->range[k]-128)*incr;
895 memmove(nsig_twob, &ray_p->range[2*k], 2);
896 nsig_2byte = NSIG_I2(nsig_twob);
897 if (nsig_2byte == 0 || nsig_2byte == 65535)
898 ray_data = NSIG_NO_ECHO2;
899 else ray_data = (float)(nsig_2byte-32768)/100.;
903 memmove(nsig_twob, &ray_p->range[2*k], 2);
904 nsig_2byte = NSIG_I2(nsig_twob);
905 if (nsig_2byte == 0 || nsig_2byte == 65535)
906 ray_data = NSIG_NO_ECHO2;
907 else ray_data = (float)nsig_2byte/100.;
910 case NSIG_DTB_PHIDP2:
911 memmove(nsig_twob, &ray_p->range[2*k], 2);
912 nsig_2byte = NSIG_I2(nsig_twob);
913 if (nsig_2byte == 0 || nsig_2byte == 65535)
914 ray_data = NSIG_NO_ECHO;
916 ray_data = 360.*(nsig_2byte-1)/65534.;
920 case NSIG_DTB_RHOHV2:
921 memmove(nsig_twob, &ray_p->range[2*k], 2);
922 nsig_2byte = NSIG_I2(nsig_twob);
923 if (nsig_2byte == 0 || nsig_2byte == 65535)
924 ray_data = NSIG_NO_ECHO2;
925 else ray_data = (float)(nsig_2byte-1)/65533.;
928 case NSIG_DTB_HCLASS2:
929 memmove(nsig_twob, &ray_p->range[2*k], 2);
930 nsig_2byte = NSIG_I2(nsig_twob);
931 if (nsig_2byte == 0 || nsig_2byte == 65535)
932 ray_data = NSIG_NO_ECHO2;
934 ray_data = nsig_2byte;
937 if (ray_data == NSIG_NO_ECHO || ray_data == NSIG_NO_ECHO2)
938 ray->range[k] = ray->h.invf(BADVAL);
940 ray->range[k] = ray->h.invf(ray_data);
943 if (data_type == NSIG_DTB_KDP)
944 printf("v[%d]->sweep[%d]->ray[%d]->range[%d] = %f, %d, %f\n",
945 ifield, i, j, k, ray->h.f(ray->range[k]),
946 (int)ray_p->range[k], ray_data);
951 nsig_free_sweep(nsig_sweep);
954 /* Do not reset radar->h.nvolumes. It is already set properly. */
955 if (radar_verbose_flag)
956 fprintf(stderr, "Max index of radar->v[0..%d]\n", radar->h.nvolumes);
959 /** close nsig file **/
962 radar = RSL_prune_radar(radar);
963 /** return radar pointer **/