/* hrdrad -- Convert HRD P-3 data to Candis format. author of main routines -- Dave Raymond Carlos Lopez-Carrillo November 11, 2002 -------------------------- hrdrad-exp -- modification of hrdrad, so it can be more easily mantained. Furthermore, I've decided to eliminate all kind of thresholds that can be done once we have a candis file. Unfolding is not an aoption any more. It has to be afterward. This might be a little more costly but it helps in disentangling translation from proccesing. I also have included in the output the nyquist frequency as a 1D field. Although it is normally a constant parameter thruoughout a single flight, it can be changed -- as it happend during EPIC mission 7. So, this field will serve later in selecting the appropiate unfolding procedure. Spectral width will also be included if it is in the input data. ------------------------------------------------------------------ MOD-December-5-2003 --------------------------------------------- December 5, 2003 I'm reinserting the following modification that got dropped in the second release that was developed to analyse EPIC data. "a modification MOD-July-7-1998 was added to handel situations where the field da (drift) is not generated by the translator. This situations arise when hrdrawair2 is used instead of hrdrawair to translate the insitu data from original to Candis Format." This MODification was introduced during the analysis of the TogaCoare data. ------------------------------------------------------------------ compilation sequence: gcc -O -c hrdrad-exp.c gcc -O -c getp3rec.c gcc -O -c unfoldit.c gcc -o hrdrad-exp hrdrad.o getp3rec.o unfoldit.o -lm -lcdf */ #include #include #include "cdfhdr.h" #define MAXCHAR 200 /*Max characters on the calling line*/ #define MAX 10000 /* max size of input record */ #define OUTMAX 4096 /* max size of uncompress buffer */ #define RAYHEAD 22 /* size of ray header in shorts */ #define D2R 0.01745329278 /*conversion factor to go from degrees to radians pi/180*/ #define USAGE "Usage:hrdradx -[t|l] [-fba] [-x af] [-s t1:t2] [-c c_raz:c_rel:c_heading:c_pitch:c_fgate]" /*Vars to prove headers */ short year_h,mon_h,day_h,hour_h,min_h,sec_h; short ac_id_h; char flt_id_h[9]; short d_head_size_h,r_head_size_h,tzone_h; char proj_id_h[17]; short lf_data_mode_h,lf_vargates_h,lf_num_in_h,lf_num_out_h; float lf_firstgate_h,lf_gatestep_h; short ta_data_mode_h,ta_vargates_h,ta_num_in_h,ta_num_out_h; float ta_firstgate_h,ta_gatestep_h; float nyv_h; /* input buffer */ short buff[MAX]; long length; short getshort(); /* uncompress buffer */ short ubuff[OUTMAX]; long outlen; /* stuff for header buffer */ short year,mon,day,hour,min,sec; short ac_id; char flt_id[9]; short d_head_size,r_head_size,tzone; char proj_id[17]; short lf_data_mode,lf_vargates,lf_num_in,lf_num_out; float lf_firstgate,lf_gatestep; short ta_data_mode,ta_vargates,ta_num_in,ta_num_out; float ta_firstgate,ta_gatestep; float nyv; float aux01,aux02; /* stuff for data header */ short sweep,record,radar; long bpointer; /* stuff for ray header */ short raybytes,refl,vel,width,ryear,rmon,rday,rhour,rmin,rsec,rmsec; float lon,lat,alt,ua,va,wa,u,v,w,raz,rel,pitch,roll,drift,heading; /* Candis stuff */ char hbuf[HBMAX][LINE]; char cline[LINE],pline[LINE]; float *sbuff,*vbuff; long nsbuff,nvbuff; float *range,*nyq_vel,*raycnt,*cdfdbz,*cdfvr,*cdfwidth,*cdfwid; float *cdfyear,*cdfmon,*cdfday,*cdfhour,*cdfmin,*cdfsec,*cdfmsec; float *cdflon,*cdflat,*cdfalt,*cdfua,*cdfva,*cdfwa,*cdfu,*cdfv,*cdfw; float *cdfraz,*cdfrel,*cdfpitch,*cdfroll,*cdfdrift,*cdfheading; float *cdfaz,*cdfel,*cdftime; int cdfgates; float bad,badlim; static long iraycnt; char *corr; char *afilename; static float *monotime; /* monotonic increasing time */ static float *nyqvel; /* aircraft file stuff */ char abuf[HBMAX][LINE]; float *sabuf,*vabuf; long nsabuf,nvabuf; float *aua,*ava,*awa,*au,*av,*aw,*alon,*alat,*aalt,*apitch,*aroll,*aheading, *adrift,*tptime; static long actime,bctime,maxtime; struct field *fp; FILE *afile; /*MOD-July-7-1998*/ static int drift_flag; /* additive corrections */ float c_raz,c_rel; /* raw azimuth and elevation (degrees) */ float c_heading; /* heading (degrees) */ float c_pitch; /* pitch angle (degrees) */ float c_fgate; /* distance of first gate (km) */ float c_junk; /* fake entry to detect bad format */ /* getopt stuff */ extern char *optarg; extern int optind; int c; static char argu[MAXCHAR]; char ACF_name[5]; /* function to fix time */ float mtimefix(float time_new,int day_new); /*function to get and check arguments */ void fun_cgeta(int argc,char **argv); void fun_candis_setup(); void fun_use_ac_navigation(float *lon,float *lat,float *alt, float *ua,float *va,float *wa,float *u,float *v,float *w, float *pitch,float *roll,float *drift,float *heading); int taflag,lfflag,vrangeflag,gateaverage,timeselect,flag_acf; int front,back; int gateaverage,timeselect; float time1,time2; float crange,drange,adrange; main(int argc, char **argv) { int value,nvars,loop,poop; int headflag; int p3file; int unc; float reqheading,tolerance,rolltol,difference; float vscale; float dbzvalue; fun_cgeta(argc,argv); /* only one header allowed */ headflag = 0; p3file = 0; /* read P-3 radar records one at a time */ while ((length = getp3rec(MAX,buff)) >= 0) { /*loop on records -- both header and data*/ if (length > 0) { /* if record length zero, skip it */ if (getshort(buff,0) == 0) { /* check if this is a header record */ /* skip header if we get a second one -- Dave Jorgensen says that number of gates and gate spacing never changes -- thus, we can do this. However, other parameters may change, for instance the nyquist velocity was changed during EPIC mission 7 */ if (headflag == 0) { /* first header. decode it to working vars*/ p3header(buff,length,&year,&mon,&day,&hour,&min,&sec, &ac_id,flt_id,&d_head_size,&r_head_size,&tzone,proj_id, &lf_data_mode,&lf_vargates,&lf_num_in,&lf_num_out, &lf_firstgate,&lf_gatestep, &ta_data_mode,&ta_vargates,&ta_num_in,&ta_num_out, &ta_firstgate,&ta_gatestep,&nyv); /* only do processed data */ if ((taflag && ta_data_mode != 1) || (!taflag && lf_data_mode != 1)) { fprintf(stderr,"hrdrad: time series data not accepted\n"); exit(1); } fun_candis_setup(); headflag++; } else { /*other headers */ fprintf(stderr,"hrdrad: warning: another volume header\n"); /* decode the header to probe vars*/ p3header(buff,length,&year_h,&mon_h,&day_h,&hour_h,&min_h,&sec_h, &ac_id,flt_id_h,&d_head_size_h,&r_head_size_h,&tzone_h,proj_id_h, &lf_data_mode_h,&lf_vargates_h,&lf_num_in_h,&lf_num_out_h, &lf_firstgate_h,&lf_gatestep_h, &ta_data_mode_h,&ta_vargates_h,&ta_num_in_h,&ta_num_out_h, &ta_firstgate_h,&ta_gatestep_h,&nyv_h); if (nyv_h != nyv ){ fprintf(stderr,"hrdrad: warning: nyquist velocity has change\n"); fprintf(stderr,"from nyv =%f \t to nyv_n = %f ; last ray time %f\n",nyv,nyv_h,*monotime); /*stop if requested on input params */ /* fprintf(stderr," .. stopping!!!!\n"); */ /* stop(); */ /* otherwise update working var */ nyv = nyv_h; } else /* print out some info */ fprintf(stderr,"nyv =%f \t last ray time %f\n",nyv,*monotime); /* continue; anyway with record reading*/ } } else /*check if this is a data record */ if (getshort(buff,0) == 1) { /* data record */ if (headflag > 0) { /* procced only if a header record has been read */ /* get data record header */ p3dhead(buff,length,&bpointer,&sweep,&record,&radar); /* only keep data from selected radar */ /* old condition if ((taflag && radar == 1) || (!taflag && radar == 2)) continue; */ if ((taflag && radar == 2) || (!taflag && radar == 1)) { while (bpointer + RAYHEAD <= length) { /* loop on rays */ /* get ray header; p3rhead output: lon, lat (degrees) -- longitude and latitude alt (km) -- altitude ua,va,wa (m/s) -- aircraft velocity u,v,w (m/s) -- wind velocity raz,rel (degrees) -- raw azimuth and elevation of ray pitch,roll,heading,drift (degrees) -- aircraft orientation and drift */ p3rhead(buff,length,&bpointer,&raybytes,&refl,&vel,&width, &ryear,&rmon,&rday,&rhour,&rmin,&rsec,&rmsec,&lon,&lat,&alt, &ua,&va,&wa,&u,&v,&w,&raz,&rel,&pitch,&roll,&drift,&heading); /* get ray data */ /*fprintf(stderr,"\n\n\nrayBYTES = %i \n\n",raybytes); */ /*fprintf(stderr,"\n\n\nrefl = %i, vel =%i, width =%i\n\n",refl,vel,width); HERE T1*/ if ( (unc = uncomp(buff,length,&bpointer,raybytes,&outlen,OUTMAX,ubuff)) == 0) { /* fill scalar fields */ *raycnt = iraycnt/1000.; *cdfyear = ryear; *cdfmon = rmon; *cdfday = rday; *cdfhour = rhour; *cdfmin = rmin; *cdfsec = rsec; *cdfmsec = rmsec; /* compute time of day in ksec */ *cdftime = *cdfhour*3.6 + *cdfmin*0.06 + (*cdfsec + *cdfmsec*0.001)*0.001; /* fix time to be monotonic and Modify condition about time interval selection to be in the monotonic time */ *monotime = mtimefix(*cdftime,*cdfday); /*do only data within the selected times*/ if (timeselect && (*monotime > time1 && *monotime < time2)) { /* test front-back conditions for tail radar -- skip rays not desired */ if (taflag && !back && rel <= 0.) continue; if (taflag && !front && rel > 0.) continue; *cdfraz = raz + c_raz; *cdfrel = rel + c_rel; if (afilename) { /*use navigation from in situ tape if desired */ fun_use_ac_navigation(&lon,&lat,&alt,&ua,&va,&wa,&u,&v,&w,&pitch,&roll,&drift,&heading); } /* <-- AFILE */ /* fill rest of scalar fields */ *cdflon = lon; *cdflat = lat; *cdfalt = alt; *cdfua = ua; *cdfva = va; *cdfwa = wa; *cdfu = u; *cdfv = v; *cdfw = w; *cdfpitch = pitch + c_pitch; *cdfroll = roll; *cdfheading = heading + c_heading; /* MOD-Jul-7-1998 */ if (drift_flag == 0 ) *cdfdrift = drift - c_heading; else *cdfdrift = bad; /* cdfdrift = drift - c_heading; */ /* compute earth-relative elevation and azimuth */ if (taflag) { /*angles are expected in degrees and results are also in degrees*/ ta_angles(*cdfroll,*cdfpitch,*cdfheading,*cdfraz,*cdfrel,cdfaz,cdfel); } else { lf_angles(*cdfroll,*cdfpitch,*cdfheading,*cdfraz,*cdfrel,cdfaz,cdfel); } /* fill radar data fields */ /* log reflectivity */ nvars = refl + vel + width; vscale = nyv/128.; if (nvars < 1 || nvars > 3) { fprintf(stderr,"hrdrad: nvars = %d is out of range\n",nvars); exit(1); } if (refl) { /* reflectivity range is encoded as follow: value meaning 0 = no reflectivity data available 1 = -31.5 dbz 64 = 0 dbz 255 = +95.5 dbz with 0.5 dbz steps. Therefore, the reflectivity in decibels is given by DBZ = (value - 64)/2 */ poop = 0; for (loop = 0; loop < cdfgates; loop++) { if (gateaverage && loop < 64) poop += 2; value = ubuff[poop*nvars]; if (gateaverage && loop < 128) poop++; poop++; /* convert dbz to reflectivity/1000 */ dbzvalue = (value - 64)/2.; cdfdbz[loop] = (value == 0 ? bad : dbzvalue); } } else { for (loop = 0; loop < cdfgates; loop++) cdfdbz[loop] = bad; } if (taflag) { /*Tail radar*/ if (vel) { /* velocity is scaled respect to nyquist velocity and it is encoded as follow: value meaning 0 = no velocity data available 1 = max velocity towards radar 128 = zero velocity 255 = max velocity away from radar */ poop = 0; for (loop = 0; loop < cdfgates; loop++) { if (gateaverage && loop < 64) poop += 2; value = ubuff[1 + poop*nvars]; if (gateaverage && loop < 128) poop++; poop++; cdfvr[loop] = (value == 0 ? bad : vscale*(value - 128)); } } else { for (loop = 0; loop < cdfgates; loop++) cdfvr[loop] = bad; } if (width) { /* spectral width range from 1/256 to 255/256 0 means width data was not available at this range. To convert from width to m/s multyply by the nyquist velocity */ poop = 0; for (loop = 0; loop < cdfgates; loop++) { if (gateaverage && loop < 64) poop += 2; value = ubuff[2 + poop*nvars]; if (gateaverage && loop < 128) poop++; poop++; /* cdfwidth[loop] = (value == 0 ? bad : nyv*(value/256)); */ cdfwidth[loop] = (nyv*(value/256)); } } else { for (loop = 0; loop < cdfgates; loop++) cdfwidth[loop] = bad; } } /*if tail radar*/ /* write a variable slice */ /* but include nyquist velocity*/ *nyqvel=nyv; putslice(stdout,nvbuff,vbuff); /* increment ray count */ iraycnt++; } /* <-- TIMESELECT */ else { /*check if we have completed time interval*/ if (timeselect && (*monotime > time2)) break; /*out from loop on raysss*/ } } else { /* break out from loop on raysss and go to the next record */ fprintf(stderr,"UNC= %i\n",unc); fprintf(stderr,"hrdrad: bad ray at raycnt = %g -- skip record\n", *raycnt); break; } } /*loop in rays*/ } /*do only data from selected radar */ } else { /* continue; record reading*/ fprintf(stderr,"hrdrad: data record before header record -- skip!\n"); } } /*process data record*/ } else { /* do as follows for zero length records */ p3file++; if (headflag) { fprintf(stderr, "hrdrad: zero length record in file %d at time %.3f, raycnt %.3f\n", p3file,*cdftime,*raycnt); } } /*check if we have completed time interval*/ if (timeselect && (*monotime > time2)) break; /*out of record reading loop*/ } /*read records loop*/ /* end */ fprintf(stderr,"hrdrad: ending on record of length %d\n",length); fflush(stderr); exit(0); } /*EO_main*/ /*--------------------------------------------------------------------*/ /* Function to calculate the elevation and azimuth angles of a vector with respect to a horizontal plane over the Earth from the pitch, roll, heading, elevation and azimuth of the same vector with respect to the WP-3 aircraft reference frame. The azimuth and elevation are for the tail radar. Author -- Lucio Lopez */ /* This function expects the angles in degrees and it also calculates az and el in degrees. */ ta_angles(roll,pitch,hdg,raz,rel,az,el) float roll,pitch,hdg,raz,rel; /* input angles */ float *az,*el; /* output angles */ { float r,p,h,a,e; /* components of unit vector relative to horizontal reference frame */ float vex,vey,vez; /* constant to change angles from degrees to radians, 2pi/360. */ float k; k = 2.*3.14159/360.; /* transformation of angles units. */ r = k*roll; p = k*pitch; h = k*hdg; a = k*raz; e = k*rel; /* It was assumed that the order in the angle rotations is heading, pitch and roll. So, in order to transform the aircraft coordinate system to a horizontal coordinate system, we first rotate back the roll angle then the pitch angle and finally the heading angle. It is also assumed that positive pitch angles are for the aircraft nose up the horizontal, positive roll angles are for the right wing up the horizontal and the positive heading is measured from the north in a clockwise sense. */ vex = cos(e)*cos(h)*sin(a+r) - cos(e)*sin(h)*sin(p)*cos(a+r) + sin(e)*sin(h)*cos(p); vey = sin(e)*cos(h)*cos(p) - cos(e)*sin(h)*sin(a+r) - cos(e)*cos(h)*sin(p)*cos(a+r); vez = sin(p)*sin(e) + cos(p)*cos(e)*cos(a+r); /* now get the elevation and azimuth angles with respect to the horizontal coordinate system. Then azimuth is positive from the north in a clockwise sense, and elevation is from the horizontal where positive is up the horizontal plane. */ *el = asin(vez)/k; *az = atan2(vex,vey)/k; if (*az < 0.) *az += 360.; } /*--------------------------------------------------------------------*/ /* Function to calculate the elevation and azimuth angles of a vector with respect to a horizontal plane over the Earth from the pitch, roll, heading, elevation and azimuth of the same vector with respect to the WP-3 aircraft reference frame. The azimuth and elevation are for the lower fuselage radar. Author -- Lucio Lopez */ /* This function expects the angles in degrees and it also calculates az and el in degrees. */ lf_angles(roll,pitch,hdg,raz,rel,az,el) float roll,pitch,hdg,raz,rel; /* input angles */ float *az,*el; /* output angles */ { float r,p,h,a,e; /* components of unit vector relative to horizontal reference frame */ float vex,vey,vez; /* constant to change angles from degrees to radians, 2pi/360. */ float k; k = 2.*3.14159/360.; /* transformation of angles units. */ r = k*roll; p = k*pitch; h = k*hdg; a = k*raz; e = k*rel; /* It was assumed that the order in the angle rotations is heading, pitch and roll. So, in order to transform the aircraft coordinate system to a horizontal coordinate system, we first rotate back the roll angle then the pitch angle and finally the heading angle. It is also assumed that positive pitch angles are for the aircraft nose up the horizontal, positive roll angles are for the right wing up the horizontal and the positive heading is measured from the north in a clockwise sense. */ vex = cos(e)*sin(a)*(cos(h)*cos(r)+sin(h)*sin(p)*sin(r)) + sin(e)*(cos(h)*sin(r) - sin(h)*sin(p)*cos(r)) + cos(e)*cos(a)*sin(h)*cos(p); vey = cos(e)*sin(a)*(cos(h)*sin(p)*sin(r)-sin(h)*cos(r)) - sin(e)*(sin(h)*sin(r) + cos(h)*sin(p)*cos(r)) + cos(e)*cos(a)*cos(h)*cos(p); vez = cos(e)*cos(a)*sin(p) - cos(e)*sin(a)*cos(p)*sin(r) + cos(p)*cos(r)*sin(e); /* now get the elevation and azimuth angles with respect to the horizontal coordinate system. Then azimuth is positive from the north in a clockwise sense, and elevation is from the horizontal where positive is up the horizontal plane. */ *el = asin(vez)/k; *az = atan2(vex,vey)/k; if (*az < 0.) *az += 360.; } /*--------------------------------------------------------------------*/ float mtimefix(float time_new,int day_new) { /* Carlos Lopez-Carrillo -- December 15,95 mtimefix -- fix time to be monotonic Increasing (in Kiloseconds) Note: (day_new != day_old) is the "new day test" date: December 21, 1995 status EXP. */ #define KSEC_DAY 86.4 #define SKIPTIMECOND -1.0 static float time_old = -1.0; /* this alway accept the 1st time and allows to set day_old on the 1st call */ static int day_old; static float shift = 0.0; /* shift in time as the day changes */ if(time_old <= time_new) { /* is time increasing */ time_old = time_new; /* actualize time_old */ day_old = day_new; /* actualize hour_old */ } else { if (day_new != day_old) { /* is it a new day */ time_old = time_new; /* actualize time_old */ day_old = day_new; /* actualize day_old */ shift += KSEC_DAY; /* set shift */ }/* Note:The actualization in time_old, MUST BE before the time fixing.Because time_old must be compared to cdftime not to monotime.(monotime is the result of the fixing). */ else return(SKIPTIMECOND); /* set skip condition */ } time_new += shift; /* fixing time */ return(time_new); } /* EOmtimefix */ /*--------------------------------------------------------------------*/ void fun_cgeta(int argc,char **argv){ char c; int l,k,j; /* initialize vars */ taflag = 0; lfflag = 0; cdfgates = 0; front = back = 0; timeselect = 0; gateaverage = 0; corr = (char*)NULL; afilename = (char*)NULL; c_raz = c_rel = c_heading = c_pitch = c_fgate = 0.; flag_acf = 0; /*Get Arguments */ while ((c = getopt(argc,argv,"tlc:x:fbas:")) != -1) { switch (c) { case 't' : taflag = 1; break; case 'l' : lfflag = 1; break; case 'x' : afilename = optarg; flag_acf=optind - 1; break; case 'f' : front = 1; break; case 'b' : back = 1; break; case 'a' : gateaverage = 1; break; case 's' : timeselect = 1; if (sscanf(optarg,"%f:%f:%f",&time1,&time2,&c_junk) != 2) { fprintf(stderr, "hrdrad: format of time selection string '%s' incorrect\n", optarg); fprintf(stderr,USAGE); exit(1); } break; case 'c' : corr = optarg; if (sscanf(optarg,"%f:%f:%f:%f:%f:%f",&c_raz,&c_rel,&c_heading, &c_pitch,&c_fgate,&c_junk) != 5) { fprintf(stderr,"hrdrad: format of correction string '%s' incorrect\n",optarg); fprintf(stderr,USAGE); exit(1); } break; default : fprintf(stderr,USAGE); exit(1); } } /*check arguments */ if (taflag + lfflag != 1) { fprintf(stderr,"hrdrad: select either -t or -l\n"); fprintf(stderr,USAGE); exit(1); } if (taflag && !front && !back) { fprintf(stderr,"hrdrad: select either -f or -b or both when -t selected\n"); fprintf(stderr,USAGE); exit(1); } if (!taflag && gateaverage) { fprintf(stderr,"hrdrad: no gate averaging on lower fuselage radar\n"); fprintf(stderr,USAGE); exit(1); } if (timeselect && (time2 < time1)) { fprintf(stderr,"time2=%f must be bigger than time1=%f",time2,time1); fprintf(stderr,USAGE); exit(1); } /* Try open aircraft in situ file -- assumed to be created by hrdair * and cdfcated into a single variable slice */ if (afilename) { if ((afile = fopen(afilename,"r")) == NULL) { fprintf(stderr,"hrdrad: can't open in situ aircraft file %s\n",afilename); exit(1); } /* read the file */ gethdr(afile,abuf,HBMAX); nsabuf = elemcnt(abuf,HBMAX,'s'); nvabuf = elemcnt(abuf,HBMAX,'v'); sabuf = getbuff(nsabuf); vabuf = getbuff(nvabuf); getslice(afile,nsabuf,sabuf); getslice(afile,nvabuf,vabuf); fclose(afile); /* make sure file has been cdfcated and has tptime as index field */ tptime = getptr2(abuf,HBMAX,sabuf,'s','r',"tptime",&fp); if (tptime == NULL || fp->dim != 1 || strcmp(fp->dname1,"tptime")) { fprintf(stderr,"hrdrad: %s doesn't have 'tptime' as index field\n", afilename); fprintf(stderr," did you forget to cdfcat this file?\n"); exit(1); } maxtime = fp->dsize1; /* get pointers to other fields of interest */ aua = getptr(abuf,HBMAX,vabuf,'v','d',"gsx"); ava = getptr(abuf,HBMAX,vabuf,'v','d',"gsy"); awa = getptr(abuf,HBMAX,vabuf,'v','d',"wgs"); au = getptr(abuf,HBMAX,vabuf,'v','d',"wx"); av = getptr(abuf,HBMAX,vabuf,'v','d',"wy"); aw = getptr(abuf,HBMAX,vabuf,'v','d',"wz"); alon = getptr(abuf,HBMAX,vabuf,'v','d',"lon"); alat = getptr(abuf,HBMAX,vabuf,'v','d',"lat"); /* get altitude either from radar altitude or geopotential altitude -- * the difference should not be important in this context */ if ((aalt = getptr(abuf,HBMAX,vabuf,'v','r',"ra")) == NULL) { aalt = getptr(abuf,HBMAX,vabuf,'v','d',"ga"); } apitch = getptr(abuf,HBMAX,vabuf,'v','d',"pc"); aroll = getptr(abuf,HBMAX,vabuf,'v','d',"rl"); aheading = getptr(abuf,HBMAX,vabuf,'v','d',"hd"); /*MOD-July-7-1998*/ if ((adrift = getptr(abuf,HBMAX,vabuf,'v','r',"da")) == NULL) drift_flag=1; else drift_flag=0; /* adrift = getptr(abuf,HBMAX,vabuf,'v','d',"da"); */ } /* end of in situ file calcs */ /* if so far so good, make a list of arguments to display in header*/ /*get string with all called arguments in it*/ if (flag_acf) { /*repalce actual file name, which is already in afilename, with "acf" */ strcpy(ACF_name,"acf"); argv[flag_acf] = ACF_name; } l = 0; for ( k=0 ; k < argc; k++) { j=0; while ( (c = argv[k][j++]) != '\0' ) { if (l < MAXCHAR) argu[l++] = c; else { fprintf(stderr,"Not enough room to hold line arguments\n"); exit(1);} } argu[l++] = ' '; } argu[l] = '\n'; fprintf(stderr,"%s\n",argu); } /*--------------------------------------------------------------------*/ void fun_candis_setup() { int loop_f; /* Assume that the number of gates with data = number of output range bins. * If cdfgates is nonzero and within the possible range, assume it already * contains the number of desired gates. For gate averaging with variable * spacing on the tail radar, the total number of gates at 300 spacing * will be 256/4 + 128/2 + 128 = 256. (Hack, hack, hack!) */ /*as it is, at this point cdfgates is zero so, cdfgates should be set to ta_num_out*/ /*fiund out how many gates there are*/ if (taflag && !gateaverage) { cdfgates = (cdfgates > 0 && cdfgates <= ta_num_out ? cdfgates : ta_num_out); } else if (taflag && gateaverage) { if (ta_vargates) { cdfgates = (cdfgates > 0 && cdfgates <= 256 ? cdfgates : 256); } else { fprintf(stderr,"hrdrad: no fixed tail gates with averaging yet\n"); exit(1); } } else { cdfgates = (cdfgates > 0 && cdfgates <= lf_num_out ? cdfgates : lf_num_out); } /* Now Make a Null header */ nullhdr(hbuf,HBMAX,"float"); /*write out calling line*/ sprintf(cline,"%s",argu); addcline(hbuf,HBMAX,cline); if (afilename) { if (strlen(afilename) < 50) sprintf(cline," acf = %s\n",afilename); else sprintf(cline,"af = toolong_to_display\n"); addcline(hbuf,HBMAX,cline); } /*write out useful comments*/ if (taflag) addcline(hbuf,HBMAX," not unfolded\n"); sprintf(cline," project = %s\n",proj_id); addcline(hbuf,HBMAX,cline); sprintf(cline," flight = %s\n",flt_id); addcline(hbuf,HBMAX,cline); sprintf(cline," aircraft = %d\n",ac_id); addcline(hbuf,HBMAX,cline); sprintf(cline," %d-%d-%d %d:%d:%d\n",year,mon,day,hour,min,sec); addcline(hbuf,HBMAX,cline); sprintf(cline," time zone - GMT = %d min\n",tzone); addcline(hbuf,HBMAX,cline); sprintf(cline," gates = %d\n",cdfgates); addcline(hbuf,HBMAX,cline); sprintf(cline," Nyquist Velocity = %4.2f \n",nyv); addcline(hbuf,HBMAX,cline); /* now include parameters */ badlim = BADLIM; bad = BAD; sprintf(pline,"%g",badlim); addpar(hbuf,HBMAX,"badlim",pline); sprintf(pline,"%g",bad); addpar(hbuf,HBMAX,"bad",pline); /* static fields */ addfld(hbuf,HBMAX,'s',"range",1000.,0.,'l',1,"range",cdfgates); addfld(hbuf,HBMAX,'s',"nyq_vel",100.,0.,'l',0); /* variable fields */ addfld(hbuf,HBMAX,'v',"raycnt",1000.,0.,'l',0); addfld(hbuf,HBMAX,'v',"year",1.,0.,'l',0); addfld(hbuf,HBMAX,'v',"mon",1.,0.,'l',0); addfld(hbuf,HBMAX,'v',"day",1.,0.,'l',0); addfld(hbuf,HBMAX,'v',"hour",1.,0.,'l',0); addfld(hbuf,HBMAX,'v',"min",1.,0.,'l',0); addfld(hbuf,HBMAX,'v',"sec",1.,0.,'l',0); addfld(hbuf,HBMAX,'v',"msec",1.,0.,'l',0); addfld(hbuf,HBMAX,'v',"lon",180.,0.,'l',0); addfld(hbuf,HBMAX,'v',"lat",180.,0.,'l',0); addfld(hbuf,HBMAX,'v',"alt",1000.,0.,'l',0); addfld(hbuf,HBMAX,'v',"ua",10.,0.,'l',0); addfld(hbuf,HBMAX,'v',"va",10.,0.,'l',0); addfld(hbuf,HBMAX,'v',"wa",10.,0.,'l',0); addfld(hbuf,HBMAX,'v',"u",10.,0.,'l',0); addfld(hbuf,HBMAX,'v',"v",10.,0.,'l',0); addfld(hbuf,HBMAX,'v',"w",10.,0.,'l',0); addfld(hbuf,HBMAX,'v',"raz",180.,0.,'l',0); addfld(hbuf,HBMAX,'v',"rel",180.,0.,'l',0); addfld(hbuf,HBMAX,'v',"pitch",180.,0.,'l',0); addfld(hbuf,HBMAX,'v',"roll",180.,0.,'l',0); addfld(hbuf,HBMAX,'v',"drift",180.,0.,'l',0); addfld(hbuf,HBMAX,'v',"heading",180.,0.,'l',0); addfld(hbuf,HBMAX,'v',"az",180.,0.,'l',0); addfld(hbuf,HBMAX,'v',"el",180.,0.,'l',0); addfld(hbuf,HBMAX,'v',"time",100.,0.,'l',0); /* adding field for monotonic increasing time */ addfld(hbuf,HBMAX,'v',"monotime",100.,0.,'l',0); addfld(hbuf,HBMAX,'v',"nyqvel",100.,0.,'l',0); addfld(hbuf,HBMAX,'v',"dbz",1.,0.,'s',1,"range",cdfgates); if (taflag) { addfld(hbuf,HBMAX,'v',"vr",1.,0.,'s',1,"range",cdfgates); addfld(hbuf,HBMAX,'v',"sw",1.,0.,'s',1,"range",cdfgates); } /* allocate space */ nsbuff = elemcnt(hbuf,HBMAX,'s'); sbuff = getbuff(nsbuff); nvbuff = elemcnt(hbuf,HBMAX,'v'); vbuff = getbuff(nvbuff); /* allocate pointers */ range = getptr(hbuf,HBMAX,sbuff,'s','d',"range"); nyq_vel = getptr(hbuf,HBMAX,sbuff,'s','d',"nyq_vel"); raycnt = getptr(hbuf,HBMAX,vbuff,'v','d',"raycnt"); cdfyear = getptr(hbuf,HBMAX,vbuff,'v','d',"year"); cdfmon = getptr(hbuf,HBMAX,vbuff,'v','d',"mon"); cdfday = getptr(hbuf,HBMAX,vbuff,'v','d',"day"); cdfhour = getptr(hbuf,HBMAX,vbuff,'v','d',"hour"); cdfmin = getptr(hbuf,HBMAX,vbuff,'v','d',"min"); cdfsec = getptr(hbuf,HBMAX,vbuff,'v','d',"sec"); cdfmsec = getptr(hbuf,HBMAX,vbuff,'v','d',"msec"); cdflon = getptr(hbuf,HBMAX,vbuff,'v','d',"lon"); cdflat = getptr(hbuf,HBMAX,vbuff,'v','d',"lat"); cdfalt = getptr(hbuf,HBMAX,vbuff,'v','d',"alt"); cdfua = getptr(hbuf,HBMAX,vbuff,'v','d',"ua"); cdfva = getptr(hbuf,HBMAX,vbuff,'v','d',"va"); cdfwa = getptr(hbuf,HBMAX,vbuff,'v','d',"wa"); cdfu = getptr(hbuf,HBMAX,vbuff,'v','d',"u"); cdfv = getptr(hbuf,HBMAX,vbuff,'v','d',"v"); cdfw = getptr(hbuf,HBMAX,vbuff,'v','d',"w"); cdfpitch = getptr(hbuf,HBMAX,vbuff,'v','d',"pitch"); cdfroll = getptr(hbuf,HBMAX,vbuff,'v','d',"roll"); cdfdrift = getptr(hbuf,HBMAX,vbuff,'v','d',"drift"); cdfheading = getptr(hbuf,HBMAX,vbuff,'v','d',"heading"); cdfraz = getptr(hbuf,HBMAX,vbuff,'v','d',"raz"); cdfrel = getptr(hbuf,HBMAX,vbuff,'v','d',"rel"); cdfaz = getptr(hbuf,HBMAX,vbuff,'v','d',"az"); cdfel = getptr(hbuf,HBMAX,vbuff,'v','d',"el"); cdftime = getptr(hbuf,HBMAX,vbuff,'v','d',"time"); /* monotonic time stuff*/ monotime = getptr(hbuf,HBMAX,vbuff,'v','d',"monotime"); nyqvel = getptr(hbuf,HBMAX,vbuff,'v','d',"nyqvel"); cdfdbz = getptr(hbuf,HBMAX,vbuff,'v','d',"dbz"); if (taflag) { cdfvr = getptr(hbuf,HBMAX,vbuff,'v','d',"vr"); cdfwidth = getptr(hbuf,HBMAX,vbuff,'v','d',"sw"); } /* zero time and raycnt */ *cdftime = *raycnt = 0.; /* compute range values in km */ /* this is a hack, since the input file doesn't contain variable range specs */ /* Dave Jorgensen assures me that this always works */ /* Gate averaging is limited to the tail radar, and results in drange = .3 km */ crange = .001*(taflag ? ta_firstgate : lf_firstgate) + c_fgate; if (gateaverage) drange = 0.3; else drange = .001*(taflag ? ta_gatestep : lf_gatestep); vrangeflag = (taflag ? ta_vargates : lf_vargates); /* add range index parameters when gate spacing is uniform */ if (!taflag || gateaverage || cdfgates <= 256 || !vrangeflag) { sprintf(cline,"%g",crange); addpar(hbuf,HBMAX,"range0",cline); sprintf(cline,"%g",drange); addpar(hbuf,HBMAX,"drange",cline); } /* do range comps */ for (loop_f = 0; loop_f < cdfgates; loop_f++) { range[loop_f] = crange; adrange = drange; if (!gateaverage && vrangeflag && loop_f >= 256) adrange = 2.* drange; if (!gateaverage && vrangeflag && loop_f >= 384) adrange = 4.* drange; crange += adrange; } /* unambiguous velocity */ *nyq_vel = nyv; /* write header and static slice to standard output */ puthdr(stdout,hbuf,HBMAX); putslice(stdout,nsbuff,sbuff); }/*EO_fun_candis_setup*/ /*--------------------------------------------------------------------*/ void fun_use_ac_navigation(float *lon,float *lat,float *alt, float *ua,float *va,float *wa,float *u,float *v,float *w, float *pitch,float *roll,float *drift,float *heading) { /*globals required but not changed here : maxtime (set in fun_cgeta), monotime */ long actime,bctime; float radartime,ahi,alo,denom; /* set time to initial time plus one */ actime = 1; /* first, find the correct time bracket */ radartime = 1000.*(*monotime); while (radartime < tptime[actime] && actime > 1) actime--; while (radartime > tptime[actime] && actime < maxtime - 1) actime++; bctime = actime - 1; /* next, find the correct interpolation factors */ denom = tptime[actime] - tptime[bctime]; if (denom < 1.e-5) { fprintf(stderr,"hrdrad: non-monotonic time in aircraft data at %f\n", tptime[actime]); exit(1); } ahi = (radartime - tptime[bctime])/denom; if (ahi > 1.01 || ahi < -0.01) { fprintf(stderr,"hrdrad: radar time %f not bracketed by %f @ %f\n", radartime,tptime[bctime],tptime[actime]); exit(1); } alo = 1. - ahi; /* do the interpolations */ *lon = ahi*alon[actime] + alo*alon[bctime]; *lat = ahi*alat[actime] + alo*alat[bctime]; *alt = (ahi*aalt[actime] + alo*aalt[bctime])/1000.; *ua = ahi*aua[actime] + alo*aua[bctime]; *va = ahi*ava[actime] + alo*ava[bctime]; *wa = ahi*awa[actime] + alo*awa[bctime]; *u = ahi*au[actime] + alo*au[bctime]; *v = ahi*av[actime] + alo*av[bctime]; *w = ahi*aw[actime] + alo*aw[bctime]; *pitch = ahi*apitch[actime] + alo*apitch[bctime]; *roll = ahi*aroll[actime] + alo*aroll[bctime]; /* drift = ahi*adrift[actime] + alo*adrift[bctime]; */ *heading = ahi*aheading[actime] + alo*aheading[bctime]; /* MOD-Jul-7-1998 */ if (drift_flag == 0 ) *drift = ahi*adrift[actime] + alo*adrift[bctime]; else *drift = bad; } /*EO_fun_use_ac_navigation*/