#include /* Structure definitions */ #define RADIUS (7.5/60.) /* Default target radius */ #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 adeg(x) ((x)*3600*1000) /* Degree into mas */ /* Options for GSC2 queries....... */ int gsc2_options = 0 ; /* Linked records, to sort them... */ typedef struct linked_GSC2 { struct linked_GSC2 *prev ; GSC2rec rec ; } LGSC2 ; static int mrec = 100 ; static int irec, truncated ; static int input_id ; /* ID(1) Zone (2) */ static char *prompt[] = { "Center", /* input_id=0 */ "GSC2-ID", /* input_id=1 */ "GSC2 HTM", /* input_id=2 */ } ; static char whole_gsc2 ; /* -whole option */ LGSC2 *last, *arec, *just_added; /* 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 ; /* Other options */ static int opted = 6 ; /* Default edit Decimal, Epoch, ID */ /* Definitions of the record fields */ struct s_fields { char letter ; char selected ; short order ; double factor ; double lim[2] ; }; static struct s_fields thefields[] = { {'r',0,0,1.}, /* MUST BE #0 */ {'i',0,0,1.}, /* MUST BE #1 */ {'x',0,0,(180./M_PI)}, /* MUST BE #2 */ {'y',0,0,(180./M_PI)}, /* MUST BE #3 */ {'a',0,0,(M_PI/180.)}, /* Must be #4 */ {'d',0,0,(M_PI/180.)}, /* Must be #5 */ {'c',0,0,1.}, /* Class*/ {'j',0,0,1.}, /* a */ {'e',0,0,1.}, /* Elli */ {'o',0,0,1.}, /* Ep */ {'F',0,0,1.}, /* Mag. */ {'J',0,0,1.}, /* Mag. */ {'V',0,0,1.}, /* Mag. */ {'N',0,0,1.}, /* Mag. */ {0} } ; static struct s_fields *compare_fields[16] ; static struct s_fields *check_fields[16] ; static double _factor_ = 1 ; static long matched ; /* Matching records */ static int stopid; /* 1 = Stop after last ID read */ static char usage[] = "\ Usage: gsc22 [-HELP] [-R root_name] [-r[s] [min,]radius] [-b[s] x[,y]]\n\ [-e edit_opt] [-f input_file] [-m max_records]\n\ [-c center | -g GSC2 | -z htm] [-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 GSC2-ID, e=edit ObsDate, x=edit x,y\n\ -f: specifies an input file (default stdin)\n\ -g: query from an GSC2.2 identification (e.g.N01230121234)\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\ -zxx: search on a HTM zone (e.g. N0123012)\n\ -whole: search on the whole sky\n\ -R: Root (directory) name where the GSC2 files are located ($GSC2root)\n\ ====The abbreviations of the parameters (symbolized !)are:\n\ a=Alpha d=Delta o=Epoch g=GSC2.2 r=distance x=proj.E y=proj.N \n\ [FJVN]=mag s=status c=class j=MajAxis e=Eccentr\n\ (The GSC.2. Catalogue, Version August 2001)\n\ "; static char HELP[] = "\ Access to GSC2.2 at CDS\n\ =============================================================================\n\ \n\ The fields have the following meaning:\n\ \n\ GSC2-ID = GSC2.2 Name (8-character for HTM zone + number up to 5 digits)\n\ RA+Dec = Position (Alpha, Delta) in degrees at Epoch\n\ Epoch = Epoch of positions\n\ e_RA,e_Dec= Mean errors on RA*cos(dec) and Dec (arcsec)\n\ Fmag ħmag = Magnitude in F photographic band (red), and its mean error\n\ Jmag ħmag = Magnitude in J photographic band (blue), and its mean error\n\ Vmag ħmag = Magnitude in V photographic band (green), and its mean error\n\ Nmag ħmag = Magnitude in N photographic band (near-IR), and its mean error\n\ Class = 0=star 1=Galaxy, 2=Blend, 3=Non-star, 4=?, 5=Defect\n\ Axis = Semi-major axis (pixels)\n\ ecc = Eccentricity of ellipse\n\ PA = Position angle (N thru E) of ellipse\n\ Status = Source status flag (see GSC2 documentation)\n\ r(\") = distance from specified center, in arc seconds.\n\ =============================================================================\n\ For details concerning GSC2.2, see http://ad.usno.navy.mil/gsc2/\n\ =============================================================================\n\ "; /*================================================================== Find the GSC2_ROOT parameter *==================================================================*/ static void set_root(char *pgm) /*++++++++++++++++ .PURPOSE Define the GSC2root .RETURNS .REMARKS Take care of soft links... -----------------*/ { struct stat buf; char name[1024]; char *a, *p ; int i; if (!pgm) return; if (lstat(pgm, &buf)) return ; /* Error */ for (i=0; S_ISLNK(buf.st_mode); i++) { if (gsc2_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 ; pgm = name ; if (gsc2_options&1) fprintf(stderr, " %s\n", pgm) ; if (lstat(pgm, &buf)) return ; /* Error */ } a = malloc(12+strlen(pgm)) ; sprintf(a, "GSC2root=%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) { strcpy(p, "/bindat") ; putenv(a) ; } else free(a) ; } /*================================================================== Get the parameters *==================================================================*/ static struct s_fields *get_field(struct s_fields *f, char *p) /*++++++++++++++++ .PURPOSE Interpret the meaning of a parameter .RETURNS The relevant field / NULL -----------------*/ { while (f->letter) { if (f->letter == *p) return(f) ; 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 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 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) ; } /*================================================================== Check routine (return 0 = don't keep, 1 = keep) *==================================================================*/ static int check(GSC2rec *rec) { struct s_fields **f; double o[2], u[3], cosdec, du, dr2 ; float *amag ; int i, flag ; /* Verify first the Position */ flag = 0 ; rec->xy[0] = rec->xy[1] = NULLxy ; if (du2max >= 0) { cosdec = cos(rec->de) ; /* RADIANS */ u[0] = cosdec * cos(rec->ra) ; u[1] = cosdec * sin(rec->ra) ; u[2] = sin(rec->de); 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)->letter) { case_xy: if (!flag) { cosdec = cos(rec->de) ; /* RADIANS */ u[0] = cosdec * cos(rec->ra) ; u[1] = cosdec * sin(rec->ra) ; u[2] = sin(rec->de); 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_mag: if (*amag < (*f)->lim[0]) return(0) ; if (*amag > (*f)->lim[1]) return(0) ; continue ; case 'F': amag = &(rec->Fmag); goto case_mag; case 'J': amag = &(rec->Jmag); goto case_mag; case 'V': amag = &(rec->Vmag); goto case_mag; case 'N': amag = &(rec->Nmag); goto case_mag; case 'o': if (rec->ep < (*f)->lim[0]) return(0) ; if (rec->ep > (*f)->lim[1]) return(0) ; continue; case 'c': if (rec->class < (*f)->lim[0]) return(0) ; if (rec->class > (*f)->lim[1]) return(0) ; continue ; case 'd': if (rec->de < (*f)->lim[0]) return(0) ; if (rec->de > (*f)->lim[1]) return(0) ; continue ; case 'e': if (rec->e < (*f)->lim[0]) return(0) ; if (rec->e > (*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 ; case 'j': if (rec->a < (*f)->lim[0]) return(0) ; if (rec->a > (*f)->lim[1]) return(0) ; continue ; } return(1) ; } /*================================================================== Comparison routines *==================================================================*/ static int compare(GSC2rec *a, GSC2rec *b) { struct s_fields **f; float *af, *bf ; double dif ; int r ; for (f=compare_fields ; *f; f++) { af = bf = (float *)0 ; switch((*f)->letter) { case 'x': dif = a->xy[0] - b->xy[0]; break ; case 'y': dif = a->xy[1] - b->xy[1]; break ; case 'a': dif = a->ra - b->ra; break; case 'd': dif = a->de - b->de; break; case 'o': dif = a->ep - b->ep; break; case 'r': dif = a->rho - b->rho; break; case 'c': return((a->class - b->class) * (*f)->order) ; case 'e': af = &(a->e); bf = &(b->e) ; break ; case 'j': af = &(a->a); bf = &(b->a) ; break ; case 'F': af = &(a->Fmag); bf = &(b->Fmag) ; break ; case 'J': af = &(a->Jmag); bf = &(b->Jmag) ; break ; case 'V': af = &(a->Vmag); bf = &(b->Vmag) ; break ; case 'N': af = &(a->Nmag); bf = &(b->Nmag) ; break ; } if (af) { dif = *af - *bf ; if (__isnan(dif)) { if (__isnan(*af)) { if (__isnan(*bf)) return(0) ; return(100) ; } return(-100) ; } } if (dif > 0) return((*f)->order) ; else if (dif) return(-(*f)->order) ; } return(0) ; } /*================================================================== Sorting the Records *==================================================================*/ static int add_rec(GSC2rec *new) /*++++++++++++++++ .PURPOSE Add a new record, link it. .RETURNS 0 (not kept) / 1 -----------------*/ { LGSC2 *curr, *prec, *next ; int k, n, dif, pdif ; /* (0) Check if the new record has any use... */ if ((irec == mrec) && (compare(new, &(last->rec)) >= 0)) { truncated++ ; return(0) ; } curr = prec = (LGSC2 *)0 ; /* Have to insert the new record between prec and curr... */ /* Compare with record added just before -- can be a good starting point */ if (just_added) { pdif = compare(new, &(just_added->rec)); if (pdif <= 0) curr = just_added ; if ((pdif == 0) && (just_added->prev)) prec = just_added, curr = just_added->prev; } else pdif = 9999 ; /* When there is a large number of records, use random checks to find out where the new record is to be inserted. The number of random checks is log2(Nrecs). This choice of the starting point (curr) requires that the compare function returns a value which is a scalable difference between two records. In this test, prec is set ONLY when when two records are identical. This means that no further displacement in the linked list is needed. */ if (pdif && (irec & (~255))) for (k=(irec>>3); k; k >>= 1) { n = rand()%irec ; next = &(arec[n]) ; dif = compare(new, &(next->rec)) ; /* printf("....Random #%d, dif=%d\n", n, dif) ; */ if (dif > 0) continue ; if (dif == 0) { prec = next; curr = next->prev; break; } if (curr && (dif <= pdif)) continue ; /* I'm away ! */ curr = next ; pdif = dif ; } if (!curr) curr=last ; /* (1) Find where in the list we've to insert the new record */ n = 0 ; if (!prec) while (curr && (compare(new, &(curr->rec)) < 0)) prec=curr, curr=curr->prev, n++ ; /* printf("....(tested %d links)\n", n) ; */ /* (2) If max attained, put the new in place of last */ if (irec == mrec) { truncated++ ; next = last ; if (prec == last) prec = (LGSC2*)0; else last = last->prev ; } else next = arec + irec++ ; /* (3) Insert the new record, and set the links */ just_added = next ; *(&(next->rec)) = *new ; next->prev = curr; if (prec) prec->prev = next; else last = next ; if (gsc2_options&1) { #if 1 fprintf(stderr, "----irec = %3d %8d r=%7.2f\n", irec, just_added->rec.id, just_added->rec.rho); for (prec=last; prec; prec=prec->prev) fprintf(stderr, "%8d", prec->rec.id) ; fprintf(stderr, "\n") ; #endif } return(1) ; } /*================================================================== Search in a Circle *==================================================================*/ static int digest(GSC2rec *rec) /*++++++++++++++++ .PURPOSE Digest routine: check and display if necessary .RETURNS 0 -----------------*/ { int st ; st = check(rec) ; if (st <= 0) return(st) ; matched += 1 ; if (compare_fields[0]) add_rec(rec) ; else { if (irec >= mrec) { truncated++; return(-1) ; } irec++ ; puts(gsc2a(rec, opted)); } return(0) ; } /*================================================================== 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; { long tested = 0; char *p, *a, *pgm; GSC2rec rec ; /* Keep program name to define the default GSC2 root name */ pgm = argv[0]; /* *argv[0] == '/' ? argv[0] : (char *)0 ; */ if ((!getenv("GSC2root")) && (p = getenv("GSC2_ROOT"))) { a = malloc(10+strlen(p)) ; strcpy(a, "GSC2root=") ; strcat(a, p) ; putenv(a) ; } if (!getenv("GSC2root")) set_root(pgm) ; /* Just loop on read */ if (gsc2_open(0)) exit(1); while(1) { if (gsc2_read(&rec) <= 0) break ; tested++; puts(gsc2a(&rec, opted)); } fprintf(stderr, "====Number of tested records: %ld\n", tested) ; }