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
44 extern int radar_verbose_flag;
46 /* Some parameter headaches are prevented when vol is declared global. */
48 extern int read_entire_lassen_file(FILE *f, Lassen_volume *vol);
50 /**********************************************************************/
52 /* lassen_load_sweep */
54 /* By: John Merritt */
55 /* Space Applications Corporation */
57 /**********************************************************************/
58 void lassen_load_sweep(Sweep *s, int isweep_num, int ifield, int period, Lassen_sweep *ptr)
60 float c = RSL_SPEED_OF_LIGHT;
63 unsigned char *ray_data;
67 Range (*invf)(float x);
70 /* calibration period flags m.whimpey */
80 if (s == NULL) return;
81 f = (float (*)(Range x))NULL; /* This quiets the pedantic warning. */
82 invf = (Range (*)(float x))NULL;
83 if (ifield == ZT_INDEX) {kk = OFF_UZ; invf = ZT_INVF; f = ZT_F;}
84 if (ifield == DZ_INDEX) {kk = OFF_CZ; invf = DZ_INVF; f = DZ_F;}
85 if (ifield == VR_INDEX) {kk = OFF_VEL; invf = VR_INVF; f = VR_F;}
86 if (ifield == SW_INDEX) {kk = OFF_WID; invf = SW_INVF; f = SW_F;}
87 if (ifield == ZD_INDEX) {kk = OFF_ZDR; invf = ZD_INVF; f = ZD_F;}
88 if (ifield == PH_INDEX) {kk = OFF_PHI; invf = PH_INVF; f = PH_F;}
89 if (ifield == RH_INDEX) {kk = OFF_RHO; invf = RH_INVF; f = RH_F;}
90 if (ifield == LR_INDEX) {kk = OFF_LDR; invf = LR_INVF; f = LR_F;}
91 if (ifield == KD_INDEX) {kk = OFF_KDP; invf = KD_INVF; f = KD_F;}
92 if (ifield == TI_INDEX) {kk = OFF_TIME; invf = TI_INVF; f = TI_F;}
94 elev = (float)ptr->fangle*360.0/16384.0;
95 Vu = c*((float)vol.prf/10.)/(4.*(float)vol.freq*100000.0);
97 s->h.sweep_num = ptr->sweep;
99 s->h.nrays = ptr->numrays;
100 s->h.beam_width = 1.0; /* What is it really? */
101 s->h.horz_half_bw = .5;
102 s->h.vert_half_bw = .5;
106 for(i=0;i<(int)ptr->numrays;i++) {
108 if (aray == NULL) continue;
109 s->ray[i] = RSL_new_ray((int)aray->numgates);
110 s->ray[i]->h.month = aray->month;
111 s->ray[i]->h.day = aray->day;
112 s->ray[i]->h.year = (int)aray->year + 1900;
113 if (s->ray[i]->h.year < 1980) s->ray[i]->h.year += 100; /* Year > 2000. */
114 s->ray[i]->h.hour = aray->hour;
115 s->ray[i]->h.minute = aray->minute;
116 s->ray[i]->h.sec = aray->second;
117 s->ray[i]->h.azimuth = (float)aray->vangle*360.0/16384.0;
119 s->ray[i]->h.ray_num = i;
120 s->ray[i]->h.elev = elev;
121 s->ray[i]->h.elev_num = s->h.sweep_num;
122 s->ray[i]->h.range_bin1 = aray->rangeg1;
123 s->ray[i]->h.gate_size = aray->gatewid;
124 s->ray[i]->h.vel_res = 0.5; /* What is this really? */
126 s->ray[i]->h.fix_angle = s->h.elev;
127 s->ray[i]->h.frequency = (float)vol.freq*1.0e-4; /* GHz */
128 s->ray[i]->h.wavelength = c / s->ray[i]->h.frequency * 1.0e-9;
129 s->ray[i]->h.prf = (int)aray->prf/10;
130 s->ray[i]->h.nyq_vel = s->ray[i]->h.prf * s->ray[i]->h.wavelength / 4.0;
131 if (s->ray[i]->h.prf != 0)
132 s->ray[i]->h.unam_rng = c / (2.0 * s->ray[i]->h.prf * 1000.0); /* km */
134 s->ray[i]->h.unam_rng = 0.0;
135 s->ray[i]->h.pulse_width = (float)(aray->p_width * 0.05);
136 s->ray[i]->h.pulse_count = aray->n_pulses;
138 s->ray[i]->h.beam_width = 1.0; /* What is it really? */
140 s->ray[i]->h.invf = invf;
142 ray_data = (unsigned char *)aray;
144 m = aray->offset[kk];
147 /* conversion changes by m.whimpey */
148 for(j=0; j<s->ray[i]->h.nbins; j++) {
150 case OFF_UZ: /* UZ field. */
151 case OFF_CZ: /* CZ field. */
152 /* Apply a 1.4 dB correction. Ken Glasson did this. */
153 /* Removed 1.4 dB correction 09/27/2006--BMRC doesn't use it. */
154 if (period == BERRIMAH1)
155 x[j] = ((float)ray_data[m+j] - 56.0)/2.0; /* Removed +1.4 */
157 x[j] = ((float)ray_data[m+j] - 64.0)/2.0; /* Removed +1.4 */
159 case OFF_VEL: /* VR field */
160 if (period == BERRIMAH1)
161 x[j] = (float)(Vu*((double)ray_data[m+j]-128.0)/128.);
162 else if (period == SCSMEX)
163 x[j] = (float)(Vu*((double)(ray_data[m+j]^0x80)-128.0)/127.);
165 x[j] = (float)(Vu*((double)ray_data[m+j]-128.0)/127.);
166 /* fprintf(stderr,"Velocity for ray[%d] at x[%d] = %f, nyquist = %f\n",i,j,x[j],s->ray[i]->h.nyq_vel); */
168 case OFF_WID: /* SW field */
169 if (period == BERRIMAH1)
170 x[j] = (float)(Vu*(double)ray_data[m+j]/100.);
172 x[j] = (float)(Vu*(double)ray_data[m+j]/256.);
174 case OFF_ZDR: /* ZD field */
175 if (period < MCTEX_EARLY) break;
176 if (period <= BERRIMAH2)
177 x[j] = ((float)ray_data[m+j] - 64.0)/21.25;
179 x[j] = ((float)ray_data[m+j] - 128.)*18./254.;
181 case OFF_PHI: /* PH field */
182 if (period < MCTEX_EARLY) break;
183 if (period <= BERRIMAH2)
184 x[j] = ((float)ray_data[m+j] - 128.0)*32.0/127.0;
185 else if (period == GUNN_PT)
186 x[j] = ((float)ray_data[m+j] - 64.5)*360.0/254.0;
187 else if (period > GUNN_PT)
188 x[j] = 90.0+((float)ray_data[m+j] - 128.0)*180.0/254.0;
189 /*Extra 90 degrees added by Scott Collis, CAWCR/BMRC 2008*/
191 case OFF_RHO: /* RH field */
192 if (period < MCTEX_EARLY) break;
193 if (period <= BERRIMAH2)
194 x[j] = sqrt((float)ray_data[m+j]/256.822 + 0.3108);
196 x[j] = (((float)ray_data[m+j]-1.)*1.14/254.) + 0.01;
198 case OFF_LDR: /* LR field -- no 'period' here. */
199 x[j] = ((float)ray_data[m+j] - 250.0)/6;
201 case OFF_KDP: /* KP field -- no 'period' here. */
202 x[j] = ((float)ray_data[m+j]); /* This conversion is taken
203 from SIGMET. Is this right? */
205 case OFF_TIME: /* TIME field -- no 'period' here. */
206 x[j] = ((float)ray_data[m+j]); /* This conversion is taken
207 from SIGMET. Is this right? */
210 if (ray_data[m+j] == 0) x[j] = BADVAL;
212 for (j=0; j<s->ray[i]->h.nbins; j++)
213 s->ray[i]->range[j] = invf(x[j]);
219 /**********************************************************************/
221 /* RSL_lassen_to_radar */
223 /* By: John Merritt */
224 /* Space Applications Corporation */
226 /**********************************************************************/
227 Radar *RSL_lassen_to_radar(char *infile)
229 /* Lassen specific. */
232 int period; /* m.whimpey changed early variable to period */
234 int q[MAX_RADAR_VOLUMES];
235 extern int rsl_qfield[];
236 extern int *rsl_qsweep; /* See RSL_read_these_sweeps in volume.c */
237 extern int rsl_qsweep_max;
242 /* Listing of the RSL field type indexes, in the order that
243 * LASSEN stores them.
245 int rsl_index[] = {ZT_INDEX, DZ_INDEX, VR_INDEX, SW_INDEX, ZD_INDEX,
246 PH_INDEX, RH_INDEX, LR_INDEX, KD_INDEX, TI_INDEX};
248 char *ltype[] = {"UZ", "CZ", "Vel", "Wid", "Zdr", "Phi",
249 "Rho", "Ldr", "Kdp", "Time"};
252 struct compare_date { /* introduced by m.whimpey */
260 /* M.Whimpey changes made 19990422 for extra periods */
262 struct compare_date cvrt_date[NUM_DATES] =
263 { {1992, 1, 1, 0, 0, 0},
264 {1995, 11, 1, 0, 0, 0},
265 {1995, 11, 25, 20, 5, 0},
266 {1996, 1, 1, 0, 0, 0},
267 {1997, 10, 1, 0, 0, 0},
268 {1997, 11, 20, 3, 40, 0},
269 {1998, 05, 04, 0, 0, 0},
270 {1998, 05, 17, 03, 47, 0} };
272 int d; /* date counter */
274 unsigned long vt; /* vol time */
275 unsigned long dt; /* date time */
278 /* open Lassen file */
279 if (infile == NULL) {
282 f = fdopen(save_fd, "r");
284 if((f=fopen(infile, "r"))==(FILE *)NULL) {
288 f = uncompress_pipe(f); /* Transparently, use gunzip. */
289 #define NEW_BUFSIZ 16384
290 setvbuf(f,NULL,_IOFBF,(size_t)NEW_BUFSIZ); /* Faster i/o? */
292 if((read_entire_lassen_file(f, &vol)) == 0)
294 perror("RSL_lassen_to_radar ... read_entire_lassen_file");
300 if (radar_verbose_flag) {
301 fprintf(stderr,"\n Version = %d",vol.version);
302 fprintf(stderr,"\n Volume = %d",vol.volume);
303 fprintf(stderr,"\n Numsweeps = %d",vol.numsweeps);
304 fprintf(stderr,"\n Time = %2.2d/%2.2d/%2.2d %2.2d:%2.2d:%2.2d - %2.2d:%2.2d:%2.2d",
305 (int)vol.month, (int)vol.day, (int)vol.year,
306 (int)vol.shour, (int)vol.sminute, (int)vol.ssecond,
307 (int)vol.ehour, (int)vol.eminute, (int)vol.esecond);
308 fprintf(stderr,"\n Angle: start %d, stop %d", (int)vol.a_start, (int)vol.a_stop);
309 fprintf(stderr,"\n");
312 /* determine which period the lassen volume belongs m.whimpey */
313 vt = (vol.year-90) * 32140800;
314 vt += vol.month * 2678400;
315 vt += vol.day * 86400;
316 vt += vol.shour * 3600;
317 vt += vol.sminute * 60;
320 for(d=0; d<NUM_DATES; d++) {
321 dt = (cvrt_date[d].year-1990) * 32140800;
322 dt += cvrt_date[d].month * 2678400;
323 dt += cvrt_date[d].day * 86400;
324 dt += cvrt_date[d].hour * 3600;
325 dt += cvrt_date[d].minute * 60;
326 dt += cvrt_date[d].second;
334 fprintf(stderr, "%s: Error Vol date before first known!\n",
340 /* Max. expected volumes. */
341 radar = RSL_new_radar(MAX_RADAR_VOLUMES);
343 radar->h.month = vol.month;
344 radar->h.day = vol.day;
345 radar->h.year = vol.year + 1900;
346 radar->h.hour = vol.shour;
347 radar->h.minute = vol.sminute;
348 radar->h.sec = vol.ssecond;
349 strcpy(radar->h.radar_type, "lassen");
350 radar->h.nvolumes = MAX_RADAR_VOLUMES;
351 memcpy(&radar->h.radar_name, vol.radinfo.radar_name, 8);
352 memcpy(&radar->h.name, vol.radinfo.site_name, 8);
353 memcpy(&radar->h.city, "????", 4);
354 memcpy(&radar->h.state,"AU", 2);
355 radar->h.latd = vol.radinfo.latitude.degree;
356 radar->h.latm = vol.radinfo.latitude.minute;
357 radar->h.lats = vol.radinfo.latitude.second;
358 /* Is there a problem with the minutes/seconds when negative?
359 * The degree/minute/sec all should have the same sign.
361 if (radar->h.latd < 0) {
362 if (radar->h.latm > 0) radar->h.latm *= -1;
363 if (radar->h.lats > 0) radar->h.lats *= -1;
365 radar->h.lond = vol.radinfo.longitude.degree;
366 radar->h.lonm = vol.radinfo.longitude.minute;
367 radar->h.lons = vol.radinfo.longitude.second;
368 if (radar->h.lond < 0) {
369 if (radar->h.lonm > 0) radar->h.lonm *= -1;
370 if (radar->h.lons > 0) radar->h.lons *= -1;
372 radar->h.height = vol.radinfo.antenna_height;
377 /* iterate for each sweep in radar volume */
379 /* Determine which field types exist. The array 'q' is a boolean to
380 * force one print messages, if requested.
382 * Well, I tried looping through all the rays and examining aray->flags.<?>,
383 * but, it turns out that these flags bounce around, ie. toggle, for
384 * fields that don't really exist. Therefore, I cannot just logically OR
385 * the field flags together to determine which fields exist -- is this
386 * a bug when the file was created?
388 * What I've seen in lass2uf is to examine the first ray of the first sweep,
389 * but, I think that is unreliable too, due to the above finding.
391 * The solution, now, is to examine the existance of OFFSET values.
392 * This seems to be consistant, throughout the volume.
393 * The OFFSET is really what is used, anyway, to extract the data.
395 memset(q, 0, sizeof(q));
396 for (i=0; i<vol.numsweeps; i++) {
398 for (j=0; j<ptr->numrays; j++) {
400 for (k=0; k<NUMOFFSETS; k++) {
401 if (aray->offset[k] != 0 && !q[rsl_index[k]]) {
402 /* From RSL_select_fields */
403 if (rsl_qfield[rsl_index[k]] == 1)
409 if (radar_verbose_flag)
410 fprintf(stderr,"\n Fields are (Lassen nomenclature):");
411 for (k=0; k<NUMOFFSETS; k++) {
412 i = rsl_index[k]; /* Lassen index order to RSL index order translation. */
414 if (radar_verbose_flag) fprintf(stderr," %s", ltype[k]);
417 if (radar_verbose_flag) fprintf(stderr,"\n");
419 if (radar_verbose_flag) fprintf(stderr," Fields are (RSL nomenclature):");
420 for (k=0; k<NUMOFFSETS; k++) {
422 /* BTW, it doesn't matter if we allocate volumes, sweeps, or rays that
423 * have no data. The radar library will check for NULL pointers.
425 i = rsl_index[k]; /* Lassen index order to RSL index order translation. */
427 radar->v[i] = RSL_new_volume(vol.numsweeps);
428 radar->v[i]->h.f = RSL_f_list[i];
429 radar->v[i]->h.invf = RSL_invf_list[i];
430 if (radar_verbose_flag) fprintf(stderr," %s", RSL_ftype[i]);
431 if (k >= 2 && radar_verbose_flag) fprintf(stderr," "); /* Alignment. */
434 if (radar_verbose_flag) fprintf(stderr,"\n");
436 for(j=0;j<radar->h.nvolumes; j++) {
437 for(i=0;i<(int)vol.numsweeps;i++) {
438 if (rsl_qsweep != NULL) {
439 if (i > rsl_qsweep_max) break;
440 if (rsl_qsweep[i] == 0) continue;
444 radar->v[j]->sweep[i] = RSL_new_sweep(ptr->numrays);
446 /* 'period' is a flag for different calibrations */
447 lassen_load_sweep(radar->v[j]->sweep[i], i, j, period, ptr);
451 radar = RSL_prune_radar(radar);
455 Radar *RSL_lassen_to_radar(char *infile)
457 fprintf(stderr, "LASSEN is not installed in this version of RSL.\n");
458 fprintf(stderr, "Reinstall RSL w/ -DHAVE_LASSEN in the Makefile.\n");