Index: /tags/fixcat-1-0/Ohana/src/fixcat/Makefile =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/Makefile (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/Makefile (revision 2512) @@ -0,0 +1,63 @@ +include ../../Configure +CC = gcc -g # override this for now for minimal test +HOME = $(ROOT)/src/fixcat +CONFIG = $(ROOT)/config +PROGRAM = fixcat + +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) + +FIXCAT = \ +$(SRC)/fixcat.$(ARCH).o $(SRC)/gcatalog.$(ARCH).o \ +$(SRC)/coordops.$(ARCH).o $(SRC)/sorts.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o $(SRC)/config.$(ARCH).o \ +$(SRC)/Fread.$(ARCH).o $(SRC)/check_lockfile.$(ARCH).o \ +$(SRC)/find_images.$(ARCH).o $(SRC)/match_images.$(ARCH).o \ +$(SRC)/gregions.$(ARCH).o $(SRC)/aregion.$(ARCH).o \ +$(SRC)/find_funnymags.$(ARCH).o $(SRC)/gcatstats.$(ARCH).o\ +$(SRC)/find_matches.$(ARCH).o $(SRC)/check_permissions.$(ARCH).o + + +OBJ = $(FIXCAT) + +default: $(PROGRAM) + +$(MARKSTAR): $(INC)/markstar.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 $(BIN) ]; then mkdir -p $(BIN); 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/fixcat-1-0/Ohana/src/fixcat/doc/example.txt =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/doc/example.txt (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/doc/example.txt (revision 2512) @@ -0,0 +1,13 @@ + +rats: foreach f ( */*.cpt ) +foreach? echo $f +foreach? echo $f >> bright.log +foreach? (markstar -v $f >> bright.dat) >>& bright.log +foreach? end + + + + echo $f + echo $f >> bright.log + (markstar -v $f >> bright.dat) >>& bright.log + end Index: /tags/fixcat-1-0/Ohana/src/fixcat/include/markstar.h =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/include/markstar.h (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/include/markstar.h (revision 2512) @@ -0,0 +1,42 @@ +# include +# include + +int VERBOSE; +int RESET; +int FORCE_RUN; + +/* global variables set in parameter file */ +char CATDIR[256]; +char GSCDIR[256]; +char ImageCat[256]; +char GSCFILE[256]; + +double RADIUS; +double TRAIL_WIDTH; +int NBINS; +int NPTSINLINE; +double MIN_DENSITY; +double NSIGMA; + +double BRIGHT_HALO_MAG; +double BRIGHT_HALO_SLOPE; +double BRIGHT_XTRAIL_WIDTH; +double BRIGHT_XTRAIL_MAG; +double BRIGHT_XTRAIL_SLOPE; +double BRIGHT_YTRAIL_WIDTH; +double BRIGHT_YTRAIL_MAG; +double BRIGHT_YTRAIL_SLOPE; + +double GHOST_MAG; +double GHOST_RADIUS; +double OPTICAL_AXIS1; +double OPTICAL_AXIS2; + +typedef struct { + Coords coords; + float *X, *Y; + int *N; + double RA[2], DEC[2]; + double Area, density, spacing; +} CatStats; + Index: /tags/fixcat-1-0/Ohana/src/fixcat/src/ConfigInit.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/ConfigInit.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/ConfigInit.c (revision 2512) @@ -0,0 +1,45 @@ +# include "markstar.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); + + /* used in other pipeline functions */ + ScanConfig (config, "IMAGE_CATALOG", "%s", 0, ImageCat); + ScanConfig (config, "CATDIR", "%s", 0, CATDIR); + ScanConfig (config, "GSCDIR", "%s", 0, GSCDIR); + ScanConfig (config, "GSCFILE", "%s", 0, GSCFILE); + /* unique to markstar */ + ScanConfig (config, "SEARCH_RADIUS", "%lf", 0, &RADIUS); + ScanConfig (config, "TRAIL_WIDTH", "%lf", 0, &TRAIL_WIDTH); + ScanConfig (config, "NANGLE_BINS", "%d", 0, &NBINS); + ScanConfig (config, "NPTSINLINE", "%d", 0, &NPTSINLINE); + ScanConfig (config, "MIN_DENSITY", "%lf", 0, &MIN_DENSITY); + ScanConfig (config, "SPACE_SIGMA", "%lf", 0, &NSIGMA); + + ScanConfig (config, "BRIGHT_HALO_MAG", "%lf", 0, &BRIGHT_HALO_MAG); + ScanConfig (config, "BRIGHT_HALO_SLOPE", "%lf", 0, &BRIGHT_HALO_SLOPE); + ScanConfig (config, "BRIGHT_XTRAIL_WIDTH", "%lf", 0, &BRIGHT_XTRAIL_WIDTH); + ScanConfig (config, "BRIGHT_XTRAIL_MAG", "%lf", 0, &BRIGHT_XTRAIL_MAG); + ScanConfig (config, "BRIGHT_XTRAIL_SLOPE", "%lf", 0, &BRIGHT_XTRAIL_SLOPE); + ScanConfig (config, "BRIGHT_YTRAIL_WIDTH", "%lf", 0, &BRIGHT_YTRAIL_WIDTH); + ScanConfig (config, "BRIGHT_YTRAIL_MAG", "%lf", 0, &BRIGHT_YTRAIL_MAG); + ScanConfig (config, "BRIGHT_YTRAIL_SLOPE", "%lf", 0, &BRIGHT_YTRAIL_SLOPE); + + ScanConfig (config, "GHOST_MAG", "%lf", 0, &GHOST_MAG); + ScanConfig (config, "GHOST_RADIUS", "%lf", 0, &GHOST_RADIUS); + ScanConfig (config, "OPTICAL_AXIS1", "%lf", 0, &OPTICAL_AXIS1); + ScanConfig (config, "OPTICAL_AXIS2", "%lf", 0, &OPTICAL_AXIS2); + + free (config); +} Index: /tags/fixcat-1-0/Ohana/src/fixcat/src/aregion.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/aregion.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/aregion.c (revision 2512) @@ -0,0 +1,186 @@ +# include "markstar.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; } + while (ra >= 360.0) { ra -= 360.0; } + + if (dec >= 86.25) { + sprintf (file, "%s/n8230/pole.cpt\0", GSCDIR); + 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); + exit (0); + } + temp[5] = 0; + sprintf (file, "%s/%s/%s.cpt\0", GSCDIR, 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/fixcat-1-0/Ohana/src/fixcat/src/check_lockfile.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/check_lockfile.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/check_lockfile.c (revision 2512) @@ -0,0 +1,39 @@ +# include "markstar.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/fixcat-1-0/Ohana/src/fixcat/src/check_permissions.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/check_permissions.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/check_permissions.c (revision 2512) @@ -0,0 +1,63 @@ +# include "markstar.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/fixcat-1-0/Ohana/src/fixcat/src/config.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/config.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/config.c (revision 2512) @@ -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/fixcat-1-0/Ohana/src/fixcat/src/coordops.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/coordops.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/coordops.c (revision 2512) @@ -0,0 +1,342 @@ +# include "markstar.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/fixcat-1-0/Ohana/src/fixcat/src/find_bright_stars.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/find_bright_stars.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/find_bright_stars.c (revision 2512) @@ -0,0 +1,142 @@ +# include "markstar.h" + +find_bright_stars (catalog, catstats) +Catalog catalog[]; +CatStats catstats[]; +{ + + int i, j, n, m, first_j, Nave, Ngsc; + Catalog GSCdata; + double MinRA, MinDEC, MaxRA, MaxDEC, RaCenter, DecCenter; + Coords *tcoords; + float *X1, *Y1, *X2, *Y2; + int *N1, *N2; + char *mark; + double dX, dY, dR, MaxDist, MaxDist1, MaxRadius, radius, radius2; + + load_gsc_data (&GSCdata, catstats); + + Nave = catalog[0].Naverage; + Ngsc = GSCdata.Naverage; + + /* in the function below, it is better to have the catalog with + more stars associated with index 2 (j) */ + X2 = catstats[0].X; + Y2 = catstats[0].Y; + N2 = catstats[0].N; + ALLOCATE (X1, float, Ngsc); + ALLOCATE (Y1, float, Ngsc); + ALLOCATE (N1, int, Ngsc); + ALLOCATE (mark, char, Nave); + bzero (mark, Nave); + + tcoords = &catstats[0].coords; + for (i = 0; i < Ngsc; i++) { + fRD_to_XY (&X1[i], &Y1[i], GSCdata.average[i].R, GSCdata.average[i].D, tcoords); + N1[i] = i; + } + if (Ngsc > 1) sort_lists (X1, Y1, N1, Ngsc); + + /* first find stellar halos */ + /* max radius (mag = -1) */ + /** the j index moves quickly and is better associated with the catalog with more stars */ + MaxRadius = BRIGHT_HALO_SLOPE * (-1.0 - BRIGHT_HALO_MAG); + /** find catalog stars near GSC stars **/ + for (i = j = 0; (i < Ngsc) && (j < Nave); ) { + dX = X1[i] - X2[j]; + if (dX <= -2*MaxRadius) { + i++; + continue; + } + if (dX >= 2*MaxRadius) { + j++; + continue; + } + /* negative dX: j is too large, positive dX, i is too large */ + first_j = j; + radius = MAX (BRIGHT_HALO_SLOPE * (0.001*GSCdata.average[N1[i]].M - BRIGHT_HALO_MAG), 0.0); + if (radius == 0) { + i++; + continue; + } + radius2 = radius*radius; + for (; (dX > -2*radius) && (j < Nave); j++) { + dX = X1[i] - X2[j]; + dY = Y1[i] - Y2[j]; + dR = dX*dX + dY*dY; + if (dR < radius2) { /* new measurement of this star */ + mark[j] = TRUE; + } + } + j = first_j; + i++; + } + + /* next find y spikes */ + /** find catalog stars near GSC stars **/ + for (i = j = 0; (i < Ngsc) && (j < Nave); ) { + dX = X1[i] - X2[j]; + if (dX <= -BRIGHT_YTRAIL_WIDTH) { + i++; + continue; + } + if (dX >= BRIGHT_YTRAIL_WIDTH) { + j++; + continue; + } + /* negative dX: j is too large, positive dX, i is too large */ + first_j = j; + MaxDist = MAX (BRIGHT_YTRAIL_SLOPE * (0.001*GSCdata.average[N1[i]].M - BRIGHT_YTRAIL_MAG), 0.0); + for (; (dX > -BRIGHT_YTRAIL_WIDTH) && (j < Nave); j++) { + dX = X1[i] - X2[j]; + dY = Y1[i] - Y2[j]; + if ((fabs(dX) < BRIGHT_YTRAIL_WIDTH) && (fabs(dY) < MaxDist)) { /* star on spike */ + mark[j] = TRUE; + } + } + j = first_j; + i++; + } + + /* next find x spikes */ + /** find catalog stars near GSC stars **/ + MaxDist = MAX (BRIGHT_XTRAIL_SLOPE * (-1.0 - BRIGHT_XTRAIL_MAG), 0.0); + for (i = j = 0; (i < Ngsc) && (j < Nave); ) { + dX = X1[i] - X2[j]; + if (dX <= -MaxDist) { + i++; + continue; + } + if (dX >= MaxDist) { + j++; + continue; + } + /* negative dX: j is too large, positive dX, i is too large */ + first_j = j; + MaxDist1 = MAX (BRIGHT_XTRAIL_SLOPE * (0.001*GSCdata.average[N1[i]].M - BRIGHT_XTRAIL_MAG), 0.0); + for (; (dX > -MaxDist1) && (j < Nave); j++) { + dX = X1[i] - X2[j]; + dY = Y1[i] - Y2[j]; + if ((fabs(dY) < BRIGHT_XTRAIL_WIDTH) && (fabs(dX) < MaxDist1)) { /* star on spike */ + mark[j] = TRUE; + } + } + j = first_j; + i++; + } + + /* done with search, mark selected stars */ + for (i = 0; i < Nave; i++) { + if (mark[i]) { + catalog[0].average[N2[i]].code = ID_BLEED; + } + } + + free (X1); + free (Y1); + free (N1); + free (mark); + free (GSCdata.average); + +} + Index: /tags/fixcat-1-0/Ohana/src/fixcat/src/find_funnymags.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/find_funnymags.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/find_funnymags.c (revision 2512) @@ -0,0 +1,27 @@ +# include "markstar.h" + +find_funnymags (catalog, images, Nimage) +Catalog catalog[]; +Image images[]; +int Nimage; +{ + + int i; + int *bad; + + ALLOCATE (bad, int, Nimage); + bzero (bad, Nimage*sizeof(int)); + + for (i = 0; i < catalog[0].Nmeasure; i++) { + if (catalog[0].measure[i].M < 9000) { + bad[catalog[0].image[i]] ++; + } + } + + for (i = 0; i < Nimage; i++) { + if (bad[i] > 0) { + fprintf (stdout, "%s %d %d\n", images[i].name, images[i].tzero, bad[i]); + } + } + +} Index: /tags/fixcat-1-0/Ohana/src/fixcat/src/find_ghosts.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/find_ghosts.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/find_ghosts.c (revision 2512) @@ -0,0 +1,70 @@ +# include "markstar.h" + +find_ghosts (catalog, catstats, filename, image, Nimage) +Catalog catalog[]; +CatStats catstats[]; +char *filename; +Image image[]; +int Nimage; +{ + + int i, j, first_j, Nave, Ngsc, Nregions; + Catalog GSCdata, Ghostdata; + double MinRA, MinDEC, MaxRA, MaxDEC, RaCenter, DecCenter; + Coords *tcoords; + float *X1, *Y1, *X2, *Y2; + int *N1, *N2, *match; + char *mark; + double dX, dY, dR, RADIUS2, BRIGHT_RADIUS, radius, X, Y, R, D; + Image timage; + GSCRegion *region, *gregions2(); + + for (i = 0; i < Nimage; i++) { + + /* coords structure for ghost image */ + timage = image[i]; + timage.coords.cdelt1 = -1*image[i].coords.cdelt1; + timage.coords.cdelt2 = -1*image[i].coords.cdelt2; + timage.coords.crpix1 = 2*OPTICAL_AXIS1 - image[i].coords.crpix1; + timage.coords.crpix2 = 2*OPTICAL_AXIS2 - image[i].coords.crpix2; + + region = gregions2 (&timage, &Nregions); + + /* find ghost stars */ + load_gsc_data_ghost (&GSCdata, region, Nregions, &timage); + + /* refect ghost stars to find locations of ghosts */ + ALLOCATE (Ghostdata.average, Average, MAX (GSCdata.Naverage, 1)); + for (j = 0; j < GSCdata.Naverage; j++) { + RD_to_XY (&X, &Y, GSCdata.average[j].R, GSCdata.average[j].D, &timage); + fXY_to_RD (&Ghostdata.average[j].R, &Ghostdata.average[j].D, X, Y, &image[i]); + Ghostdata.average[j].M = GSCdata.average[j].M; + } + Ghostdata.Naverage = GSCdata.Naverage; + + /* match ghosts and image stars, mark ghosts */ + find_matches (catalog, catstats, &Ghostdata, i); + + free (Ghostdata.average); + free (GSCdata.average); + } + +} + +/* + + just to be clear on some of the terms: + + the "ghost" is the fuzzy patch of light on an image caused by + a reflected star + + the "ghost star" is the star in the sky which causes a ghost + + the "ghost image" is the image location on the sky + where ghost stars may come from. + + the "image stars" are observed stars at the location of the + ghost - these are detections of the ghost + + +*/ Index: /tags/fixcat-1-0/Ohana/src/fixcat/src/find_group.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/find_group.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/find_group.c (revision 2512) @@ -0,0 +1,69 @@ +# include "markstar.h" + +find_group (catstats, mark, Npts, i, ANGLE) + CatStats catstats[]; + double *ANGLE; + int i, Npts; + char *mark; +{ + + float *R, *D; + int j, N, bin, Ndegbin; + unsigned short int *A; + double Ra, De, dR, dD, angle, limit; + + /* assign some parameter values */ + + R = catstats[0].X; + D = catstats[0].Y; + limit = 10.0*sqrt((double)(M_PI*RADIUS*RADIUS*catstats[0].density)/(double)(NBINS)); + Ndegbin = NBINS / 180.0; + ALLOCATE (A, unsigned short int, NBINS) + bzero (A, NBINS*sizeof(short)); + + /* look for points concentrated in an angle bin */ + if (mark[i]) return (FALSE); + Ra = R[i]; + De = D[i]; + N = 0; + /* points east */ + for (j = i + 1; (j < Npts) && (R[j] - Ra < RADIUS); j++) { + if (mark[j]) continue; + dD = D[j] - De; + if (fabs(dD) < RADIUS) { + dR = R[j] - Ra; + angle = atan2 (dD,dR); + if (!finite(angle)) + continue; + if (angle < 0) angle += M_PI; + bin = 1 + angle*DEG_RAD*Ndegbin; + A[bin] ++; + if (A[bin] > limit) { + *ANGLE = (bin - 1.0) / (DEG_RAD*Ndegbin); + return (TRUE); + } + N ++; + } + } + /* points west */ + for (j = i - 1; (j >= 0) && (Ra - R[j] < RADIUS); j--) { + if (mark[j]) continue; + dD = D[j] - De; + if (fabs(dD) < RADIUS) { + dR = R[j] - Ra; + angle = atan2 (dD,dR); + if (!finite(angle)) + continue; + if (angle < 0) angle += M_PI; + bin = 1 + angle*DEG_RAD*Ndegbin; + A[bin] ++; + if (A[bin] > limit) { + *ANGLE = (bin - 1.0) / (DEG_RAD*Ndegbin); + return (TRUE); + } + N ++; + } + } + return (FALSE); +} + Index: /tags/fixcat-1-0/Ohana/src/fixcat/src/find_images.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/find_images.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/find_images.c (revision 2512) @@ -0,0 +1,261 @@ +# include "markstar.h" +# define OLDSTYLE 1 +# if (OLDSTYLE) +double opening_angle (); +# endif + +Image *find_images (catstats, Nimages) +CatStats catstats[]; +int *Nimages; +{ + + Header header; + Image *timage, *image; + int i, j, k, found, nimage, Nimage, NIMAGE, NTIMAGE, Nloop, Nlast; + int n, Nim, status; + FILE *f; + double Xc[6], Yc[6], Xi[6], Yi[6], r, d, x, y, dx, dy; + Coords *tcoords; + + /* we make positional comparisons in the projection of catalog */ + tcoords = &catstats[0].coords; + /* define catalog corners */ + Xc[0] = catstats[0].RA[0]; Yc[0] = catstats[0].DEC[0]; + Xc[1] = catstats[0].RA[1]; Yc[1] = catstats[0].DEC[0]; + Xc[2] = catstats[0].RA[1]; Yc[2] = catstats[0].DEC[1]; + Xc[3] = catstats[0].RA[0]; Yc[3] = catstats[0].DEC[1]; + Xc[4] = catstats[0].RA[0]; Yc[4] = catstats[0].DEC[0]; + Xc[5] = catstats[0].RA[1]; Yc[5] = catstats[0].DEC[1]; + for (j = 0; j < 6; j++) { + r = Xc[j]; d = Yc[j]; + RD_to_XY (&Xc[j], &Yc[j], r, d, tcoords); + } + /* find Y positions of RA center */ + r = 0.5*(catstats[0].RA[0] + catstats[0].RA[1]); + d = catstats[0].DEC[0]; + RD_to_XY (&x, &y, r, d, tcoords); + Yc[0] = MIN (y, Yc[0]); + Yc[1] = MIN (y, Yc[1]); + Yc[4] = MIN (y, Yc[4]); + /* find Y positions of RA center */ + r = 0.5*(catstats[0].RA[0] + catstats[0].RA[1]); + d = catstats[0].DEC[1]; + RD_to_XY (&x, &y, r, d, tcoords); + Yc[2] = MAX (y, Yc[2]); + Yc[3] = MAX (y, Yc[3]); + Yc[5] = MAX (y, Yc[5]); + + dx = 0.02*(Xc[2] - Xc[0]); + dy = 0.02*(Yc[2] - Yc[0]); + Xc[0] -= dx; Yc[0] -= dy; + Xc[1] += dx; Yc[1] -= dy; + Xc[2] += dx; Yc[2] += dy; + Xc[3] -= dx; Yc[3] += dy; + Xc[4] -= dx; Yc[4] -= dy; + Xc[5] -= dx; Yc[5] -= dy; + + /* check if image datafile exists, get header, number of images */ + if (!fits_read_header (ImageCat, &header)) { + fprintf (stderr, "ERROR: No images in catalog %s (1)\n", ImageCat); + exit (0); + } + Nimage = 0; + fits_scan (&header, "NIMAGES", "%d", 1, &Nimage); + if (Nimage == 0) { + fprintf (stderr, "ERROR: No images in catalog %s (1)\n", ImageCat); + exit (0); + } + + /* get ready to read data on images */ + f = fopen (ImageCat, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: trouble opening Image catalog: %s (2)\n", ImageCat); + exit (0); + } + fseek (f, header.size, SEEK_SET); + + /* set up buffers for images, temporary storage */ + NTIMAGE = 100; + ALLOCATE (timage, Image, NTIMAGE); + NIMAGE = 100; + ALLOCATE (image, Image, NIMAGE); + nimage = 0; + Nloop = Nimage / NTIMAGE + 1; + Nlast = Nimage % NTIMAGE; + + /* read in images in groups of NTIMAGE (100) */ + for (n = 0; n < Nloop; n++) { + Nim = (n == Nloop - 1) ? Nlast : NTIMAGE; + status = Fread (timage, sizeof(Image), Nim, f, "image"); + if (status != Nim) { + fprintf (stderr, "ERROR: couldn't read images from image catalog: %s\n", ImageCat); + exit (0); + } + /* test each image in block */ + for (i = 0; i < Nim; i++) { + /* 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; + Xi[5] = timage[i].NX; Yi[5] = timage[i].NY; + found = FALSE; + /* transform to tcoords */ + if (catstats[0].DEC[1] > 86.25) { /* pole */ + for (j = 0; j < 6; j++) { + XY_to_RD (&r, &d, Xi[j], Yi[j], &timage[i].coords); + if (d > catstats[0].DEC[0] - 0.5) found = TRUE; + } + } else { + for (j = 0; j < 6; 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 edges cross */ + for (j = 0; (j < 5) && !found; j++) { + for (k = 0; (k < 5) && !found; k++) { + found |= edge_check (&Xi[j], &Yi[j], &Xc[k], &Yc[k]); + } + } + } + if (found) { + image[nimage] = timage[i]; + image[nimage].code = 0; + nimage ++; + if (nimage == NIMAGE) { + NIMAGE += 100; + REALLOCATE (image, Image, NIMAGE); + } + } + } + } + + if (VERBOSE) { + for (i = 0; i < nimage; i++) { + XY_to_RD (&r, &d, 0.5*image[i].NX, 0.5*image[i].NY, &image[i].coords); + fprintf (stderr, "associated images: %d %8.4f %8.4f %10d %6d %5.3f %6.3f %6.3f\n", + i, r, d, image[i].tzero, image[i].nstar, 0.001*image[i].secz, + 0.001*image[i].Mcal, 0.001*image[i].dMcal); + } + } + + REALLOCATE (image, Image, MAX (nimage, 1)); + free (timage); + *Nimages = nimage; + fclose (f); + return (image); +} + +# if (OLDSTYLE) +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); + +} +#else + +int edge_check (x1, y1, x2, y2) +double *x1, *y1, *x2, *y2; +{ + + double dot; + double dx1, dx2, dy1, dy2, dx3, dy3, x, y; + double cross, L1, L3, d1, d2; + + dx1 = x1[1] - x1[0]; dy1 = y1[1] - y1[0]; + dx2 = x2[1] - x2[0]; dy2 = y2[1] - y2[0]; + + cross = dx2*dy1 - dx1*dy2; + + if (cross == 0) { /* lines are parallel, are they inline? */ + dx3 = x2[1] - x1[0]; dy2 = y2[1] - y1[0]; + L1 = hypot (dx1,dy1); + L3 = hypot (dx3,dy3); + + dot = fabs (dx1*dx3 + dy1*dy3) / (L1*L3); + if (dot == 1.0) { /* lines are inline, do they overlap? */ + d1 = (x1[1]-x2[1])*(x1[1]-x2[0]) + (y1[1]-y2[1])*(y1[1]-y2[0]); + d2 = (x1[0]-x2[1])*(x1[0]-x2[0]) + (y1[0]-y2[1])*(y1[0]-y2[0]); + if (d1*d2 < 0) { /* lines overlap */ + return (TRUE); + } else { + return (FALSE); + } + } else { + return (FALSE); + } + } + + x = (dx1*dx2*(y2[0] - y1[0]) + (x1[0]*dy1*dx2 - x2[0]*dy2*dx1)) / cross; + + if (dx1 != 0) { + y = y1[0] + (x - x1[0])*dy1/dx1; + } else { + y = y2[0] + (x - x2[0])*dy2/dx2; + } + + d1 = (x - x1[0])*(x - x1[1]) + (y - y1[0])*(y - y1[1]); + d2 = (x - x2[0])*(x - x2[1]) + (y - y2[0])*(y - y2[1]); + if ((d1 > 0) || (d2 > 0)) { + return (FALSE); + } + + return (TRUE); + +} + + +#endif + Index: /tags/fixcat-1-0/Ohana/src/fixcat/src/find_line.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/find_line.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/find_line.c (revision 2512) @@ -0,0 +1,161 @@ +# include "markstar.h" + +find_line (catstats, mark, Npts, i, M, B, Angle) + CatStats catstats[]; + double *M, *B, Angle; + int i, Npts; + char *mark; +{ + + float *R, *D; + int j, N; + double X, Y, X2, Y2, XY, m, b, det; + double dR, dD, Ra, De, angle; + char Flipped; + + R = catstats[0].X; + D = catstats[0].Y; + + /* fit a line to points near line */ + Ra = R[i]; + De = D[i]; + X = Ra; + Y = De; + X2 = Ra*Ra; + Y2 = De*De; + XY = Ra*De; + N = 1; + /* points to the east */ + for (j = i + 1; (j < Npts) && (R[j] - Ra < RADIUS); j++) { + dD = D[j] - De; + if (fabs(dD) < RADIUS) { + dR = R[j] - Ra; + angle = atan2 (dD,dR); + if (!finite(angle)) continue; + if (angle < 0) angle += M_PI; + if (fabs(angle - Angle) < RAD_DEG) { + X += R[j]; + Y += D[j]; + X2 += R[j]*R[j]; + Y2 += D[j]*D[j]; + N ++; + XY += R[j]*D[j]; + } + } + } + /* points to the west */ + for (j = i - 1; !mark[i] && (j >= 0) && (Ra - R[j] < RADIUS); j--) { + if (mark[j]) continue; + dD = D[j] - De; + if (fabs(dD) < RADIUS) { + dR = R[j] - Ra; + angle = atan2 (dD,dR); + if (!finite(angle)) continue; + if (angle < 0) angle += M_PI; + if (fabs(angle - Angle) < RAD_DEG) { + X += R[j]; + Y += D[j]; + X2 += R[j]*R[j]; + Y2 += D[j]*D[j]; + N ++; + XY += R[j]*D[j]; + } + } + } + /* determine coeffs */ + Flipped = 0; + det = 1.0 / (X2*N - X*X); + m = det * (XY*N - X*Y); + b = det * (X2*Y - XY*X); + if (fabs(m) > 1.1) { /* use a line of R = m*D + b instead */ + /* fprintf (stderr, "high slope object: %f %f -> ", m, b); */ + det = 1.0 / (Y2*N - Y*Y); + m = det * (XY*N - X*Y); + b = det * (Y2*X - XY*Y); + Flipped = 1; + /* fprintf (stderr, "%f %f\n", m, b); */ + } + + *M = m; + *B = b; + return (Flipped); + +} + + +find_better_line (catstats, mark, Npts, i, M, B, axis) + CatStats catstats[]; + double *M, *B; + int i, Npts, axis; + char *mark; +{ + + float *R, *D; + int j, N; + double X, Y, X2, Y2, XY, m, b, det; + double dR, dD, Ra, De, delta; + + R = catstats[0].X; + D = catstats[0].Y; + + /* fit a line to points near line */ + Ra = R[i]; + De = D[i]; + X = Y = X2 = Y2 = XY = N = 0; + m = *M; b = *B; + + /* points to the east */ + for (j = i; (j < Npts) && (R[j] - Ra < RADIUS); j++) { + dD = D[j] - De; + if (fabs(dD) < RADIUS) { + dR = R[j] - Ra; + if (axis == 1) + delta = R[j] - m*D[j] - b; + else + delta = D[j] - m*R[j] - b; + if (fabs(delta) < 2*TRAIL_WIDTH) { + X += R[j]; + Y += D[j]; + X2 += R[j]*R[j]; + Y2 += D[j]*D[j]; + N ++; + XY += R[j]*D[j]; + } + } + } + /* points to the west */ + for (j = i - 1; (j >= 0) && (Ra - R[j] < RADIUS); j--) { + dD = D[j] - De; + if (fabs(dD) < RADIUS) { + dR = R[j] - Ra; + if (axis == 1) + delta = R[j] - m*D[j] - b; + else + delta = D[j] - m*R[j] - b; + if (fabs(delta) < 2*TRAIL_WIDTH) { + X += R[j]; + Y += D[j]; + X2 += R[j]*R[j]; + Y2 += D[j]*D[j]; + N ++; + XY += R[j]*D[j]; + } + } + } + /* determine coeffs */ + if (axis == 0) { + det = 1.0 / (X2*N - X*X); + m = det * (XY*N - X*Y); + b = det * (X2*Y - XY*X); + } else { + det = 1.0 / (Y2*N - Y*Y); + m = det * (XY*N - X*Y); + b = det * (Y2*X - XY*Y); + } + + /* fprintf (stderr, "%f %f %d\n", m, b, N); */ + + *M = m; + *B = b; + +} Index: /tags/fixcat-1-0/Ohana/src/fixcat/src/find_matches.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/find_matches.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/find_matches.c (revision 2512) @@ -0,0 +1,87 @@ +# include "markstar.h" + +find_matches (catalog, catstats, Ghostdata, Nimage) +CatStats catstats[]; +Catalog catalog[], Ghostdata[]; +int Nimage; +{ + + int i, j, k, n, m, N, M, first_j; + double X, Y, RADIUS, RADIUS2; + float *X1, *Y1, *X2, *Y2; + float dX, dY, dR; + int *N1, *N2; + int Nstar, Nghost, Ng; + unsigned int flags; + Coords *tcoords; + + Nstar = catalog[0].Naverage; + Nghost = Ghostdata[0].Naverage; + + if (Nghost < 1) return (0); + + /* it is better to have the catalog with fewer stars + assigned to the X1 set */ + X2 = catstats[0].X; + Y2 = catstats[0].Y; + N2 = catstats[0].N; + + ALLOCATE (X1, float, Nghost); + ALLOCATE (Y1, float, Nghost); + ALLOCATE (N1, int, Nghost); + + /* project ghosts to the frame of the catalog */ + tcoords = &catstats[0].coords; + for (i = 0; i < Nghost; i++) { + fRD_to_XY (&X1[i], &Y1[i], Ghostdata[0].average[i].R, Ghostdata[0].average[i].D, tcoords); + N1[i] = i; + } + if (Nghost > 1) sort_lists (X1, Y1, N1, Nghost); + + /* choose a radius for matches */ + RADIUS2 = GHOST_RADIUS*GHOST_RADIUS; + + /** find matched stars **/ + for (i = j = 0; (i < Nghost) && (j < Nstar); ) { + dX = X1[i] - X2[j]; + if (dX <= -2*GHOST_RADIUS) { + i++; + continue; + } + if (dX >= 2*GHOST_RADIUS) { + j++; + continue; + } + /* negative dX: j is too large, positive dX, i is too large */ + first_j = j; + for (; (dX > -2*GHOST_RADIUS) && (j < Nstar); j++) { + dX = X1[i] - X2[j]; + dY = Y1[i] - Y2[j]; + dR = dX*dX + dY*dY; + if (dR < RADIUS2) { + /* this object may be a ghost star, + but only mark those measurements on the correct image */ + M = N2[j]; + m = catalog[0].average[M].offset; + Ng = 0; + for (k = 0; k < catalog[0].average[M].Nm; k++) { + if (catalog[0].image[m+k] == Nimage) { + catalog[0].measure[m+k].average |= GHOST_DATA; + Ng ++; + } + } + /* all measurements are ghosts */ + if (catalog[0].average[M].Nm == Ng) { + catalog[0].average[M].code = ID_GHOST; + } + } + } + j = first_j; + i++; + } + + free (X1); + free (Y1); + free (N1); + +} Index: /tags/fixcat-1-0/Ohana/src/fixcat/src/find_trails.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/find_trails.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/find_trails.c (revision 2512) @@ -0,0 +1,41 @@ +# include "markstar.h" + +find_trails (catalog, catstats) +Catalog catalog[]; +CatStats catstats[]; +{ + + int i, j, N, Nave, axis, marked; + double density, spacing, Area, Angle, m, b, RaCenter, DecCenter; + double MinRA, MaxRA, MinDEC, MaxDEC; + float *X1, *Y1; + char *mark; + int *N1; + Coords tcoords; + + Nave = catalog[0].Naverage; + N1 = catstats[0].N; + + ALLOCATE (mark, char, Nave); + bzero (mark, Nave); + + for (i = 0; i < Nave; i++) { + /* already marked, ignore */ + if ((mark[i]) || + (catalog[0].average[N1[i]].code == ID_BLEED) || + (catalog[0].average[N1[i]].code == ID_GHOST)) + continue; + /* a good star, ignore */ + if ((catalog[0].average[N1[i]].Nm > 3) && + (catalog[0].average[N1[i]].Nm > 4*catalog[0].average[N1[i]].Nn)) { + continue; + } + if (find_group (catstats, mark, Nave, i, &Angle)) { + /* this point has an excess nearby concentration, find the line */ + axis = find_line (catstats, mark, Nave, i, &m, &b, Angle); + find_better_line (catstats, mark, Nave, i, &m, &b, axis); + mark_trail (catstats, mark, Nave, i, m, b, axis, catalog); + } + } +} + Index: /tags/fixcat-1-0/Ohana/src/fixcat/src/fixcat.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/fixcat.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/fixcat.c (revision 2512) @@ -0,0 +1,49 @@ +# include "markstar.h" + +main (argc, argv) +int argc; +char **argv; +{ + + FILE *f; + int i, Nimages; + Image *image, *find_images(); + Catalog catalog; + CatStats catstats; + 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); + } + if (argc < 2) { + fprintf (stderr, "ERROR: Usage: fixcat (catalog)\n"); + exit (0); + } + + /* if lockfile exists, program will complain and quit */ + check_lockfile (); + + gcatalog (argv[1], &catalog); + + gcatstats (&catalog, &catstats); + + image = find_images (&catstats, &Nimages); + + match_images (&catalog, image, Nimages); + + find_funnymags (&catalog, image, Nimages); + + 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/fixcat-1-0/Ohana/src/fixcat/src/gcatalog.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/gcatalog.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/gcatalog.c (revision 2512) @@ -0,0 +1,108 @@ +# include "markstar.h" + +gcatalog (catname, catalog) +char *catname; +Catalog catalog[]; +{ + + char filename[128], line[64]; + int Nitems, nitems; + int i, Nmeas, Nmiss, done; + FILE *f; + struct tm *local; + struct timeval now; + unsigned int NotTrail; + unsigned short NotBad; + + sprintf (filename, "%s/%s\0", CATDIR, catname); + + /* read catalog header */ + if (!fits_read_header (filename, &catalog[0].header)) { + fprintf (stderr, "ERROR: file doesn't exist %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)); + if (catalog[0].Naverage == 0) { + /* this is a valid ending state: no stars */ + if (VERBOSE) fprintf (stderr, "SUCCESS: no stars yet in catalog %s\n", filename); + fclose (f); + clear_lockfile (); + 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); + } + + if (RESET) { + NotBad = (0xffff ^ 0xc003); + for (i = 0; i < catalog[0].Naverage; i++) { + if ((catalog[0].average[i].code & ID_MOVING) && + (catalog[0].average[i].code & ID_BAD_OBJECT) && + (catalog[0].average[i].code & 0x0003)) { + /* this will set the correct bit for each to 0 */ + catalog[0].average[i].code &= NotBad; + } + } + for (i = 0; i < catalog[0].Nmeasure; i++) { + catalog[0].measure[i].average &= ~PART_OF_TRAIL; + } + } + + + fclose (f); + +} Index: /tags/fixcat-1-0/Ohana/src/fixcat/src/gcatstats.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/gcatstats.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/gcatstats.c (revision 2512) @@ -0,0 +1,66 @@ +# include "markstar.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); + + /* double check on region RA and DEC ranges */ + DecCenter = 0.5*(MinDEC + MaxDEC); + RaCenter = 0.5*(MinRA + MaxRA); + if (MaxDEC > 86.25) { /* we are on the pole */ + DecCenter = 90.0; + RaCenter = 0.0; + } + + catstats[0].RA[0] = MinRA; + catstats[0].RA[1] = MaxRA; + catstats[0].DEC[0] = MinDEC; + catstats[0].DEC[1] = MaxDEC; + catstats[0].Area = (MaxDEC - MinDEC)*(MaxRA - MinRA) / cos (RAD_DEG*DecCenter); + if (MaxDEC > 86.25) { /* we are on the pole */ + catstats[0].Area = 44.2; + } + catstats[0].density = catalog[0].Naverage / (3600*3600*catstats[0].Area); + catstats[0].spacing = 1.0 / (NSIGMA * catstats[0].density * 2*TRAIL_WIDTH); + fprintf (stderr, "Area, density, spacing: %f %f %f\n", catstats[0].Area, + catstats[0].density, catstats[0].spacing); + /* number of stars per square arcsec */ + + /** 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; + 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, &catstats[0].coords); + 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/fixcat-1-0/Ohana/src/fixcat/src/gregions.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/gregions.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/gregions.c (revision 2512) @@ -0,0 +1,134 @@ +# include "markstar.h" + +GSCRegion *gregions1 (catstats, Nregions) +CatStats catstats[]; +int *Nregions; +{ + + GSCRegion *region; + int i, j, k, x, y, done, nregion, nregion2, NREGION; + double ra, dec, dx, dy, Xo[4], Yo[4]; + FILE *f; + + f = fopen (GSCFILE, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: could not open GSC region file %s\n", GSCFILE); + exit (0); + } + nregion = 0; + NREGION = 10; + ALLOCATE (region, GSCRegion, NREGION); + + if (catstats[0].DEC[0] == 86.25) { /* pole region */ + dec = 86.0; + for (ra = 0.1; ra < 370; ra+= 30.0) { + aregion (®ion[nregion], f, ra, dec); + done = FALSE; + for (j = 0; (j < nregion - 1) && !done; j++) { + if (!strcmp (region[nregion].filename, region[j].filename)) { + nregion --; + done = TRUE; + } + } + nregion ++; + if (nregion == NREGION) { + NREGION += 10; + REALLOCATE (region, GSCRegion, NREGION); + } + } + } else { + for (x = 0; x < 2; x++) { + for (y = 0; y < 2; y++) { + for (dx = -0.1; dx <= 0.1; dx += 0.2) { + for (dy = -0.1; dy <= 0.1; dy += 0.2) { + ra = catstats[0].RA[x] + dx; + dec = catstats[0].DEC[y] + dy; + aregion (®ion[nregion], f, ra, dec); + done = FALSE; + for (j = 0; (j < nregion) && !done; j++) { + if (!strcmp (region[nregion].filename, region[j].filename)) { + done = TRUE; + } + } + if (!done) { + nregion ++; + } + if (nregion == NREGION) { + NREGION += 10; + REALLOCATE (region, GSCRegion, NREGION); + } + } + } + } + } + } + + if (VERBOSE) { + fprintf (stderr, "using %d regions\n", nregion); + for (i = 0; i < nregion; i++) { + fprintf (stderr, "region %d: %f %f %f %f\n", i, region[i].RA[0], region[i].RA[1], region[i].DEC[0], region[i].DEC[1]); + } + } + + REALLOCATE (region, GSCRegion, MAX (nregion, 1)); + *Nregions = nregion; + + fclose (f); + return (region); + +} + +GSCRegion *gregions2 (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); + 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/fixcat-1-0/Ohana/src/fixcat/src/load_gsc_data.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/load_gsc_data.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/load_gsc_data.c (revision 2512) @@ -0,0 +1,61 @@ +# include "markstar.h" +# define NBYTES 160000 +# define BYTES_STAR 23 +# define BLOCK 1000 +# define DNSTARS 1000 + +load_gsc_data (catalog, catstats) +Catalog catalog[]; +CatStats catstats[]; +{ + + /* load data from the GSC files */ + char filename[128]; + char *tbuffer; + int nstar, NSTARS; + int i, j, Nbytes, nbytes, Nregions; + double R, D, M, MagLimit; + FILE *f; + GSCRegion *region, *gregions1(), *gregions2(); + + region = gregions1 (catstats, &Nregions); + + MagLimit = MAX (MAX (BRIGHT_HALO_MAG, BRIGHT_XTRAIL_MAG), BRIGHT_YTRAIL_MAG); + nstar = 0; + NSTARS = DNSTARS; + ALLOCATE (tbuffer, char, (BLOCK*BYTES_STAR)); + ALLOCATE (catalog[0].average, Average, NSTARS); + Nbytes = BLOCK*BYTES_STAR; + + for (j = 0; j < Nregions; j++) { + f = fopen (region[j].filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "GSC file for region %s missing\n", region[j].filename); + exit (0); + } + + while ((nbytes = fread (tbuffer, 1, Nbytes, f)) > 0) { + for (i = 0; i < nbytes / BYTES_STAR; i++) { + dparse (&M, 3, &tbuffer[i*BYTES_STAR]); + if (M > MagLimit) continue; + dparse (&R, 1, &tbuffer[i*BYTES_STAR]); + dparse (&D, 2, &tbuffer[i*BYTES_STAR]); + catalog[0].average[nstar].R = R; + catalog[0].average[nstar].D = D; + catalog[0].average[nstar].M = M * 1000.0; + nstar++; + if (nstar == NSTARS - 1) { + NSTARS += DNSTARS; + REALLOCATE (catalog[0].average, Average, NSTARS); + } + } + } + fclose (f); + } + + free (tbuffer); + REALLOCATE (catalog[0].average, Average, nstar); + catalog[0].Naverage = nstar; + +} + Index: /tags/fixcat-1-0/Ohana/src/fixcat/src/load_gsc_data_ghost.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/load_gsc_data_ghost.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/load_gsc_data_ghost.c (revision 2512) @@ -0,0 +1,75 @@ +# include "markstar.h" +# define NBYTES 160000 +# define BYTES_STAR 23 +# define BLOCK 1000 +# define DNSTARS 1000 + +/* this routine is basically identical to load_gsc_data, but + it is not limited to a single region file, and it merges the + results */ + +load_gsc_data_ghost (catalog, region, Nregion, image) +Catalog catalog[]; +GSCRegion *region; +int Nregion; +Image image[]; +{ + + /* load data from the GSC files */ + char filename[128]; + char *tbuffer; + int nstar, NSTARS; + int i, j, Nbytes, nbytes; + double R, D, M, X, Y; + FILE *f; + Coords *tcoords; + int MinX, MinY, MaxX, MaxY; + + nstar = 0; + NSTARS = DNSTARS; + ALLOCATE (tbuffer, char, (BLOCK*BYTES_STAR)); + ALLOCATE (catalog[0].average, Average, NSTARS); + tcoords = &image[0].coords; + MinX = 0; + MinY = 0; + MaxX = image[0].NX; + MaxY = image[0].NY; + + for (j = 0; j < Nregion; j++) { + + f = fopen (region[j].filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "no GSC file for region %s (2)\n", region[j].filename); + exit (0); + } + + Nbytes = BLOCK*BYTES_STAR; + while ((nbytes = fread (tbuffer, 1, Nbytes, f)) > 0) { + for (i = 0; i < nbytes / BYTES_STAR; i++) { + dparse (&M, 3, &tbuffer[i*BYTES_STAR]); + if (M > GHOST_MAG) continue; + dparse (&R, 1, &tbuffer[i*BYTES_STAR]); + dparse (&D, 2, &tbuffer[i*BYTES_STAR]); + RD_to_XY (&X, &Y, R, D, tcoords); + if ((X < MinX) || (X > MaxX) || (Y < MinY) || (Y > MaxY)) continue; + catalog[0].average[nstar].R = R; + catalog[0].average[nstar].D = D; + catalog[0].average[nstar].M = M * 1000.0; + nstar++; + if (nstar == NSTARS - 1) { + NSTARS += DNSTARS; + REALLOCATE (catalog[0].average, Average, NSTARS); + } + } + } + fclose (f); + + } + + free (tbuffer); + REALLOCATE (catalog[0].average, Average, MAX (nstar, 1)); + catalog[0].Naverage = nstar; + +} + + Index: /tags/fixcat-1-0/Ohana/src/fixcat/src/mark_trail.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/mark_trail.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/mark_trail.c (revision 2512) @@ -0,0 +1,272 @@ +# include "markstar.h" +int *make_common_list (); +int *check_common_list (); + +mark_trail (catstats, mark, Nave, i, m, b, axis, catalog) + CatStats catstats[]; + char *mark; + int i, Nave, axis; + double m, b; + Catalog catalog[]; +{ + + float *R, *D; + int j, k, jj, kk, N, NPTS, marked, start, end, Nm; + int *good, *seq, *N1; + double *dist, *dist2; + double d2, di, Di, dD, n, Dist, scale, norm; + double spacing; + int *list, Nlist; + int M, Nmeas, thisimage; + + R = catstats[0].X; + D = catstats[0].Y; + N1 = catstats[0].N; + spacing = catstats[0].spacing; + + NPTS = 200; + ALLOCATE (good, int, NPTS); + ALLOCATE (dist, double, NPTS); + ALLOCATE (dist2, double, NPTS); + ALLOCATE (seq, int, NPTS); + + /* Find all points which lie near line */ + /* save the entry number and distance along line */ + N = 0; + scale = sqrt (1.0 + m*m); + + if (axis == 1) { + for (j = 0; j < Nave; j++) { + norm = scale * fabs(R[j] - m*D[j] - b); + if (!mark[j] && (norm < TRAIL_WIDTH)) { + good[N] = j; + dist[N] = scale * (D[j] - D[i] + m*(R[j] - R[i])); + seq[N] = N; + N++; + if (N == NPTS - 1) { + NPTS += 200; + REALLOCATE (good, int, NPTS); + REALLOCATE (dist, double, NPTS); + REALLOCATE (dist2, double, NPTS); + REALLOCATE (seq, int, NPTS); + } + } + } + } else { + for (j = 0; j < Nave; j++) { + norm = scale * fabs(D[j] - m*R[j] - b); + if (!mark[j] && (norm < TRAIL_WIDTH)) { + good[N] = j; + dist[N] = scale * (R[j] - R[i] + m*(D[j] - D[i])); + seq[N] = N; + N++; + if (N == NPTS - 1) { + NPTS += 200; + REALLOCATE (good, int, NPTS); + REALLOCATE (dist, double, NPTS); + REALLOCATE (dist2, double, NPTS); + REALLOCATE (seq, int, NPTS); + } + } + } + } + + if (N < NPTSINLINE) + return (0); + + sort_seq (dist, seq, N); + + start = -1; end = -1; + for (j = 0; j < N-1; j++) { + /* if we have part of a line, and next point is in the line, check for common images */ + if ((start != -1) && (fabs(dist[j] - dist[j+1]) < spacing)) { + list = check_common_list (list, N1[good[seq[j+1]]], &Nlist, catalog); + if (Nlist == 0) { /* if no common images, dump list and continue */ + end = j + 1; + if (end - start > NPTSINLINE) { + if (axis == 0) + fprintf (stderr, "marking line %f %f %d pts %d %d %d, %f %f\n", m, b, N, end-start, start, end, dist[start], dist[end-1]); + else + fprintf (stderr, "marking line %f %f %d pts %d %d %d, %f %f\n", 1.0/m, -1.0*b/m, N, end-start, start, end, dist[start], dist[end-1]); + for (k = start; k < end; k++) { + M = N1[good[seq[k]]]; + mark[good[seq[k]]] = TRUE; + /* we need to mark measurements from all images in common on the line */ + Nm = 0; + for (jj = 0; jj < catalog[0].average[M].Nm; jj++) { + Nmeas = catalog[0].average[M].offset + jj; + thisimage = catalog[0].image[Nmeas]; + for (kk = 0; kk < Nlist; kk++) { + if (thisimage == list[kk]) { + catalog[0].measure[Nmeas].average |= PART_OF_TRAIL; + Nm ++; + } + } + } + /* if there is only 1 measurement, mark object as bad */ + if (catalog[0].average[M].Nm == Nm) { + catalog[0].average[M].code = ID_TRAIL; + } + catalog[0].average[M].code = ID_TRAIL; + } + } + start = -1; + end = -1; + } + } + /* if we haven't yet found a line segment, check for the beginning */ + if ((start < 0) && (fabs(dist[j] - dist[j+1]) < spacing)) { + start = j; + list = make_common_list (N1[good[seq[j]]], N1[good[seq[j+1]]], &Nlist, catalog); + if (Nlist == 0) { /* if no common images, move on */ + start = -1; + free (list); + } + } + /* if we have a complete line, check for validity. if it has enough members, + mark them and continue searching for lines */ + if ((start != -1) && ((fabs(dist[j] - dist[j+1]) >= spacing) || (j == N-2))) { + end = j + 1; + if (end - start > NPTSINLINE) { + if (axis == 0) + fprintf (stderr, "marking line %f %f %d pts %d %d %d, %f %f\n", m, b, N, end-start, start, end, dist[start], dist[end-1]); + else + fprintf (stderr, "marking line %f %f %d pts %d %d %d, %f %f\n", 1.0/m, -1.0*b/m, N, end-start, start, end, dist[start], dist[end-1]); + for (k = start; k < end; k++) { + M = N1[good[seq[k]]]; + mark[good[seq[k]]] = TRUE; + if (catalog[0].average[M].code == ID_BLEED) continue; + /* we need to mark measurements from all images in common on the line */ + Nm = 0; + for (jj = 0; jj < catalog[0].average[M].Nm; jj++) { + Nmeas = catalog[0].average[M].offset + jj; + thisimage = catalog[0].image[Nmeas]; + for (kk = 0; kk < Nlist; kk++) { + if (thisimage == list[kk]) { + catalog[0].measure[Nmeas].average |= PART_OF_TRAIL; + Nm ++; + } + } + } + /* if there is only 1 measurement, mark object as bad */ + if (catalog[0].average[M].Nm == Nm) { + catalog[0].average[M].code = ID_TRAIL; + } + catalog[0].average[M].code = ID_TRAIL; + } + } + free (list); + start = -1; + end = -1; + } + } +} + + +/* I is Average seq number for star 1, J for star 2 */ +/* make a list of images in common between two measurements */ +int *make_common_list (I, J, Nlist, catalog) +int I, J, *Nlist; +Catalog catalog[]; +{ + + int i, j, N1, N2, nlist, NLIST; + int *list; + int k, already; + + NLIST = 50; + ALLOCATE (list, int, NLIST); + nlist = 0; + + for (i = 0; i < catalog[0].average[I].Nm; i++) { + N1 = catalog[0].average[I].offset + i; + if (catalog[0].image[N1] == -1) + continue; /* not a real measurement */ + for (j = 0; j < catalog[0].average[J].Nm; j++) { + N2 = catalog[0].average[J].offset + j; + if (catalog[0].image[N2] == -1) + continue; /* not a real measurement */ + if (catalog[0].image[N1] == catalog[0].image[N2]) { + already = FALSE; + for (k = 0; !already && (k < nlist); k++) { + if (catalog[0].image[N1] == list[k]) { + already = TRUE; + } + } + if (!already) { + list[nlist] = catalog[0].image[N1]; + nlist ++; + if (nlist == NLIST - 1) { + NLIST += 50; + REALLOCATE (list, int, NLIST); + } + } + } + } + } + + REALLOCATE (list, int, MAX(nlist, 1)); + *Nlist = nlist; + return (list); + +} + + + +/* J is Average seq number for star */ + +int *check_common_list (inlist, J, Nlist, catalog) +int *inlist, J, *Nlist; +Catalog catalog[]; +{ + + int i, j, N2, Ninlist; + int *list, NLIST, nlist; + int already, found, k; + + NLIST = 50; + ALLOCATE (list, int, NLIST); + nlist = 0; + Ninlist = *Nlist; + + for (j = 0; j < catalog[0].average[J].Nm; j++) { + N2 = catalog[0].average[J].offset + j; + found = FALSE; + for (i = 0; !found && (i < Ninlist); i++) { + if (catalog[0].image[N2] == -1) + continue; /* not a real measurement */ + if (inlist[i] == catalog[0].image[N2]) { + found = TRUE; + already = FALSE; + for (k = 0; !already && (k < nlist); k++) { + if (inlist[i] == list[k]) { + already = TRUE; + } + } + if (!already) { + list[nlist] = inlist[i]; + nlist ++; + if (nlist == NLIST - 1) { + NLIST += 50; + REALLOCATE (list, int, NLIST); + } + } + } + } + } + + /* if there are no common images, return the input list */ + if (nlist != 0) { + free (inlist); + REALLOCATE (list, int, MAX(nlist, 1)); + *Nlist = nlist; + return (list); + } else { + return (inlist); + } + +} + + +/* measurements associated with image number -1 are not + measurements we made, but are added, like USNO */ Index: /tags/fixcat-1-0/Ohana/src/fixcat/src/markstar.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/markstar.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/markstar.c (revision 2512) @@ -0,0 +1,66 @@ +# include "markstar.h" + +main (argc, argv) +int argc; +char **argv; +{ + + FILE *f; + int i, Nstars, Nimage, Nregions, Nmissed; + Image *image, *find_images(); + Catalog catalog; + CatStats catstats; + 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); + } + 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: markstar (catalog)\n"); + exit (0); + } + + /* if lockfile exists, program will complain and quit */ + check_lockfile (); + check_permissions (argv[1]); + + gcatalog (argv[1], &catalog); + + gcatstats (&catalog, &catstats); + + image = find_images (&catstats, &Nimage); + + match_images (&catalog, image, Nimage); + + find_bright_stars (&catalog, &catstats); + + find_ghosts (&catalog, &catstats, argv[1], image, Nimage); + + find_trails (&catalog, &catstats); + + 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/fixcat-1-0/Ohana/src/fixcat/src/match_images.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/match_images.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/match_images.c (revision 2512) @@ -0,0 +1,57 @@ +# include "markstar.h" + +match_images (catalog, image, Nimage) +Catalog catalog[]; +Image image[]; +int Nimage; +{ + + int j, k, found; + unsigned int *start, *stop; + + /* this must be allocated so future free will not fail */ + ALLOCATE (catalog[0].image, int, MAX (catalog[0].Nmeasure, 1)); + if (catalog[0].Naverage == 0) { + if (VERBOSE) fprintf (stderr, "no stars in catalog, skipping\n"); + return (FALSE); + } + + ALLOCATE (start, unsigned int, Nimage); + ALLOCATE (stop, unsigned int, Nimage); + for (j = 0; j < Nimage; j++) { + start[j] = image[j].tzero - MAX(0.05*image[j].trate*image[j].NY, 1); + stop[j] = image[j].tzero + MAX(1.05*image[j].trate*image[j].NY, 1); + } + + for (j = 0; j < catalog[0].Nmeasure; j++) { + found = FALSE; + if (catalog[0].measure[j].t == 0) { + catalog[0].image[j] = -1; + found = TRUE; + } + for (k = 0; (k < Nimage) && !found; k++) { + if ((catalog[0].measure[j].t >= start[k]) && + (catalog[0].measure[j].t <= stop[k])) { + catalog[0].image[j] = k; + found = TRUE; + } + } + if (!found) { + fprintf (stderr, "ERROR: can't find source image for this measurement: %d\n", + catalog[0].measure[j].t); + exit (0); + } + } + free (start); + free (stop); +} + + /* this routine uses the time of each measurement to match the +measurement with an image. Since the measurement is only store to 1 +sec accuracy, which corresponds to roughly 30 rows at nominal speed, +we can't tell exactly which image the star come from. However, this +doesn't matter, and in fact this helps a bit: a measurement from the +top of one image is the same as from the bottom of the next. +Therefore, we intentionally blur the edges of the images by 5%, which +will help to tie together neighboring images... */ + Index: /tags/fixcat-1-0/Ohana/src/fixcat/src/sorts.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/sorts.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/sorts.c (revision 2512) @@ -0,0 +1,175 @@ + + +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; + } +} + +sort_entries (X, Y, N) + double *X, *Y; + int N; +{ + int l,j,ir,i; + double tX, tY; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tX = X[l]; + tY = Y[l]; + } + else { + tX = X[ir]; + X[ir] = X[0]; + tY = Y[ir]; + Y[ir] = Y[0]; + if (--ir == 0) { + X[0] = tX; + Y[0] = tY; + 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]; + j += (i=j) + 1; + } + else j = ir + 1; + } + X[i] = tX; + Y[i] = tY; + } +} + + + +sort (X, N) + double *X; + int N; +{ + int l,j,ir,i; + double tX; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tX = X[l]; + } + else { + tX = X[ir]; + X[ir] = X[0]; + if (--ir == 0) { + X[0] = tX; + 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]; + j += (i=j) + 1; + } + else j = ir + 1; + } + X[i] = tX; + } +} + + +sort_seq (X, S, N) + double *X; + int *S; + int N; +{ + int l,j,ir,i; + int tS; + double tX; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tX = X[l]; + tS = S[l]; + } + else { + tX = X[ir]; + X[ir] = X[0]; + tS = S[ir]; + S[ir] = S[0]; + if (--ir == 0) { + X[0] = tX; + 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]; + S[i] = S[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + X[i] = tX; + S[i] = tS; + } +} + Index: /tags/fixcat-1-0/Ohana/src/fixcat/src/wcatalog.c =================================================================== --- /tags/fixcat-1-0/Ohana/src/fixcat/src/wcatalog.c (revision 2512) +++ /tags/fixcat-1-0/Ohana/src/fixcat/src/wcatalog.c (revision 2512) @@ -0,0 +1,74 @@ +# include "markstar.h" + +wcatalog (catname, catalog) +char *catname; +Catalog catalog[]; +{ + + int i, Nitems, nitems, status, mode; + char filename[128], line[256]; + 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, "MARKSTAR", "%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); + +} +