/*++++++++++++++ .IDENTIFICATION usnob1.c .LANGUAGE C .AUTHOR Francois Ochsenbein [CDS] .ENVIRONMENT USNO-B1.0 Catalogue .KEYWORDS CDS Catalogue Server .VERSION 1.0 16-Nov-2002 .VERSION 1.1 04-Dec-2002: Arrange the values with a binary tree .VERSION 1.2 08-Dec-2002: Use tree with 3 nodes .VERSION 1.3 15-Jan-2003: Accept just RA/Dec limits. Accept soft link .VERSION 1.4 28-Jan-2003: Accept RELATIVE soft links .VERSION 1.5 13-Sep-2003: Verify RA/Dec in correct range .VERSION 1.6 31-Jan-2004: A few errors... .COMMENTS Access to USNO-B1.0 catalogue Remember that the environment variable USNOBroot provides the root of USNOB compressed catalogue; its default is derived from the program name IF the program resides in a .../bin or .../src directory. If everything fails, the absolute default of USNOBroot is /USNOB. ---------------*/ #include /* Structure definitions */ #define RADIUS (7.5/60.) /* Default target radius */ #ifndef COMPUTE_EpPOS /* 0 by +/- 1 by Matrix 2 by Cartesian */ #define COMPUTE_EpPOS 0 #endif #include #include /* Malloc */ #include #include #include #include #define ITEMS(a) (sizeof(a)/sizeof(a[0])) #define MIN(a,b) ((a)<(b) ? (a) : (b)) #define MAX(a,b) ((a)>(b) ? (a) : (b)) #define issign(c) (((c)=='+')||((c)=='-')) #define sind(x) sin(x*(M_PI/180.)) #define cosd(x) cos(x*(M_PI/180.)) #define tand(x) tan(x*(M_PI/180.)) #define asind(x) (180./M_PI)*asin(x) #define atand(x) (180./M_PI)*atan(x) #define atan2d(x,y) (180./M_PI)*atan2(x,y) #define adeg(x) ((x)*3600*1000) /* Degree into mas */ /* Options for USNOB queries....... */ int usnob_options = 0 ; /* Linked records, to sort them... */ typedef struct linked_USNOB { struct linked_USNOB *lt ; struct linked_USNOB *eq ; struct linked_USNOB *gt ; USNOBrec rec ; } LUSNOB ; static USNOBrec rec0; /* NULL Record */ static long compare_calls; static int mrec = 100 ; static int irec, truncated ; static int input_id ; /* ID(1) Zone (2) */ static char *prompt[] = { "Center", /* input_id=0 */ "USNOB-ID", /* input_id=1 */ "USNOB Zone", /* input_id=2 */ "RA+Dec Limits", /* input_id=3 */ } ; static char whole_usnob ; /* -whole option */ LUSNOB *last, *arec, *root; /* Definitions related to Center of Field */ static double radius = 0 ; static double boxy[2] ; static int (*digest_routine)() ; static double localR[3][3] ; /* Matrix, [0] = dir.cos. of center */ static double du2max = -1 ; static char optE; /* Other options */ static int opted = 6 ; /* Default edit Decimal, Epoch, ID */ /* Definitions of the record fields */ struct s_fields { char name[4] ; /* letter */ char matching; /* 0=exact 1=1letter, 2=2letters */ char selected ; short order ; double factor ; long lim[2] ; }; static struct s_fields thefields[] = { {"r",0,0,0,3.6e6}, /* MUST BE #0 */ {"i",0,0,0,1.}, /* MUST BE #1 */ {"x",0,0,0,(180./M_PI)*3.6e6}, /* MUST BE #2 */ {"y",0,0,0,(180./M_PI)*3.6e6}, /* MUST BE #3 */ {"a",0,0,0,3.6e6}, /* Must be #4 */ {"d",0,0,0,3.6e6}, /* Must be #5 */ {"p",0,0,0,1.}, /* pmt #6 */ {"e",0,0,0,10.}, /* Mean Epoch */ {"mb1",0,0,0,100.}, /* Mag. blue 1 */ {"mb2",0,0,0,100.}, /* Mag. blue 2 */ {"mr1",0,0,0,100.}, /* Mag. red 1 */ {"mr2",0,0,0,100.}, /* Mag. red 2 */ {"mi",0,0,0,100.}, /* Mag. infrared*/ {"o",0,0,0,1.}, /* Nobs (ndet) */ {"cb1",0,0,0,1.}, /* Class cb #14 */ {"cb2",0,0,0,1.}, /* Class cb1... */ {"cr1",0,0,0,1.}, /* Class cb1... */ {"cr2",0,0,0,1.}, /* Class cb1... */ {"ci",0,0,0,1.}, /* Class cb1... */ {"t",0,0,0,USNOB_TYC}, /* Tycho-2 Star */ {""} } ; #define FIELDS_m 8 #define FIELDS_mb 8 #define FIELDS_mr 10 #define FIELDS_c 14 #define FIELDS_cb 14 #define FIELDS_cr 16 static struct s_fields *compare_fields[32] ; static struct s_fields *check_fields[32] ; static double _factor_ = 1 ; static long matched ; /* Matching records */ static int stopid; /* 1 = Stop after last ID read */ static char optEmsg[] = "\ #---------Position recomputed to Epoch."; static char usage[] = "\ Usage: usnob1 [-HELP] [-R root_name] [-E] [-r[s] [min,]radius] [-b[s] x[,y]]\n\ [-e edit_opt] [-f input_file] [-m max_records]\n\ [-c center | -i ID | -z zone] [-l! min,max] [-s !] [-whole]\n\ "; static char help[] = "\ -HELP: display column explanations\n\ -b: target box in arcmin ; -bs = target box in arcsec\n\ -r: target radius in arcmin ; -rs = target radius in arcsec\n\ -c: target center in decimal or sexagesimal (default in stdin)\n\ -e: s=edit Sexa; m=edit mas; i=edit USNOB-ID, x=edit x,y\n\ -f: specifies an input file (default stdin)\n\ -id: query from an USNOB1 (8-digit) number\n\ -m: max number of stars to retrieve\n\ -l!: Set the limits (range) on one of the parameters (below)\n\ -s!: Sort the result by the parameter ! (list below)\n\ -zxxxx: search on a zone (between 0000 and 1799)\n\ -whole: search on the whole sky\n\ -E: Compute the position of mean Epoch (apply the proper motion back)\n\ -R: Root (directory) name where the USNOB files are located ($USNOBroot)\n\ ====The abbreviations of the parameters (symbolized !)are:\n\ a=Alpha c=anyClass d=Delta e=Epoch i=USNOB1 m=anyMag o=Nobs\n\ p=pmotion r=distance s=sigmaPos x=proj.E y=proj.N z=Zone\n\ 'c' and 'm' may be follwed by the color b r, b1 b2 r1 r2 i or O J E F N\n\ "; static char HELP[] = "\ Access to USNO-B1.0 Catalog at CDS\n\ =============================================================================\n\ \n\ The fields have the following meaning:\n\ \n\ USNO-B1.0 = Sequential number in original catalogue, as ZZZZ-NNNNNNN, \n\ where ZZZZ is the zone number from 0000 to 1799\n\ Tycho-2 = Name in the 'Tycho-2 Catalog' (Hog et al., 2000, Cat. )\n\ Fir Tycho-2 stars, all parameters are derived from Tycho-2\n\ RA,Dec = Position (Alpha, Delta) at Epoch and Equinox J2000\n\ sRA,sDE = Mean Error on (RAcos(DE)) and (DE) at mean Epoch\n\ Epoch = Mean Epoch of the observations\n\ pmRA,pmDE = relative proper motions, in arcsec/yr, of the object\n\ P = probability, in units of 0.1, of the proper motion value\n\ spA,spD = mean error on proper motions pmRA and pmDE resp.\n\ Fit = fit value, on RA and DE, in units of 0.1arcsec\n\ N = Number of detections, in the range 2 to 5\n\ MsY = flags M=Motion confirmed in another catalog, \n\ s=object in diffraction spike\n\ Y=object in YS4 (in preparation)\n\ -----------------------------------------------------------------------------\n\ 2 to 5 photometries are given in blue (O or J), red (E or F), and near-IR (N)\n\ For each photometry, the following parametres are given:\n\ Xmag = magnitude derived from the photograohic material (accuracy: 0.3mag)\n\ C = calibration 0 to 9, 0=bright photometric standard on plate, \n\ 1=faint photometric standard on plate, 2=one plate away, etc.\n\ Surv. = Survey number and field number as S-FFFC:\n\ 0- and 1- are POSS-I, 2-, 3- and 7- are POSS-II, 4- is SERC-J,\n\ 5- is ESO-R, 6- is AAO-R, 8- and 9- are SERC-I.\n\ cl = classification star/galaxy, where 0 is non-stellar and 11 is stellar\n\ ++++ WARNING, this value may be unreliable ++++\n\ xi,eta= distance of photocentre from the mean position of all images\n\ -----------------------------------------------------------------------------\n\ r(\") = distance from specified center, in arc seconds.\n\ x,y(\") = position from center toward East (x) and North (y)\n\ =============================================================================\n\ For details concerning USNO-B1.0, see http://www.nofs.navy.mil/data/fchpix/\n\ =============================================================================\n\ "; /*================================================================== Find the USNOBroot parameter *==================================================================*/ static void set_root(char *pgm) /*++++++++++++++++ .PURPOSE Define the USNOBroot .RETURNS .REMARKS Take care of soft links... -----------------*/ { struct stat buf; char name[1024]; char *a, *p ; int i, r; if (!pgm) return; if (lstat(pgm, &buf)) return ; /* Error */ for (i=0; S_ISLNK(buf.st_mode); i++) { if (usnob_options&1) fprintf(stderr, "....Pgm %s ->", pgm) ; if (i>20) return; /* Give up */ i = readlink(pgm, name, sizeof(name)) ; if (i<=0) return; /* Give up */ if (i >= sizeof(name)-1) return; name[i] = 0 ; if (name[0] != '/') { /* RELATIVE NAME */ for (r=strlen(pgm)-1; (r>=0) && (pgm[r] != '/'); r--) ; ++r; if ((i+r) >= sizeof(name)-1) { fprintf(stderr, "****Too long filename %s(..)%s\n", pgm, name); return; } for (a=name+strlen(name); a>=name; a--) a[r] = a[0]; while(--r >= 0) name[r] = pgm[r]; } pgm = name ; if (usnob_options&1) fprintf(stderr, " %s\n", pgm) ; if (lstat(pgm, &buf)) return ; /* Error */ } a = malloc(12+strlen(pgm)) ; sprintf(a, "USNOBroot=%s", pgm) ; for (p = a + strlen(a) -1 ; (p>a) && (*p != '/'); p--) ; if (p>a) for (--p; (p>a) && (*p != '/'); p--) ; if ((strncmp(p, "/bin/", 5) == 0) ||(strncmp(p, "/src/", 5) == 0)) { *p = 0; putenv(a) ; } else free(a) ; } /*================================================================== Get the parameters *==================================================================*/ static struct s_fields *dup_field(struct s_fields *field) /*++++++++++++++++ .PURPOSE Clone a field .RETURNS The copy -----------------*/ { struct s_fields *f; f = (struct s_fields *)malloc(sizeof(struct s_fields)); memcpy(f, field, sizeof(struct s_fields)) ; return(f); } static struct s_fields *get_field(struct s_fields *fields, char *p) /*++++++++++++++++ .PURPOSE Interpret the meaning of a parameter .RETURNS The relevant field / NULL -----------------*/ { static char *equiv[] = { "Ob1", "Er1", "Jb2", "Fr2", "Ni", "Ii", (char *)0}; char **a, letter; struct s_fields *f; for (f=fields; f->name[0]; f++) { if (f->name[0] != *p) continue; if (f->name[1] == 0) return(f) ; if (p[1] == 0) return(f) ; letter = *p; for (a=equiv; *a; a++) { if (p[1] == **a) break; } if (*a) p = *a; if (f->name[1] != p[1]) continue; if (p[2] == 0) return(f); if (f->name[2] == p[2]) return(f); } return((struct s_fields *)0) ; } /*============================================================================*/ static void tr_ou(double o[2], double u[3]) /*++++++++++++++++ .PURPOSE Compute direction cosines .RETURNS --- -----------------*/ { u[0] = u[1] = cosd(o[1]) ; u[0] *= cosd(o[0]) ; u[1] *= sind(o[0]) ; u[2] = sind(o[1]) ; } static int tr_uo(double u[3], double o[2]) /*++++++++++++++++ .PURPOSE Compute RA+Dec from direction cosines .RETURNS 1 / 0 (x=y=z=0) -----------------*/ { double r2; /* sqrt(x*x+y*y) */ r2 = u[0]*u[0] + u[1]*u[1]; o[0] = 0.e0; if (r2 == 0.e0) { /* in case of poles */ o[0] = 999.; o[1] = 99.; if (u[2] == 0.e0) return(0); o[1] = ( u[2]>0.e0 ? 90.e0 : -90.e0); o[0] = 0.; return(1); } o[1] = atand ( u[2] / sqrt(r2)); o[0] = atan2d (u[1] , u[0] ); if (o[0] < 0.e0) o[0] += 360.e0; return (1); } static void tr_oR ( o , R ) /*++++++++++++++++ .PURPOSE Creates the rotation matrix R[3][3] defined as R[0] (first row) = unit vector towards Zenith R[1] (second row) = unit vector towards East R[2] (third row) = unit vector towards North .RETURNS --- -----------------*/ double o[2]; /* IN: original angles */ double R[3][3]; /* OUT: rotation matrix */ { R[2][2] = cosd(o[1]); R[0][2] = sind(o[1]); R[1][1] = cosd(o[0]); R[1][0] = -sind(o[0]); R[1][2] = 0.e0; R[0][0] = R[2][2] * R[1][1]; R[0][1] = -R[2][2] * R[1][0]; R[2][0] = -R[0][2] * R[1][1]; R[2][1] = R[0][2] * R[1][0]; } static void tr_uu( u1 , u2 , R ) /*++++++++++++++++ .PURPOSE Rotates the unit vector u1 to u2, as (u2) = (R) * (u1) .RETURNS --- -----------------*/ double u1[3]; /* IN: Unit vector */ double u2[3]; /* OUT: Resulting unit vector after rotation */ double R[3][3]; /* IN: rotation matrix (e.g. created by tr_oR)*/ { register int i; double u_stack[3]; /* allows same address for input/output */ for (i=0; i<3; i++) u_stack[i] = R[i][0]*u1[0] +R[i][1]*u1[1] + R[i][2]*u1[2] ; u2[0] = u_stack[0]; /* copies to output */ u2[1] = u_stack[1]; /* copies to output */ u2[2] = u_stack[2]; /* copies to output */ } static void tr_uu1( u1 , u2 , R ) /*++++++++++++++++ .PURPOSE Rotates the unit vector u1 to u2, as (u2) = (R)^-1^ * (u1) .RETURNS --- -----------------*/ double u1[3]; /* IN: Unit vector */ double u2[3]; /* OUT: Resulting unit vector after rotation */ double R[3][3]; /* IN: rotation matrix (e.g. created by tr_oR)*/ { register int i; double u_stack[3]; /* allows same address for input/output */ for (i=0; i<3; i++) u_stack[i] = R[0][i]*u1[0] +R[1][i]*u1[1] + R[2][i]*u1[2] ; u2[0] = u_stack[0]; /* copies to output */ u2[1] = u_stack[1]; /* copies to output */ u2[2] = u_stack[2]; /* copies to output */ } /*===========================================================================*/ static int match_unum (char *string, double *value) /*++++++++++++++++ .PURPOSE Analyze the incoming and get unsigned number (no E notation) .RETURNS How many bytes matched .REMARKS Initial blanks return 0 -----------------*/ { char *p; for (p=string; isdigit(*p); p++) ; if (*p == '.') for (++p; isdigit(*p); p++) ; if (p == string) *value = 0 ; else *value = _factor_ * atof(string) ; return(p-string); } static int match_num (char *string, double *value) /*++++++++++++++++ .PURPOSE Analyze the incoming and get signed number (no E notation) .RETURNS How many bytes matched .REMARKS The internal _value_ is set -----------------*/ { char *p, *s ; for (s=string; *s == ' '; s++) ; p = s ; if ((*p == '+') || (*p == '-')) p++ ; p += match_unum(p, value) ; if (*s == '-') *value = -*value ; return(p-string); } static int get_lim (char *string, double values[2]) /*++++++++++++++++ .PURPOSE Find out the limits (2 numbers separated by ,) .RETURNS 0 / 1 / 2 / 3 .REMARKS Accept p=min,max=min+max ==> sqrt(min1^2+min2^2) / -----------------*/ { int found, round, n ; double val ; char *p ; found = 0 ; p = string ; values[0] = values[1] = 0 ; round = 0 ; NextRound: round ^= 1; n = match_num(p, &val) ; if (n) { found |= 1, p += n ; if (round) values[0] = val ; else values[0] = sqrt(values[0]*values[0] + val*val) ; } if (*p) { p++ ; n = match_num(p, &val) ; if (n) { found |= 2, p += n ; } if (round) values[1] = val ; else values[1] = sqrt(values[1]*values[1] + val*val) ; } if (*p && round) goto NextRound; return(found) ; } static int get_usnob (char *string, long id[2]) /*++++++++++++++++ .PURPOSE Interpret a string for USNOB1 Numbers (zone-number) .RETURNS Number of bytes matched .REMARKS Accept ZONE-* -----------------*/ { double v ; char *p ; p = string ; p += match_unum(p, &v) ; id[0] = v ; if (*p) { /* Zone is specified */ ++p; p += match_unum(p, &v) ; if (*p) v=0; /* There are non-numeric chars, e.g. 12* */ id[1] = v ; } else id[1] = id[0] ; return(p-string) ; } static int get_center (char *string, double center[2]) /*++++++++++++++++ .PURPOSE Interpret a string for RA + DEC .RETURNS Number of bytes matched -----------------*/ { double v, factor ; char *p, *s ; factor = _factor_ ; center[0] = center[1] = 0 ; p = string ; while (*p == ' ') p++ ; p += match_unum(p, center) ; while (*p == ' ') p++ ; if (*p == ',') /* For those HEASARC guys ! */ for (++p; *p == ' '; p++) ; if ((*p == '+') || (*p == '-')) ; else { /* Sexagesimal */ while ((*p == ':') || isdigit(*p)) { _factor_ /= 60.; if (*p == ':') p++ ; p += match_unum(p, &v); while (*p == ' ') p++ ; center[0] += v ; } center[0] *= 15. ; /* Right Ascension */ } s = p ; if ((*s == '+') || (*s == '-')) p++ ; _factor_ = 1 ; while ((*p == ':') || isdigit(*p)) { if (*p == ':') p++ ; p += match_unum(p, &v); while (*p == ' ') p++ ; center[1] += v ; _factor_ /= 60. ; } if (*s == '-') center[1] = -center[1] ; _factor_ = factor ; /* Restore */ return(p-string) ; } /*================================================================== Check routine (return 0 = don't keep, 1 = keep) *==================================================================*/ static int matching_color(char *col, USNOBrec *rec) /*++++++++++++++++ .PURPOSE Match the color 0=b1, 1=b2, etc... .RETURNS 0..4 => the required index -----------------*/ { int i; if (!col[0]) { for (i=0; i<5 && (rec->mag[i]>=9999); i++) ; if (i<5) return(i); return(0); } if (col[0] == 'r') i = 1; else if (col[0] == 'b') i = 0; else return(4); if (col[1]) /* Which color is specified */ return(i|((col[1]-'1')<<1)); if (rec->mag[i]>=9999) i += 2; return(i); } static int remove_pm(USNOBrec *rec) /*++++++++++++++++ .PURPOSE Match the color 0=b1, 1=b2, etc... .RETURNS 0 (no proper motion) / 1 .REMARKS Use matrix dx -sin(a) -sin(d)cos(a) dy = cos(a) -sin(d)sin(a) . mu1.t dz 0 +cos(d) mu2t -----------------*/ { double R[3][3], u[3], v[3], o[2], t; int ra, sd; if ((rec->pmra | rec->pmsd) == 0) /* Zero proper motion */ return(0); t = (rec->epoch - 20000)/10.; if (t == 0) return(0); /* E.g. for Tycho-2 */ o[0] = rec->ra; o[1] = rec->sd - adeg(90); o[0] /= 3.6e6; o[1] /= 3.6e6; v[1] = rec->pmra*t/3.6e6; v[2] = rec->pmsd*t/3.6e6; #if COMPUTE_EpPOS == 0 if (fabs(o[1]) > 84) { /* Close to Pole */ v[0] = 1; v[1] *= (M_PI/180.); v[2] *= (M_PI/180.); tr_oR(o, R); tr_uu1(v, u, R); tr_uo(u, o); } else { o[0] += v[1]/cosd(o[1]); if (o[0] < 0.) o[0] += 360.; if (o[0] >= 360.) o[0] -= 360.; o[1] += v[2]; } #else # if COMPUTE_EpPOS == 1 /* Use Full Matrix expression */ v[0] = 1; v[1] *= (M_PI/180.); v[2] *= (M_PI/180.); tr_oR(o, R); tr_uu1(v, u, R); # else /* Compute Cartesian diffarence */ tr_ou(o,u); R[0][1] = -sind(o[0]); R[1][1] = cosd(o[0]); R[2][1] = 0; R[0][2] = sind(o[1]); R[1][2] = R[0][2]*R[0][1]; R[2][2] = cosd(o[1]); R[0][2]*= -R[1][1]; u[0] += R[0][1]*v[1] + R[0][2]*v[2]; u[1] += R[1][1]*v[1] + R[1][2]*v[2]; u[2] += R[2][2]*v[2]; # endif tr_uo(u, o); #endif /* Convert to mas */ ra = adeg(o[0]); sd = adeg(o[1]) + adeg(90); #if 0 /* How large is the change ?? */ fprintf(stderr, "....Change for %4d-%07d (%+05.1fyr) %6d %6d\n", rec->zone, rec->id, t, ra-rec->ra, sd-rec->sd); #endif return(1); } static int check(USNOBrec *rec) /*++++++++++++++++ .PURPOSE Verify the various constraints .RETURNS 1 (OK) / 0 (does not fit) -----------------*/ { struct s_fields **f; double o[2], u[3], du, dr2 ; int i, flag ; /* Verify first the Position */ flag = 0 ; rec->xy[0] = rec->xy[1] = NULLxy ; if (du2max >= 0) { o[0] = rec->ra / 3.6e6 ; o[1] = rec->sd / 3.6e6 - 90. ; tr_ou(o, u); flag |= 1 ; du = u[0] - localR[0][0] ; dr2 = du * du ; du = u[1] - localR[0][1] ; dr2 += du * du ; du = u[2] - localR[0][2] ; dr2 += du * du ; if (dr2 > du2max) return(0) ; rec->rho = thefields[0].factor * 2.*asind(0.5*sqrt(dr2)) ; if (opted&8) { /* Compute x,y */ tr_uu(u, u, localR) ; rec->xy[0] = thefields[2].factor*u[1]/u[0] ; rec->xy[1] = thefields[3].factor*u[2]/u[0] ; flag |= 2 ; } } else rec->rho = -1 ; /* printf("#id=%4d,%8d, rho=%8d, mags=%4d,%4d,%4d\n", rec->zone,rec->id,rec->rho, rec->mB, rec->mR, rec->ci) ; */ for (f=check_fields ; *f; f++) switch((*f)->name[0]) { case_xy: if (!flag) { o[0] = rec->ra / 3.6e6 ; o[1] = rec->sd / 3.6e6 - 90. ; tr_ou(o, u); flag |= 1 ; } if (!(flag&2)) { tr_uu(u, u, localR) ; rec->xy[0] = thefields[2].factor*u[1]/u[0] ; rec->xy[1] = thefields[3].factor*u[2]/u[0] ; flag |= 2 ; } if (rec->xy[i] < (*f)->lim[0]) return(0) ; if (rec->xy[i] > (*f)->lim[1]) return(0) ; break ; case 'x': i = 0 ; goto case_xy ; case 'y': i = 1 ; goto case_xy ; case 'a': if ((*f)->lim[0] <= (*f)->lim[1]) { if (rec->ra < (*f)->lim[0]) return(0) ; if (rec->ra > (*f)->lim[1]) return(0) ; } else { if ((rec->ra < (*f)->lim[0]) && (rec->ra > (*f)->lim[1])) return(0) ; } continue ; case 'm': i = matching_color((*f)->name+1, rec); if (rec->mag[i] < (*f)->lim[0]) return(0) ; if (rec->mag[i] > (*f)->lim[1]) return(0) ; continue ; case 'c': if (rec->flags&USNOB_TYC) return(0); i = matching_color((*f)->name+1, rec); if (i>=5) return(0); if (rec->phot[i].field == 0) return(0); if (rec->phot[i].stargal < (*f)->lim[0]) return(0); if (rec->phot[i].stargal > (*f)->lim[1]) return(0); continue ; case 'd': if (rec->sd < (*f)->lim[0]) return(0) ; if (rec->sd > (*f)->lim[1]) return(0) ; continue ; case 'e': if (rec->epoch < (*f)->lim[0]) return(0) ; if (rec->epoch > (*f)->lim[1]) return(0) ; continue ; case 'i': /* lim[0] = Zone, lim[1] = ID */ if (rec->zone != (*f)->lim[0]) return(stopid) ; if ((*f)->lim[1] == 0) /* Whole Zone */ continue; if (rec->id < (*f)->lim[1]) return(0) ; if (rec->id > (*f)->lim[1]) return(stopid) ; continue ; case 'o': if (rec->ndet < (*f)->lim[0]) return(0); if (rec->ndet > (*f)->lim[1]) return(0); continue ; case 't': if ((rec->flags&USNOB_TYC) != (*f)->lim[0]) return(0); case 'p': if (rec->pmtot < 0) { u[0] = rec->pmra ; u[1] = rec->pmsd ; rec->pmtot = 0.5 + sqrt(u[0]*u[0] + u[1]*u[1]) ; } if (rec->pmtot < (*f)->lim[0]) return(0) ; if (rec->pmtot > (*f)->lim[1]) return(0) ; continue ; case 'r': if (rec->rho < (*f)->lim[0]) return(0) ; if (rec->rho > (*f)->lim[1]) return(0) ; continue ; } return(1) ; } /*================================================================== Comparison routines *==================================================================*/ static int compare(USNOBrec *a, USNOBrec *b) /*++++++++++++++++ .PURPOSE Compare two records according to the sort options .RETURNS Difference (a-b) -----------------*/ { struct s_fields **f; double u[3] ; int i, j; compare_calls++; for (f=compare_fields ; *f; f++) switch((*f)->name[0]) { case 'x': return((a->xy[0] - b->xy[0]) * (*f)->order) ; case 'y': return((a->xy[1] - b->xy[1]) * (*f)->order) ; case 'a': return((a->ra - b->ra) * (*f)->order) ; case 'c': i = matching_color((*f)->name+1, a); j = matching_color((*f)->name+1, b); return((a->phot[i].stargal - b->phot[j].stargal) * (*f)->order) ; case 'i': return((a->id - b->id) * (*f)->order) ; case 'o': return((a->ndet - b->ndet) * (*f)->order) ; case 'd': return((a->sd - b->sd) * (*f)->order) ; case 'm': i = matching_color((*f)->name+1, a); j = matching_color((*f)->name+1, b); return((a->mag[i] - b->mag[j]) * (*f)->order) ; case 'p': if (a->pmtot < 0) { u[0] = a->pmra ; u[1] = a->pmsd ; a->pmtot = 0.5 + sqrt(u[0]*u[0] + u[1]*u[1]) ; } if (b->pmtot < 0) { u[0] = b->pmra ; u[1] = b->pmsd ; b->pmtot = 0.5 + sqrt(u[0]*u[0] + u[1]*u[1]) ; } return((a->pmtot - b->pmtot) * (*f)->order) ; case 'e': return((a->epoch - b->epoch) * (*f)->order) ; case 'r': return((a->rho - b->rho) * (*f)->order) ; } return(0) ; } /*================================================================== Sorting the Records *==================================================================*/ static int add_rec(USNOBrec *new) /*++++++++++++++++ .PURPOSE Add a new record, link it. .RETURNS 0 (not kept) / 1 -----------------*/ { LUSNOB *node, *prev, *n; long comp1; int diff; comp1 = compare_calls; /* (0) Check if the new record has any use... */ n = (LUSNOB *)0; /* n -> last */ if (irec == mrec) { if (!last) for (last=root; last->gt; last=last->gt) n = last; diff = compare(new, &(last->rec)); if (diff >= 0) { truncated++; return(0); } } /* (1) Find where in the list we've to insert the new record */ node = root; prev = (LUSNOB *)0; while(node) { diff = compare(new, &(node->rec)); prev = node; if (diff == 0) break; node = diff < 0 ? node->lt : node->gt ; } /* (2) If max attained, put the new in place of last */ if (irec == mrec) { truncated++ ; node = last ; /* Where to store new record */ if (last->eq) { /* Several last records... */ node = last->eq; last->eq = node->eq ; } else { /* n->gt must give LAST record */ if (!n) for (last=root; last->gt; last=last->gt) n = last; if (!n) { /* Happens when root == last ! */ root = root->lt; last = root; } else { /* The last is now given by n */ last = n; last->gt = node->lt; } } while(last->gt) last = last->gt; } else node = arec + irec++ ; /* (3) Insert the new record, and set the links */ *(&(node->rec)) = *new ; node->lt = node->gt = node->eq = (LUSNOB *)0; if (!prev) { root = node; last = (LUSNOB *)0; } else { if (diff < 0) prev->lt = node; else if (diff > 0) { /* I may have prev == last ... */ prev->gt = node; if (last) while(last->gt) last = last->gt; } else { node->eq = prev->eq; prev->eq = node; } } #if 0 fprintf(stderr, "....Adding#%d: comparisons=%ld\n", irec, compare_calls-comp1); #endif return(1) ; } static void print_nodes(LUSNOB *node) /*++++++++++++++++ .PURPOSE Print all nodes in order .RETURNS --- .REMARKS Recursivity = simplicity !! -----------------*/ { LUSNOB *n; if (!node) return; if (node->lt) print_nodes(node->lt); for (n=node; n; n=n->eq) puts(usnob2a(&(n->rec), opted)); if (node->gt) print_nodes(node->gt); } /*================================================================== Search in a Circle *==================================================================*/ static int digest(USNOBrec *rec) /*++++++++++++++++ .PURPOSE Digest routine: check and display if necessary .RETURNS 0 -----------------*/ { int st ; /* --- Apply the Position Transformation if needed */ if (optE) remove_pm(rec); st = check(rec) ; if (st <= 0) return(st) ; matched += 1 ; if (compare_fields[0]) add_rec(rec) ; else { if (irec >= mrec) { truncated++; matched = 0; return(-1) ; } irec++ ; puts(usnob2a(rec, opted)); } return(0) ; } long usnob_loop() /*++++++++++++++++ .PURPOSE Loop on read & test of USNOB records .RETURNS Number of tested records. -----------------*/ { long tested = 0 ; USNOBrec rec ; while(1) { if (usnob_read(&rec) <= 0) break ; tested++; if (digest(&rec) < 0) break ; } return(tested) ; } long usnob_pos(long gra[2], long gsd[2]) /*++++++++++++++++ .PURPOSE Launch Search on USNOB stars from limits in RA + DE .RETURNS Number of tested records. .REMARKS Merge with specified RA / DE limits -----------------*/ { long ra[2], sd[2], sra[4], tra[4], tested ; int i, j ; tested = 0 ; if ((!thefields[5].selected) && (!thefields[4].selected)) return(usnob_search(gra, gsd, digest)) ; if (gsd) sd[0] = gsd[0], sd[1] = gsd[1] ; else sd[0] = 0, sd[1] = adeg(180) ; if (thefields[5].selected) { /* Limits in DE */ sd[0] = MAX(sd[0], thefields[5].lim[0]) ; sd[1] = MIN(sd[1], thefields[5].lim[1]) ; if (sd[0] > sd[1]) return(0); /* Mismatch Dec */ } sra[0] = tra[0] = 0; sra[1] = tra[1] = adeg(360)-1 ; /* Selected RA */ sra[2] = tra[2] = sra[3] = tra[3] = -1 ; if (gra) { tra[0] = gra[0]; if (gra[0] <= gra[1]) tra[1] = gra[1] ; else tra[3] = gra[1], tra[2] = 0 ; } if (thefields[4].selected) { /* Limits in RA */ sra[0] = thefields[4].lim[0] ; if (thefields[4].lim[0] <= thefields[4].lim[1]) sra[1] = thefields[4].lim[1] ; else sra[3] = thefields[4].lim[1], sra[2] = 0 ; } for (i=0; i<4; i += 2) for (j=0; j<4; j += 2) { if (sra[j] < 0) continue ; if (tra[i] < 0) continue ; ra[0] = MAX(tra[i], sra[j]) ; ra[1] = MIN(tra[i+1], sra[j+1]) ; if (ra[0] > ra[1]) continue ; tested += usnob_search(ra, sd, digest) ; } return(tested) ; } long usnob_center(double center[2]) /*++++++++++++++++ .PURPOSE Search USNOB stars from a central position: - either within circle of radius (degrees) - or within box of half-dimensions boxy .RETURNS Number of tested records. -----------------*/ { long ra[2], sd[2]; unsigned char *p ; double value, sr, cosdec, da, dd, maxrad ; ra[0] = 0, ra[1] = adeg(360) - 1 ; dd = 180. ; if ((radius==0) && (boxy[1]==0)) radius = RADIUS ; if (radius > 0) dd = maxrad = radius ; if (boxy[1] > 0) /* Search in a Box */ dd = maxrad = sqrt(boxy[0]*boxy[0] + boxy[1]*boxy[1]) ; /* Don't accept too large searches... */ if (dd > 45.) { fprintf(stderr, "****Radius [%.5f]deg or Box [%.5fx%.5f] too large\n", maxrad, boxy[0], boxy[1]) ; exit(1) ; } /* Derive first the limits in Dec and RA */ value = center[1] - dd ; sd[0] = (90.+value) * 3.6e6 ; value = center[1] + dd ; sd[1] = (90.+value) * 3.6e6 ; if (sd[0] < 0) sd[0] = 0 ; if (sd[1] > adeg(180)) sd[1] = adeg(180) ; cosdec = cosd(center[1]) ; sr = sind(maxrad) ; du2max = 2. * sind(maxrad/2.) ; /* Compute limits on RA */ if (sr < cosdec) { /* Pole not included in circle */ da = asind(sr/cosdec) ; ra[0] = (center[0] - da)*3.6e6 ; if (ra[0] < 0) ra[0] += adeg(360) ; ra[1] = (center[0] + da)*3.6e6 ; ra[1] %= adeg(360) ; } /* Define the constants */ tr_oR(center, localR) ; du2max *= du2max ; /* Display the Constraints */ if (usnob_options&1) { printf("#....usnob_center(%.5f, %.5f) maxrad=%.5f, limits_mas:\n", center[0], center[1], maxrad) ; printf("#.... ra=[%ld,%ld]", ra[0], ra[1]) ; printf(" sd=[%ld,%ld]\n", sd[0], sd[1]) ; } /* Merge the -la & -ld limits, and Launch Search */ return(usnob_pos(ra, sd)) ; } /*================================================================== Main Program *==================================================================*/ /* Signal: if interrupted, return the signal number */ #include void OnSignal(signo) { if (isatty(2)) fprintf(stderr, "\n****Signal #%d received, STOP\n", signo) ; exit(signo) ; } main (argc, argv) int argc; char **argv; { char line[BUFSIZ], *p, *a ; double center[2] ; double flims[2] ; char *the_center, *pgm, argline[BUFSIZ] ; struct s_fields *f, *fc ; int goon = 1 ; long tested, along, id ; int i, sign ; int non_flagged_arg = 0 ; int limradec = 0; USNOBrec rec ; FILE *input_file ; if (argc < 2) { fprintf(stderr, "%s%s", usage, help) ; exit(1) ; } for (i=1; i<35; i++) signal(i, OnSignal) ; /* Keep program name to define the default USNOB root name */ pgm = argv[0] ; /* fprintf(stderr, "....pgm=%s\n", pgm ? pgm : "(nil)") ; */ /* NOTE: The options can also be RA-DEC Radius */ the_center = line ; input_file = stdin ; while (--argc > 0) { p = *++argv; if (*p == '-') switch(p[1]) { case 'H': /* HELP */ printf("%s", HELP); exit(0) ; case 'h': /* Help */ fprintf(stderr, "%s%s", usage, help) ; exit(0) ; case 'f': /* Next parameter = input file */ argc--, argv++; input_file = fopen(*argv, "r") ; if (!input_file) { perror(*argv) ; exit(1) ; } continue ; case 'b': /* Next parameter = Box limits */ argc--, argv++; switch(p[2]) { case 'd': _factor_ = 1; break ; case 's': _factor_ = 1./3600. ; break ; default: _factor_ = 1./60. ; break ; } if (get_lim(*argv, flims) < 2) flims[1] = flims[0] ; _factor_ = 1.; flims[0] /= 2.0 ; flims[1] /= 2.0 ; boxy[0] = flims[0] ; boxy[1] = flims[1] ; thefields[2].lim[1] = thefields[2].factor*tand(flims[0]) ; thefields[3].lim[1] = thefields[3].factor*tand(flims[1]) ; thefields[2].lim[0] = -thefields[2].lim[1] ; thefields[3].lim[0] = -thefields[3].lim[1] ; thefields[2].selected = thefields[3].selected = 1 ; continue ; case 'E': /* Remove the proper motion effect */ optE = 1; continue ; case 'R': /* Root Name of Catalogue */ if (p[2]) p += 2; else p = *++argv, argc-- ; a = malloc(12+strlen(p)) ; sprintf(a, "USNOBroot=%s", p) ; putenv(a) ; pgm = (char *)0 ; /* Forget about implied from program name */ continue ; case 'r': /* Next parameter = Radius */ argc--, argv++; switch(p[2]) { case 'd': _factor_ = 1; break ; case 's': _factor_ = 1./3600. ; break ; default: _factor_ = 1./60. ; break ; } if (get_lim(*argv, flims) < 2) { flims[1] = flims[0] ; flims[0] = 0 ; } radius = flims[1] ; _factor_ = 1.; thefields[0].lim[0] = thefields[0].factor*flims[0] ; thefields[0].lim[1] = thefields[0].factor*flims[1] ; continue ; case 'm': /* Max number of records */ argc--, argv++; mrec = atoi(*argv) ; if (mrec < 1) mrec = 1 ; continue ; case 'e': /* Edit option */ p += 2; /* Is the argument of sort glued with -e ? */ if (!*p) p = *++argv, argc-- ; opted = 0 ; if (isdigit(*p)) opted = atoi(p) ; else while(*p) switch(*(p++)) { case '.': case ' ': case 'd': continue ; case 's': opted |= 1 ; continue ; case 'i': opted |= 2 ; continue ; case 'e': opted |= 4 ; continue ; case 'x': opted |= 8 ; continue ; case 'm': opted |= 16; continue ; default: fprintf(stderr, "****Bad argument: -e '%s'\n%s", --p, usage); exit(1) ; } continue ; case 'i': /* Input ID */ input_id = 1 ; thefields[1].selected = 1; if (argc > 1) { /* ID as argument */ stopid = -1 ; /* Indicates to stop after last ID read */ argc--, argv++; the_center = strdup(*argv) ; } continue ; case 'z': /* A Zone */ input_id = 2 ; if (argc > 1) { /* Zone as argument */ usnob_stop(1) ; /* Must stop at EOF */ argc--, argv++; the_center = strdup(*argv) ; } continue ; case 'l': /* Set up the limits */ argc--, argv++; i = get_lim(*argv, flims); if (i==1) flims[1] = flims[0] ; if (i==2) flims[0] = -0x7fffffff ; if (!(f = get_field(thefields, p+2))) { fprintf(stderr, "****Unknown field in %s\n%s", p, usage); exit(1) ; } /* Exact name match, or just first few letters ? matching=1 when only 'm', matching=2 when 'mb' */ if (f->matching = strcmp(f->name, p+2)) f->matching = strlen(p)-2; f->selected = 1 ; f->lim[0] = flims[0] * f->factor ; f->lim[1] = flims[1] * f->factor ; if (p[2] == 'a') limradec |= 1; else if (p[2] == 'd') { /* Limits in Declination */ f->lim[0] += adeg(90), f->lim[1] += adeg(90) ; limradec |= 2; } if ((f->lim[0] > f->lim[1]) && (p[2] != 'a')) along = f->lim[0], f->lim[0] = f->lim[1], f->lim[1] = along ; continue ; case 's': /* Sort Order */ p += 2; /* Is the argument of sort glued with -s ? */ if (!*p) p = *++argv, argc-- ; for (i = 0 ; *p && (i < ITEMS(compare_fields)); i++, p++) { sign = 1 ; if (*p == '-') sign = -1, p++ ; else if (*p == '+') p++ ; f = get_field(thefields, p) ; if (!f) { fprintf(stderr, "****Unknown field in -s %s\n%s", *argv, usage) ; exit(1) ; } f = dup_field(f); f->order = sign ; /* Skip the rest of the word, until next sort */ for (sign=1; p[1] && f->name[sign] && isalnum(*p); p++) sign++; f->name[sign] = 0; /* e.g. magnitudes truncated */ compare_fields[i] = f ; } continue ; case 'c': /* Get Center */ argc--, argv++; if (argc < 1) { fprintf(stderr, "****Unspecified center in -c argument\n%s", usage) ; exit(1) ; } strcpy(argline, *argv); the_center = argline ; while (argc > 1) { /* Position as several arguments ? */ if (issign(argv[1][0])) { if (! isdigit(argv[1][1])) break ; } else if (! isdigit(argv[1][0])) break ; argc--, argv++ ; strcat(argline, " ") ; strcat(argline, *argv) ; } continue ; case 'v': /* Verbose */ usnob_options |= 1 ; continue ; case 'w': /* Whole sky */ if (strcmp(p, "-whole") == 0) { whole_usnob = 1 ; mrec = 1999999999 ; continue ; } /* NO BREAK */ default: fprintf(stderr, "****Unknown argument: %s\n%s", p, usage); exit(1) ; } /* Non-Option: if */ if (isdigit(*p) && (non_flagged_arg < 2)) { non_flagged_arg++ ; if (the_center == line) { the_center = *argv; continue ; } get_lim(*argv, flims); radius = flims[0]/60. ; continue ; } fprintf(stderr, "****Unknown argument: %s\n%s", p, usage); exit(1) ; } /* Limits in RA and DEC can be ok */ if ((input_id == 0) && (limradec == 3)) { input_id = 3; the_center = (char *)0; } /* Define USNOB from the program name */ if (!getenv("USNOBroot")) set_root(pgm); #if 0 if (pgm && (!getenv("USNOBroot"))) { a = malloc(12+strlen(pgm)) ; sprintf(a, "USNOBroot=%s", pgm) ; for (p = a + strlen(a) -1 ; (p>a) && (*p != '/'); p--) ; if (p>a) for (--p; (p>a) && (*p != '/'); p--) ; if ((strncmp(p, "/bin/", 5) == 0) ||(strncmp(p, "/src/", 5) == 0)) { *p = 0; putenv(a) ; } else free(a) ; } #endif /* Set up the check_fields */ for (f=thefields, i=0; f->name[0]; f++) { if (!f->selected) continue; check_fields[i++] = fc = f; /* dup_field(f); */ if (fc->matching) fc->name[fc->matching] = 0; } /* Set up the array of results */ if (compare_fields[0]) arec = (LUSNOB *)malloc(mrec*sizeof(LUSNOB)) ; while (goon) { compare_calls = 0; if (whole_usnob) goon = 0 ; else if (the_center == line) { if (isatty(0)) printf("----Give %s: ", prompt[input_id]) ; if (!fgets(line, sizeof(line), input_file)) break ; p = line + strlen(line) ; while ((p>line) && (iscntrl(p[-1]))) p-- ; *p = 0 ; if (strncasecmp(line, "quit", 4) == 0) break ; if (ispunct(line[0])) { puts(line); continue ; } } else goon = 0 ; irec = truncated = tested = matched = 0 ; if (whole_usnob) { printf("#USNOB %s\n", "(whole)") ; puts(usnob_head(opted)) ; tested = usnob_pos((long *)0, (long *)0); } else switch(input_id) { case 1: /* Get from USNOB-IDs */ usnob_options &= ~0x10 ; get_usnob(the_center, thefields[1].lim); printf("#USNOB %s\n", the_center) ; if (optE) printf("%s\n", optEmsg); puts(usnob_head(opted)) ; if (usnob_set(thefields[1].lim[0], thefields[1].lim[1]) == 0) tested = usnob_loop() ; break ; case 2: /* Get in Zone */ usnob_options &= ~0x10 ; printf("#USNOB Zone %s\n", the_center) ; puts(usnob_head(opted)) ; if (usnob_zopen(atoi(the_center)) == 0) tested = usnob_loop() ; break ; case 0: /* Get from Center */ usnob_options &= ~0x10 ; if (compare_fields[0]) { if (((compare_fields[0])->name[0] == 'r') &&((compare_fields[0])->order > 0)) usnob_options |= 0x10 ; } if (get_center(the_center, center) < 0) continue ; printf("#Center: %s\n", the_center) ; /* v1.5: Verify RA/Dec in correct range */ if ((center[0]>=360.0) || (center[1]<-90.) || (center[1]>90.)) { printf("#***bad: %s\t***(center out of limits)\n", the_center); break; } if (optE) printf("%s\n", optEmsg); puts(usnob_head(opted)) ; tested = usnob_center(center) ; break; case 3: /* Use RA+Dec Lim. */ usnob_options |= 0x10; printf("#Limits (mas) in RA[%d,%d] and SPD[%d,%d]\n", thefields[4].lim[0], thefields[4].lim[1], thefields[5].lim[0], thefields[5].lim[1]); if (optE) printf("%s\n", optEmsg); puts(usnob_head(opted)); tested = usnob_pos(thefields[4].lim, thefields[5].lim); break; } /*---- List saved records in order (replace backward by forward link) */ if (compare_fields[0]) { print_nodes(root); root = last = (LUSNOB *)0; } if (truncated) { printf("#+++Truncated to %d out of ", mrec) ; if (matched) printf("%ld", matched) ; else printf("(not-computed)") ; printf(" matches (%ld tested)\n", tested) ; } else printf("#--- %ld matches (%ld tested)\n", matched, tested); printf("#...sort required %d total comparisons\n", compare_calls); } usnob_close() ; exit(0); }