Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/Makefile =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/Makefile (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/Makefile (revision 21890) @@ -0,0 +1,67 @@ +include ../../Configure +CC = gcc -g # override this for now for minimal test +HOME = $(ROOT)/src/addusno +CONFIG = $(ROOT)/config +PROGRAM = addusno + +fix: addusno needs to be brought up-to-date with other addstar/addrefs functions +install: fix + +BIN = $(HOME)/bin +INC = $(HOME)/include +SRC = $(HOME)/src +MAN = $(HOME)/doc +DESTBIN = $(LBIN) +DESTLIB = $(LLIB) +DESTINC = $(LINC) +DESTMAN = $(LMAN) + +# +INCS = -I$(INC) -I$(LINC) -I$(XINC) +LIBS = -L$(LLIB) -lFITS -lohana -lm +CFLAGS = -o $*.$(ARCH).o $(INCS) +CCFLAGS = $(INCS) $(LIBS) + +ADDUSNO = \ +$(SRC)/addusno.$(ARCH).o $(SRC)/gcatalog.$(ARCH).o \ +$(SRC)/gcatstats.$(ARCH).o $(SRC)/getusno.$(ARCH).o \ +$(SRC)/find_matches.$(ARCH).o $(SRC)/wcatalog.$(ARCH).o \ +$(SRC)/sort_lists.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/check_lockfile.$(ARCH).o \ +$(SRC)/find_proper.$(ARCH).o $(SRC)/reset_catalog.$(ARCH).o \ +$(SRC)/check_permissions.$(ARCH).o + +OBJ = $(ADDUSNO) + +default: $(PROGRAM) + +$(ADDUSNO): $(INC)/addusno.h + +# dependancy rules for binary code ########################## +$(PROGRAM): $(BIN)/$(PROGRAM).$(ARCH) + +$(BIN)/$(PROGRAM).$(ARCH): $(OBJ) + @if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi + $(CC) $(OBJ) -o $(BIN)/$(PROGRAM).$(ARCH) $(CCFLAGS) + +install: $(DESTBIN)/$(PROGRAM) + +$(DESTBIN)/$(PROGRAM): $(BIN)/$(PROGRAM).$(ARCH) + @if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN); fi + rm -f $(DESTBIN)/$(PROGRAM) + cp $(BIN)/$(PROGRAM).$(ARCH) $(DESTBIN)/$(PROGRAM) + +# utilities ################################################# +clean: + rm -f $(BIN)/*.$(ARCH) + rm -f `find . -name "*.o"` + rm -f `find . -name "*~"` + rm -f `find . -name "#*"` + +.SUFFIXES: .$(ARCH).o + +.c.$(ARCH).o: + $(CC) $(CFLAGS) -c $< + + Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/include/addstar.h =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/include/addstar.h (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/include/addstar.h (revision 21890) @@ -0,0 +1,14 @@ +# include +# include + +/* global variables set in parameter file */ +char ImageCat[256]; +char ImageTemplate[256]; +char CatTemplate[256]; +char GSCFILE[256]; +char CATDIR[256]; +double NSIGMA; +double ALPHA; +int VERBOSE; +int XOVERSCAN, YOVERSCAN; +int NEWPHOTCODE, PHOTCODE; Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/include/addusno.h =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/include/addusno.h (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/include/addusno.h (revision 21890) @@ -0,0 +1,46 @@ +# include +# include + +/* global variables set in parameter file */ +char CATDIR[256]; +int VERBOSE; +int RESET; +double RADIUS; +double PROPER; +int USNO_RED; +int USNO_BLUE; +char CDROM[256]; +int FORCE_RUN; +char PhotCodeFile[256]; + +typedef struct { + double X; + double Y; + double R; + double D; + double M, dM; + char dophot; + double sky; + double fx, fy, df; + double Mgal, Map; + int found; +} Stars; + +typedef struct { + Coords coords; + float *X, *Y; + int *N; + double RA[2], DEC[2]; + double Area, density, spacing; +} CatStats; + +typedef struct { + double R, D; + double r, b; +} USNOdata; + +typedef struct { + int *match; + float *X, *Y; + int *N; +} USNOstats; Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/ConfigInit.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/ConfigInit.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/ConfigInit.c (revision 21890) @@ -0,0 +1,25 @@ +# include "addusno.h" + +ConfigInit (int argc, char **argv) { + + char *config, *file; + + /*** load configuration info ***/ + file = SelectConfigFile (&argc, argv, "ptolemy"); + config = LoadConfigFile (file); + if (config == (char *) NULL) { + fprintf (stderr, "ERROR: can't find configuration file %s\n", file); + if (file != (char *) NULL) free (file); + exit (0); + } + if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file); + + ScanConfig (config, "PHOTCODE_FILE", "%s", 0, PhotCodeFile); + ScanConfig (config, "USNO_RADIUS", "%lf", 0, &RADIUS); + ScanConfig (config, "USNO_PROPER", "%lf", 0, &PROPER); + ScanConfig (config, "USNO_CDROM", "%s", 0, CDROM); + ScanConfig (config, "CATDIR", "%s", 0, CATDIR); + + free (config); + free (file); +} Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/addusno.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/addusno.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/addusno.c (revision 21890) @@ -0,0 +1,77 @@ +# include "addusno.h" + +main (int argc, char **argv) { + + int i, N, Nusno; + Catalog catalog; + CatStats catstats; + USNOstats usnostats; + USNOdata *usno, *getusno(); + struct timeval now, then; + PhotCodeData photcodes; + + gettimeofday (&then, (void *) NULL); + ConfigInit (argc, argv); + + VERBOSE = FALSE; + if ((i = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (i, &argc, argv); + } + FORCE_RUN = FALSE; + if ((i = get_argument (argc, argv, "-f"))) { + FORCE_RUN = TRUE; + remove_argument (i, &argc, argv); + } + RESET = FALSE; + if ((i = get_argument (argc, argv, "-reset"))) { + RESET = TRUE; + remove_argument (i, &argc, argv); + } + if (argc < 2) { + fprintf (stderr, "ERROR: Usage: addusno (region)\n"); + exit (0); + } + + if (!LoadPhotcodes (PhotCodeFile)) { + fprintf (stderr, "error loading photcodes\n"); + exit (0); + } + + USNO_RED = GetPhotcodebyName ("USNO_RED"); + if (!USNO_RED) { + fprintf (stderr, "ERROR: USNO_RED photcode not found in photcode table\n"); + exit (0); + } + USNO_BLUE = GetPhotcodebyName ("USNO_BLUE"); + if (!USNO_BLUE) { + fprintf (stderr, "ERROR: USNO_BLUE photcode not found in photcode table\n"); + exit (0); + } + + /* if lockfile exists, program will complain and quit */ + check_lockfile (); + check_permissions (argv[1]); + + gcatalog (argv[1], &catalog); + + if (RESET) reset_catalog (&catalog); + + gcatstats (&catalog, &catstats); + + usno = getusno (&usnostats, &catstats, &Nusno); + + find_matches (&catstats, &catalog, &usnostats, usno, Nusno); + + find_proper (&catstats, &catalog, &usnostats, usno, Nusno); + + wcatalog (argv[1], &catalog); + + if (VERBOSE) { + gettimeofday (&now, (void *) NULL); + fprintf (stderr, "%s: elapsed time = %.2f sec\n", argv[1], + (now.tv_sec - then.tv_sec) + 1e-6*(now.tv_usec - then.tv_usec)); + } + clear_lockfile (); + fprintf (stderr, "SUCCESS\n"); +} Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/check_lockfile.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/check_lockfile.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/check_lockfile.c (revision 21890) @@ -0,0 +1,39 @@ +# include "addusno.h" + +check_lockfile () +{ + + FILE *f; + char filename[128]; + struct stat filestat; + + sprintf (filename, "%s/lock\0", CATDIR); + if (stat (filename, &filestat) != -1) { + fprintf (stderr, "ERROR: catalog %s is locked, try again later\n", CATDIR); + exit (0); + } + + f = fopen (filename, "w"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't set lock file %s\n", filename); + exit (0); + } + fclose (f); + +} + +clear_lockfile () +{ + + char filename[128], line[256]; + struct stat filestat; + + sprintf (filename, "%s/lock\0", CATDIR); + if (stat (filename, &filestat) != -1) { + sprintf (line, "rm %s\0", filename); + system (line); + } else { + fprintf (stderr, "can't remove lockfile, why not?\n"); + } + +} Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/check_permissions.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/check_permissions.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/check_permissions.c (revision 21890) @@ -0,0 +1,63 @@ +# include "addusno.h" + +check_permissions (char *basefile) { + + FILE *f; + char *c, dir[256], filename[256]; + struct stat filestat; + uid_t fuid, uid; + gid_t fgid, gid; + int status; + + uid = getuid(); + gid = getgid(); + + /* check permission to write to directory */ + sprintf (filename, "%s\0", basefile); + c = strrchr (filename, '/'); + if (c == (char *) NULL) { + strcpy (dir, "."); + } else { + *c = 0; + strcpy (dir, filename); + } + status = stat (dir, &filestat); + if (status == -1) { + fprintf (stderr, "ERROR: can't write to %s\n", dir); + exit (0); + } + if (((uid == filestat.st_uid) && (filestat.st_mode & S_IRWXU)) || + ((gid == filestat.st_gid) && (filestat.st_mode & S_IRWXG)) || + (filestat.st_mode & S_IRWXO)) { + } else { + fprintf (stderr, "ERROR: can't write to %s\n", dir); + exit (0); + } + + /* check permission to write to file */ + sprintf (filename, "%s\0", basefile); + status = stat (filename, &filestat); + if (status == 0) { /* file exists, are permissions OK? */ + if (((uid == filestat.st_uid) && (filestat.st_mode & S_IRUSR) && (filestat.st_mode & S_IWUSR)) || + ((gid == filestat.st_gid) && (filestat.st_mode & S_IRGRP) && (filestat.st_mode & S_IWGRP)) || + ((filestat.st_mode & S_IROTH) && (filestat.st_mode & S_IWOTH))) { + } else { + fprintf (stderr, "ERROR: can't write to %s\n", filename); + exit (0); + } + } + + /* check permission to write to backup file */ + sprintf (filename, "%s~\0", basefile); + status = stat (filename, &filestat); + if (status == 0) { /* file exists, are permissions OK? */ + if (((uid == filestat.st_uid) && (filestat.st_mode & S_IRUSR) && (filestat.st_mode & S_IWUSR)) || + ((gid == filestat.st_gid) && (filestat.st_mode & S_IRGRP) && (filestat.st_mode & S_IWGRP)) || + ((filestat.st_mode & S_IROTH) && (filestat.st_mode & S_IWOTH))) { + } else { + fprintf (stderr, "ERROR: can't write to %s\n", filename); + exit (0); + } + } + +} Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/config.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/config.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/config.c (revision 21890) @@ -0,0 +1,110 @@ +# include +# define D_NBYTES 4096 + +char *LoadConfigFile (filename) +char *filename; +{ + + FILE *f; + int i, done, Nbytes, NBYTES, size; + char *config; + + f = fopen (filename, "r"); + if (f == NULL) { + fprintf (stderr, "couldn't find %s\n", filename); + return ((char *) NULL); + } + + NBYTES = D_NBYTES; + ALLOCATE (config, char, NBYTES); + + size = 0; + done = FALSE; + for (i = 0; !done; i++) { + Nbytes = Fread (&config[i*D_NBYTES], sizeof(char), D_NBYTES, f, "char"); + size += Nbytes; + if (Nbytes < D_NBYTES) + done = TRUE; + else { + NBYTES += D_NBYTES; + REALLOCATE (config, char, NBYTES); + } + } + + config[size] = '\n'; + config[size+1] = 0; + REALLOCATE (config, char, size + 2); + + fclose (f); + + return (config); +} + +int ScanConfig /* we expect one more field: the pointer to the value requested */ +# ifndef ANSI +(config, field, mode, va_alist) + char *config; char *field, *mode; va_dcl +# else +(char config[], + char field[], + char mode[],...) +# endif +{ + + int i, j; + char *p, *p2, tmp[256]; + va_list argp; + double value; + +# ifndef ANSI + va_start (argp); +# else + va_start (argp, N); +# endif + + /* find the correct line with field */ + p2 = config; + for (i = 0; ; i++) { + if (!strncmp (field, p2, strlen(field))) { + p = p2 + strlen (field); + break; + } + else { + p2 = strchr (p2, '\n'); + if (p2 == (char *) NULL) { + fprintf (stderr, "entry %s not found in config file\n", field); + return (FALSE); + } + else p2++; + } + } + if (!strcmp (mode, "%s")) { + p2 = strchr (p, '\n'); + if (p2 == (char *) NULL) + p2 = config + strlen(config); + bcopy (p, tmp, (p2-p)); + tmp[(p2-p)] = 0; + stripwhite (tmp); + p2 = va_arg (argp, char *); + strcpy (p2, tmp); + } + else { + + /* try to get a numerical value from the field */ + value = strtod (p, &p2); + if ((*p2 == 'd') || (*p2 == 'D')) + value *= pow (10.0, atof (p2 + 1)); + + if (!strcmp (mode, "%d")) *va_arg (argp, int *) = value; + if (!strcmp (mode, "%u")) *va_arg (argp, unsigned *) = value; + if (!strcmp (mode, "%ld")) *va_arg (argp, long *) = value; + if (!strcmp (mode, "%hd")) *va_arg (argp, short *) = value; + if (!strcmp (mode, "%f")) *va_arg (argp, float *) = value; + if (!strcmp (mode, "%lf")) *va_arg (argp, double *) = value; + + } + + va_end (argp); + return (TRUE); + +} Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/coordops.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/coordops.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/coordops.c (revision 21890) @@ -0,0 +1,342 @@ +# include "addstar.h" + +XY_to_RD (ra, dec, x, y, coords) +double *ra, *dec; +double x, y; +Coords coords[]; +{ + + double L, M, X, Y, T; + double R, sphi, cphi, stht, ctht; + double alpha, delta, salp, calp, sdel, sdp, cdp; + + *ra = 0; + *dec = 0; + +# if 1 + if (!strcmp(&coords[0].ctype[4], "-PLY")) { + X = coords[0].cdelt1*(x - coords[0].crpix1 + x*x*coords[0].polyterms[0][0] + x*y*coords[0].polyterms[1][0] + y*y*coords[0].polyterms[2][0] + + x*x*x*coords[0].polyterms[3][0] + x*x*y*coords[0].polyterms[4][0] + x*y*y*coords[0].polyterms[5][0] + y*y*y*coords[0].polyterms[6][0]); + Y = coords[0].cdelt2*(y - coords[0].crpix2 + x*x*coords[0].polyterms[0][1] + x*y*coords[0].polyterms[1][1] + y*y*coords[0].polyterms[2][1] + + x*x*x*coords[0].polyterms[3][1] + x*x*y*coords[0].polyterms[4][1] + x*y*y*coords[0].polyterms[5][1] + y*y*y*coords[0].polyterms[6][1]); + } else { +# else + { +# endif + X = coords[0].cdelt1*(x - coords[0].crpix1); + Y = coords[0].cdelt2*(y - coords[0].crpix2); + } + + L = (X*coords[0].pc1_1 + Y*coords[0].pc1_2); + M = (X*coords[0].pc2_1 + Y*coords[0].pc2_2); + + if (!strcmp(&coords[0].ctype[4], "-PLY") || !strcmp(&coords[0].ctype[4], "-TAN") || !strcmp(&coords[0].ctype[4], "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) { + R = hypot (L,M); + if ((L == 0) && (M == 0)) { + sphi = 0; + cphi = 1; + } + else { + sphi = L / R; + cphi = -M / R; + } + + if (!strcmp(&coords[0].ctype[4], "-PLY") || !strcmp(&coords[0].ctype[4], "-TAN")) { + if (R == 0) { + stht = 1.0; + ctht = 0.0; + } + else { + T = DEG_RAD / R; + stht = T / sqrt ( 1.0 + T*T); + ctht = 1.0 / sqrt ( 1.0 + T*T); + } + } + if (!strcmp(&coords[0].ctype[4], "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) { + ctht = RAD_DEG * R; + stht = sqrt (1 - ctht*ctht); + } + + sdp = sin(RAD_DEG*coords[0].crval2); + cdp = cos(RAD_DEG*coords[0].crval2); + + sdel = stht*sdp - ctht*cphi*cdp; + salp = ctht*sphi; + calp = stht*cdp + ctht*cphi*sdp; + alpha = atan2 (salp, calp); + delta = asin (sdel); + + *ra = DEG_RAD*alpha + coords[0].crval1; + *dec = DEG_RAD*delta; + } + if (!strcmp(&coords[0].ctype[4], "-LIN") || !strcmp(&coords[0].ctype[0], "GENE")) { + *ra = L + coords[0].crval1; + *dec = M + coords[0].crval2; + } + +} + + +fXY_to_RD (ra, dec, x, y, coords) +float *ra, *dec; +double x, y; +Coords coords[]; +{ + + double L, M, X, Y, T; + double R, sphi, cphi, stht, ctht; + double alpha, delta, salp, calp, sdel, sdp, cdp; + + *ra = 0; + *dec = 0; + + if (!strcmp(&coords[0].ctype[4], "-PLY")) { + X = coords[0].cdelt1*(x - coords[0].crpix1 + x*x*coords[0].polyterms[0][0] + x*y*coords[0].polyterms[1][0] + y*y*coords[0].polyterms[2][0] + + x*x*x*coords[0].polyterms[3][0] + x*x*y*coords[0].polyterms[4][0] + x*y*y*coords[0].polyterms[5][0] + y*y*y*coords[0].polyterms[6][0]); + Y = coords[0].cdelt2*(y - coords[0].crpix2 + x*x*coords[0].polyterms[0][1] + x*y*coords[0].polyterms[1][1] + y*y*coords[0].polyterms[2][1] + + x*x*x*coords[0].polyterms[3][1] + x*x*y*coords[0].polyterms[4][1] + x*y*y*coords[0].polyterms[5][1] + y*y*y*coords[0].polyterms[6][1]); + } else { + X = coords[0].cdelt1*(x - coords[0].crpix1); + Y = coords[0].cdelt2*(y - coords[0].crpix2); + } + L = (X*coords[0].pc1_1 + Y*coords[0].pc1_2); + M = (X*coords[0].pc2_1 + Y*coords[0].pc2_2); + + if (!strcmp(&coords[0].ctype[4], "-PLY") || !strcmp(&coords[0].ctype[4], "-TAN") || !strcmp(&coords[0].ctype[4], "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) { + R = hypot (L,M); + if ((L == 0) && (M == 0)) { + sphi = 0; + cphi = 1; + } + else { + sphi = L / R; + cphi = -M / R; + } + + if (!strcmp(&coords[0].ctype[4], "-PLY") || !strcmp(&coords[0].ctype[4], "-TAN") ) { + if (R == 0) { + stht = 1.0; + ctht = 0.0; + } + else { + T = DEG_RAD / R; + stht = T / sqrt ( 1.0 + T*T); + ctht = 1.0 / sqrt ( 1.0 + T*T); + } + } + if (!strcmp(&coords[0].ctype[4], "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) { + ctht = RAD_DEG * R; + stht = sqrt (1 - ctht*ctht); + } + + sdp = sin(RAD_DEG*coords[0].crval2); + cdp = cos(RAD_DEG*coords[0].crval2); + + sdel = stht*sdp - ctht*cphi*cdp; + salp = ctht*sphi; + calp = stht*cdp + ctht*cphi*sdp; + alpha = atan2 (salp, calp); + delta = asin (sdel); + + *ra = DEG_RAD*alpha + coords[0].crval1; + *dec = DEG_RAD*delta; + } + if (!strcmp(&coords[0].ctype[4], "-LIN") || !strcmp(&coords[0].ctype[0], "GENE")) { + *ra = L + coords[0].crval1; + *dec = M + coords[0].crval2; + } + +} + + +RD_to_XY (x, y, ra, dec, coords) +double *x, *y; +double ra, dec; +Coords coords[]; +{ + + double tmp_d; + double X, Y, sphi, cphi, stht; + double salp, calp, sdel, cdel, sdp, cdp; + + if (!strcmp(&coords[0].ctype[4], "-PLY")) { + /* fprintf (stderr, "approximate to polynomial TAN plane fit\n"); */ + } + + *x = 0; + *y = 0; + if (!strcmp(&coords[0].ctype[0], "GENE") || !strcmp(&coords[0].ctype[4], "-LIN")) { + X = (ra - coords[0].crval1); + Y = (dec - coords[0].crval2); + } + + if (!strcmp(&coords[0].ctype[4], "-PLY") || !strcmp(&coords[0].ctype[4], "-TAN") || !strcmp(&coords[0].ctype[4], "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) { + sdp = sin(RAD_DEG*coords[0].crval2); + cdp = cos(RAD_DEG*coords[0].crval2); + salp = sin(RAD_DEG*(ra - coords[0].crval1)); + calp = cos(RAD_DEG*(ra - coords[0].crval1)); + sdel = sin(RAD_DEG*dec); + cdel = cos(RAD_DEG*dec); + + stht = sdel*sdp + cdel*cdp*calp; /* sin(theta) */ + sphi = cdel*salp; /* = cos(theta)*sin(phi) */ + cphi = cdel*sdp*calp - sdel*cdp; /* = cos(theta)*cos(phi) */ + + if (!strcmp(&coords[0].ctype[4], "-PLY") || !strcmp(&coords[0].ctype[4], "-TAN") ) { + X = DEG_RAD * sphi / stht; + Y = -DEG_RAD * cphi / stht; + } + if (!strcmp(&coords[0].ctype[4], "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) { + X = DEG_RAD * sphi; + Y = -DEG_RAD * cphi; + } + X = X; + Y = Y; + } + + tmp_d = 1.0 / (coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1); + *x = tmp_d * (coords[0].pc2_2*X - coords[0].pc1_2*Y) / coords[0].cdelt1 + coords[0].crpix1; + *y = tmp_d * (coords[0].pc1_1*Y - coords[0].pc2_1*X) / coords[0].cdelt2 + coords[0].crpix2; + +} + + +fRD_to_XY (x, y, ra, dec, coords) +float *x, *y; +double ra, dec; +Coords coords[]; +{ + + double tmp_d; + double X, Y, sphi, cphi, stht; + double salp, calp, sdel, cdel, sdp, cdp; + + *x = 0; + *y = 0; + if (!strcmp(&coords[0].ctype[0], "GENE") || !strcmp(&coords[0].ctype[4], "-LIN")) { + X = (ra - coords[0].crval1); + Y = (dec - coords[0].crval2); + } + + if (!strcmp(&coords[0].ctype[4], "-PLY") || !strcmp(&coords[0].ctype[4], "-TAN") || !strcmp(&coords[0].ctype[4], "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) { + sdp = sin(RAD_DEG*coords[0].crval2); + cdp = cos(RAD_DEG*coords[0].crval2); + salp = sin(RAD_DEG*(ra - coords[0].crval1)); + calp = cos(RAD_DEG*(ra - coords[0].crval1)); + sdel = sin(RAD_DEG*dec); + cdel = cos(RAD_DEG*dec); + + stht = sdel*sdp + cdel*cdp*calp; /* sin(theta) */ + sphi = cdel*salp; /* = cos(theta)*sin(phi) */ + cphi = cdel*sdp*calp - sdel*cdp; /* = cos(theta)*cos(phi) */ + + if (!strcmp(&coords[0].ctype[4], "-PLY") || !strcmp(&coords[0].ctype[4], "-TAN") ) { + X = DEG_RAD * sphi / stht; + Y = -DEG_RAD * cphi / stht; + } + if (!strcmp(&coords[0].ctype[4], "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) { + X = DEG_RAD * sphi; + Y = -DEG_RAD * cphi; + } + X = X; + Y = Y; + } + + tmp_d = 1.0 / (coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1); + *x = tmp_d * (coords[0].pc2_2*X - coords[0].pc1_2*Y) / coords[0].cdelt1 + coords[0].crpix1; + *y = tmp_d * (coords[0].pc1_1*Y - coords[0].pc2_1*X) / coords[0].cdelt2 + coords[0].crpix2; + +} + + + +GetCoords (coords, header) +Coords coords[]; +Header header[]; +{ + + int status; + double rotate; + + rotate = 0.0; + coords[0].crval1 = coords[0].crpix1 = coords[0].cdelt1 = 0.0; + coords[0].crval2 = coords[0].crpix2 = coords[0].cdelt2 = 0.0; + coords[0].pc1_1 = coords[0].pc2_2 = 1.0; + coords[0].pc2_1 = coords[0].pc1_2 = 0.0; + strcpy (coords[0].ctype, "NONE"); + + status = TRUE; + if (fits_scan (header, "CTYPE2", "%s", 1, coords[0].ctype)) { + status = fits_scan (header, "CRVAL1", "%lf", 1, &coords[0].crval1); + status &= fits_scan (header, "CRPIX1", "%lf", 1, &coords[0].crpix1); + status &= fits_scan (header, "CRVAL2", "%lf", 1, &coords[0].crval2); + status &= fits_scan (header, "CRPIX2", "%lf", 1, &coords[0].crpix2); + if (fits_scan (header, "CDELT1", "%lf", 1, &coords[0].cdelt1)) { + status &= fits_scan (header, "CDELT2", "%lf", 1, &coords[0].cdelt2); + if (fits_scan (header, "CROTA2", "%lf", 1, &rotate)) { + coords[0].pc1_1 = cos(rotate*RAD_DEG); + coords[0].pc1_2 = -sin(rotate*RAD_DEG); + coords[0].pc2_1 = sin(rotate*RAD_DEG); + coords[0].pc2_2 = cos(rotate*RAD_DEG); + } + if (fits_scan (header, "PC001001", "%lf", 1, &coords[0].pc1_1)) { + status &= fits_scan (header, "PC001002", "%lf", 1, &coords[0].pc1_2); + status &= fits_scan (header, "PC002001", "%lf", 1, &coords[0].pc2_1); + status &= fits_scan (header, "PC002002", "%lf", 1, &coords[0].pc2_2); + } + } + else { + if (fits_scan (header, "CD1_1", "%lf", 1, &coords[0].pc1_1)) { + status &= fits_scan (header, "CD1_2", "%lf", 1, &coords[0].pc1_2); + status &= fits_scan (header, "CD2_1", "%lf", 1, &coords[0].pc2_1); + status &= fits_scan (header, "CD2_2", "%lf", 1, &coords[0].pc2_2); + coords[0].cdelt1 = coords[0].cdelt2 = 1.0; + } + else { + status = FALSE; + } + } + } + else { + if (fits_scan (header, "RA_O", "%lf", 1, &coords[0].crval1)) { + status = fits_scan (header, "RA_X", "%lf", 1, &coords[0].pc1_1); + status &= fits_scan (header, "RA_Y", "%lf", 1, &coords[0].pc1_2); + status &= fits_scan (header, "DEC_O", "%lf", 1, &coords[0].crval2); + status &= fits_scan (header, "DEC_X", "%lf", 1, &coords[0].pc2_1); + status &= fits_scan (header, "DEC_Y", "%lf", 1, &coords[0].pc2_2); + coords[0].crpix1 = coords[0].crpix2 = 0.0; + coords[0].cdelt1 = coords[0].cdelt2 = 1.0; + strcpy (coords[0].ctype, "GENE"); + } + } + if (!status) { + fprintf (stderr, "error getting all elements for coordinate mode %d\n", coords[0].ctype); + coords[0].crval1 = coords[0].crpix1 = coords[0].cdelt1 = 0.0; + coords[0].crval2 = coords[0].crpix2 = coords[0].cdelt2 = 0.0; + coords[0].pc1_1 = coords[0].pc2_2 = 1.0; + coords[0].pc2_1 = coords[0].pc1_2 = 0.0; + strcpy (coords[0].ctype, "NONE"); + } +} + + + /* -PLY projection is an extrapolation of the -TAN projection. + In addition to the usual linear terms of CRPIXi, PC00i00j, there are + higher order polynomial terms (up to 3rd order): + Axis 1 terms: + PCA1X2Y0 = coords.polyterm[0][0] = x^2 + PCA1X1Y1 = coords.polyterm[1][0] = xy + PCA1X0Y2 = coords.polyterm[2][0] = y^2 + PCA1X3Y0 = coords.polyterm[3][0] = x^3 + PCA1X2Y1 = coords.polyterm[4][0] = x^2 y + PCA1X1Y2 = coords.polyterm[5][0] = x y^2 + PCA1X0Y3 = coords.polyterm[6][0] = y^2 + Axis 2 terms: + PCA2X2Y0 = coords.polyterm[0][1] = x^2 + PCA2X1Y1 = coords.polyterm[1][1] = xy + PCA2X0Y2 = coords.polyterm[2][1] = y^2 + PCA2X3Y0 = coords.polyterm[3][1] = x^3 + PCA2X2Y1 = coords.polyterm[4][1] = x^2 y + PCA2X1Y2 = coords.polyterm[5][1] = x y^2 + PCA2X0Y3 = coords.polyterm[6][1] = y^2 + */ Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/dumpcat.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/dumpcat.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/dumpcat.c (revision 21890) @@ -0,0 +1,60 @@ +# include "addstar.h" + +main (argc, argv) +int argc; +char **argv; +{ + + int i, j, Nstars, Nimage, Nregions; + Average *average; + Measure *measure; + Missing *missing; + Catalog catalog; + GSCRegion region; + FILE *f; + + if (argc < 4) { + fprintf (stderr, "USAGE: dumpcat catalog output mode(0,1,2,3)\n"); + exit (0); + } + + strcpy (region.filename, argv[1]); + + gcatalog (®ion, &catalog); + + f = fopen (argv[2], "w"); + if (f == (FILE *) NULL) { + fprintf (stderr, "failed to open file %s for write\n", argv[2]); + clear_lockfile (); + exit (0); + } + + average = catalog.average; + for (i = 0; i < catalog.Naverage; i++) { + fprintf (f, "%10.6f %10.6f %7.3f %3d %3d %8d %d\n", average[i].R, average[i].D, + 0.001*average[i].M, average[i].Nm, average[i].Nn, average[i].offset, catalog.measure[average[i].offset].t); + if ((atof(argv[3]) == 1) || (atof(argv[3]) == 3)) { + measure = &catalog.measure[average[i].offset]; + for (j = 0; j < average[i].Nm; j++) { + fprintf (f, " %d %d %f %f %d\n", measure[j].dR, measure[j].dD, 0.001*measure[j].M, + measure[j].t, measure[j].average); + } + } + if ((atof(argv[3]) == 2) || (atof(argv[3]) == 3)) { + if (average[i].Nn > 0) { + missing = &catalog.missing[average[i].missing]; + for (j = 0; j < average[i].Nn; j++) { + fprintf (f, " %f\n", missing[j].t); + } + } else { + fprintf (f, "no missing obs\n"); + } + } + } + + + + fclose (f); + + +} Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/dumpusno.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/dumpusno.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/dumpusno.c (revision 21890) @@ -0,0 +1,144 @@ +# include +# include + +# define NZONE 24 + +int SPDzone[] = { + 0, 75, 450, 375, 1500, 1650, 300, 1425, 1725, 525, 1275, 225, + 675, 150, 600, 1575, 750, 975, 900, 1050, 1125, 1200, 825, 1350}; + +int USNOdisk[] = { + 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 5, 5, 6, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10}; + +# ifndef ALLOCATE +# define ALLOCATE(X,T,S) \ + X=(T *)malloc((unsigned) ((S)*sizeof(T)));\ + if(X==NULL) \ + { \ + fprintf(stderr,"failed to malloc X\n");\ + exit(0);\ + } +# define REALLOCATE(X,T,S) \ + X=(T *)realloc(X,(unsigned) ((S)*sizeof(T))); \ + if(X==NULL) \ + { \ + fprintf(stderr,"failed to realloc X\n"); \ + exit(0); \ + } +# endif /* ALLOCATE */ + +main (argc, argv) +int argc; +char **argv; +{ + + long int offset; + int i, bin, first, last, Nstars, nitems, Nitems, Nbins; + float hours[100]; + int start[100], number[100], *buffer, *buf; + char filename[128], c; + FILE *f; + float r, b; + double RA0, DEC0, RA1, DEC1, tmp; + int iDEC0, iDEC1, iRA0, iRA1; + int spd, spd_start, spd_end, disk; + double atof(); + + if (argc != 6) { + fprintf (stderr, "USAGE: addusno RA DEC dRA dDEC (mnt)\n"); + exit (0); + } + + RA0 = atof (argv[1]); + DEC0 = atof (argv[2]); + RA1 = RA0 + atof (argv[3]); + DEC1 = DEC0 + atof (argv[4]); + iRA0 = RA0 * 360000.0; + iRA1 = RA1 * 360000.0; + iDEC0 = (DEC0 + 90.0) * 360000.0; + iDEC1 = (DEC1 + 90.0) * 360000.0; + + spd_start = (int)((DEC0 + 90) / 7.5) * 75.0; + tmp = (DEC1 + 90) / 7.5; + if (tmp > (int)(tmp)) { + spd_end = (int)(1 + (DEC1 + 90) / 7.5) * 75.0; + } else { + spd_end = (int)(0 + (DEC1 + 90) / 7.5) * 75.0; + } + + for (spd = spd_start; spd < spd_end; spd += 75) { + disk = -1; + for (i = 0; i < NZONE; i++) { + if (spd == SPDzone[i]) + disk = USNOdisk[i]; + } + if (disk < 0) { + fprintf (stderr, "ERROR: can't find cdrom for spd %d\n", spd); + exit (0); + } + + /* load accelerator file */ + sprintf (filename, "%s/zone%04d.acc\0", argv[5], spd); + fprintf (stderr, "reading from %s\n", filename); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "can't open file %s, is cdrom %d in drive?\n", filename, disk); + fprintf (stderr, "press return when ready to continue: "); + fscanf (stdin, "%c", &c); + fprintf (stderr, "trying again...\n"); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "still can't open file %s, is cdrom %d in drive?\n", filename, disk); + exit (0); + } + } + for (i = 0; fscanf (f, "%f %d %d", &hours[i], &start[i], &number[i]) != EOF; i++); + Nbins = i; + fclose (f); + + first = RA0 / 3.75; + last = RA1 / 3.75; + if ((first >= Nbins) || (last >= Nbins)) { + fprintf (stderr, "RA out of range\n"); + exit (0); + } + + /* open data file */ + sprintf (filename, "%s/zone%04d.cat\0", argv[5], spd); + fprintf (stderr, "reading from %s\n", filename); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "can't open file %s\n", filename); + exit (0); + } + /* advance file pointer to first slice */ + offset = 3*sizeof(int)*(start[first] - 1); + fseek (f, offset, SEEK_SET); + /* on each loop, load data from an RA slice of the catalog */ + for (bin = first; bin < last + 1; bin++) { + Nitems = 3*number[bin]; + ALLOCATE (buffer, int, Nitems); + nitems = Fread (buffer, sizeof(int), Nitems, f, "int"); + if (nitems != Nitems) { + fprintf (stderr, "error reading data from file %s (%d, %d, %d)\n", + filename, start[bin], number[bin], nitems); + exit (0); + } + buf = buffer; + /* print out data from slice within RA and DEC range */ + for (i = 0; i < number[bin]; i++, buf+=3) { + if ((buf[0] > iRA0) && (buf[0] < iRA1) && + (buf[1] > iDEC0) && (buf[1] < iDEC1)) { + r = 0.1*(buf[2] - 1000*((int)(buf[2]/1000))); + b = 0.1*((int)(buf[2] - 1000000*((int)(buf[2]/1000000))) / 1000); + fprintf (stdout, "%f %f %d %f %f\n", + buf[0]/360000.0, buf[1]/360000.0 - 90.0, buf[2], r, b); + } + } + free (buffer); + } + fclose (f); + } +} + + Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/find_matches.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/find_matches.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/find_matches.c (revision 21890) @@ -0,0 +1,157 @@ +# include "addusno.h" + +find_matches (catstats, catalog, usnostats, usno, Nusno) +CatStats catstats[]; +Catalog catalog[]; +USNOstats usnostats[]; +USNOdata usno[]; +int Nusno; +{ + + int i, j, k, n, m, N, first_j; + double RADIUS2; + float *X1, *Y1, *X2, *Y2, *usnodist; + float dX, dY, dR; + int *N1, *N2, *next, last, *catmatch; + int Nave, Nmeas, NMEAS, Nmatch; + unsigned int flags; + Measure *tmpmeasure; + Coords *tcoords; + + X1 = catstats[0].X; + Y1 = catstats[0].Y; + N1 = catstats[0].N; + Nave = catalog[0].Naverage; + + X2 = usnostats[0].X; + Y2 = usnostats[0].Y; + N2 = usnostats[0].N; + + ALLOCATE (usnodist, float, Nave); + ALLOCATE (catmatch, int, Nave); + bzero (usnodist, Nave*sizeof(float)); + bzero (catmatch, Nave*sizeof(int)); + + /* set up link listed pointers for new measurements */ + Nmatch = 0; + Nmeas = catalog[0].Nmeasure; + NMEAS = Nmeas + 1000; + ALLOCATE (next, int, NMEAS); + REALLOCATE (catalog[0].measure, Measure, NMEAS); + /* set up pointers for linked list of measurements */ + for (i = 0; i < Nmeas - 1; i++) { + next[i] = i+1; + } + next[i] = -1; + last = i; + + /* choose a radius for matches */ + /* 0.02 corrects cerror to arcsec from storage units */ + RADIUS2 = RADIUS*RADIUS; + + /** find matched stars **/ + for (i = j = 0; (i < Nave) && (j < Nusno); ) { + if (catalog[0].average[N1[i]].code & ID_MOVING) { + /* this is not a star, skip */ + i++; + continue; + } + dX = X1[i] - X2[j]; + if (dX <= -2*RADIUS) { + i++; + continue; + } + if (dX >= 2*RADIUS) { + j++; + continue; + } + /* negative dX: j is too large, positive dX, i is too large */ + first_j = j; + for (; (dX > -2*RADIUS) && (j < Nusno); j++) { + dX = X1[i] - X2[j]; + dY = Y1[i] - Y2[j]; + dR = dX*dX + dY*dY; + if (dR < RADIUS2) { /* new measurement of this star */ + /** insert star in measurement list */ + n = N1[i]; /* N1 refers to the average[] list */ + N = N2[j]; /* N2 refers to the usno[] list */ + if ((catalog[0].average[n].code & ID_USNO) && (usnodist[i] < dR)) { + /* existing USNO match is closer than this new one, skip this one */ + continue; + } + if (catalog[0].average[n].code & ID_USNO) { + /* this USNO match is closer than previous one, fix usno[N].match */ + /* old one isn't a match, unset it */ + usnostats[0].match[catmatch[i]] = -1; + continue; + } + usnodist[i] = dR; + m = catalog[0].average[n].offset; /* first measurement of this star */ + for (k = 0; k < catalog[0].average[n].Nm - 1; k++) + m = next[m]; + next[Nmeas+1] = next[m]; /* insert 2 measurements in linked list */ + next[Nmeas] = Nmeas + 1; + next[m] = Nmeas; + if (next[Nmeas+1] == -1) { /* last just was moved */ + last = Nmeas+1; + } + Nmatch ++; + + /** add measurements for this star **/ + catalog[0].measure[Nmeas].dR = 360000.0*(catalog[0].average[n].R - usno[N].R); + catalog[0].measure[Nmeas].dD = 360000.0*(catalog[0].average[n].D - usno[N].D); + catalog[0].measure[Nmeas].M = 1000.0*fabs(usno[N].r); + catalog[0].measure[Nmeas].Mcal= 0; /* above measurement is exact */ + catalog[0].measure[Nmeas].dM = 100; /* error in input files stored in thousandths of mag */ + catalog[0].measure[Nmeas].t = 0; /* a flag: if 0, image is not in database */ + catalog[0].measure[Nmeas].averef = n; + catalog[0].measure[Nmeas].source = USNO_RED; + catalog[0].measure[Nmeas].dophot = 0; + catalog[0].measure[Nmeas+1].dR = catalog[0].measure[Nmeas].dR; + catalog[0].measure[Nmeas+1].dD = catalog[0].measure[Nmeas].dD; + catalog[0].measure[Nmeas+1].M = 1000.0*fabs(usno[N].b); + catalog[0].measure[Nmeas+1].Mcal= 0; /* above measurement is exact */ + catalog[0].measure[Nmeas+1].dM = 100; /* error in input files stored in thousandths of mag */ + catalog[0].measure[Nmeas+1].t = 0; /* a flag: if 0, image is not in database */ + catalog[0].measure[Nmeas+1].averef = n; + catalog[0].measure[Nmeas+1].source = USNO_BLUE; + catalog[0].measure[Nmeas+1].dophot = 0; + /* add flag in average to mark as matched with the USNO catalog */ + catalog[0].average[n].code |= ID_USNO; + /* add matched value to usno list, to keep track for next round */ + usnostats[0].match[N] = i; + catmatch[i] = N; + /* we add two measurement for each star: 1 in r, 1 in b */ + catalog[0].average[n].Nm +=2; + Nmeas +=2; + if (Nmeas == NMEAS - 2) { /* just to be safe... */ + NMEAS = Nmeas + 1000; + REALLOCATE (next, int, NMEAS); + REALLOCATE (catalog[0].measure, Measure, NMEAS); + } + } + } + j = first_j; + i++; + } + + REALLOCATE (catalog[0].measure, Measure, Nmeas); + + /* fix order of Measure (memory intensive, but fast) */ + N = 0; + ALLOCATE (tmpmeasure, Measure, Nmeas); + for (i = 0; i < Nave; i++) { + n = catalog[0].average[i].offset; + catalog[0].average[i].offset = N; + for (k = 0; k < catalog[0].average[i].Nm; k++, N++) { + tmpmeasure[N] = catalog[0].measure[n]; + n = next[n]; + } + } + free (catalog[0].measure); + catalog[0].measure = tmpmeasure; + + catalog[0].Nmeasure = Nmeas; + if (VERBOSE) fprintf (stderr, "Nusno, Nave, Nmeas: %d %d %d, (%d matches)\n", Nusno, Nave, Nmeas, Nmatch); +} + Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/find_proper.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/find_proper.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/find_proper.c (revision 21890) @@ -0,0 +1,158 @@ +# include "addusno.h" + +find_proper (catstats, catalog, usnostats, usno, Nusno) +CatStats catstats[]; +Catalog catalog[]; +USNOstats usnostats[]; +USNOdata usno[]; +int Nusno; +{ + + int i, j, k, n, m, N, first_j; + double RADIUS2, PROPER2; + float *X1, *Y1, *X2, *Y2, *usnodist; + float dX, dY, dR, dR2; + int *N1, *N2, *next, last; + int Nave, Nmeas, NMEAS, Nmatch; + unsigned int flags; + Measure *tmpmeasure; + Coords *tcoords; + int already_matched, far_enough; + + X1 = catstats[0].X; + Y1 = catstats[0].Y; + N1 = catstats[0].N; + Nave = catalog[0].Naverage; + + /* no need to do this twice!! */ + ALLOCATE (usnodist, float, Nave); + bzero (usnodist, Nave*sizeof(float)); + + X2 = usnostats[0].X; + Y2 = usnostats[0].Y; + N2 = usnostats[0].N; + + /* set up link listed pointers for new measurements */ + Nmatch = 0; + Nmeas = catalog[0].Nmeasure; + NMEAS = Nmeas + 1000; + ALLOCATE (next, int, NMEAS); + REALLOCATE (catalog[0].measure, Measure, NMEAS); + /* set up pointers for linked list of measurements */ + for (i = 0; i < Nmeas - 1; i++) { + next[i] = i+1; + } + next[i] = -1; + last = i; + + /* choose a radius for matches */ + PROPER2 = PROPER*PROPER; + RADIUS2 = RADIUS*RADIUS; + + /** find matched stars **/ + for (i = j = 0; (i < Nave) && (j < Nusno); ) { + if (catalog[0].average[N1[i]].code & ID_MOVING) { + /* this is not a star, skip */ + i++; + continue; + } + if (catalog[0].average[N1[i]].Nm < 3) { + /* may just be a noise spike, skip */ + i++; + continue; + } + if (catalog[0].average[N1[i]].code & ID_USNO) { + /* already matched with USNO, skip this one */ + i++; + continue; + } + dX = X1[i] - X2[j]; + if (dX <= -2*PROPER) { + i++; + continue; + } + if (dX >= 2*PROPER) { + j++; + continue; + } + /* negative dX: j is too large, positive dX, i is too large */ + first_j = j; + for (; (dX > -2*PROPER) && (j < Nusno); j++) { + dX = X1[i] - X2[j]; + dY = Y1[i] - Y2[j]; + dR = dX*dX + dY*dY; + if (dR < PROPER2) { + n = N1[i]; /* N1 refers to the average[] list */ + N = N2[j]; /* N2 refers to the usno[] list */ + if (usnostats[0].match[N] > -1) + continue; + if ((catalog[0].average[n].code & ID_USNO) && (usnodist[i] < dR)) { + /* existing USNO match is closer than this new one, skip this one */ + continue; + } + usnodist[i] = dR; + m = catalog[0].average[n].offset; /* first measurement of this star */ + for (k = 0; k < catalog[0].average[n].Nm - 1; k++) + m = next[m]; + next[Nmeas+1] = next[m]; /* insert 2 measurements in linked list */ + next[Nmeas] = Nmeas + 1; + next[m] = Nmeas; + if (next[Nmeas+1] == -1) { /* last just was moved */ + last = Nmeas+1; + } + Nmatch ++; + + /** add measurements for this star **/ + catalog[0].measure[Nmeas].dR = 360000.0*(catalog[0].average[n].R - usno[N].R); + catalog[0].measure[Nmeas].dD = 360000.0*(catalog[0].average[n].D - usno[N].D); + catalog[0].measure[Nmeas].M = 1000.0*fabs(usno[N].r); + catalog[0].measure[Nmeas].Mcal= 0; /* above measurement is exact */ + catalog[0].measure[Nmeas].dM = 100; /* error in input files stored in thousandths of mag */ + catalog[0].measure[Nmeas].t = 0; /* a flag: if 0, image is not in database */ + catalog[0].measure[Nmeas].averef = n; + catalog[0].measure[Nmeas].source = USNO_RED; + catalog[0].measure[Nmeas+1].dR = catalog[0].measure[Nmeas].dR; + catalog[0].measure[Nmeas+1].dD = catalog[0].measure[Nmeas].dD; + catalog[0].measure[Nmeas+1].M = 1000.0*fabs(usno[N].b); + catalog[0].measure[Nmeas+1].Mcal= 0; /* above measurement is exact */ + catalog[0].measure[Nmeas+1].dM = 100; /* error in input files stored in thousandths of mag */ + catalog[0].measure[Nmeas+1].t = 0; /* a flag: if 0, image is not in database */ + catalog[0].measure[Nmeas+1].averef = n; + catalog[0].measure[Nmeas+1].source = USNO_BLUE; + /* add flag in average to mark as matched with the USNO catalog */ + catalog[0].average[n].code |= (ID_PROPER | ID_USNO); + + /* we add two measurement for each star: 1 in r, 1 in b */ + catalog[0].average[n].Nm +=2; + Nmeas +=2; + if (Nmeas == NMEAS - 2) { /* just to be safe... */ + NMEAS = Nmeas + 1000; + REALLOCATE (next, int, NMEAS); + REALLOCATE (catalog[0].measure, Measure, NMEAS); + } + } + } + j = first_j; + i++; + } + + REALLOCATE (catalog[0].measure, Measure, Nmeas); + + /* fix order of Measure (memory intensive, but fast) */ + N = 0; + ALLOCATE (tmpmeasure, Measure, Nmeas); + for (i = 0; i < Nave; i++) { + n = catalog[0].average[i].offset; + catalog[0].average[i].offset = N; + for (k = 0; k < catalog[0].average[i].Nm; k++, N++) { + tmpmeasure[N] = catalog[0].measure[n]; + n = next[n]; + } + } + free (catalog[0].measure); + catalog[0].measure = tmpmeasure; + + catalog[0].Nmeasure = Nmeas; + if (VERBOSE) fprintf (stderr, "Nusno, Nave, Nmeas: %d %d %d, (%d matches)\n", Nusno, Nave, Nmeas, Nmatch); +} + Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/gcatalog.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/gcatalog.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/gcatalog.c (revision 21890) @@ -0,0 +1,100 @@ +# include "addusno.h" + +gcatalog (catname, catalog) +char *catname; +Catalog catalog[]; +{ + + char filename[128]; + int Nitems, nitems; + int i, Nmeas, Nmiss, done; + FILE *f; + struct tm *local; + struct timeval now; + + + sprintf (filename, "%s/%s\0", CATDIR, catname); + + /* read catalog header */ + if (!fits_read_header (filename, &catalog[0].header)) { + fprintf (stderr, "ERROR: can't open catalog file: %s\n", filename); + exit (0); + } + + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't open catalog file: %s\n", filename); + exit (0); + } + fseek (f, catalog[0].header.size, SEEK_SET); + + /** find number of stars, measurements **/ + catalog[0].Nmissing = catalog[0].Naverage = catalog[0].Nmeasure = 0; + fits_scan (&catalog[0].header, "NSTARS", "%d", 1, &catalog[0].Naverage); + fits_scan (&catalog[0].header, "NMEAS", "%d", 1, &catalog[0].Nmeasure); + fits_scan (&catalog[0].header, "NMISS", "%d", 1, &catalog[0].Nmissing); + ALLOCATE (catalog[0].average, Average, MAX (catalog[0].Naverage, 1)); + ALLOCATE (catalog[0].measure, Measure, MAX (catalog[0].Nmeasure, 1)); + ALLOCATE (catalog[0].missing, Missing, MAX (catalog[0].Nmissing, 1)); + + done = FALSE; + fits_scan (&catalog[0].header, "ADDUSNO", "%t", 1, &done); + if (done && !FORCE_RUN) { + if (VERBOSE) fprintf (stderr, "SUCCESS: file already up-to-date %s\n", filename); + fclose (f); + clear_lockfile (); + exit (0); + } + + if (catalog[0].Naverage == 0) { + if (VERBOSE) fprintf (stderr, "no stars yet in catalog %s\n", filename); + clear_lockfile (); + fclose (f); + exit (0); + } + + /* read average values */ + Nitems = catalog[0].Naverage; + nitems = Fread (catalog[0].average, sizeof(Average), Nitems, f, "average"); + if (nitems != Nitems) { + fprintf (stderr, "ERROR: failed to read data from catalog file %s (1)\n", filename); + fclose (f); + exit (0); + } + + /* read measurements */ + Nitems = catalog[0].Nmeasure; + nitems = Fread (catalog[0].measure, sizeof(Measure), Nitems, f, "measure"); + if (nitems != Nitems) { + fprintf (stderr, "ERROR: failed to read data from catalog file %s (2)\n", filename); + fclose (f); + exit (0); + } + + /* read missing */ + Nitems = catalog[0].Nmissing; + nitems = Fread (catalog[0].missing, sizeof(Missing), Nitems, f, "missing"); + if (nitems != Nitems) { + fprintf (stderr, "ERROR: failed to read data from catalog file %s (3)\n", filename); + fclose (f); + exit (0); + } + + if (VERBOSE) fprintf (stderr, "read %d stars from catalog file %s (%d measurements, %d missing)\n", + catalog[0].Naverage, filename, catalog[0].Nmeasure, catalog[0].Nmissing); + + for (i = Nmeas = Nmiss = 0; i < catalog[0].Naverage; i++) { + Nmeas += catalog[0].average[i].Nm; + Nmiss += catalog[0].average[i].Nn; + } + if ((Nmeas != catalog[0].Nmeasure) || (Nmiss != catalog[0].Nmissing)) { + fprintf (stderr, "ERROR: data in catalog %s is corrupt, sums don't check\n"); + fprintf (stderr, "ERROR: Nmeas: %d, %d\n", Nmeas, catalog[0].Nmeasure); + fprintf (stderr, "ERROR: Nmiss: %d, %d\n", Nmiss, catalog[0].Nmissing); + exit (0); + } + + fclose (f); + +} + Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/gcatstats.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/gcatstats.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/gcatstats.c (revision 21890) @@ -0,0 +1,64 @@ +# include "addusno.h" + +gcatstats (catalog, catstats) +Catalog catalog[]; +CatStats catstats[]; +{ + + int i; + double RaCenter, DecCenter; + double MinRA, MaxRA, MinDEC, MaxDEC; + float *X1, *Y1; + int *N1; + Coords *tcoords; + + fits_scan (&catalog[0].header, "RA0", "%lf", 1, &MinRA); + fits_scan (&catalog[0].header, "RA1", "%lf", 1, &MaxRA); + fits_scan (&catalog[0].header, "DEC0", "%lf", 1, &MinDEC); + fits_scan (&catalog[0].header, "DEC1", "%lf", 1, &MaxDEC); + + /* + MinRA = MaxRA = catalog[0].average[0].R; + MinDEC = MaxDEC = catalog[0].average[0].D; + for (i = 0; i < catalog[0].Naverage; i++) { + MinRA = MIN (catalog[0].average[i].R, MinRA); + MaxRA = MAX (catalog[0].average[i].R, MaxRA); + MinDEC = MIN (catalog[0].average[i].D, MinDEC); + MaxDEC = MAX (catalog[0].average[i].D, MaxDEC); + } + */ + /* double check on region RA and DEC ranges */ + DecCenter = 0.5*(MinDEC + MaxDEC); + RaCenter = 0.5*(MinRA + MaxRA); + + catstats[0].RA[0] = MinRA; + catstats[0].RA[1] = MaxRA; + catstats[0].DEC[0] = MinDEC; + catstats[0].DEC[1] = MaxDEC; + + /** allocate local arrays **/ + ALLOCATE (catstats[0].X, float, catalog[0].Naverage); + ALLOCATE (catstats[0].Y, float, catalog[0].Naverage); + ALLOCATE (catstats[0].N, int, catalog[0].Naverage); + + /* project onto rectilinear grid with 1 arcsec pixels, sort by X */ + /* reference for coords is center of field */ + catstats[0].coords.crval1 = RaCenter; + catstats[0].coords.crval2 = DecCenter; + catstats[0].coords.crpix1 = catstats[0].coords.crpix2 = 0.0; + catstats[0].coords.cdelt1 = catstats[0].coords.cdelt2 = 1.0 / 3600.0; + catstats[0].coords.pc1_1 = catstats[0].coords.pc2_2 = 1.0; + catstats[0].coords.pc1_2 = catstats[0].coords.pc2_1 = 0.0; + strcpy (catstats[0].coords.ctype, "RA---TAN"); + + X1 = catstats[0].X; + Y1 = catstats[0].Y; + tcoords = &catstats[0].coords; + for (i = 0; i < catalog[0].Naverage; i++, X1++, Y1++) { + fRD_to_XY (X1, Y1, catalog[0].average[i].R, catalog[0].average[i].D, tcoords); + catstats[0].N[i] = i; + } + if (catalog[0].Naverage > 1) sort_lists (catstats[0].X, catstats[0].Y, catstats[0].N, catalog[0].Naverage); + +} + Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/getusno.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/getusno.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/getusno.c (revision 21890) @@ -0,0 +1,156 @@ +# include "addusno.h" +# define NZONE 24 + +int SPDzone[] = { + 0, 75, 450, 375, 1500, 1650, 300, 1425, 1725, 525, 1275, 225, + 675, 150, 600, 1575, 750, 975, 900, 1050, 1125, 1200, 825, 1350}; + +int USNOdisk[] = { + 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 5, 5, 6, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10}; + +USNOdata *getusno (usnostats, catstats, Nusno) +CatStats catstats[]; +USNOstats usnostats[]; +int *Nusno; +{ + + long int offset; + int i, bin, first, last, Nstars, nitems, Nitems, Nbins, *N1; + float hours[100]; + int start[100], number[100], *buffer, *buf; + char filename[128], c; + FILE *f; + float r, b, *X1, *Y1; + double RA0, DEC0, RA1, DEC1; + int iDEC0, iDEC1, iRA0, iRA1; + int spd, spd_start, spd_end, disk; + int NUSNO, nusno; + USNOdata *usno; + Coords *tcoords; + + iRA0 = catstats[0].RA[0] * 360000.0; + iRA1 = catstats[0].RA[1] * 360000.0; + iDEC0 = (catstats[0].DEC[0] + 90.0) * 360000.0; + iDEC1 = (catstats[0].DEC[1] + 90.0) * 360000.0; + + spd_start = (int)( (catstats[0].DEC[0] + 90) / 7.5) * 75.0; + DEC1 = (catstats[0].DEC[1] + 90) / 7.5; + if (DEC1 > (int)(DEC1)) { + spd_end = (int)(1 + (catstats[0].DEC[1] + 90) / 7.5) * 75.0; + } else { + spd_end = (int)(0 + (catstats[0].DEC[1] + 90) / 7.5) * 75.0; + } + + NUSNO = 5000; + ALLOCATE (usno, USNOdata, NUSNO); + nusno = 0; + + for (spd = spd_start; spd < spd_end; spd += 75) { + disk = -1; + for (i = 0; i < NZONE; i++) { + if (spd == SPDzone[i]) + disk = USNOdisk[i]; + } + if (disk < 0) { + fprintf (stderr, "ERROR: can't find cdrom for spd %d\n", spd); + exit (0); + } + + /* load accelerator file */ + sprintf (filename, "%s/zone%04d.acc\0", CDROM, spd); + fprintf (stderr, "reading from %s\n", filename); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "can't open file %s, is cdrom %d in drive?\n", filename, disk); + fprintf (stderr, "press return when ready to continue: "); + fscanf (stdin, "%c", &c); + fprintf (stderr, "trying again...\n"); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "still can't open file %s, is cdrom %d in drive?\n", filename, disk); + exit (0); + } + } + for (i = 0; fscanf (f, "%f %d %d", &hours[i], &start[i], &number[i]) != EOF; i++); + Nbins = i; + fclose (f); + + first = catstats[0].RA[0] / 3.75; + if ((catstats[0].RA[1] / 3.75) == (int) (catstats[0].RA[1] / 3.75)) + last = catstats[0].RA[1] / 3.75; + else + last = 1 + catstats[0].RA[1] / 3.75; + + if ((first > Nbins) || (last > Nbins)) { + fprintf (stderr, "RA out of range\n"); + exit (0); + } + + /* open data file */ + sprintf (filename, "%s/zone%04d.cat\0", CDROM, spd); + fprintf (stderr, "reading from %s\n", filename); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "can't open file %s\n", filename); + exit (0); + } + /* advance file pointer to first slice */ + offset = 3*sizeof(int)*(start[first] - 1); + fseek (f, offset, SEEK_SET); + /* on each loop, load data from an RA slice of the catalog */ + for (bin = first; bin < last; bin++) { + Nitems = 3*number[bin]; + ALLOCATE (buffer, int, Nitems); + nitems = Fread (buffer, sizeof(int), Nitems, f, "int"); + if (nitems != Nitems) { + fprintf (stderr, "error reading data from file %s (%d, %d, %d)\n", + filename, start[bin], number[bin], nitems); + exit (0); + } + buf = buffer; + /* print out data from slice within RA and DEC range */ + for (i = 0; i < number[bin]; i++, buf+=3) { + if ((buf[0] > iRA0) && (buf[0] < iRA1) && + (buf[1] > iDEC0) && (buf[1] < iDEC1)) { + usno[nusno].R = buf[0]/360000.0; + usno[nusno].D = buf[1]/360000.0 - 90.0; + usno[nusno].r = 0.1*(buf[2] - 1000*((int)(buf[2]/1000))); + usno[nusno].b = 0.1*((int)(buf[2] - 1000000*((int)(buf[2]/1000000))) / 1000); + nusno ++; + if (nusno == NUSNO - 1) { + NUSNO += 5000; + REALLOCATE (usno, USNOdata, NUSNO); + } + } + } + free (buffer); + } + fclose (f); + } + + ALLOCATE (usnostats[0].X, float, nusno); + ALLOCATE (usnostats[0].Y, float, nusno); + ALLOCATE (usnostats[0].N, int, nusno); + ALLOCATE (usnostats[0].match, int, nusno); + + X1 = usnostats[0].X; + Y1 = usnostats[0].Y; + N1 = usnostats[0].N; + + /* project onto rectilinear grid with 1 arcsec pixels, sort by X */ + /* reference for coords is center of field */ + tcoords = &catstats[0].coords; + for (i = 0; i < nusno; i++) { + fRD_to_XY (&X1[i], &Y1[i], usno[i].R, usno[i].D, tcoords); + N1[i] = i; + usnostats[0].match[i] = -1; + } + if (nusno > 1) sort_lists (X1, Y1, N1, nusno); + + REALLOCATE (usno, USNOdata, MAX (nusno, 1)); + *Nusno = nusno; + return (usno); + +} + + Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/add_unfound.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/add_unfound.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/add_unfound.c (revision 21890) @@ -0,0 +1,35 @@ + /* incorporate unmatched image stars -- after all catalogs searched... */ + for (i = 0; i < Nstars; i++) { + if (!found[i]) { + catalog[0].average[Nave].R = stars[N1[i]].R; + catalog[0].average[Nave].D = stars[N1[i]].D; + catalog[0].average[Nave].M = 1000.0*stars[N1[i]].M; + catalog[0].average[Nave].Nm = 1; + catalog[0].average[Nave].Xp = catalog[0].average[Nave].Xm = 0.0; + catalog[0].average[Nave].offset = Nmeas; + catalog[0].measure[Nmeas].dR = 0.0; + catalog[0].measure[Nmeas].dD = 0.0; + catalog[0].measure[Nmeas].M = 1000.0*stars[N1[i]].M; + catalog[0].measure[Nmeas].dM = stars[N1[i]].dM; + catalog[0].measure[Nmeas].t = 3600.0 + stars[N1[i]].Y*0.1; + catalog[0].measure[Nmeas].average = (0x00ffffff & Nave); + /* time of star = time of image + offset + Y*rate */ + Nave ++; + if (Nave == 0x1000000) { + fprintf (stderr, "too many stars in catalog\n"); + clear_lockfile (); + exit (0); + } + if (Nave == NAVE) { + NAVE = Nave + 1000; + REALLOCATE (next, int, NAVE); + REALLOCATE (catalog[0].average, Average, NAVE); + } + Nmeas ++; + if (Nmeas == NMEAS) { + NMEAS = Nmeas + 1000; + REALLOCATE (next, int, NMEAS); + REALLOCATE (catalog[0].measure, Measure, NMEAS); + } + } + } Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/addstar.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/addstar.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/addstar.c (revision 21890) @@ -0,0 +1,77 @@ +# include "addstar.h" + +main (argc, argv) +int argc; +char **argv; +{ + + FILE *f; + int i, N, Nstars, Nimage, Nregions, Nmissed; + Stars *stars, *gstars(); + Image image, *pimage, *gimages(); + Catalog catalog; + GSCRegion *region, *gregions(); + struct timeval now, then; + + gettimeofday (&then, (void *) NULL); + ConfigInit (argc, argv); + + VERBOSE = FALSE; + if ((i = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (i, &argc, argv); + } + NEWPHOTCODE = FALSE; + if ((N = get_argument (argc, argv, "-p"))) { + NEWPHOTCODE = TRUE; + remove_argument (N, &argc, argv); + PHOTCODE = atof (argv[N]); + remove_argument (N, &argc, argv); + } + if (argc < 2) { + fprintf (stderr, "ERROR: Usage: addstar filename\n"); + exit (0); + } + + /* if lockfile exists, program will complain and quit */ + check_lockfile (); + + stars = gstars (argv[1], &Nstars, &image); + + pimage = gimages (&image, &Nimage); + region = gregions (&image, &Nregions); + + for (i = 0; i < Nregions; i++) { + gcatalog (®ion[i], &catalog); + find_matches (®ion[i], stars, Nstars, &catalog, &image, pimage, Nimage); + wcatalog (®ion[i], &catalog); + } + + for (Nmissed = i = 0; i < Nstars; i++) { + if (stars[i].found == -1) { + fprintf (stderr, "%d %f %f %f %f\n", i, stars[i].R, stars[i].D, stars[i].M, stars[i].dM); + Nmissed ++; + } + } + if (Nmissed) fprintf (stderr, "WARNING: %d stars in image were missed!\n", Nmissed); + + wimage (argv[1], &image, Nstars); + + if (VERBOSE) { + gettimeofday (&now, (void *) NULL); + fprintf (stderr, "%s: elapsed time = %.2f sec\n", argv[1], + (now.tv_sec - then.tv_sec) + 1e-6*(now.tv_usec - then.tv_usec)); + } + clear_lockfile (); + fprintf (stderr, "SUCCESS\n"); +} + + + /* stars within a file are not all within 1 region file: + 1) check if unfound stars are in current region + 2) if not, keep unfound for next pass */ + + +/* input star lists should be adjusted so 1 e/sec is mag 25.000 */ + + Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/aregion.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/aregion.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/aregion.c (revision 21890) @@ -0,0 +1,188 @@ +# include "addstar.h" + +double BigDecBounds[] = {0.0, 7.5, 15.0, 22.5, 30.0, 37.5, 45.0, + 52.5, 60.0, 67.5, 75.0, 82.5, 90.0, + 0.0, -7.5, -15.0, -22.5, -30.0, -37.5, -45.0, + -52.5, -60.0, -67.5, -75.0, -82.5, -90.0}; +char *DecSections[] = {"N0000", "N0730", "N1500", "N2230", "N3000", "N3730", "N4500", + "N5230", "N6000", "N6730", "N7500", "N8230", "weirdness", + "S0000", "S0730", "S1500", "S2230", "S3000", "S3730", "S4500", + "S5230", "S6000", "S6730", "S7500", "S8230", "weirdness"}; + +char *Dec2Sections[] = {"n0000", "n0730", "n1500", "n2230", "n3000", "n3730", "n4500", + "n5230", "n6000", "n6730", "n7500", "n8230", "weirdness", + "s0000", "s0730", "s1500", "s2230", "s3000", "s3730", "s4500", + "s5230", "s6000", "s6730", "s7500", "s8230", "weirdness"}; + +char *disk[] = {"disk 1", "disk 1", "disk 1", "disk 1", "disk 1", "disk 1", "disk 1", + "disk 1", "disk 1", "disk 1", "disk 1", "disk 1", "weirdness", + "disk 1", "disk 2", "disk 2", "disk 2", "disk 2", "disk 2", "disk 2", + "disk 2", "disk 2", "disk 2", "disk 2", "disk 2", "weirdness"}; + +int NBigRASections [] = {48, 47, 45, 43, 40, 36, 32, 28, 21, 15, 9, 3, 3, 48, 47, 45, 43, 40, 36, 32, 28, 21, 15, 9, 3, 3}; + +int NDecLines[] = {593, 584, 551, 530, 522, 465, 406, 362, 280, 198, 123, 24, + 0, 597, 578, 574, 577, 534, 499, 442, 376, 294, 212, 144, 48}; + +/* find region file which contains ra, dec */ +aregion (region, f, ra, dec) +GSCRegion region[]; +FILE *f; +double ra, dec; +{ + + + char buffer[28800], temp[50], file[50]; + double dr, dd; + double RA0, RA1, DEC0, DEC1; + int i, NBigDec, NLINES, done, nregion; + + while (ra < 0) { + ra += 360.0; + } + + if (dec >= 86.25) { + sprintf (file, "%s/n8230/pole.cpt\0", CATDIR); + region[0].DEC[0] = 86.25; + region[0].DEC[1] = 93.75; + region[0].RA[0] = -180.0; + region[0].RA[1] = 540.0; + strcpy (region[0].filename, file); + return; + } + + NBigDec = -1; + for (i = 0; i < 12; i++) { +# ifdef DEBUG + fprintf (stderr, "%d %f %f %f\n", i, dec, BigDecBounds[i], BigDecBounds[i+1]); +# endif + if ((dec >= BigDecBounds[i]) && (dec < BigDecBounds[i+1])) { + NBigDec = i; + break; + } + } + if (NBigDec < 0) { + for (i = 13; i < 24; i++) { +# ifdef DEBUG + fprintf (stderr, "%d %f %f %f\n", i, dec, BigDecBounds[i], BigDecBounds[i+1]); +# endif + if ((dec < BigDecBounds[i]) && (dec >= BigDecBounds[i+1])) { + NBigDec = i; + break; + } + } + } + if (NBigDec < 0) { + fprintf (stderr, "dec out of range: %f\n", dec); + } + + NLINES = 0; + for (i = 0; i < NBigDec; i++) { + NLINES += NDecLines[i]; + } + fseek (f, 5*2880 + 48*NLINES, SEEK_SET); + + done = FALSE; + Fread (buffer, 1, 48*NDecLines[NBigDec], f, "char"); + for (i = 0; !done && (i < NDecLines[NBigDec]); i++) { + strncpy (temp, &buffer[i*48], 48); + temp[49] = 0; + hms_to_deg (&RA0, &RA1, &DEC0, &DEC1, &temp[7]); + if (RA1 < RA0) RA1 += 360.0; +# ifdef DEBUG + fprintf (stderr, "%f %f %f %f %f %f %s\n", DEC0, dec, DEC1, RA0, ra, RA1, temp); +# endif + if ((dec >= 0) && (dec >= DEC0) && (dec < DEC1) && (ra >= RA0) && (ra < RA1)) { + done = TRUE; + } + if ((dec < 0) && (dec < DEC0) && (dec >= DEC1) && (ra >= RA0) && (ra < RA1)) { + done = TRUE; + } + } + if (!done) { + fprintf (stderr, "error in search: %f %f\n", ra, dec); + clear_lockfile (); + exit (0); + } + temp[5] = 0; + sprintf (file, "%s/%s/%s.cpt\0", CATDIR, Dec2Sections[NBigDec],&temp[1]); + if (DEC0 < DEC1) { + region[0].DEC[0] = DEC0; + region[0].DEC[1] = DEC1; + } else { + region[0].DEC[0] = DEC1; + region[0].DEC[1] = DEC0; + } + region[0].RA[0] = RA0; + region[0].RA[1] = RA1; + strcpy (region[0].filename, file); + return; +} + + +/**********/ +int hms_to_deg (h0, h1, d0, d1, string) + char *string; + double *h0, *h1, *d0, *d1; +{ + + int flag_d0, flag_d1, flag_h0, flag_h1; + double tmp; + + *d0 = *h0 = *d1 = *h1 = 0; + + flag_h0 = dparse (h0, 1, string); + flag_h1 = dparse (h1, 4, string); + flag_d0 = dparse (d0, 7, string); + flag_d1 = dparse (d1, 9, string); + *h0 *= flag_h0; + *h1 *= flag_h1; + *d0 *= flag_d0; + *d1 *= flag_d1; + + dparse (&tmp, 2, string); + *h0 += tmp/60.0; + dparse (&tmp, 3, string); + *h0 += tmp/3600.0; + + dparse (&tmp, 5, string); + *h1 += tmp/60.0; + dparse (&tmp, 6, string); + *h1 += tmp/3600.0; + + dparse (&tmp, 8, string); + *d0 += tmp/60.0; + + dparse (&tmp, 10, string); + *d1 += tmp/60.0; + + *h0 *= 15*flag_h0; + *h1 *= 15*flag_h1; + *d0 *= flag_d0; + *d1 *= flag_d1; + + return (TRUE); +} + + +/* + if (buffer[i*48 + 19] == ' ') continue; + strncpy (temp, &buffer[i*48 + 19], 20); temp[20] = 0; + hms_to_deg (&RA, &DEC, temp, ' ', 3); + + + dBigRA = 360.0 / (int)(0.5 + 48*cos(3.1415927*0.5*(BigDecBounds[NBigDec] + BigDecBounds[NBigDec + 1])/180.0)); + + NBigRA = ra / dBigRA; + + NBig = NBigRA; + for (i = 0; (i < 12) && (i < NBigDec); i++) { + NBig += NBigRASections[i]; + } + for (i = 13; (i < 24) && (i < NBigDec); i++) { + NBig += NBigRASections[i]; + } + + + fprintf (stderr, "%d %d %d %d %f %f -> %f %f %f\n", NBigDec, NBigRA, NBigRASections[NBigDec], NBig, ra, dec, BigDecBounds[NBigDec], BigDecBounds[NBigDec + 1], dBigRA); +*/ Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/gimages.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/gimages.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/gimages.c (revision 21890) @@ -0,0 +1,214 @@ +# include "addstar.h" +double opening_angle (); + +Image *gimages (image, Npimage) +Image image[]; +int *Npimage; +{ + + int i, j, k, NIMAGE, Nimage, addtolist; + int NTIMAGE, Ntimage, ntimage; + int NPIMAGE, npimage; + FILE *f; + Image *timage, *pimage; + Coords tcoords; + Header header; + double r, d; + double Xi[5], Yi[5], Xo[5], Yo[5]; /* image and original corners */ + + /* check if image datafile exists, get header */ + if (!fits_read_header (ImageCat, &header)) { + if (VERBOSE) fprintf (stderr, "New image catalog %s (1)\n", ImageCat); + *Npimage = 0; + return ((Image *) NULL); + } + + /* project onto rectilinear grid with 1 arcsec pixels */ + /* we keep the original crpix1,2 and crref1,2 */ + tcoords = image[0].coords; + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + tcoords.pc1_1 = tcoords.pc2_2 = 1.0; + tcoords.pc1_2 = tcoords.pc2_1 = 0.0; + strcpy (tcoords.ctype, "RA---TAN"); + + /* define original corners */ + Xo[0] = 0; Yo[0] = 0; + Xo[1] = image[0].NX; Yo[1] = 0; + Xo[2] = image[0].NX; Yo[2] = image[0].NY; + Xo[3] = 0; Yo[3] = image[0].NY; + Xo[4] = 0; Yo[4] = 0; /* so we can make a loop easily */ + for (j = 0; j < 5; j++) { + XY_to_RD (&r, &d, Xo[j], Yo[j], &image[0].coords); + RD_to_XY (&Xo[j], &Yo[j], r, d, &tcoords); + } + + /* get ready to read data on images */ + f = fopen (ImageCat, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "New image catalog: %s (2)\n", ImageCat); + free (header.buffer); + *Npimage = 0; + return ((Image *) NULL); + } + fseek (f, header.size, SEEK_SET); + + Nimage = 0; + fits_scan (&header, "NIMAGES", "%d", 1, &Nimage); + + NTIMAGE = 100; + ALLOCATE (timage, Image, NTIMAGE); + + npimage = 0; + NPIMAGE = 20; + ALLOCATE (pimage, Image, NPIMAGE); + + for (Ntimage = 0; Ntimage < Nimage; Ntimage += ntimage) { + ntimage = Fread (timage, sizeof(Image), NTIMAGE, f, "image"); + for (i = 0; i < ntimage; i++) { + addtolist = FALSE; + /* define image corners */ + Xi[0] = 0; Yi[0] = 0; + Xi[1] = timage[i].NX; Yi[1] = 0; + Xi[2] = timage[i].NX; Yi[2] = timage[i].NY; + Xi[3] = 0; Yi[3] = timage[i].NY; + Xi[4] = 0; Yi[4] = 0; /* so we can make a loop easily */ + /* transform to tcoords */ + for (j = 0; j < 5; j++) { + XY_to_RD (&r, &d, Xi[j], Yi[j], &timage[i].coords); + RD_to_XY (&Xi[j], &Yi[j], r, d, &tcoords); + } + /* check if any edges cross */ + for (j = 0; j < 4; j++) { + for (k = 0; k < 4; k++) { + addtolist |= edge_check (&Xi[j], &Yi[j], &Xo[k], &Yo[k]); + } + } + + /* + fprintf (stderr, "%f %f %f %d %f %f %d %d %d %f %f %s: %d\n", 0.1*timage[i].detection_limit, 0.1*timage[i].saturation_limit, 0.02*timage[i].cerror, + timage[i].tzero, 1e-4*timage[i].trate, + 0.001*timage[i].secz, timage[i].nstar, timage[i].NX, timage[i].NY, 0.001*timage[i].Mcal, 0.1*timage[i].Xm, timage[i].name, addtolist); + + fprintf (stderr, "%s\n", timage[i].coords.ctype); + fprintf (stderr, "%f %f\n", timage[i].coords.crval1, timage[i].coords.crval2); + fprintf (stderr, "%f %f\n", timage[i].coords.crpix1, timage[i].coords.crpix2); + fprintf (stderr, "%f %f\n", timage[i].coords.pc1_1, timage[i].coords.pc1_2); + fprintf (stderr, "%f %f\n", timage[i].coords.pc2_1, timage[i].coords.pc2_2); + fprintf (stderr, "%f %f\n", timage[i].coords.cdelt1, timage[i].coords.cdelt2); + */ + + if (addtolist) { + pimage[npimage] = timage[i]; + npimage ++; + if (npimage == NPIMAGE) { + NPIMAGE += 20; + REALLOCATE (pimage, Image, NPIMAGE); + } + } + } + + if (ntimage == 0) { + fprintf (stderr, "ERROR: expected %d images, only found %d\n", Nimage, Ntimage); + clear_lockfile (); + exit (0); + } + } + + if (VERBOSE) fprintf (stderr, "found %d overlapping images\n", npimage); + + fclose (f); + *Npimage = npimage; + return (pimage); + +} + +int edge_check (x1, y1, x2, y2) +double *x1, *y1, *x2, *y2; +{ + + double theta1, theta2; + double Theta1, Theta2; + + theta1 = opening_angle (x1[0], y1[0], x2[0], y2[0], x1[1], y1[1]); + theta2 = opening_angle (x1[0], y1[0], x2[0], y2[0], x2[1], y2[1]); + + if (theta1*theta2 < 0.0) { + return (FALSE); + } + + if (fabs(theta1) < fabs(theta2)) { + return (FALSE); + } + + Theta1 = theta1; + Theta2 = theta2; + theta1 = opening_angle (x2[0], y2[0], x1[1], y1[1], x2[1], y2[1]); + theta2 = opening_angle (x2[0], y2[0], x1[1], y1[1], x1[0], y1[0]); + + + if (theta1*theta2 < 0.0) { + return (FALSE); + } + + if (fabs(theta1) < fabs(theta2)) { + return (FALSE); + } + + return (TRUE); + +} + +double opening_angle (x1, y1, x2, y2, x3, y3) +double x1, y1, x2, y2, x3, y3; +{ + + double dx1, dy1, dx2, dy2, ct, st, theta; + + dx1 = x1 - x2; + dy1 = y1 - y2; + + dx2 = x3 - x2; + dy2 = y3 - y2; + + ct = (dx1*dx2 + dy1*dy2); + st = (dx1*dy2 - dx2*dy1); + + theta = atan2 (st, ct); + + return (theta); + +} + + + + +/* probably not needed now +char *_parse_nextword(); + + +int fparse (X, NX, line) +float *X; +int NX; +char *line; +{ + + int i; + char *word; + char *ptr; + + word = line; + for (i = 0; i < NX - 1; i++) + word = _parse_nextword (word); + + *X = strtod (word, &ptr); + if (ptr == word) + return (FALSE); + else { + if (word[0] == '-') + return (-1); + else + return (1); + } +} + +*/ Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/gregions.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/gregions.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/gregions.c (revision 21890) @@ -0,0 +1,57 @@ +# include "addstar.h" + +GSCRegion *gregions (image, Nregions) +Image *image; +int *Nregions; +{ + + GSCRegion *region; + FILE *f; + double x, y; + double dr, dd, dec, ra; + int i, j, done, nregion, NREGION; + + f = fopen (GSCFILE, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't find GSC region file %s\n", GSCFILE); + clear_lockfile (); + exit (0); + } + + /* find regions at image corners */ + NREGION = 10; + ALLOCATE (region, GSCRegion, NREGION); + nregion = 0; + + /* look for new regions on grid across image */ + for (x = 0.0; x <= 1.0; x+=0.25) { + for (y = 0.0; y <= 1.0; y+=0.25) { + XY_to_RD (&ra, &dec, image[0].NX*(1.1*x - 0.05), image[0].NY*(1.1*y - 0.05), &image[0].coords); + aregion (®ion[nregion], f, ra, dec); + done = FALSE; + for (j = 0; (j < nregion) && !done; j++) { + if (!strcmp (region[nregion].filename, region[j].filename)) { + nregion --; + done = TRUE; + } + } + nregion ++; + if (nregion == NREGION) { + NREGION += 10; + REALLOCATE (region, GSCRegion, NREGION); + } + } + } + + if (VERBOSE) { + fprintf (stderr, "found %d region files:\n", nregion); + for (i = 0; i < nregion; i++) { + fprintf (stderr, " %d %s\n", i, region[i].filename); + } + } + *Nregions = nregion; + + fclose (f); + return (region); + +} Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/gstars.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/gstars.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/gstars.c (revision 21890) @@ -0,0 +1,189 @@ +# include "addstar.h" +# define D_NSTARS 1000 +# define BYTES_STAR 31 +# define BLOCK 1000 +#include +#include + +Stars *gstars (file, NSTARS, image) +char file[]; +int *NSTARS; +Image *image; +{ + + FILE *f; + Header header; + int i, j, nstar; + int nbytes, Nbytes; + Stars *stars; + char *buffer, line[64], *c; + struct tm timeptr; + double tmp; + + if (!fits_read_header (file, &header)) { + fprintf (stderr, "ERROR: can't find image file %s\n", file); + clear_lockfile (); + exit(0); + } + + /* open file */ + f = fopen (file, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't find data file %s\n", file); + clear_lockfile (); + exit(0); + } + fseek (f, header.size, SEEK_SET); + + c = strrchr (file, 0x2f); + if (c == (char *) NULL) { + strcpy (image[0].name, file); + } else { + strcpy (image[0].name, (c+1)); + } + /* get astrometry information */ + strcpy (image[0].coords.ctype, "NONE"); + fits_scan (&header, "CTYPE1", "%s", 1, image[0].coords.ctype); + if (strcmp (image[0].coords.ctype, "RA---PLY")) { + fprintf (stderr, "ERROR: wrong astrometric info in header\n"); + clear_lockfile (); + exit (0); + } + fits_scan (&header, "CDELT1", "%f", 1, &image[0].coords.cdelt1); + fits_scan (&header, "CDELT2", "%f", 1, &image[0].coords.cdelt2); + fits_scan (&header, "CRVAL1", "%lf", 1, &image[0].coords.crval1); + fits_scan (&header, "CRVAL2", "%lf", 1, &image[0].coords.crval2); + fits_scan (&header, "CRPIX1", "%f", 1, &image[0].coords.crpix1); + fits_scan (&header, "CRPIX2", "%f", 1, &image[0].coords.crpix2); + fits_scan (&header, "PC001001", "%f", 1, &image[0].coords.pc1_1); + fits_scan (&header, "PC001002", "%f", 1, &image[0].coords.pc1_2); + fits_scan (&header, "PC002001", "%f", 1, &image[0].coords.pc2_1); + fits_scan (&header, "PC002002", "%f", 1, &image[0].coords.pc2_2); + fits_scan (&header, "CERROR", "%lf", 1, &tmp); + /* RA Terms */ + fits_scan (&header, "PCA1X2Y0", "%f", 1, &image[0].coords.polyterms[0][0]); + fits_scan (&header, "PCA1X1Y1", "%f", 1, &image[0].coords.polyterms[1][0]); + fits_scan (&header, "PCA1X0Y2", "%f", 1, &image[0].coords.polyterms[2][0]); + fits_scan (&header, "PCA1X3Y0", "%f", 1, &image[0].coords.polyterms[3][0]); + fits_scan (&header, "PCA1X2Y1", "%f", 1, &image[0].coords.polyterms[4][0]); + fits_scan (&header, "PCA1X1Y2", "%f", 1, &image[0].coords.polyterms[5][0]); + fits_scan (&header, "PCA1X0Y3", "%f", 1, &image[0].coords.polyterms[6][0]); + /* Dec Terms */ + fits_scan (&header, "PCA2X2Y0", "%f", 1, &image[0].coords.polyterms[0][1]); + fits_scan (&header, "PCA2X1Y1", "%f", 1, &image[0].coords.polyterms[1][1]); + fits_scan (&header, "PCA2X0Y2", "%f", 1, &image[0].coords.polyterms[2][1]); + fits_scan (&header, "PCA2X3Y0", "%f", 1, &image[0].coords.polyterms[3][1]); + fits_scan (&header, "PCA2X2Y1", "%f", 1, &image[0].coords.polyterms[4][1]); + fits_scan (&header, "PCA2X1Y2", "%f", 1, &image[0].coords.polyterms[5][1]); + fits_scan (&header, "PCA2X0Y3", "%f", 1, &image[0].coords.polyterms[6][1]); + image[0].coords.Npolyterms = 2; /* how many do we use? */ + /* CERROR in data file is in pixels, convert to 20*arcsec */ + image[0].cerror = tmp * 50.0 * image[0].coords.cdelt1 * 3600.0; + while (image[0].coords.crval1 < 0) image[0].coords.crval1 += 360.0; + while (image[0].coords.crval1 > 360.0) image[0].coords.crval1 -= 360.0; + + fits_scan (&header, "NAXIS1", "%hd", 1, &image[0].NX); + fits_scan (&header, "NAXIS2", "%hd", 1, &image[0].NY); + fits_scan (&header, "PHOTCODE", "%hd", 1, &image[0].source); + if (NEWPHOTCODE) { + image[0].source = PHOTCODE; + } + image[0].NX -= XOVERSCAN; + image[0].NY -= YOVERSCAN; + + tmp = 0; + fits_scan (&header, "FLIMIT", "%lf", 1, &tmp); + image[0].detection_limit = tmp * 10.0; + + tmp = 0; + fits_scan (&header, "FSATUR", "%lf", 1, &tmp); + image[0].saturation_limit = tmp * 10.0; + + if (!fits_scan (&header, "TZERO", "%d", 1, &image[0].tzero)) { + fits_scan (&header, "UT", "%s", 1, line); + /* remove : characters */ + for (c = strchr (line, 0x3a); c != (char *) NULL; c = strchr (c, 0x3a)) + *c = ' '; + sscanf (line, "%d %d %d", &timeptr.tm_hour, &timeptr.tm_min, &timeptr.tm_sec); + fits_scan (&header, "DATE-UT", "%s", 1, line); + for (c = strchr (line, 0x2f); c != (char *) NULL; c = strchr (c, 0x2f)) + *c = ' '; + sscanf (line, "%d %d %d", &timeptr.tm_year, &timeptr.tm_mon, &timeptr.tm_mday); + if (timeptr.tm_year < 90) { /* correct for year 2000 turnover */ + timeptr.tm_year += 100; + } + image[0].tzero = mktime (&timeptr); + } + + tmp = 0; + fits_scan (&header, "TRATE", "%lf", 1, &tmp); + image[0].trate = 10000 * tmp; + + tmp = 0; + fits_scan (&header, "AIRMASS", "%lf", 1, &tmp); + image[0].secz = 1000*tmp; + + /* secz is in units milli-airmass */ + image[0].Mcal = ALPHA*(image[0].secz - 1000); + image[0].Xm = 30.0; + image[0].code = 0; + bzero (image[0].dummy, sizeof(image[0].dummy)); + + /* find number of stars */ + fits_scan (&header, "NSTARS", "%d", 1, &image[0].nstar); + if (image[0].nstar == 0) { + fprintf (stderr, "ERROR: can't get NSTARS from header\n"); + clear_lockfile (); + exit (0); + } + ALLOCATE (stars, Stars, image[0].nstar); + Nbytes = image[0].nstar*BYTES_STAR; + + /* load in stars by blocks of 1000 */ + nstar = 0; + ALLOCATE (buffer, char, (BLOCK*BYTES_STAR)); + for (i = 0; i < (int)(Nbytes / (BLOCK*BYTES_STAR)); i++) { + nbytes = Fread (buffer, 1, (BLOCK*BYTES_STAR), f, "char"); + if (nbytes != BLOCK*BYTES_STAR) { + fprintf (stderr, "ERROR: failed to read in stars (1)\n"); + clear_lockfile (); + exit (0); + } + for (j = 0; j < BLOCK; j++, nstar++) { + dparse (&stars[nstar].X, 1, &buffer[j*BYTES_STAR]); + dparse (&stars[nstar].Y, 2, &buffer[j*BYTES_STAR]); + dparse (&stars[nstar].M, 3, &buffer[j*BYTES_STAR]); + dparse (&stars[nstar].dM, 4, &buffer[j*BYTES_STAR]); + XY_to_RD (&stars[nstar].R, &stars[nstar].D, stars[nstar].X, stars[nstar].Y, &image[0].coords); + stars[nstar].found = -1; + } + } + /* left over fraction of a block */ + nbytes = Fread (buffer, 1, (Nbytes % (BLOCK*BYTES_STAR)), f, "char"); + if (nbytes != (Nbytes % (BLOCK*BYTES_STAR))) { + fprintf (stderr, "ERROR: failed to read in stars (2)\n"); + clear_lockfile (); + exit (0); + } + for (j = 0; j < nbytes / BYTES_STAR; j++, nstar++) { + dparse (&stars[nstar].X, 1, &buffer[j*BYTES_STAR]); + dparse (&stars[nstar].Y, 2, &buffer[j*BYTES_STAR]); + dparse (&stars[nstar].M, 3, &buffer[j*BYTES_STAR]); + dparse (&stars[nstar].dM, 4, &buffer[j*BYTES_STAR]); + XY_to_RD (&stars[nstar].R, &stars[nstar].D, stars[nstar].X, stars[nstar].Y, &image[0].coords); + stars[nstar].found = -1; + } + + if (image[0].nstar != nstar) { + fprintf (stderr, "ERROR: failed to read in all stars (%d != %d)\n", image[0].nstar, nstar); + clear_lockfile (); + exit (0); + } + free (header.buffer); + free (buffer); + + if (VERBOSE) fprintf (stderr, "read %d stars from target file\n", image[0].nstar); + *NSTARS = image[0].nstar; + + return (stars); +} Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/sort_stars.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/sort_stars.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/sort_stars.c (revision 21890) @@ -0,0 +1,61 @@ +# include "relphot.h" + +void sort_stars (radec, stars, Nstars) +int **radec; +Star *stars; +int Nstars; +{ + + int i; + double *RAs; + + ALLOCATE (radec[0], int, Nstars); + ALLOCATE (RAs, double, Nstars); + for (i = 0; i < Nstars; i++) { + radec[0][i] = i; + RAs[i] = stars[i].RA; + } + + sort_seq (radec[0], RAs, Nstars); + free (RAs); + +} + + +void sort_seq (seq, value, N) +int *seq; +double *value; +int N; +{ + int l,j,ir,i; + int temp; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + temp = seq[--l]; + } + else { + temp = seq[ir]; + seq[ir] = seq[0]; + if (--ir == 0) { + seq[0] = temp; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && value[seq[j]] < value[seq[j+1]]) ++j; + if (value[temp] < value[seq[j]]) { + seq[i]=seq[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + seq[i] = temp; + } +} + + Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/string.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/string.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/string.c (revision 21890) @@ -0,0 +1,64 @@ +# include +# include +# define TRUE 1 +# define FALSE 0 + +int scan_line (f, line) +FILE *f; +char line[]; +{ + + int i, status; + char c; + + status = EOF + 1; + + for (i = 0, c = 0; (c != '\n') && (status != EOF); i++) { + status = fscanf (f, "%c", &c); + line[i] = c; + } + line[i - 1] = 0; /* this could make things crash! */ + + if (i > 1) { + status = EOF + 1; + } + + return (status); + +} + +char *_parse_nextword(string) +char *string; +{ + if (string == (char *) NULL) return ((char *) NULL); + + for (; isspace (*string); string++); + for (; (*string != 0) && !isspace (*string); string++); + for (; isspace (*string); string++); + return (string); +} + +int dparse (X, NX, line) +double *X; +int NX; +char *line; +{ + + int i; + char *word; + char *ptr; + + word = line; + for (i = 0; i < NX - 1; i++) + word = _parse_nextword (word); + + *X = strtod (word, &ptr); + if (ptr == word) + return (FALSE); + else { + if (word[0] == '-') + return (-1); + else + return (1); + } +} Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/wimage.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/wimage.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/old/wimage.c (revision 21890) @@ -0,0 +1,67 @@ +# include "addstar.h" + +wimage (filename, image, Nstars) +char *filename; +Image image[]; +int Nstars; +{ + + int Nimages, status; + FILE *f; + Header header; + + if (!fits_read_header (ImageCat, &header)) { + if (VERBOSE) fprintf (stderr, "can't find %s, creating a new one\n", ImageCat); + if (!fits_read_header (ImageTemplate, &header)) { + fprintf (stderr, "ERROR: can't find template header %s\n", ImageTemplate); + clear_lockfile (); + exit (0); + } + f = fopen (ImageCat, "w"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't create/open image catalog file: %s\n", ImageCat); + clear_lockfile (); + exit (0); + } + } + + Nimages = 0; + fits_scan (&header, "NIMAGES", "%d", 1, &Nimages); + Nimages ++; + fits_modify (&header, "NIMAGES", "%d", 1, Nimages); + + f = fopen (ImageCat, "r+"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't create/open image catalog file: %s\n", ImageCat); + clear_lockfile (); + exit (0); + } + + /* position to begining of file to write header */ + fseek (f, 0, SEEK_SET); + status = Fwrite (header.buffer, 1, header.size, f, "char"); + if (status != header.size) { + fprintf (stderr, "ERRPR: failed writing data to image header\n"); + clear_lockfile (); + exit (0); + } + + fseek (f, 0, SEEK_END); + fprintf (stderr, "Image: %d\n", sizeof(Image)); + status = Fwrite (image, sizeof(Image), 1, f, "image"); + if (status != 1) { + fprintf (stderr, "ERRPR: failed writing data to image catalog\n"); + clear_lockfile (); + exit (0); + } + + fclose (f); + + +} + +/* the image we add in this routine is always a new image. + We only need to do two things: + 1) append the data for the image to the end of the file + 2) update the NIMAGES field in the header (which means, read in header, re-write) + */ Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/reset_catalog.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/reset_catalog.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/reset_catalog.c (revision 21890) @@ -0,0 +1,86 @@ +# include "addusno.h" + +reset_catalog (catalog) +Catalog catalog[]; +{ + + int i, j, k, n, m, N, M, found; + int *next, last, this; + int Nave, Nmeas, Ndel; + unsigned int flags; + Measure *tmpmeasure; + unsigned short NotUSNO, NotProper; + + Nave = catalog[0].Naverage; + + /* set up link listed pointers for new measurements */ + Nmeas = catalog[0].Nmeasure; + ALLOCATE (next, int, Nmeas); + /* set up pointers for linked list of measurements */ + for (i = 0; i < Nmeas - 1; i++) { + next[i] = i+1; + } + next[i] = -1; + last = i; + + /** find matched stars **/ + NotUSNO = 0xffff ^ ID_USNO; + NotProper = 0xffff ^ ID_PROPER; + for (i = 0; i < Nave; i++) { + found = FALSE; + if (catalog[0].average[i].code & ID_USNO) { + catalog[0].average[i].code &= NotUSNO; + found = TRUE; + } + if (catalog[0].average[i].code & ID_PROPER) { + catalog[0].average[i].code &= NotProper; + found = TRUE; + } + if (found) { + M = catalog[0].average[i].offset; + Ndel = 0; + for (k = 0; k < catalog[0].average[i].Nm; k++) { + if (catalog[0].measure[M+k].t == 0) { + if (M+k > 0) { + this = next[M+k]; + if (next[M+k] == -1) { + next[M+k-1] = -1; + } else { + for (j = M+k+1; (j < Nmeas) && (next[j] == -2); j++); + if (j == Nmeas) { + next[M+k-1] = -1; + fprintf (stderr, "warning, why didn't I find last link to start with?\n"); + } + else { + next[M+k-1] = this; + } + } + } + next[M+k] = -2; /* delete link */ + Ndel ++; + } + } + catalog[0].average[i].Nm -= Ndel; + Nmeas -= Ndel; + } + } + + /* fix order of Measure (memory intensive, but fast) */ + N = 0; + ALLOCATE (tmpmeasure, Measure, Nmeas); + for (i = 0; i < Nave; i++) { + n = catalog[0].average[i].offset; + catalog[0].average[i].offset = N; + for (k = 0; k < catalog[0].average[i].Nm; k++, N++) { + tmpmeasure[N] = catalog[0].measure[n]; + n = next[n]; + } + } + free (catalog[0].measure); + catalog[0].measure = tmpmeasure; + + if (VERBOSE) fprintf (stderr, "deleted %d reference %d -> %d)\n", + catalog[0].Nmeasure - Nmeas, catalog[0].Nmeasure, Nmeas); + catalog[0].Nmeasure = Nmeas; +} + Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/sort_lists.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/sort_lists.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/sort_lists.c (revision 21890) @@ -0,0 +1,48 @@ + +sort_lists (X, Y, S, N) +float *X, *Y; +int *S, N; +{ + int l,j,ir,i; + double tX, tY, tS; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tX = X[l]; + tY = Y[l]; + tS = S[l]; + } + else { + tX = X[ir]; + X[ir] = X[0]; + tY = Y[ir]; + Y[ir] = Y[0]; + tS = S[ir]; + S[ir] = S[0]; + if (--ir == 0) { + X[0] = tX; + Y[0] = tY; + S[0] = tS; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && X[j] < X[j+1]) j++; + if (tX < X[j]) { + X[i] = X[j]; + Y[i] = Y[j]; + S[i] = S[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + X[i] = tX; + Y[i] = tY; + S[i] = tS; + } +} Index: /tags/ohana/addusno-1-0/Ohana/src/addusno/src/wcatalog.c =================================================================== --- /tags/ohana/addusno-1-0/Ohana/src/addusno/src/wcatalog.c (revision 21890) +++ /tags/ohana/addusno-1-0/Ohana/src/addusno/src/wcatalog.c (revision 21890) @@ -0,0 +1,74 @@ +# include "addstar.h" + +wcatalog (catname, catalog) +char *catname; +Catalog catalog[]; +{ + + char filename[256], line[256]; + int i, Nitems, nitems, status, mode; + FILE *f; + struct stat filestat; + + sprintf (filename, "%s/%s\0", CATDIR, catname); + + status = stat (filename, &filestat); + if (status == 0) { /* file exists, make backup copy */ + sprintf (line, "mv %s %s~\0", filename, filename); + status = system (line); + if (status) { + fprintf (stderr, "ERROR: unable to create %s~, exiting\n", filename); + exit (0); + } + } + + if (catalog[0].Naverage == 0) { + if (VERBOSE) fprintf (stderr, "no stars in catalog, skipping\n"); + return (FALSE); + } + + f = fopen (filename, "w"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't create new catalog file: %s\n", filename); + exit (0); + } + mode = S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH; + chmod (filename, mode); + + fits_modify (&catalog[0].header, "NSTARS", "%d", 1, catalog[0].Naverage); + fits_modify (&catalog[0].header, "NMEAS", "%d", 1, catalog[0].Nmeasure); + fits_modify (&catalog[0].header, "NMISS", "%d", 1, catalog[0].Nmissing); + + fits_modify (&catalog[0].header, "ADDUSNO", "%t", 1, TRUE); + + nitems = Fwrite (catalog[0].header.buffer, 1, catalog[0].header.size, f, "char"); + if (nitems != catalog[0].header.size) { + fprintf (stderr, "ERROR: failed to write header\n"); + exit (0); + } + + Nitems = catalog[0].Naverage; + nitems = Fwrite (catalog[0].average, sizeof(Average), Nitems, f, "average"); + if (nitems != Nitems) { + fprintf (stderr, "ERROR: failed to write catalog file aves %s\n", filename); + exit (0); + } + + Nitems = catalog[0].Nmeasure; + nitems = Fwrite (catalog[0].measure, sizeof(Measure), Nitems, f, "measure"); + if (nitems != Nitems) { + fprintf (stderr, "ERROR: failed to write catalog file meas %s\n", filename); + exit (0); + } + + Nitems = catalog[0].Nmissing; + nitems = Fwrite (catalog[0].missing, sizeof(Missing), Nitems, f, "missing"); + if (nitems != Nitems) { + fprintf (stderr, "ERROR: failed to write catalog file miss %s\n", filename); + exit (0); + } + + fclose (f); + +} +