3 This is the TRMM Office Radar Software Library.
4 Copyright (C) 1996, 1997
6 Space Applications Corporation
9 This library is free software; you can redistribute it and/or
10 modify it under the terms of the GNU Library General Public
11 License as published by the Free Software Foundation; either
12 version 2 of the License, or (at your option) any later version.
14 This library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 Library General Public License for more details.
19 You should have received a copy of the GNU Library General Public
20 License along with this library; if not, write to the Free
21 Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
25 Michael Whimpey from BMRC, Australia changed code, for the
26 different callibrations, between, Pre_mctex, mctex, Gunn_Pt
29 Michael Whimpey added more callibration periods over berrimah and scsmex.
30 please see code where there is m.whimpey
43 extern int radar_verbose_flag;
45 /* Some parameter headaches are prevented when vol is declared global. */
47 extern int read_entire_lassen_file(FILE *f, Lassen_volume *vol);
49 /**********************************************************************/
51 /* lassen_load_sweep */
53 /* By: John Merritt */
54 /* Space Applications Corporation */
56 /**********************************************************************/
57 void lassen_load_sweep(Sweep *s, int isweep_num, int ifield, int period, Lassen_sweep *ptr)
59 float c = RSL_SPEED_OF_LIGHT;
62 unsigned char *ray_data;
66 Range (*invf)(float x);
69 /* calibration period flags m.whimpey */
79 if (s == NULL) return;
80 f = (float (*)(Range x))NULL; /* This quiets the pedantic warning. */
81 invf = (Range (*)(float x))NULL;
82 if (ifield == ZT_INDEX) {kk = OFF_UZ; invf = ZT_INVF; f = ZT_F;}
83 if (ifield == DZ_INDEX) {kk = OFF_CZ; invf = DZ_INVF; f = DZ_F;}
84 if (ifield == VR_INDEX) {kk = OFF_VEL; invf = VR_INVF; f = VR_F;}
85 if (ifield == SW_INDEX) {kk = OFF_WID; invf = SW_INVF; f = SW_F;}
86 if (ifield == ZD_INDEX) {kk = OFF_ZDR; invf = ZD_INVF; f = ZD_F;}
87 if (ifield == PH_INDEX) {kk = OFF_PHI; invf = PH_INVF; f = PH_F;}
88 if (ifield == RH_INDEX) {kk = OFF_RHO; invf = RH_INVF; f = RH_F;}
89 if (ifield == LR_INDEX) {kk = OFF_LDR; invf = LR_INVF; f = LR_F;}
90 if (ifield == KD_INDEX) {kk = OFF_KDP; invf = KD_INVF; f = KD_F;}
91 if (ifield == TI_INDEX) {kk = OFF_TIME; invf = TI_INVF; f = TI_F;}
93 elev = (float)ptr->fangle*360.0/16384.0;
94 Vu = c*((float)vol.prf/10.)/(4.*(float)vol.freq*100000.0);
96 s->h.sweep_num = ptr->sweep;
98 s->h.nrays = ptr->numrays;
99 s->h.beam_width = 1.0; /* What is it really? */
100 s->h.horz_half_bw = .5;
101 s->h.vert_half_bw = .5;
105 for(i=0;i<(int)ptr->numrays;i++) {
107 if (aray == NULL) continue;
108 s->ray[i] = RSL_new_ray((int)aray->numgates);
109 s->ray[i]->h.month = aray->month;
110 s->ray[i]->h.day = aray->day;
111 s->ray[i]->h.year = (int)aray->year + 1900;
112 if (s->ray[i]->h.year < 1980) s->ray[i]->h.year += 100; /* Year > 2000. */
113 s->ray[i]->h.hour = aray->hour;
114 s->ray[i]->h.minute = aray->minute;
115 s->ray[i]->h.sec = aray->second;
116 s->ray[i]->h.azimuth = (float)aray->vangle*360.0/16384.0;
118 s->ray[i]->h.ray_num = i;
119 s->ray[i]->h.elev = elev;
120 s->ray[i]->h.elev_num = s->h.sweep_num;
121 s->ray[i]->h.range_bin1 = aray->rangeg1;
122 s->ray[i]->h.gate_size = aray->gatewid;
123 s->ray[i]->h.vel_res = 0.5; /* What is this really? */
125 s->ray[i]->h.fix_angle = s->h.elev;
126 s->ray[i]->h.frequency = (float)vol.freq*1.0e-4; /* GHz */
127 s->ray[i]->h.wavelength = c / s->ray[i]->h.frequency * 1.0e-9;
128 s->ray[i]->h.prf = (int)aray->prf/10;
129 s->ray[i]->h.nyq_vel = s->ray[i]->h.prf * s->ray[i]->h.wavelength / 4.0;
130 if (s->ray[i]->h.prf != 0)
131 s->ray[i]->h.unam_rng = c / (2.0 * s->ray[i]->h.prf * 1000.0); /* km */
133 s->ray[i]->h.unam_rng = 0.0;
134 s->ray[i]->h.pulse_width = (float)(aray->p_width * 0.05);
135 s->ray[i]->h.pulse_count = aray->n_pulses;
137 s->ray[i]->h.beam_width = 1.0; /* What is it really? */
139 s->ray[i]->h.invf = invf;
141 ray_data = (unsigned char *)aray;
143 m = aray->offset[kk];
146 /* conversion changes by m.whimpey */
147 for(j=0; j<s->ray[i]->h.nbins; j++) {
149 case OFF_UZ: /* UZ field. */
150 case OFF_CZ: /* CZ field. */
151 /* Apply a 1.4 dB correction. Ken Glasson did this. */
152 /* Removed 1.4 dB correction 09/27/2006--BMRC doesn't use it. */
153 if (period == BERRIMAH1)
154 x[j] = ((float)ray_data[m+j] - 56.0)/2.0; /* Removed +1.4 */
156 x[j] = ((float)ray_data[m+j] - 64.0)/2.0; /* Removed +1.4 */
158 case OFF_VEL: /* VR field */
159 if (period == BERRIMAH1)
160 x[j] = (float)(Vu*((double)ray_data[m+j]-128.0)/128.);
161 else if (period == SCSMEX)
162 x[j] = (float)(Vu*((double)(ray_data[m+j]^0x80)-128.0)/127.);
164 x[j] = (float)(Vu*((double)ray_data[m+j]-128.0)/127.);
165 /* fprintf(stderr,"Velocity for ray[%d] at x[%d] = %f, nyquist = %f\n",i,j,x[j],s->ray[i]->h.nyq_vel); */
167 case OFF_WID: /* SW field */
168 if (period == BERRIMAH1)
169 x[j] = (float)(Vu*(double)ray_data[m+j]/100.);
171 x[j] = (float)(Vu*(double)ray_data[m+j]/256.);
173 case OFF_ZDR: /* ZD field */
174 if (period < MCTEX_EARLY) break;
175 if (period <= BERRIMAH2)
176 x[j] = ((float)ray_data[m+j] - 64.0)/21.25;
178 x[j] = ((float)ray_data[m+j] - 128.)*18./254.;
180 case OFF_PHI: /* PH field */
181 if (period < MCTEX_EARLY) break;
182 if (period <= BERRIMAH2)
183 x[j] = ((float)ray_data[m+j] - 128.0)*32.0/127.0;
184 else if (period == GUNN_PT)
185 x[j] = ((float)ray_data[m+j] - 64.5)*360.0/254.0;
186 else if (period > GUNN_PT)
187 x[j] = 90.0+((float)ray_data[m+j] - 128.0)*180.0/254.0;
188 /*Extra 90 degrees added by Scott Collis, CAWCR/BMRC 2008*/
190 case OFF_RHO: /* RH field */
191 if (period < MCTEX_EARLY) break;
192 if (period <= BERRIMAH2)
193 x[j] = sqrt((float)ray_data[m+j]/256.822 + 0.3108);
195 x[j] = (((float)ray_data[m+j]-1.)*1.14/254.) + 0.01;
197 case OFF_LDR: /* LR field -- no 'period' here. */
198 x[j] = ((float)ray_data[m+j] - 250.0)/6;
200 case OFF_KDP: /* KP field -- no 'period' here. */
201 x[j] = ((float)ray_data[m+j]); /* This conversion is taken
202 from SIGMET. Is this right? */
204 case OFF_TIME: /* TIME field -- no 'period' here. */
205 x[j] = ((float)ray_data[m+j]); /* This conversion is taken
206 from SIGMET. Is this right? */
209 if (ray_data[m+j] == 0) x[j] = BADVAL;
211 for (j=0; j<s->ray[i]->h.nbins; j++)
212 s->ray[i]->range[j] = invf(x[j]);
218 /**********************************************************************/
220 /* RSL_lassen_to_radar */
222 /* By: John Merritt */
223 /* Space Applications Corporation */
225 /**********************************************************************/
226 Radar *RSL_lassen_to_radar(char *infile)
228 /* Lassen specific. */
231 int period; /* m.whimpey changed early variable to period */
233 int q[MAX_RADAR_VOLUMES];
234 extern int rsl_qfield[];
235 extern int *rsl_qsweep; /* See RSL_read_these_sweeps in volume.c */
236 extern int rsl_qsweep_max;
241 /* Listing of the RSL field type indexes, in the order that
242 * LASSEN stores them.
244 int rsl_index[] = {ZT_INDEX, DZ_INDEX, VR_INDEX, SW_INDEX, ZD_INDEX,
245 PH_INDEX, RH_INDEX, LR_INDEX, KD_INDEX, TI_INDEX};
247 char *ltype[] = {"UZ", "CZ", "Vel", "Wid", "Zdr", "Phi",
248 "Rho", "Ldr", "Kdp", "Time"};
251 struct compare_date { /* introduced by m.whimpey */
259 /* M.Whimpey changes made 19990422 for extra periods */
261 struct compare_date cvrt_date[NUM_DATES] =
262 { {1992, 1, 1, 0, 0, 0},
263 {1995, 11, 1, 0, 0, 0},
264 {1995, 11, 25, 20, 5, 0},
265 {1996, 1, 1, 0, 0, 0},
266 {1997, 10, 1, 0, 0, 0},
267 {1997, 11, 20, 3, 40, 0},
268 {1998, 05, 04, 0, 0, 0},
269 {1998, 05, 17, 03, 47, 0} };
271 int d; /* date counter */
273 unsigned long vt; /* vol time */
274 unsigned long dt; /* date time */
277 /* open Lassen file */
278 if (infile == NULL) {
281 f = fdopen(save_fd, "rb");
283 if((f=fopen(infile, "rb"))==(FILE *)NULL) {
287 f = uncompress_pipe(f); /* Transparently, use gunzip. */
288 #define NEW_BUFSIZ 16384
289 setvbuf(f,NULL,_IOFBF,(size_t)NEW_BUFSIZ); /* Faster i/o? */
291 if((read_entire_lassen_file(f, &vol)) == 0)
293 perror("RSL_lassen_to_radar ... read_entire_lassen_file");
299 if (radar_verbose_flag) {
300 fprintf(stderr,"\n Version = %d",vol.version);
301 fprintf(stderr,"\n Volume = %d",vol.volume);
302 fprintf(stderr,"\n Numsweeps = %d",vol.numsweeps);
303 fprintf(stderr,"\n Time = %2.2d/%2.2d/%2.2d %2.2d:%2.2d:%2.2d - %2.2d:%2.2d:%2.2d",
304 (int)vol.month, (int)vol.day, (int)vol.year,
305 (int)vol.shour, (int)vol.sminute, (int)vol.ssecond,
306 (int)vol.ehour, (int)vol.eminute, (int)vol.esecond);
307 fprintf(stderr,"\n Angle: start %d, stop %d", (int)vol.a_start, (int)vol.a_stop);
308 fprintf(stderr,"\n");
311 /* determine which period the lassen volume belongs m.whimpey */
312 vt = (vol.year-90) * 32140800;
313 vt += vol.month * 2678400;
314 vt += vol.day * 86400;
315 vt += vol.shour * 3600;
316 vt += vol.sminute * 60;
319 for(d=0; d<NUM_DATES; d++) {
320 dt = (cvrt_date[d].year-1990) * 32140800;
321 dt += cvrt_date[d].month * 2678400;
322 dt += cvrt_date[d].day * 86400;
323 dt += cvrt_date[d].hour * 3600;
324 dt += cvrt_date[d].minute * 60;
325 dt += cvrt_date[d].second;
333 fprintf(stderr, "%s: Error Vol date before first known!\n",
339 /* Max. expected volumes. */
340 radar = RSL_new_radar(MAX_RADAR_VOLUMES);
342 radar->h.month = vol.month;
343 radar->h.day = vol.day;
344 radar->h.year = vol.year + 1900;
345 radar->h.hour = vol.shour;
346 radar->h.minute = vol.sminute;
347 radar->h.sec = vol.ssecond;
348 strcpy(radar->h.radar_type, "lassen");
349 radar->h.nvolumes = MAX_RADAR_VOLUMES;
350 memcpy(&radar->h.radar_name, vol.radinfo.radar_name, 8);
351 memcpy(&radar->h.name, vol.radinfo.site_name, 8);
352 memcpy(&radar->h.city, "????", 4);
353 memcpy(&radar->h.state,"AU", 2);
354 radar->h.latd = vol.radinfo.latitude.degree;
355 radar->h.latm = vol.radinfo.latitude.minute;
356 radar->h.lats = vol.radinfo.latitude.second;
357 /* Is there a problem with the minutes/seconds when negative?
358 * The degree/minute/sec all should have the same sign.
360 if (radar->h.latd < 0) {
361 if (radar->h.latm > 0) radar->h.latm *= -1;
362 if (radar->h.lats > 0) radar->h.lats *= -1;
364 radar->h.lond = vol.radinfo.longitude.degree;
365 radar->h.lonm = vol.radinfo.longitude.minute;
366 radar->h.lons = vol.radinfo.longitude.second;
367 if (radar->h.lond < 0) {
368 if (radar->h.lonm > 0) radar->h.lonm *= -1;
369 if (radar->h.lons > 0) radar->h.lons *= -1;
371 radar->h.height = vol.radinfo.antenna_height;
376 /* iterate for each sweep in radar volume */
378 /* Determine which field types exist. The array 'q' is a boolean to
379 * force one print messages, if requested.
381 * Well, I tried looping through all the rays and examining aray->flags.<?>,
382 * but, it turns out that these flags bounce around, ie. toggle, for
383 * fields that don't really exist. Therefore, I cannot just logically OR
384 * the field flags together to determine which fields exist -- is this
385 * a bug when the file was created?
387 * What I've seen in lass2uf is to examine the first ray of the first sweep,
388 * but, I think that is unreliable too, due to the above finding.
390 * The solution, now, is to examine the existance of OFFSET values.
391 * This seems to be consistant, throughout the volume.
392 * The OFFSET is really what is used, anyway, to extract the data.
394 memset(q, 0, sizeof(q));
395 for (i=0; i<vol.numsweeps; i++) {
397 for (j=0; j<ptr->numrays; j++) {
399 for (k=0; k<NUMOFFSETS; k++) {
400 if (aray->offset[k] != 0 && !q[rsl_index[k]]) {
401 /* From RSL_select_fields */
402 if (rsl_qfield[rsl_index[k]] == 1)
408 if (radar_verbose_flag)
409 fprintf(stderr,"\n Fields are (Lassen nomenclature):");
410 for (k=0; k<NUMOFFSETS; k++) {
411 i = rsl_index[k]; /* Lassen index order to RSL index order translation. */
413 if (radar_verbose_flag) fprintf(stderr," %s", ltype[k]);
416 if (radar_verbose_flag) fprintf(stderr,"\n");
418 if (radar_verbose_flag) fprintf(stderr," Fields are (RSL nomenclature):");
419 for (k=0; k<NUMOFFSETS; k++) {
421 /* BTW, it doesn't matter if we allocate volumes, sweeps, or rays that
422 * have no data. The radar library will check for NULL pointers.
424 i = rsl_index[k]; /* Lassen index order to RSL index order translation. */
426 radar->v[i] = RSL_new_volume(vol.numsweeps);
427 radar->v[i]->h.f = RSL_f_list[i];
428 radar->v[i]->h.invf = RSL_invf_list[i];
429 if (radar_verbose_flag) fprintf(stderr," %s", RSL_ftype[i]);
430 if (k >= 2 && radar_verbose_flag) fprintf(stderr," "); /* Alignment. */
433 if (radar_verbose_flag) fprintf(stderr,"\n");
435 for(j=0;j<radar->h.nvolumes; j++) {
436 for(i=0;i<(int)vol.numsweeps;i++) {
437 if (rsl_qsweep != NULL) {
438 if (i > rsl_qsweep_max) break;
439 if (rsl_qsweep[i] == 0) continue;
443 radar->v[j]->sweep[i] = RSL_new_sweep(ptr->numrays);
445 /* 'period' is a flag for different calibrations */
446 lassen_load_sweep(radar->v[j]->sweep[i], i, j, period, ptr);
450 radar = RSL_prune_radar(radar);
454 Radar *RSL_lassen_to_radar(char *infile)
456 fprintf(stderr, "LASSEN is not installed in this version of RSL.\n");
457 fprintf(stderr, "Reinstall RSL w/ -DHAVE_LASSEN in the Makefile.\n");